We present the Flagship galaxy mock a simulated catalogue of billions of galaxies designed to support the scientific exploitation of the Euclid mission. Euclid is a mediumclass mission of the European Space Agency optimised to determine the properties of dark matter and dark energy on the largest scales of the Universe. It probes structure formation over more than 10 billion years primarily from the combination of weak gravitational lensing and galaxy clustering data. The breath of Euclids data will also foster a wide variety of scientific analyses. The Flagship simulation was developed to provide a realistic approximation to the galaxies that will be observed by Euclid and used in its scientific analyses. We ran a stateoftheart Nbody simulation with four trillion particles producing a lightcone on the fly. From the dark matter particles we produced a catalogue of 16 billion haloes in one octant of the sky in the lightcone up to redshift z3. We then populated these haloes with mock galaxies using a halo occupation distribution and abundance matching approach calibrating the free parameters of the galaxy mock against observed correlations and other basic galaxy properties. Modelled galaxy properties include luminosity and flux in several bands redshifts positions and velocities spectral energy distributions shapes and sizes stellar masses star formation rates metallicities emission line fluxes and lensing properties. We selected a final sample of 3.4 billion galaxies with a magnitude cut of HElt26 where we are complete. We have performed a comprehensive set of validation tests to check the similarity to observational data and theoretical models. In particular our catalogue is able to closely reproduce the main characteristics of the weak lensing and galaxy clustering samples to be used in the missions main cosmological analysis. abridged
['Astrophysics - Cosmology and Nongalactic Astrophysics']
Euclid. V. The Flagship galaxy mock catalogue a comprehensive simulation for the Euclid mission
{'Euclid . V. The Flagship galaxy mock catalogue: a comprehensive simulation for the Euclid mission': "V. Kansal 157 , 158 , D. Karagiannis 159 , 160 , C. C. Kirkpatrick 94 , J. Le Graet 82 , L. Legrand 161 , J. Lesgourgues 70 , T. I. Liaudat 136 , A. Loureiro 162 , 163 , J. Macias-Perez 164 , M. Magliocchetti 79 , C. Mancini 69 , F. Mannucci 165 , R. Maoli 166 , 71 , C. J. A. P. Martins 167 , 65 , L. Maurin 52 , R. B. Metcalf 103 , 6 , M. Migliaccio 168 , 169 , M. Miluzio 53 , A. Mora 170 , C. Moretti 30 , 132 , 16 , 18 , 17 , G. Morgante 6 , S. Nadathur 148 , L. Nicastro 6 , Nicholas A. Walton 125 , M. Oguri 171 , 172 , L. Patrizii 57 , V. Popa 112 , A. Pourtsidou 35 , 173 , P. Reimberg 8 , I. Risso 121 , P.-F. Rocci 52 , R. P. Rollins 35 , B. Rusholme 174 , M. Sahlén 175 , A. G. Sánchez 37 , J. Schaye 24 , J. A. Schewtschenko 35 , A. Schneider 3 , M. Schultheis 106 , M. Sereno 6 , 57 , F. Shankar 176 , A. Shulevski 177 , 178 , 179 , A. Silvestri 143 , P. Simon 102 , A. Spurio Mancini 180 , 39 , S. A. Stanford 181 , K. Tanidis 34 , C. Tao 82 , N. Tessore 93 , G. Testera 62 , M. Tewes 102 , S. Toft 90 , 45 , 44 , S. Tosi 61 , 62 , A. Troja 114 , 80 , M. Tucci 43 , C. Valieri 57 , J. Valiviita 96 , 97 , D. Vergani 6 , F. Vernizzi 122 , G. Verza 182 , 127 , P. Vielzeuf 82 , J. R. Weaver 183 , L. Zalesky 144 , P. Dimauro 71 , 184 P.-A. Duc 139 , Y. Fang 59 , A. M. N. Ferguson 35 , C. M. Gutierrez 185 , I. Kovaˇci'c 186 , S. Kruk 53 , A. M. C. Le Brun 10 , A. Montoro 1 , 2 , C. Murray 107 , L. Pagano 137 , 138 , D. Paoletti 6 , 120 , E. Sarpa 30 , 132 , 17 , A. Viitanen 94 , 71 , J. Martín-Fleitas 170 , and L. Y. A. Yung 92 , 187  \n(A ffi liations can be found after the references)  \nMay 24, 2024", 'ABSTRACT': "We present the Flagship galaxy mock, a simulated catalogue of billions of galaxies designed to support the scientific exploitation of the Euclid mission. Euclid is a medium-class mission of the European Space Agency optimised to determine the properties of dark matter and dark energy on the largest scales of the Universe. It probes structure formation over more than 10 billion years primarily from the combination of weak gravitational lensing and galaxy clustering data. The breath of Euclid 's data will also foster a wide variety of scientific analyses. The Flagship simulation was developed to provide a realistic approximation to the galaxies that will be observed by Euclid and used in its scientific analyses. We ran a state-of-the-art N -body simulation with four trillion particles, producing a lightcone on the fly. From the dark matter particles, we produced a catalogue of 16 billion haloes in one octant of the sky in the lightcone up to redshift z = 3. We then populated these haloes with mock galaxies using a halo occupation distribution and abundance matching approach, calibrating the free parameters of the galaxy mock against observed correlations and other basic galaxy properties. Modelled galaxy properties include luminosity and flux in several bands, redshifts, positions and velocities, spectral energy distributions, shapes and sizes, stellar masses, star formation rates, metallicities, emission line fluxes, and lensing properties. We selected a final sample of 3.4 billion galaxies with a magnitude cut of H E < 26, where we are complete. We have performed a comprehensive set of validation tests to check the similarity to observational data and theoretical models. In particular, our catalogue is able to closely reproduce the main characteristics of the weak lensing and galaxy clustering samples to be used in the mission's main cosmological analysis. Moreover, given its depth and completeness, this new galaxy mock also provides the community with a powerful tool for developing a wide range of scientific analyses beyond the Euclid mission.  \nKey words. Cosmology: observations - large-scale structure of Universe - Gravitational lensing: weak - Galaxies: evolution - Catalogues", '1. Introduction': "The discovery of the accelerated expansion of the Universe has driven a large observational e ff ort to study its cause and nature (Albrecht et al. 2006; Amendola et al. 2018). This phenomenon, usually referred to as dark energy, can be the result of a hypothesised fluid with negative pressure or the inadequacy of our current understanding of gravitation. Large observational surveys are needed to sample enough volume and a high number density of sources to properly characterize the Universe's evolution from its expansion rate and the growth rate of its structures. Current large surveys such as the Dark Energy Survey (DES; Abbott et al. 2023), the Hyper Suprime Cam Subaru Strategic Program (HSC-SSP; Aihara et al. 2018), the Kilo-Degree Survey (KiDS; Heymans et al. 2021) and the Dark Energy Spectroscopic Instrument (DESI; Dey et al. 2019) are providing data that in combination with cosmic microwave background (CMB) data place strong constraints on cosmological parameters (e.g., Planck Collaboration: Aghanim et al. 2020). Forthcoming surveys from the ground, such as the Vera C. Rubin Observatory (LSST; Ivezi'c et al. 2019), or from space, like Euclid (Laureijs et al. 2011; Euclid Collaboration: Mellier et al. 2024) and the Nancy Grace Roman Space Telescope (Akeson et al. 2019), will collect more and higher-quality data that will allow us to determine the cosmolog-  \ncal parameters, and in particular the equation of state parameter of dark energy, to unprecedented precision. With the gain in statistical precision in the measurements from these surveys, the control of systematic errors from the combination of di ff erent cosmological probes has become key to achieving the expected accuracy.  \nWithin this framework, the European Space Agency approved the Euclid mission to carry out a comprehensive survey of most of the extragalactic sky from space. The Euclid mission is thoroughly described in Euclid Collaboration: Mellier et al. (2024). Euclid will use gravitational lensing and galaxy clustering as the main probes to study cosmology. It will carry out a wide and a deep survey of approximately 14 000 and 50 deg 2 , respectively. The wide survey (EWS) will reach a magnitude limit of I E ∼ 24 . 5, and the deep survey (EDS) will push approximately two magnitudes fainter (Euclid Collaboration: Scaramella et al. 2022). In its wide survey (EWS), Euclid will take images of billions of galaxies to determine their shapes and also obtain slitless spectra of tens of millions of galaxies to determine their redshifts, using its two main science instruments: the visible imaging instrument (VIS, Euclid Collaboration: Cropper et al. 2024) and the Near In- frared Spectrometer and Photometer (NISP, Euclid Collaboration: Jahnke et al. 2024). The combination of weak gravitational lensing and galaxy clustering will provide stringent  \n,  \ncosmological constraints (Euclid Collaboration: Blanchard et al. 2020).  \nCurrent and future cosmological surveys need simulations. In their definition stages, simulations are needed to define the survey requirements, to optimize its design and to plan the survey. Once a survey is running, simulations are needed to analyse the data and interpret the results. In the case of Euclid , there is a strong simulation e ff ort to prepare the science exploitation of the data. From a programmatic point of view, the e ff ort has focused on the simulations needed to support the mission reviews and explore the optimisation of the mission's scientific reach. For that purpose, Euclid has undertaken science performance verification exercises in which comprehensive analyses of the mission are performed to check the compliance with the scientific requirements. With the Euclid launch on 1 July 2023, the simulation focus has now turned to enabling the science exploitation of the first data releases.  \nThe optimal science exploitation of the new generation of galaxy surveys, such as Euclid , demands the development of large-volume and high-mass resolution numerical simulations that reproduce the large-scale galaxy distribution that these new surveys will observe with high fidelity. Not only do these help to assess the performance with a realism that cannot be achieved otherwise, but such simulations are also an essential tool for the development of the data processing and science analysis pipelines. Given the computational cost, so far, most synthetic galaxy catalogues have been developed out of N -body simulations where only gravity is used to follow the evolution of structure (e.g., Bertschinger & Gelb 1991; Couchman et al. 1995; Stadel 2001; Harnois-Déraps et al. 2013; Menon et al. 2015; Habib et al. 2016; Potter et al. 2017; Ishiyama et al. 2021; Garrison et al. 2021; Springel et al. 2021). For a recent review on cosmological N -body simulations, see Angulo & Hahn (2022). Galaxies are introduced in these simulations populating the dark matter haloes using di ff erent methods including semi-analytical models (SAM; e.g., White & Rees 1978; White & Frenk 1991; Kau ff mann et al. 1993, 1999; Somerville & Primack 1999; Benson et al. 2000; Cole et al. 2000; Hatton et al. 2003; Springel 2005; Hirschmann et al. 2016; Lagos et al. 2018; De Lucia et al. 2024), other empirical and phenomenological models like halo occupation distribution models (HOD; e.g., Cooray & Sheth 2002; Jing et al. 1998; Benson et al. 2000; Seljak 2000; Peacock & Smith 2000; Scoccimarro et al. 2001; Berlind & Weinberg 2002; Bullock et al. 2002), abundance matching (AM; e.g., Klypin et al. 1999; Kravtsov et al. 2004; Tasitsiomi et al. 2004), and sub-halo abundance matching (SHAM; e.g., Gu et al. 2024). With the improvement of computing capabilities, hydrodynamical simulations (for a recent review, see Vogelsberger et al. 2020) are now starting to be feasible for simulating cosmologically relevant volumes (Dolag et al. 2016; Pillepich et al. 2018; Schaye et al. 2023; see also the CAMELS project, the largest compilation of hydrodynamic simulations to date: Villaescusa-Navarro et al. 2023), and are sometimes used to train and inform phenomenological methods to populate N -body simulations.  \nThe production of simulated galaxy catalogues is a prolific line of development. Several cosmological surveys have produced simulations tailored to their observed samples but there are also more general-purpose simulated mocks. These galaxy catalogues include those based on the Uchuu simulation (e.g., Ereza et al. 2023; Dong-Páez et al. 2024) used for the SDSS data analysis, catalogues produced for the DESI survey (e.g., Smith et al. 2022; Balaguera-Antolínez et al. 2023; Smith et al. 2023), catalogues developed within the Rubin-LSST DESC collaboration (e.g., LSST Dark Energy Science Collaboration (LSST  \nDESC) et al. 2021; Kovacs et al. 2022), catalogues produced for the Chinese Space-station Survey Telescope (CSST, Gu et al. 2024), and general purpose galaxy catalogues (e.g., To et al. 2024; Behroozi et al. 2019).  \nWithin Euclid , we have developed the Flagship simulation, a comprehensive simulation, in terms of including a vast number of consistent galaxy properties, to help optimise the mission and prepare its scientific analysis and exploitation. The scientific goals of the mission from the main cosmological probes, weak gravitational lensing and galaxy clustering, set the requirements of the simulation in term of mass resolution, volume, and redshift coverage. Given that modelling Euclid with full hydrodynamic simulations over the volume surveyed by the mission is computationally not possible at the moment, the approach we followed to create mock surveys was to develop a state-of-the-art N -body simulation and populate the gravitationally bound dark matter structures (haloes) with galaxies in a way that best matches observational data, placing special care into simulating consistently the weak lensing and clustering properties to enable combined probes analyses.  \nThe first production was the Euclid Flagship 1 simulation (FS1 hereafter, Potter et al. 2017). The N -body simulation was run on the Piz Daint supercomputer at the Swiss National Supercomputing Centre in 2016. A lightcone of dark matter (DM) particles was generated on the fly, replicating the simulated box (with periodic boundary conditions) as a way to fill the full lightcone volume of the Euclid survey. In this scheme, we place the observer in one corner of the central box within the lightcone volume. The ROCKSTAR halo finder (Behroozi et al. 2013) was run on the DM particle distribution to generate a halo catalogue, which is now publicly available at the CosmoHub data distribution platform 1 (Carretero et al. 2017; Tallada et al. 2020). From the halo catalogue, we produced a galaxy catalogue that was used by the Euclid Collaboration to perform some early science analyses and performance assessments of the mission as a whole. To improve the scope of the Flagship simulation, a second version, called Flagship 2 (Flagship or abbreviated as FS2 hereafter) was run in 2020. There were several improvements with respect to the first version. The mass resolution and the maximum redshift covered by the lightcone output were increased in order to improve the completeness of the resulting catalogue to encompass all the galaxies expected to be detected by Euclid . In particular, the mass resolution increased by a factor of 2 in FS2, to reach a particle mass m p = 10 9 h -1 M ⊙ , which in turn allows for modelling galaxies about one magnitude fainter than with FS1 at all redshifts (see Sect. 5), and the lightcone was extended from z = 2 . 3 in FS1 to z = 3 in FS2. We also changed the way in which the spectral energy distributions were assigned in FS2 to make the resulting photometric properties closer to those observed. Similarly to the first version, a lightcone was produced on the fly (i.e., as the simulation run) and a halo catalogue was generated in post-processing with ROCKSTAR . Using HOD and AM techniques combined with relations between observational properties, we generated a galaxy catalogue containing around five billion objects covering one octant of the sky (see Sect. 5 for further details). Positions, velocities, physical properties, lensing quantities, and photometry in multiple bands were computed for all galaxies, totalling 399 parameters per galaxy generating a catalogue of 5.9 terabytes that can be accessed through the CosmoHub platform, 2 which is hosted by the Euclid mission Spanish Science Data Center. This catalogue has been the baseline input  \nsimulation for the Euclid mission pipelines and a key ingredient for its scientific preparation before the satellite launch. In this regard, additional galaxy mocks, which will not be discussed in this paper, have been constructed within the Euclid Consortium, to address more probe-specific scientific questions and account for the variance due to modelling uncertainties.  \nThis paper describes in detail the production of the second version of the Euclid Flagship galaxy catalogue (FS2). Upon publication of this paper, we expect to make a public release of the latest version of the catalogue (version 2.1.10), which will be available at the Cosmohub web portal. This publication will serve as a reference for its usage. Although designed for the Euclid mission, the catalogue can be very useful for many other studies and future galaxy surveys given its breadth in terms of number of galaxies simulated (e.g., 3.4 billion galaxies for a magnitude limited sample with H E < 26), volume covered (one octant of the sky up to z = 3), and the wide range of galaxy properties computed that, in particular, allow us to model the galaxy clustering of both photometric and spectroscopic galaxy samples along with their weak lensing observables consistently down to sub-arcminute scales. The paper is structured as follows. In Sect. 2, we describe the production and main characteristics of the N -body FS2 dark matter simulation. In Sect. 3, we present the computation of the lensing properties. In Sect. 4, we explain the production of the halo catalogue. In Sect. 5, we describe in detail the computation of the galaxy properties. The validation of the properties of the galaxy catalogue against observational constraints and theoretical models is presented in Sect. 6. Finally, we summarise our findings and present our conclusions in Sect. 7. Unless otherwise stated, all magnitudes reported in this paper are in the AB system.", '2.1. The Flagship run': "The Euclid Flagship N -body dark matter simulation features a simulation box of 3600 h -1 Mpc on a side with 16 000 3 particles, leading to a particle mass of m p = 10 9 h -1 M ⊙ . This four trillion particle simulation is the largest N -body simulation performed to date and matches the basic science requirements of the mission as it allows us to accurately resolve 10 11 h -1 M ⊙ haloes which host the faintest galaxies Euclid will observe (i.e., model a complete sample down to the Euclid flux limit) and samples a cosmological volume comparable to the one that the satellite will survey. The simulation was performed using the PKDGRAV3 code (Potter & Stadel 2016) on the Piz Daint supercomputer at the Swiss National Supercomputing Centre (CSCS). The simulation was run with a softening length of 4.5 h -1 kpc. It uses the Euclid reference cosmology, with the following values for the matter density Ω m = 0 . 319, baryon density Ω b = 0 . 049, dark energy density (in the form of a cosmological constant) ΩΛ = 0 . 681 -Ω r -Ω ν , with a radiation density Ω r = 0 . 00005509, and a contribution from massive neutrinos Ω ν = 0 . 00140343 which is derived from the minimum neutrino mass possible (0 . 0587 eV) given the measured mixing angles and assuming a normal hierarchy (de Salas et al. 2018). Besides, the values of the other cosmological parameters are: the equation of state parameter of dark energy w de = -1 . 0, the reduced Hubble parameter at redshift z = 0 (i.e., Hubble constant), h = 0 . 67, the scalar spectral index of the initial fluctuations n s = 0 . 96, and the primordial power spectrum amplitude A s = 2 . 1 × 10 -9 (i.e., σ 8 ≃ 0 . 813 derived) at k pivot = 0 . 05 Mpc -1 .  \nUsing the Euclid reference cosmology allows comparison to many other smaller simulations from N -body codes as well from approximate techniques that also use these reference values within the collaboration. The initial conditions were realised at z = 99 with first-order Lagrangian perturbation theory (1LPT) displacements from a uniform particle grid. The transfer functions for the density field and the velocity field were generated at this initial redshift by CLASS (Lesgourgues 2011) and CONCEPT (Dakin et al. 2022). As back-scaling was not used, all linear contributions from radiation, massive neutrinos, and metric perturbations (in the N -body gauge) were included via a lookup table and applied as a small corrective particle-mesh (PM) force at each time step (Fidler et al. 2019). This ensures a match to the linear evolution of the matter density field at all redshifts when including these additional linear terms (see Fig. 1).  \nThe main N -body data product was produced on-the-fly during the simulation and is a continuous full-sky particle lightcone (to z = 3), where each particle was output exactly when the shrinking light surface sweeps by it. The full-sky coverage was achieved by replicating the simulation box. The resulting catalogue contains 31 trillion particle positions and peculiar velocities (700 TB of data), and it was used to compute the roughly 125 billion ROCKSTAR main haloes and full-sky lensing maps with a HEALPix tessellation resolution N side = 8192 (Górski et al. 2005), corresponding to 0 ' . 43 per pixel. Note that the mock galaxy catalogue described below was computed only in one octant of the sky. The chosen mock area was primarily based on computational e ffi ciency, and it was deemed adequate to model the FS2 WIDE survey footprint for the first year data release (DR1), that will cover about 2500 deg 2 , for which the mock was mainly designed.", '2.2. PKDGRAV3 power spectrum validation': "Prior to performing the Flagship simulation, the PKDGRAV3 N -body code was validated against the well-established GADGET3 (Springel et al. 2008), GADGET4 (Springel et al. 2021), ABACUS (Garrison et al. 2021), and RAMSES (Teyssier 2002) codes, which each use very di ff erent methods to solve the Poisson equation as well as di ff erent methods to integrate the equations of motion. The results of these comparisons are given in Schneider et al. (2016) and more recently in Springel et al. (2021, Fig. 50) and Garrison et al. (2019, Fig. 5). All codes agree at the 1% level up to k = 10 h Mpc -1 . Convergence of the power spectrum as a function of the particle mass and simulation box size was also investigated. Conservatively, a particle mass of m p = 10 9 h -1 M ⊙ is required to ensure 1 % convergence of the power spectrum up to k = 10 h Mpc -1 . Simulation boxes larger than 1 h -1 Gpc are su ffi cient to ensure convergence in the power spectrum (e.g., Klypin & Prada 2018). However, the requirements of a light cone to z = 3, with as little replication of the volume as possible, lead to a box of 3600 h -1 Mpc on a side. This simulated box contained 4 trillion darkmatter particles, which yields the mass resolution desired and also corresponds to the limiting number allowed by the Piz Daint supercomputer given PKDGRAV3 's memory requirements (which are about 64 bytes / particle, including the tree structure and all bu ff ers for analysis and management of the Input / Output).  \nIn Fig. 1, we compare the power spectrum measured from the Euclid Flagship N -body dark matter simulation to linear theory computed by CLASS (Lesgourgues 2011) and to the Euclid Emulator nonlinear power spectrum (Euclid Collaboration: Knabenhans et al. 2021). Comparison to the BACCO emulator (Angulo  \nFig. 1. Nonlinear power spectrum at various redshifts compared to linear theory from CLASS (Lesgourgues 2011). The use of the forwardscaling method in the N -body gauge (Fidler et al. 2019) allows the N -body simulation to agree with linear theory at all redshifts including the e ff ects from weak field GR, radiation, and massive neutrinos. The narrow ± 5% band is plotted using a linear scale, making those fluctuations more visible and allowing also for negative values. Fluctuations at low k are statistical and reflect the sample variance of the particular realisation. At high k the enhanced power is due to nonlinear growth of structure not captured by the linear theory. Note that this nonlinear contribution also includes the slight -2% dip at k ≈ 10 -1 h Mpc -1 , corresponding precisely to the first BAO peak in the power spectrum. The forward-scaling method used for the Flagship simulation guarantees a match to the linear theory of CLASS at all redshifts, which was not possible to achieve with the traditional back-scaling technique (which only guarantees this at z = 0 in the presence of massive neutrinos and radiation). The bottom panel compares the power spectra to the Euclid Emulator nonlinear power spectrum (Euclid Collaboration: Knabenhans et al. 2021).  \n<!-- image -->  \net al. 2021) shows a very similar level of agreement. The 'spikes' at lower k are due to the cosmic variance present in the realisation of the Euclid Flagship N -body simulation. When comparing to other models of the nonlinear power spectrum, such as Halofit (Takahashi et al. 2012) and HMCode2020 (Mead et al. 2021), the comparison is not quite as good, with deviations over redshift at k > 0 . 1 h Mpc -1 extending to ± 5%, and notably, these models do not accurately capture the nonlinear form of the baryonic wiggles at k approximately 0 . 1-0 . 4 h Mpc -1 .", '2.3. Dark matter clustering': 'Galaxy clustering is one of the main probes of the Euclid mission. In order to validate this probe in the galaxy mock, we first  \nstudy whether the clustering of the dark matter distribution in the lightcone behaves according to theoretical expectations. Figure 2 shows the angular power spectrum of dark matter in thin redshift shells (in the HEALPix tessellation, see Sect. 3 for details), across the full depth of the Flagship N -body lightcone, i.e., 0 . 5 < z < 3. Measurements in the simulation make use of the PolSpice code (see Szapudi et al. 2001; Chon et al. 2004; Fosalba & Szapudi 2004) 3 which corrects for the e ff ect of angular masks in our finite-sky analysis.  \nResults show that measurements in the simulation agree with linear theory expectations on large scales (low multipoles) and nonlinear theory (Takahashi et al. 2012) down to very small scales (high multipoles, ℓ ∼ 10 4 ), within sample variance errors (see dotted envelopes in lower panels for each redshift bin). We note that particle shot-noise is negligible ( < 1% for all multipoles) given the high particle density, around 90 particles / ( h -1 Mpc) 3 , in the lightcone. The agreement between Flagship measurements and the Euclid Emulator2 (EE2) predictions is expected to be at a very similar level of agreement, as discussed in Euclid Collaboration: Knabenhans et al. (2021, see in particular their Figure 13) where they show that the EE2 and halofit agree within 3% up to very small (nonlinear) scales, k < 4 h Mpc -1 , for z < 2.', '3. HEALPix lensing mass maps': "Following the approach presented in Fosalba et al. (2008) and Fosalba et al. (2015), we construct a lightcone simulation by replicating the simulation box (and translating it) around the observer. Given the large box-size used for the Flagship simulation, L box = 3600 h -1 Mpc, this approach allows us to build all-sky lensing outputs without repetition up to z ∼ 2 . 0 and with one replication up to our maximum redshift, z max = 3. Then, we decompose the dark matter lightcone into a set of all-sky concentric spherical shells of given width, ∆ r , around the observer, what we call the 'onion universe'. Each dark matter 'onion shell' is then projected onto a 2D pixelised map using the HEALPix tessellation (Górski et al. 2005). For the lensing maps presented in this paper we have chosen a shell width of ∆ r ≈ 95 . 22 megayears in lookback time, up to z = 10 (and finer at higher redshifts), and an angular resolution of ∆ θ ≈ √ 3 / π 3600 / N side ≈ 0 . ' 43, for the HEALPix map resolution N side = 8192 that we use.  \nBy combining the dark matter 'onion shells' that make up the lightcone, we can easily derive lensing observables, as explained in Fosalba et al. (2008). This approach, based on approximating the observables by a discrete sum of 2D dark-matter density maps multiplied by the appropriate lensing weights, agrees with the much more complex and CPU time consuming ray-tracing technique in the Born approximation limit, i.e., in the limit where lensing deflections are calculated using unperturbed light paths (see e.g, Hilbert et al. 2020).  \nFollowing this technique we are able to produce all-sky maps of the convergence field (which is simply related to the lensing potential in harmonic space), as well as maps for other lensing fields obtained from covariant derivatives of the lensing potential, such as the deflection angle, shear, flexion, etc. Figure 3 shows the all-sky map of the convergence field, κ , for the Flagship simulation, for sources at zs = 1, with a pixel resolution of 0 . ' 43. The colour scale shown spans over the range -σ < κ < 3 σ , where σ is the root mean square (rms) fluctuation of the full-sky convergence map.  \n2  \n/  \nC  \n)  \n1  \n+  \n(  \n10  \n10  \n10  \n10  \n10  \n10  \n0  \n1  \n2  \n3  \n4  \n5  \n1.2  \n1.0  \n0.8  \nH  \nT  \n/  \nM  \nI  \nS  \n10  \n1  \nDM clustering  z = 0.5  \nFig. 2. Angular power spectrum of the dark matter field in the lightcone across di ff erent redshifts. Plots show the clustering in the following z-bins: 0 . 49 < z < 0 . 51 (top left), 0 . 99 < z < 1 . 00 (top right), 2 . 02 < z < 2 . 08 (bottom left) and 2 . 94 < z < 3 . 06 (bottom right). The clustering in the simulation (symbols) is compared against linear (dashed) and nonlinear (solid) theoretical predictions (Takahashi et al. 2012). Particle shot-noise is also shown for reference (dot-dashed line). Residuals with respect to nonlinear theory are displayed in the lower panels, along with sample variance error envelopes (dotted).  \n<!-- image -->  \nDM clustering  z = 2.0  \nSIM  \nshotnoise  \nLinearTH  \nNon  \nLin TH  \n10  \n2  \nMultipole", 'DM clustering  z = 1.0': 'SIM  \nshotnoise  \nLinearTH  \nNon  \nLin TH  \n10  \n2  \n10  \n3  \nMultipole  \nDM clustering  z = 3.0  \nSIM  \nshotnoise  \nLinearTH  \nNon  \nLin TH  \n10  \n2  \nMultipole  \nThe angular power spectrum of the convergence field in the Born approximation reads (for a flat LCDM cosmology),  \nC κ ℓ = 9 H 4 0 Ω 2 m 4 c 4 Z d r P ( ℓ/ r , r ) ( r s -r ) 2 r 2 s a 2 , (1)  \nwhere ℓ is the multipole order, H 0 = 100 h kms -1 Mpc -1 is the Hubble constant, c is the speed of light, and r s is the comoving distance to the lensing sources (we assume all sources are located at the same redshift in this approximation) where P ( ℓ/ r , r ) is the 3D density power spectrum in the simulation at a given comoving distance r from the observer.  \nIn this approach, we can take the spherical transform of the lensing potential all-sky map to obtain the corresponding maps for the other weak-lensing observables through simple relations  \nin harmonic space (see Hu 2000). In particular, the convergence field, κ , is related to the lensing potential, ϕ , through the 2D equivalent to the usual (3D) Poisson equation, which in spherical harmonic decomposition reads  \nκℓ m = -1 2 ℓ ( ℓ + 1) ϕℓ m . (2)  \nOne can thus use this expression to derive the lensing potential at each source plane (or 2D lightcone map), and obtain other lensing observables, such as deflection and shear, through their relation to the lensing potential in harmonic space (see Fosalba et al. 2015, for details). As a basic validation of the mass maps, Fig. 4 shows the measurement of the convergence angular power spectrum in the simulation compared with theory predictions, for two di ff erent source redshifts across the lightcone. Overall there  \n10  \n3  \n10  \n4  \n2  \n/  \nC  \n)  \n1  \n+  \n(  \n2  \n/  \nC  \n)  \n1  \n+  \n(  \n10  \n10  \n10  \n10  \n10  \n10  \n0  \n1  \n2  \n3  \n4  \n5  \n1.2  \n1.0  \n0.8  \n10  \n10  \n10  \n10  \n10  \n10  \n0  \n1  \n2  \n3  \n4  \n5  \n1.2  \n1.0  \n0.8  \nH  \nT  \n/  \nM  \nI  \nS  \nH  \nT  \n/  \nM  \nI  \nS  \n10  \n1  \n10  \n1  \n10  \n3  \n10  \n4  \n10  \n4  \n<!-- image -->  \nFig. 3. Lensing convergence (colour coded) and shear field (sticks) for sources at z s = 1 over a patch of 50 deg 2 (left) and in a zoom-in of the central 1 deg 2 area (right). The convergence field colour bar displays values within the range ± 3 σ , where σ is the rms value of the full-sky map. The stick at the bottom of the zoom-in image shows a reference amplitude for the shear sticks overlaid on that area of the mass map.  \n<!-- image -->  \nis good agreement between the mass map clustering compared to theory in the full range of scales (multipoles) shown, given the sample variance errors (see figure caption for details).', '4. Halo catalogue': "The dark matter haloes were identified directly on the lightcone particle data using the ROCKSTAR halo finder (Behroozi et al. 2013). ROCKSTAR is a phase space-linking friends-of-friends method that is able to find the hierarchy of substructure from parent dark matter haloes to the smallest subhaloes. ROCKSTAR is also a high-performance parallel halo finder; however, it is not capable of handling such a massive (10 trillion particle) simulation in its standard (public) version. In order to use it for finding haloes in the FS2 particle lightcone data, we had to split the data into computational 'bricks' of 375 h -1 Mpc × 375 h -1 Mpc × 625 h -1 Mpc, each with an extended 'ghost' region of 5 h -1 Mpc on each side to avoid discontinuities. The full particle lightcone comprises 3448 such computational bricks, each of which could be computed independently on a cloud of (56 core) servers at the University of Zurich. One complication in the processing is that the data in the lightcone is over a variable expansion factor, a , as a function of depth in the lightcone, changing from 1.0 at the centre to 0.25 ( z = 3 . 0) at the edge. This fact must be accounted for when converting the particle momenta to physical peculiar velocities. 4 ROCKSTAR was written to compute halo catalogues from a set of simulation snapshots, each at a fixed expansion factor, and not a lightcone and, therefore, needed to be modified to handle particles with radially dependent expansion factor.  \nROCKSTAR uses these peculiar velocities for both linking (where the linking length in velocity space is adapted from halo to halo) as well as for 'unbinding', the process of removing particles from a halo that are deemed not to be gravitationally bound to it (in isolation). Once all haloes within a brick have been found, the parent haloes with centres in the ghost region and subhaloes of such parent haloes are removed from the catalogue so that the individual bricks fit seamlessly together.  \nThe science reach of Euclid depends on how much volume it can sample and how many tracers it can use for cosmological analysis. Its design was optimised to obtain the most stringent constraints on cosmological parameters. For its weak lensing analysis, it reaches magnitude limits of I E ≃ 24 . 5. The default minimum number of particles for ROCKSTAR to define a halo is set to 20. However, we set the minimum number of particles to define a halo to 10 to be complete at the Euclid weak lensing magnitude limit (see Sect. 4.2.1). As we use such small and poorly resolved haloes, we correct the halo masses of haloes with few particles to avoid incompleteness and discreteness e ff ects in the halo mass function (Sect. 4.2). The final all-sky halo catalogue contains 126 billion main haloes. The galaxy catalogue is generated from the halo catalogue in just one octant of the lightcone that contains 15.8 billion main haloes.", '4.1. Halo mass function': 'The ROCKSTAR halo finder produces di ff erent estimates of the halo mass. These estimates include: the mass, M fof, of the particles linked together with a friends-of-friends algorithm of linking length b = 0 . 2; the mass, M vir, contained within the virial radius; the sum of the mass, M bound, of the bound particles within  \nLensing Convergence z = 1.0  \nFigure 6 shows the ratio of the cumulative M bound HMF to the T08 M vir HMF at several redshifts spanning the redshift range of the simulation lightcone. While the slopes of the HMFs in the mass range 10 11 . 5 h -1 M ⊙ ≲ M bound ≲ 10 13 . 5 h -1 M ⊙ are similar at low redshift, z ≲ 1, at higher redshift, the slope of the M bound HMF progressively gets shallower than the T08 M vir HMF. The ratio of the abundance at a given M bound halo mass compared to T08 abundance at the same halo mass also increases with redshift. Part of the di ff erence may be due to the di ff erent power spectrum transfer function used in the Flagship run compared to the input we have given to the hmf code to compute the T08 HMF, generated with CAMB (Lewis et al. 2000) for the same cosmology. We have performed the same comparison to the D16 and C17 HMFs (not shown in the figure) finding qualitatively the same result. Nevertheless, the photometric quantities we compute for the mock galaxies depend on an abundance matching technique and therefore are not a ff ected by small changes in the HMF. The observed luminosity function is recovered for the galaxies by construction in the AM technique despite any mismatch or incompleteness observed in the mass function of the dark matter haloes.  \n<!-- image -->  \nFig. 4. Angular power spectrum of the convergence field at z s = 1 (top) and z s = 3 (bottom). Plots show measurements in the simulation (symbols) compared against linear theory (dashed line) and nonlinear (solid line) fits to simulations (Takahashi et al. 2012). The lower panels show the ratio between the simulation and nonlinear theory. Sample variance error envelopes are displayed as dotted lines.  \n<!-- image -->  \nthe virial radius; the mass, M 200b, of the particles within an overdensity of 200 relative to the background density; and the mass, M 200c, of the particles within an overdensity of 200 relative to the critical density. Appendix A provides a comparison of the halo mass function (HMF) for the di ff erent halo mass definitions.  \nBased on the similar behaviour of the HMF with the di ff erent mass estimates (except for the friends-of-friends mass definition), we decided to choose the M bound definition as a default choice to build the galaxy catalogue. As our method of assigning galaxy luminosities is based on AM, the particular choice of mass estimate is not important.  \nWecompare the M bound HMF to other HMFs in the literature. We use as main reference the Tinker et al. (2008) HMF, here-  \nFig. 5. Cumulative M bound HMF of the Flagship haloes compared to the T08, D16 and C17 halo mass functions at redshift z = 0 . 1. In the top panel we show the cumulative HMFs (colours indicated in the legend). The lower panel shows the ratio of the cumulative M bound HMF to the other cumulative HMFs. The dotted line serves as a reference point when cumulative HMFs are the same.  \n<!-- image -->  \nafter T08, as it has been widely used in the literature for HMF comparison. We also compare it to the HMFs of Despali et al. (2016), D16, and of Comparat et al. (2017), C17. We use the hmf 5 (Murray et al. 2013) and COLOSSUS 6 (Diemer 2018) packages to compute the HMFs using the reference Flagship cosmological parameters. Figure 5 shows the cumulative HMFs at low redshift, z = 0 . 1 in the top panel, and the ratio of the M bound HMF to the other HMFs in the lower panel. We compute the T08, D16 and C17 HMFs for the M vir mass definition. Our M bound is almost the same as M vir for the most massive haloes but di ff ers for the lowest mass haloes due to the unbinding of particles. There is an overall o ff set of around 3-7% lower abundance in the M bound HMF compared to the other predictions for the same cosmology. For halo masses below M bound = 10 11 . 5 h -1 M ⊙ , equivalent to ∼ 300 particles, the M bound HMF starts to be incomplete. The di ff erential HMF shows the same trends as the cumulative HMF. The lightcone has little volume at z = 0 . 1 and therefore the HMF is very noisy above M bound = 10 13 . 5 h -1 M ⊙ .  \nFig. 6. Natural logarithm of the ratio of the cumulative M bound HMF to the T08 M vir HMF at several redshifts indicated in the figure legend. The normalisation between the two varies mostly within a 10% di ff erence range. The slope of the Flagship cumulative M bound HMF is shallower than the one of the T08 M vir HMF at high redshift.  \n<!-- image -->', '4.2. Mass corrections': 'As we will see later, we assign galaxy luminosities to central galaxies using a halo mass-luminosity relation derived from abundance matching between the cumulative HMF and the cumulative galaxy luminosity function. We push the detection of dark matter haloes to the very low limit of ten particles, making the e ff ects of discreteness very noticeable at the low mass end of the HMF. Furthermore, as mentioned above, below halo masses of log10[ M / ( h -1 M ⊙ )] ≲ 11 . 5 our HMF starts to be incomplete. In order to produce galaxy luminosities that are not discrete and incomplete, we need to correct the HMF for these two e ff ects.', '4.2.1. Incompleteness correction': 'In order to reach the faint absolute magnitudes that Euclid will observe, we need to detect haloes down to the corresponding low masses. For a Euclid magnitude limit of I E = 24 . 5, we need to reach absolute magnitudes around MI E -5 log10 h ≃ -13 . 0 to be complete at redshifts z ≳ 0 . 1. To reach this absolute magnitude limit, we need to reach a mass limit of M ≃ 2 × 10 10 h -1 M ⊙ . Given the resolution of the simulation, this mass corresponds to 20 particles. As we re-scale the halo masses to account for the HMF incompleteness, we need to push down to haloes identified with at least 10 particles. With the rejection of unbound particles, the M bound halo definition can have even fewer particles contributing to the halo mass. At this particle mass threshold, the halo catalogue is not complete. Nevertheless, as the two-point correlation of haloes is approximately independent of mass at low masses, the two-point correlation properties of all haloes detected will not di ff er from the one it would have had if the catalogue had been complete. That way, we can reassign the halo masses with abundance matching and assume that we are complete to the halo abundance given by the lowest number of particles and that the two-point clustering, which we use to calibrate the galaxy mock, will not change.  \nWe correct for incompleteness by reassigning the halo masses in the following way. We assume that the slope of the cumulative HMF at low masses is the same as the T08 cumulative HMF. Given that the Flagship halo abundance for the M bound definition is somewhat lower than the T08, we adjust the abundance at a mass log10[ M / ( h -1 M ⊙ )] = 11 . 5, which corresponds  \nFig. 7. Cumulative HMF at redshift z ∼ 0 . 1 for the M bound definition (blue), the T08 M vir HMF (black) and the one resulting from the mass reassignment procedure (red) that corrects for incompleteness and discreteness e ff ects.  \n<!-- image -->  \nto approximately 300 particles per halo. Above this mass threshold, there seems to be no incompleteness due to the low number of particles (see Fig. 7). We reassign the halo masses below this threshold to have the same abundance that a fiducial HMF constructed with the same faint-end slope of the T08 cumulative HMF and normalised to the Flagship cumulative HMF at log10[ M / ( h -1 M ⊙ )] = 11 . 5. The process is captured in Fig. 7 where the original cumulative HMF at redshift z ∼ 0 . 1 is shown in blue, the T08 HMF is shown in orange, and the resulting cumulative halo mass after the abundance matching procedure is shown in red. The new cumulative HMF has the same faint-end slope as the T08 HMF and the normalisation of the original HMF by construction.  \nWhile this procedure is conceptually simple, implementing it directly into the mock generation is too slow, as one needs to compute the observed cumulative HMF and to invert the T08 cumulative HMF for each galaxy. We therefore developed a faster way of implementing this correction. First, we compute the relation between the original M bound halo mass and the abundancematched halo mass. In order to be able to compute this correction e ffi ciently, we fit this relation with five parameters at all redshifts in intervals of 0.1 in redshift. We then fit each parameter as a function of redshift. We also fit the o ff sets of the cumulative HMF to the T08 values at log10 [ M / ( h -1 M ⊙ )] = 11 . 5.', '4.2.2. Discreteness correction': 'We correct for discreteness by assuming that the cumulative HMF is well defined at the mass values corresponding to a given integer number of particles. We proceed as follows. In a given volume V for each number of particles pi , we have Ni haloes with those pi particles corresponding to a halo mass Mi = pi m p. The cumulative abundance for this halo mass is ni ( ⩾ Mi ) = P ∞ i = 1 Ni / V , that is, the number of haloes with mass larger or equal than Mi divided by the volume. For all haloes with pi particles, we want to reassign their masses to distribute them according to a power-law distribution in abundance between masses Mi and Mi + 1. For low masses, it is a good approximation to assume that the cumulative HMF behaves as a power law,  \nCHMF( ⩾ Mi ) ≃ 10 α M β i , (3)  \nin a small range of masses. For each halo with pi particles and mass Mi , we draw a random number u uniformly distributed in the range [0 , 1] and use this number to obtain a halo mass M from the normalised cumulative halo mass distribution between Mi and Mi + 1,  \nu = CHMF( ⩾ M ) -CHMF( ⩾ Mi + 1) CHMF( ⩾ Mi ) -CHMF( ⩾ Mi + 1) , (4)  \nassigning a halo mass  \nM = h u ( M β i -M β i + 1 ) + M β i + 1 i 1 β . (5)  \nSo, in order to make the halo mass distribution continuous, we only need to know the slope β of the cumulative HMF as a function of mass and also as a function of redshift (as it evolves). We produce functional fits to the value of the slope to be able to have the slope at any value of the mass and redshift faster than by computing the cumulative HMF in every step.', '4.3. Halo clustering': 'Figure 8 shows the angular clustering of haloes in harmonic space at di ff erent redshifts ( z = 0 . 5 , 1 , 1 . 5 from left to right) and for two mass bins (10 13 h -1 M ⊙ > M h > 10 12 h -1 M ⊙ , and 10 12 h -1 M ⊙ > M h > 10 11 h -1 M ⊙ , top and bottom rows, respectively). The theory predictions correspond to the dark matter clustering globally re-scaled with an estimate of the linear halo bias, fitted to match the measured clustering on the lowest multipoles. Clustering corrected with a simple model for the shotnoise is also shown. The resulting (corrected) clustering is in good agreement with linear theory expectations at low multipoles, i.e, ℓ < 200 at z = 0 . 5, and ℓ < 500 at z = 1 . 5, within sample variance errors (see lower panels). These limiting multipoles are only approximate given that the scales beyond which linear matter growth and the linear halo bias model break down depend both on redshift and halo mass. We also note we have assumed that haloes are Poisson distributed which is not strictly the case (Ginzburg et al. 2017). However a proper correction of the shot-noise for halo clustering is beyond the scope of this paper.', '5. From haloes to galaxies': 'Galaxies were generated in the Flagship simulation using a prescription that includes HOD and AM techniques and observed relations between galaxy properties. The prescription follows the methodology used in populating the MICE Grand Challenge simulation with galaxies (Carretero et al. 2015, C15 hereafter). The starting point is the Flagship halo catalogue described above. We use the SciPIC algorithm, described in Appendix B, to compute the galaxy properties. We run the pipeline at the Spanish Euclid Science Data Center. In the following subsections, we describe the di ff erent steps of the galaxy catalogue production.', '5.1. Galaxy luminosities': 'Following the HOD philosophy, haloes are populated with central and satellite galaxies. Our prescription starts with a hybrid HOD and AM method, that computes the number of satellites in each halo and assigns the galaxy luminosities. Galaxy clustering measurements are used to determine the parameters and the relations implemented.  \nThe method has the following steps. First, it assumes that each halo is composed of one central galaxy and some satellites. We use a simple HOD in which the mean number of satellites (the satellite occupation) depends only on the halo mass as a power law. For all haloes with masses larger than the minimum halo mass ( M h ≥ M min), the number of central galaxies and the mean number of satellite galaxies are given by  \nN cen = 1 , ⟨ N sat ⟩ = M M 1 ! α . (6)  \nWe assign the number of satellite galaxies in each halo drawing a realisation of a Poisson distribution with mean ⟨ N sat ⟩ . We parameterise M 1 as a multiplicative factor, f times M min:  \nM 1 = f M min . (7)  \nIn MICE, we calibrated the factor f as a function of halo mass to match the SDSS two-dimensional projected clustering constraints of Zehavi et al. (2011) at low redshift. In Flagship, we adopt a constant value for the multiplicative factor in Eq. (7) and fix it to f = 15, which is approximately the mean value we used in MICE (Carretero et al. 2015). We further choose the multiplicative factor not to depend on redshift given the weak constraints on galaxy clustering at high redshift. We also set the exponent in Eq. (6) to a fixed value, α = 1, with no redshift dependence either. We will show in Sect. 6 that our galaxy mock gives results that are in good agreement with a large set of observational data.  \nWe then use AM to assign the galaxy luminosities. In order to obtain a relation between the halo mass and the central galaxy luminosity, we first calculate the cumulative density of galaxies (central and satellites), n gal, as a function of the halo mass. 7 We compute this function by integrating the CHMF, taking into account the HOD assignment,  \nn gal( > M min) = Z ∞ M min n ( M \' ) " 1 + M \' f M min ! α # d M \' , (8)  \nwhere we have used Eqs. (6) and (7) to compute the number of galaxies per halo at the mass threshold, M min. Note that both n gal and n are densities, that is, number of galaxies per unit volume. We refer to this function defined in Eq. (8) as the cumulative galaxy function (CGF).  \nThe adopted galaxy luminosity function (LF) for the AM is a variation of the function given by Dahlen et al. (2005), which is based on multi-band data taken in the Great Observatories Origins Deep Survey (GOODS; Giavalisco et al. 2004). The GOODSLFis parameterised in several optical and near-infrared bands up to redshift z = 2. We extrapolate it to higher redshifts, z = 3, and transform it to our reference band. The total luminosity function is the sum of the LFs of three populations, each with  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFig. 8. Angular power spectra of the dark matter haloes in the lightcone. Top panels show clustering for haloes of mass 10 13 h -1 M ⊙ > M h > 10 12 h -1 M ⊙ , whereas bottom panels display the corresponding results for 10 12 h -1 M ⊙ > M h > 10 11 h -1 M ⊙ . Lower panels show ratios with respect to nonlinear theory, with dotted lines displaying the envelope for sample variance errors.  \n<!-- image -->  \nFig. 9. Cumulative luminosity function of the galaxy population in the r 01 band at di ff erent redshifts used in our abundance matching procedure.  \n<!-- image -->  \nits own evolution. As most of our calibration is inherited from the MICE catalogue which was performed at low redshift using SDSS data, and in particular the New York University Value Added Galaxy Catalogue (NYU-VAGC; Blanton et al. 2005b), we choose as our reference filter the SDSS r 01 band, 8 which is the SDSS r filter redshifted to z = 0 . 1 (Blanton et al. 2003). Figure 9 shows the cumulative LF function in the r 01 band for several redshifts.  \nSeveral studies have shown that when generating a galaxy mock catalogue with AM using the observed cumulative LF, the resulting galaxy clustering amplitude for the most luminous galaxies is higher than observed. The clustering amplitude for these luminous galaxies can be lowered if scatter is applied in the luminosity assignment (e.g. Tasitsiomi et al. 2004; More et al. 2009; Behroozi et al. 2010; Reddick et al. 2013; Carretero et al. 2015). Given that the highest luminosity range is dominated by an exponential decay, the introduced scatter results in assigning higher central galaxy luminosities to haloes of lower masses, thus reducing the amplitude of the clustering at the high luminosity range. Nevertheless, at lower luminosities where the LF is mainly dominated by a power law, the inclusion of this scatter has no e ff ect in the overall shape of this function. Besides obtaining a more realistic clustering for the most luminous galaxies, introducing this scatter also reflects the fact that galaxy formation is a stochastic process. Taking this into account, we apply a scatter to the galaxy luminosities resulting from the AM step. We define an unscattered LF, Φ unscat( L ), that when convolved with a Gaussian scatter in the logarithm of the luminosity (function G in Eq. 9) gives the observed LF, Φ obs( L ),  \nΦ obs( L ) = Z ∞ -∞ Φ unscat( L \' ) G ( L -L \' ) d L \' , (9)  \nwhere the Gaussian function G has a mean of zero and a standard deviation of ∆ log10[ L / ( h -2 L ⊙ )]. To obtain the unscattered LF, we need to solve for Φ unscat( L ) in the convolution equation (9). Instead, to gain computing e ffi ciency, we approximate the e ff ect  \nFig. 10. Relation between halo mass and luminosity at various redshifts resulting from abundance matching the cumulative galaxy function and the cumulative unscattered luminosity function following Eq. (12).  \n<!-- image -->  \nof the convolution with an exponential decay factor:  \nΦ unscat( L , z , ∆ log10 L ) = Φ obs( L , z ) exp        -L L sct( z , ∆ log10 L ) ! β sct( z , ∆ log10 L )        , (10)  \nwhere we fit the parameters L sct and β sct as a function of redshift, z , and the value of the scatter, ∆ log10[ L / ( h -2 L ⊙ )]. We compute the cumulative unscattered LF as  \nn gal( > L ) = Z ∞ L Φ unscat( L \' ) d L \' , (11)  \nwhere we have omitted the dependence on redshift and the luminosity scatter. We can establish a relation between halo mass and luminosity by applying AM between the cumulative number of galaxies function (Eq. 8) and the cumulative unscattered LF (Eq. 11):  \nn gal( > M halo , z ) = n gal( > L , z , ∆ log10[ L / ( h -2 L ⊙ )]) . (12)  \nFigure 10 shows the resulting AM relations for a few redshift values from Eq. (12) for our assumed value of ∆ log10[ L / ( h -2 L ⊙ )] = 0 . 12. The relation is fitted as a double power law with parameters fitted as a function of halo mass and redshift. The luminosities of the central galaxies are assigned from their halo mass and their redshift using the relation coming from Eq. (12) and shown in Fig. 10 for a few redshift values. We then add a random realisation of the scatter in the logarithm of the luminosity to its value assuming a Gaussian distribution. We take the standard deviation of this distribution to be ∆ log10[ L / ( h -2 L ⊙ )] = 0 . 12, which is similar to the value applied in other studies and produces a consistent clustering signal in our catalogue at high luminosities and at low redshift.  \nIn order to assign luminosities to the satellites, we first compute the global LF for these galaxies by subtracting the LF for the centrals from the total LF used to compute the AM. Figure 11 provides a visual representation of this step. It shows the model LF used to compute the AM relation between mass and luminosity and the measured LF in the Flagship catalogue in a thin redshift slice centred at z = 0 . 5 for all the galaxies and for centrals and satellites only. Once we have the global LF for the satellites (orange line in Fig. 11), we assume that it is the result of the sum of the individual LFs of satellites within each halo for  \nFig. 11. Di ff erential luminosity function for central and satellite galaxies and all galaxies at redshift z = 0 . 5 in the Flagship catalogue. We also show the model LF as in Fig.9 but in its di ff erential form.  \n<!-- image -->  \nall haloes. We model the cumulative LF of the satellites within each halo using a modified four-parameter Schechter function of the form  \nN halo sat ( > L ) = A L L ∗ ( L cen) ! α exp       -L L ∗ ( L cen) ! β       , (13)  \nwhere A is the normalisation factor that ensures that the halo contains the number of satellites predicted by the realisation of the HOD with luminosities higher than the minimum luminosity assigned for that halo, L min, α is the faint slope that we fix to α = -0 . 5, and β is a parameter that controls the steepness of the bright-end exponential cut-o ff . Furthermore, L ∗ is the characteristic luminosity that depends linearly on the luminosity of the central galaxy, L ∗ = a L cen + b . Taking into account this model for each halo, we fit the three parameters, a , b , and β , to make the sum of the satellite LF of all the haloes match the global LF. We obtain the values a = 10 -0 . 5 , b = 0 and β = 1 . 5 from this process. Finally, we assign the luminosities for the satellites in each halo by randomly drawing the luminosities following Eq. (13).', '5.2. Galaxy positions and colours': "Before allocating positions and velocities to the galaxies, we assign them a ( g 01 -r 01)HOD 9 colour, defined as the colour computed using the g and r SDSS filters redshifted to z = 0 . 1 (Blanton et al. 2003). This colour will be used together with the luminosity and redshift to assign the spectral energy distribution to the galaxies (see Sect. 5.4). We define three colour populations and assign individual colours to galaxies depending on which population they belong to. The abundance of these colour populations in each halo is defined by the HOD prescription. We follow the same procedure as C15 to assign the ( g 01 -r 01)HOD colour that depends on the galaxy type. We have extended the C15 fits to fainter luminosities as the Flagship catalogue contains haloes of lower mass and thus reaches fainter luminosities. A full description of the method is given in C15. We summarise it here, noting the extensions and improvements made.  \nWe start from the colour-magnitude diagram (CMD) using the low redshift NYU-VAGC SDSS galaxy catalogue (Blanton et al. 2005a), which we show in Fig. 12. We fit the g 01 -r 01  \nFig. 12. Colour-magnitude diagram of the SDSS NYU-VAGC. The black lines show contours of the same density of galaxies in this space. The red, green, and blue solid lines show the mean of the distributions of the three Gaussian fits to the ( g 01 -r 01) colour distributions as a function of magnitude. The dotted lines show the position of the mean ± the standard deviation of the Gaussian distributions.  \n<!-- image -->  \ncolour distribution as a function of luminosity in this colourmagnitude diagram to three galaxy populations: blue, green, and red. Similarly to what was done in Baldry et al. (2004), Skibba & Sheth (2009), and Carretero et al. (2015), we assume that the three populations are characterised by Gaussian distributions. We fit their means and standard deviations as a function of luminosity. Figure 12 includes these fits as solid and dotted lines, respectively, in the colour-magnitude diagram. The red, green, and blue lines correspond to the red, green, and blue populations, respectively. Note that at high luminosities, there appears to be a small extra component redder than the main Gaussiandistributed red population. For computational convenience, we absorb this extra red population within the green population, which has no galaxies at high luminosities. This is the same procedure that was used for the MICE catalogue (Carretero et al. 2015). As an example of these fits, Fig. 13 presents these three colour distributions at the absolute magnitude Mr 01 -5 log10 h = -20 . 0. Note that the mean and standard deviations of the three populations are fit by smooth global functions that depend on the absolute magnitude and not at each value of the absolute magnitude separately.  \nWe adopt the procedure of Skibba & Sheth (2009) and followed in C15 to compute the ( g 01 -r 01)HOD colour based on the galaxy type (central or satellite). We also di ff erentiate the galaxies into three colour types defined by the Gaussian fits of the CMD. We assign the parameter color\\_kind to each galaxy in the catalogue to indicate to which colour population it belongs to: red, green, or blue. We define a function to determine the fraction of satellites that belong to the red population as a function of luminosity and another one for the green population. The HOD and the colour-magnitude diagram then determine the fraction of satellites that belong to the blue population and the fraction of centrals that belong to the three populations (see Eqs. 37-43 in C15). These functions were optimised in MICE to reproduce the clustering as a function of colour and luminosity of the SDSS sample (Zehavi et al. 2011). For Flagship we modified these functions slightly and extended them to lower luminosities. However, we did not run a proper minimisation exercise, in part due to the flexibility in the functional form that we have allowed. Figure 14 shows these functions as implemented in the Flagship catalogue. All the brightest galaxies are centrals. Then, there is a  \nFig. 13. Colour distribution of galaxies in the SDSS NYU-VAGC at absolute magnitude Mr 01 -5 log10 h = -20 . 0. The light blue histograms represents the data in the catalogue. The red, green, and blue lines are the fits to the three Gaussian populations and the black line the sum of the three.  \n<!-- image -->  \nFig. 14. Galaxy fractions as a function of absolute magnitude. The black dashed line represents the fraction of centrals, while the rest of the lines present combinations of galaxy type and galaxy colour type. The line style distinguishes centrals (dashed lines) from satellites (dotted lines) and all galaxies (solid lines). The galaxy colour types are presented with di ff erent colours, red, green, and blue.  \n<!-- image -->  \ntransition around the LF characteristic luminosity, L ⋆ , where the fraction diminishes to become approximately constant around 50% for lower luminosities. The overall fraction of red, green, and blue galaxies is constrained by the CMD, with many more red galaxies at bright luminosities and blue galaxies dominating the population at faint luminosities. At all luminosities, there is a higher fraction of red galaxies that are satellites than centrals, and the other way around for blue galaxies. The fraction of green galaxies is always low compared to the red and blue ones. Note that brighter than absolute magnitudes Mr 01 -5 log10 h ≃ -21 . 5, the galaxies tagged as green are indeed redder than the main red population. The fraction of green satellites is larger than that of green centrals at low luminosities, but the trend is reversed at bright luminosities, where all this extra population of the reddest galaxies (tagged as green) are centrals.  \nFig 15 displays the fraction of color\\_kind values as a function of luminosity at four redshifts, splitting between centrals and satellites. There are no more local minima or maxima in the fractions vs. luminosity as were seen in Fig. 14, because of the  \nFig. 15. Fraction of galaxies of di ff erent types as a function of absolute magnitude for di ff erent redshift slices ( ∆ z = ± 0 . 1), split between Red Sequence ('RS'), Green Valley ('GV'), and Blue Cloud ('BC'), and between centrals ('cens') and satellites ('sats'). The abscissae have been slightly shifted for clarity.  \n<!-- image -->  \nsplit in redshift bins. Flagship predicts that 1) luminous galaxies are redder, consistent with what is observed at low redshifts (Sandage & Visvanathan 1978); 2) satellites are more likely to be red than centrals of the same luminosity, as observed at low redshifts (van den Bosch et al. 2008); 3) centrals and satellites (in particular low-luminosity ones) are more likely to be blue at higher redshifts, in conformity with the Butcher-Oemler e ff ect (Butcher & Oemler 1978) that galaxies in higher-redshift clusters contain a higher fraction of spiral morphologies (hence bluer colours).  \nWe place the central galaxies at the centre of their haloes. Satellite galaxies are located following a triaxial Navarro, Frenk, and White (NFW; Navarro et al. 1997) profile. We use the virial radius, the concentration parameter, and the vectors of the ellipsoid semi-axes to compute the satellite positions. We obtain the concentration parameter of the haloes using the virial radius and the scale radius from the ROCKSTAR catalogue 10 for haloes with masses log10[ M halo / ( h -1 M ⊙ )] > 11 . 0, where there are enough particles per halo as to reliably determine the concentration. For this mass range, the concentration parameter agrees with a mass dependence that matches that given by Diemer & Joyce (2019) as computed with the COLOSSUS code. Below this mass thresh-  \nFig. 16. Concentration index as a function of halo mass at the redshifts marked in each panel. The blue dots correspond to the means of the distributions at each halo mass bin and are linked with a blue line. The vertical blue lines show the standard deviations of the distributions. The red line shows the Diemer & Joyce (2019) functional form values.  \n<!-- image -->  \nd, we use the relation from Diemer & Joyce (2019) to obtain the mean concentration for a given halo mass and redshift. We have computed the ratio of the standard deviation to the mean of the distribution of concentration parameters as a function of halo mass and redshift. This ratio is approximately constant with a value of 1 / 3. We therefore assign the concentration parameter for the haloes below the mass threshold M halo = 10 11 h -1 M ⊙ , drawing a realisation of a Gaussian distribution with the mean coming from the relation of Diemer & Joyce (2019) and a standard deviation that is one-third of the value of the concentration mean. Figure 16 shows with blue dots the mean of the concentration parameter as a function of halo mass at three di ff erent redshifts before correcting the values at low halo mass. The blue lines indicate the standard deviation of the concentration parameters. The red line shows the Diemer & Joyce (2019) relation. The recomputed mean concentration values at low halo mass coincide with the red line.  \nThe other ingredient for assigning positions to the satellites within their haloes is the triaxial shape of the haloes. ROCKSTAR outputs the position vectors of two of the ellipsoid semi-axes and the axis ratios. We compute the third position vector as the cross product of the other two. We take the virial radius in the semimajor axis direction as the virial radius provided by the catalogue and in the other directions we multiply the virial radius by the corresponding axis ratios. We assign the position of satellite galaxies within haloes following the procedure of Robotham  \n& Howlett (2018), 11 which we modify to compute random positions in a triaxial NFW keeping the appropriate density profile, axis ratio values, and ellipsoid orientation. We compute the positions in each halo coordinate system given by the position vectors of their semi-axes and then transform them to comoving coordinates taking into account the orientation of the ellipsoid with respect to the observer's lightcone.  \nWe implement colour segregation in the satellite galaxy distribution assuming a di ff erent concentration parameter for each of the color\\_kind samples (e.g., McDonough & Brainerd 2022). Red galaxies have the same concentrations of the dark matter in the haloes. Blue galaxies have a concentration one fourth of the concentration of the red galaxies (following Collister & Lahav 2005), while green galaxies have a concentration half of that of the red galaxies. Following C15 who placed galaxies beyond the virial radius to adjust the galaxy clustering in the MICE catalogue, we also place satellites beyond the halo virial radius up to three times its value, resulting in approximately 15% of the satellite galaxies being outside the virial radius of their haloes.", '5.3. Galaxy velocities': 'Once the galaxy positions are determined, we compute their velocities within the halo. We also use the parameter color\\_kind to produce line-of-sight velocity dispersion profiles modulated by the type of galaxy.  \nCentral galaxies are assumed to be at rest at the centre of the halo, so their velocity is the same as the centre of mass of the halo. The satellite galaxy velocities are built by solving the spherical, stationary, Jeans equation of local dynamical equilibrium,  \nd GLYPH<16> ρσ 2 r GLYPH<17> d r + 2 β ( r ) GLYPH<16> ρσ 2 r GLYPH<17> r = -ρ ( r ) G M ( r ) r 2 , (14)  \nfor given radial profiles of mass M ( r ), number density ρ ( r ), and velocity anisotropy β = 1 -σ 2 θ /σ 2 r , with σθ and σ r being the velocity dispersion in the tangential and radial directions, respectively. The velocity anisotropy is assumed to follow the Tiret et al. (2007) model:  \nβ ( r ) = β ∞ r r + r , (15)  \nβ  \nwhich is a reasonably good approximation to what is seen in Λ CDM haloes. Using the parameters found for massive clusters at z = 0 . 05 by Mamon et al. (2019), assuming that ellipticals, S0s and spirals respectively trace red sequence, green valley, and blue cloud galaxies as defined by the color\\_kind parameter, leads to the parameters given in Table 1.  \nWe proceed as follows:  \n- 1. Pre-compute the radial velocity dispersions in virial units, σ r ( r ) /v vir on a two-dimensional grid of geometrically spaced r / r vir and concentrations, where the σ r are obtained by solving (following Eq. C.1) the Jeans equation in Eq. (14) using the assumptions of Table 1;\n- 2. pre-fit three two-dimensional 5th-order polynomial approximations (one per color\\_kind ) for log10( σ r /v vir) in terms of log10( r / r vir) and log10 c dark;\n- 3. determine the full set of radial velocity dispersions using these polynomial approximations (extrapolating beyond the fit limits with the constant value at each limit);  \nTable 1. Assumptions for deriving velocity dispersionsNotes. The columns are: (1) the color\\_kind ; (2) the ratio of the concentration of the galaxy population of a given color\\_kind over that of the dark matter; (3) the asymptotic velocity anisotropy (Eq. 15) for the given color\\_kind ; and (4) the anisotropy radius (Eq. 15) in units of the halo virial radius.  \n- 4. determine the tangential components of the velocity dispersions, using σθ = σϕ = p 1 -β ( r ) σ r , where the first equality comes from the assumption of spherical symmetry ( σθ and σϕ are the azimuthal and latitudinal components respectively of the tangential velocity component);\n- 5. derive the three velocity components in spherical coordinates, assuming locally Gaussian velocity distribution functions;\n- 6. convert these physical velocities in the halo frame to peculiar velocities, by subtracting the Hubble flow, H ( z ) r , from the radial velocity component;\n- 7. convert these peculiar halo-frame velocities in spherical coordinates to Cartesian coordinates;\n- 8. convert these peculiar halo-frame Cartesian coordinates into the box frame;\n- 9. derive the line-of-sight velocities V LOS = ( R · V ) / R (where R and V are the position and velocity vectors measured relative to the observer and R the modulus of the R vector) in the box frame, assuming that the observer is at rest in the box;\n- 10. derive the redshift from the line-of-sight velocities using the relation  \n1 + z = (1 + z cos) GLYPH<18> 1 + V LOS c GLYPH<19> , (16)  \nwhere z cos is the cosmological (true) redshift, V LOS is the line-of-sight velocity component, and neglecting the peculiar motion of the observer as well as the 2nd-order transverse Doppler term which incorporates the multiplicative Lorentz factor (1 -V 2 / c 2 ) -1 / 2 (e.g., Tatum 1985).  \nThe details of steps Nos. 1 and 2 are provided in Appendix C. In a forthcoming (already coded) version of Flagship, we inserted after item No. 6 the addition of a radial infall pattern, obtained by fitting a nonlinear function of redshift and log halo mass to the radial motions inside haloes of dissipationless cosmological simulations.', '5.4. Spectral energy distributions': 'As described in Sect. 5.1, we assign a ( g 01 -r 01)HOD colour to each galaxy randomly sampling the colour distributions derived from the CMD for our three color\\_kind populations. This colour corresponds to the rest-frame colour of a galaxy in the g 01 and r 01 filters (see Sect. 5.2) or equivalently the g -r colour at redshift z = 0 . 1.  \nNext, we want to assign a spectral energy distribution (SED) to each galaxy. We take as a template basis the COSMOS SED library used in Ilbert et al. (2009) that originally comes from the SEDs of Polletta et al. (2007), complemented with templates  \nTable 2. SED template basis  \nNotes. a The reddening values sampled for the Prévot and Calzetti extinction laws are 0.0, 0.1, 0.2, 0.3, 0.4, and 0.5  \nFig. 17. Subset of six SED templates out of the 136 described in the text at redshift z = 0 . 0. The templates are normalised to have the same flux at a wavelength of 6000 Å. The figure inset indicates the colour-ordered template number (o), the original template number in Ilbert et al. (2009) ( t ), the E ( B -V ) reddening value applied to the template (ebv) and the value of the g 01 -r 01 colour ( c ).  \n<!-- image -->  \nfrom Bruzual & Charlot (2003) to expand the range of starformation histories included. In order to increase the SED coverage, each galaxy SED is computed as a linear combination of two of the templates in the SED basis.  \nFollowing Ilbert et al. (2009), our chosen SED library is composed of 31 templates. We add extinction to some of the templates to expand the colour space coverage and better represent observations. We do not apply any extinction to the first ten reddest templates. The next thirteen templates represent spectral distributions of spiral galaxies and are assumed to have the Prévot Small Magellanic Cloud extinction law (Prévot et al. 1984). The last eight bluest templates represent starburst galaxies and have the Calzetti extinction law (Calzetti et al. 2000). We build a sample of 136 templates using the 31 SEDs and the two extinction laws that we sample at six values of the E ( B -V ) colour excess (or reddening) from 0 to 0.5 in steps of 0.1. Table 2 summarises the original templates used with their order number, the extinction law applied for those templates, the reddening values sampled for each of those templates, and the total number of resulting SEDs that constitute our SED template basis.  \nFor these 136 templates, we compute the g 01 -r 01 and the COSMOS g -r colours at seven discrete redshift values covering the redshift range of the lightcone simulation, 0 < z < 3, in intervals of ∆ z = 0 . 5. Figure 17 shows a subset of six of these 136 templates at redshift z = 0 . 0.  \nWerank order the templates according to their colour at these redshift values. 12 We use the g 01 -r 01 colour at z = 0 and the  \nFig. 18. Relation between the g 01 -r 01 colour and the colour-ordered template (which is zero-ordered from 0 to 135) at redshift z = 0. The blue line connects the values of the colour and the ordered template for the 136 templates. The orange line is a smoothed version of the relation connecting the data points.  \n<!-- image -->  \nCOSMOS g -r at the rest of the redshift values for the ordering. We compute a relation between the template order number, treated as a real number but sampled at integer values, and the values of the colour at each redshift value. We smooth this relation to make the function monotonic and avoid degenerate values in the relation that will give the same template order number for di ff erent values of the colour. Figure 18 shows this relation for the g 01 -r 01 colour at redshift z = 0. We also establish relations between the g 01 -r 01 and g -r colours at all non-zero redshift steps. To obtain these relations we abundance match the ( g 01 -r 01)HOD colour distribution in the Flagship catalogue to the g -r colour distribution of the COSMOS2020 catalogue (Weaver et al. 2022) at each redshift value considered (see footnote 12). We take all galaxies in a range of ∆ z = 0 . 2 centred at the mean redshift value to build the samples for the abundance matching. The only exception is for the last redshift value, z = 3 . 0, where we take galaxies in the range 2 . 8 < z < 3 . 0 as there are no Flagship galaxies beyond that redshift. We also compute the E ( B -V ) distributions from the COSMOS2020 catalogue at the same redshift values, in this case in a redshift range of width ∆ z = 0 . 1. Figure 19 shows the normalised distribution of E ( B -V ) values in the COSMOS2020 catalogue for a sample with a magnitude limit in this catalogue of i < 24 . 5, equivalent to the 10 σ magnitude limit of the Euclid Wide Survey in the I E filter. The distribution transitions from being dominated by galaxies with no internal extinction at z = 0 to peaking at values of E ( B -V ) ∼ 0 . 2 at redshifts z ≳ 1.  \nOur procedure to assign an SED to a galaxy consists of the following steps. First, we determine the two redshift values immediately lower, z 1, and larger, z 2, than the redshift of the galaxy from the list in footnote 12, so that z 1 < z < z 2. We determine two SEDs from our 136 template basis, one at each of these two redshifts, SED ( z 1) and SED ( z 2), and compute the final SED of the galaxy as a linear combination of the two templates (Eq. 17) with weights, w 1 and w 2 proportional to their redshift di ff erence to the reference redshift as in Eq. (18),  \nSED( z ) = w 1 SED ( z 1) + w 2 SED ( z 2) , (17)  \nFig. 19. Distribution of E ( B -V ) values in the COSMOS2020 catalogue for a magnitude range i < 24 . 5. The distributions are normalised so that the total abundance for all E ( B -V ) values at each redshift adds up to 1.  \n<!-- image -->  \nw 1 = z 2 -z z 2 -z 1 , (18) w 2 = z -z 1 z 2 -z 1 .  \nWe derive the SED at each of these two redshifts from the ( g 01 -r 01)HOD colour of the galaxy. We obtain the ordered template number from the colour using the relations found previously (e.g., Fig. 18). At z = 0, we use the relation with the ( g 01 -r 01)HOD colour. At higher redshift we first obtain the COSMOS g -r colour from the ( g 01 -r 01)HOD colour using the abundance matched relations described before and then use this g -r colour to obtain the ordered template number. We take the six closest ordered templates to this ordered template number from our 136 SED template basis and compute the colour di ff erence from the original colour, either ( g 01 -r 01)HOD or g -r , depending on the redshift, and the colour of these six ordered templates. We compute a probability based on this colour di ff erence for each of these six ordered templates, P colour. This probability is assumed to be a Gaussian distribution with mean the original colour and standard deviation half the mean value of the colour di ff erences for the six ordered templates. If the standard deviation is lower than 0.015, we set it to this minimum value. We then compute the probability of the template given its E ( B -V ) value, PE ( B -V ), using the distributions of E ( B -V ) values from the COSMOS2020 catalogue obtained before (Fig. 19). We compute the final probability of each ordered template as the product of these two probabilities in colour and colour excess, P ot = P colour PE ( B -V ). We generate the cumulative distribution of this probability and draw a random number from which we obtain the ordered SED template using the inverse of the cumulative probability distribution. We perform this procedure at the redshifts z 1 and z 2 and obtain the final SED following Eqs. (17) and (18).', '5.5. Galaxy shapes and sizes': 'Galaxy shapes and sizes are assigned using phenomenological prescriptions similar to those of Miller et al. (2013) with calibration data coming from Hubble Space Telescope (HST) observations. We use two main data sets as calibrators: the Cosmic Assembly Near-Infrared Deep Extragalactic Legacy Survey 13 (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011) reductions of Dimauro et al. (2018) and the Advanced Camera  \nFig. 20. Fraction of galaxies modelled with one component model with respect to the total as a function of i -band magnitude. The green line is the resulting fraction in the Flagship catalogue when the relation in Eq. (19) is used. The orange points with error bars represent the fraction of galaxies modelled with a one-component model in the CANDELS catalogue of Dimauro et al. (2018).  \n<!-- image -->  \nfor Surveys 14 (ACS) reductions of the GOODS-South field of Niraj Welikala (hereafter referred to as HST\\_GS, private communication).  \nMost galaxies are assumed to have two components: a Sérsic profile (Sérsic 1963) component, which we will refer to as bulge, and an exponential disk component. The exception are galaxies fulfilling the following relation:  \n( g 01 -r 01)HOD > (0 . 96 + 0 . 04 z ) (19) -0 . 015 [( Mr 01 -5 log10 h ) + 24 . 0] ,  \nwhere z is the galaxy redshift and Mr 01 -5 log10 h is the absolute magnitude in the r 01 filter, that are described by a single Sérsic profile component, which we will also call bulge. We identify this modelling di ff erence with the parameter dominant\\_shape in the catalogue which is set to 0 for the one-component (bulge only) profile galaxies and to 1 for the two-component (bulge and disk) galaxies. The selection in Eq. (19) implies that the reddest galaxies at each magnitude are the ones modelled with a onecomponent (bulge only) profile. The fraction of these galaxies to the total number of galaxies varies as a function of magnitude as shown in Fig. 20, with the fraction getting smaller as magnitudes get fainter. We chose to use a simple magnitudedependent colour cut to separate the way to model the galaxies. This cut reproduces well the fraction in the CANDELS catalogue of Dimauro et al. (2018) at magnitudes brighter than i < 22, but it is lower at fainter magnitudes. For these fainter galaxies, their smaller sizes and lower signal-to-noise detections make the distinction of two profile components harder. The upward trend seen in CANDELS at i > 23 can probably be attributed to this effect. Taking this into account and as the one- and two-component models are not very di ff erent for these small sizes for the Euclid pixel size, we decided not to change this simple selection cut for faint magnitudes. As the catalogue is dominated by faint galaxies, the overall fraction of bulge-only galaxies is approximately 10% of the total galaxy population in the catalogue.  \nFig. 21. Scale length values in the HST\\_GS sample (blue dots). The black line is the fit to the median values of the distribution as a function of the i -band magnitude as given in Eq. (20)  \n<!-- image -->  \nWe start by assigning a scalelength, R h, to the galaxies based on their i -band magnitude. 15 Wecompute the median scalelength R median h in arcsec using the relation fit from the HST\\_GS catalogue, which is similar to one obtained in Miller et al. (2013),  \nR median h 1 \'\' = exp [ -1 . 675 -0 . 215 ( i -25 . 0)] . (20)  \nThe observed distribution of scalelength values in magnitude bins can be described by a power law and an exponential decay with parameters depending on magnitude. The cumulative function of such a functional form is an incomplete gamma function. Therefore, we assign the value of the scalelength as  \nR h = R median h 1 . 13 " Γ -1 1 + β α , µ !# 1 α , (21)  \nwhere Γ -1 is the inverse of the incomplete gamma function, 16 α is the exponent of the exponential decay argument, which depends on magnitude as α = 2 . 0 -0 . 08 ( i -18) with i being the i -band magnitude, β is the exponent of the power-law, which we fix at β = 1 . 5 and µ is a uniformly distributed random number in the range [0,1].  \nWe then set the luminosity fraction in the Sérsic profile, F S\'ersic , which we refer to as bulge in the catalogue, using the following relation:  \nF S\'ersic = a " 1 + 1 a ! µ -1 # , (22)  \nwhere a = 0 . 01 + 0 . 0015 (25 -i ) and µ is a uniformly distributed random number in the range [0,1]. We set a floor of 0.007 to the minimum value of a . Note that this bulge component can either be the one-component profile for galaxies with bulge \\_ fraction = 1 and dominant \\_ shape = 0 or the Sérsic profile bulge component for the two-component galaxies ( bulge \\_ fraction < 1 and dominant \\_ shape = 1).  \nFig. 22. Bulge fraction distribution for galaxies brighter than I E < 24 . 5 in the Flagship catalogue compared to the CANDELS catalogue of Dimauro et al. (2018) for the same magnitude range.  \n<!-- image -->  \nWe set the disk profile as a Sérsic profile with index n S\'ersic = 1. We set the half-light radius of the disk as R disk 50 = 1 . 678 Rh . An exponential profile contains ∼ 26% of the light within the scalelength radius. The factor 1.678 comes from the ratio of radii containing 50% and 26% of the light in a Sérsic n S\'ersic = 1 profile.  \nWe set the Sérsic half-light radius, R S\'ersic 50 , of the Sérsic profile from the disk half-light radius, R disk 50 , with  \nlog10       R S\'ersic 50 1 \'\'       = -1 . 81 + 2 . 6 1 + exp GLYPH<26> -2 . 15 GLYPH<20> log10 GLYPH<18> R disk 50 1 \'\' GLYPH<19>GLYPH<21> + 0 . 3 GLYPH<27> + 0 . 7 1 + exp GLYPH<26> -15 . 0 GLYPH<20> log10 GLYPH<18> R disk 50 1 \'\' GLYPH<19>GLYPH<21> -0 . 1 GLYPH<27> . (23)  \nThis equation tries to fit the relation between these radii in the CANDELS calibration catalogue. In Fig. 23, we show the Sérsic half-light radius, R S\'ersic 50 , distribution for composite profiles and compare it to the CANDELS. The Flagship distribution of R S\'ersic 50 values is similar to the CANDELS distribution, although somewhat narrower. The HST\\_GS sample has also measured values of the bulge component half-light radius in composite profiles. However, those values are approximately a factor of five smaller than the ones in CANDELS. The measuring methods were di ff erent in both samples. The galaxies in CANDELS were fit individually to two components, while in HST\\_GS galaxies were combined and then the bulge half-light radius of the combined galaxy was fit as a multiplicative factor of the scale length, R h. Given the individual measurement employed, we deem the CANDELS values more reliable and decided to fit to their distribution. But, one should keep in mind the systematic uncertainty in these measurements depending on the methodology employed when using these values.  \nFor one-component profiles, we assume the half-light radius of the bulge to be the same as the scalelength radius. Fig. 24 shows the distribution of these R S\'ersic 50 values for galaxies with i < 24 . 5 in the catalogue compared to the CANDELS sample.  \nWe obtain the Sérsic index for the Sérsic profile components with functions that try to reproduce the distribution of Sérsic indices in the CANDELS calibrating sample. For the composite profile, we draw Sérsic indices with  \nn S\'ersic = [( a + n S\'ersic min ) ν max] -a , (24)  \nFig. 23. Distribution of half-light radius, R S\'ersic 50 , in the Flagship catalogue for the Sérsic bulge component of galaxies simulated as two components compared to the CANDELS sample.  \n<!-- image -->  \nFig. 24. Distribution of half-light radius, R S\'ersic 50 , in the Flagship catalogue for one-component galaxies compared to the CANDELS sample.  \n<!-- image -->  \nwhere  \nν max =       a + n S\'ersic max a + n S\'ersic min       µ , (25)  \na = 0 . 2, the minimum and maximum Sérsic indices, n S\'ersic min = 0 . 5 and n S\'ersic max = 6 . 0, and µ is a uniformly distributed random number in the range [0,1].  \nFor the single component, we draw the Sérsic index from  \nn S\'ersic = n S\'ersic o ff + a " Γ -1 1 + β α , ν !# 1 α , (26)  \nwith n S\'ersic o ff = 0 . 5, a = 2 . 75, α = 1 . 9, β = 1 . 0 and Γ -1 the inverse of the incomplete gamma function as in Eq. (21). The parameter ν is given by  \nν = Γ n S\'ersic min -µ ( Γ n S\'ersic min -Γ n S\'ersic max ) , (27)  \nwith  \nΓ n S\'ersic min = Γ        1 + β α ,       n S\'ersic min -n S\'ersic o ff a       α        , (28a)  \nΓ n S\'ersic max = Γ        1 + β α ,       n S\'ersic max -n S\'ersic o ff a       α        , (28b)  \nFig. 25. Distribution of Sérsic indices, n S\'ersic , in the Flagship catalogue for the two-component galaxies (top) and the one-component galaxies (bottom) compared to the CANDELS sample for galaxies with a magnitude limit of i < 24 . 5.  \n<!-- image -->  \nwhere Γ is the incomplete gamma function; 17 n S\'ersic min = 0 . 5 and n S\'ersic max = 6 . 0, as before, and µ is a uniformly distributed random number in the range [0,1]. Figure 25 shows the distribution of Sérsic indices for the bulges of the two-component (top) and the one-component (bottom) galaxies compared to the CANDELS sample. We have restricted the comparison to a magnitude limit of i < 24 . 5. We have imposed a minimum and a maximum value to the Sérsic indices as indicated above because the Euclid pixel simulation pipeline (OU-SIM, Euclid Collaboration: Serrano et al. 2024) and subsequent Science Ground Segment pipelines were very ine ffi cient in dealing with very extended profiles. The maximum value is smaller than the maximum value allowed in the CANDELS measurements. As the distributions shown in Fig. 25 are normalised to have an integral equal to 1, the Flagship curves have larger abundances compared to CANDELS as they sample a smaller range.  \nIn the Euclid pixel simulation pipeline, we render the disk component of the two-component galaxies with the inclined disk model of Galsim 18 (Rowe et al. 2015). We compute this inclination angle, θ incl, with  \nθ incl = 90 -a " Γ -1 1 + β α , ν !# 1 α , (29)  \nwhere θ incl is in degrees and with  \nν = µ Γ θ max incl , (30)  \nand  \nΓ θ max incl = Γ " 1 + β α , θ max incl a ! α # , (31)  \n17 Implemented with the python scipy function  \nscipy.special.gammainc  \n18 https://github.com/GalSim-developers/GalSim  \nFig. 26. Distribution of inclination angles, θ incl, in the Flagship catalogue for the disk component of the two-component galaxies compared to the HST\\_GS sample for galaxies with a magnitude limit of i < 24 . 5.  \n<!-- image -->  \nwith a = 22 · , α = 1 . 39, β = 1 . 65, θ max incl = 90 · and µ a random number uniformly distributed in the range [0,1]. Figure 26 shows the distribution of inclination angles for the disk component of the galaxies simulated with two components. The distribution of inclination angles in the HST\\_GS sample is shown for comparison.  \nFor the disk component, we calculate an axis ratio, q = b / a , where b and a are the semi-minor and semi-major axis respectively, from the inclination angle with  \nq = q (cos θ incl) a + ( c sin θ incl) b , (32)  \nwith parameters a = 2 . 0153, b = 4 . 3161 and c = 0 . 35696 for the disk component of the two component galaxies.  \nWe compute the ellipticity from the axis ratio as  \nϵ = 1 -q 1 + q . (33)  \nWe compute the axis ratio of the bulge of the two components profiles with the same Eq. (32), but with parameter values that depend on the Sérsic index as  \na = 1 . 905 + 0 . 25 1 + exp GLYPH<2> 0 . 95 GLYPH<0> n S\'ersic -0 . 80 GLYPH<1>GLYPH<3> ,  \n(34a) b = 1 . 18 + 0 . 50 1 + exp GLYPH<2> -2 . 0 GLYPH<0> n S\'ersic -0 . 80 GLYPH<1>GLYPH<3> + 0 . 95 1 + exp GLYPH<2> 3 . 8 GLYPH<0> n S\'ersic -2 . 45 GLYPH<1>GLYPH<3> , (34b) c = 0 . 465 + 0 . 12 1 + exp GLYPH<2> -1 . 0 GLYPH<0> n S\'ersic -1 . 80 GLYPH<1>GLYPH<3> + 0 . 155 1 + exp GLYPH<2> 2 . 0 GLYPH<0> n S\'ersic -2 . 66 GLYPH<1>GLYPH<3> + 0 . 06 1 + exp GLYPH<2> -4 . 5 GLYPH<0> n S\'ersic -3 . 6 GLYPH<1>GLYPH<3> . (34c)  \nWe compute the ellipticity of the bulge one-component profiles with relations similar to those of Eqs. (29), (30), and (31),  \nϵ = a " Γ -1 1 + β α , ν !# 1 /α , (35)  \nArticle number, page 20 of 45  \nFig. 27. Distribution of ellipticity values, ϵ , in the Flagship catalogue for the disk component (top) and Sérsic component (middle) of the twocomponent galaxies compared to the HST\\_GS sample; and the Sérsic component of the one-component galaxies (bottom) compared to the CANDELS sample for galaxies with a magnitude limit of i < 24 . 5.  \n<!-- image -->  \nwhere  \nν = µ Γ ϵ max , (36)  \nand  \nΓ ϵ max = Γ " 1 + β α , GLYPH<18> ϵ max a GLYPH<19> α # , (37)  \nwith a = 0 . 26, α = 1 . 5, β = 0 . 6, ϵ max = 0 . 8 and µ a random number uniformly distributed in the range [0,1]. Fig 27 shows the resulting ellipticity distributions for the disk and bulge components of the galaxies simulated as two-components and onecomponent. Finally, we set the orientation angle of the galaxy randomly.', '5.6. Galaxy intrinsic alignment': "In addition to the shape modelling described in the previous subsection, we model shapes and orientations in a separate process in order to simulate the intrinsic alignment (IA) signal, as detailed in Ho ff mann et al. (in prep.).  \nOur IA model consists of two steps. Approximating each galaxy as a 3D ellipsoid, we assign in the first step two 3D  \naxis ratios to each object, taking into account its redshift, absolute magnitude, and rest-frame colour. The parameters of the shape model component have been calibrated such that the distribution of projected 2D axis ratios matches the observed distribution from the COSMOS survey. The 3D orientations of the galaxies' principal axes are modelled in a second step, using two di ff erent methods for central and satellite galaxies. Centrals are aligned with the principal axes of their host halo, while satellites are pointed towards their host halo centre. These initial orientations are then randomised depending on each galaxy's redshift, magnitude, and colour, which allows for calibrating the dependence of the resulting IA signal on galaxy properties. Once the 3D axis ratios and orientation have been assigned, the 2D intrinsic shear components are obtained via projection along the observer's line-of-sight.  \nA novel feature of the IA model, in comparison to its predecessors, is its calibration against alignment signals derived from multiple constraining datasets through a thorough exploration of the simulation parameter space. At redshifts below z = 0 . 36 we calibrate against the IA signal measured from luminous red galaxies in the LOWZ sample of the Baryon Oscillation Spectroscopic Survey (Singh & Mandelbaum 2016) as well as samples of red and blue galaxies from the Sloan Digital Sky Survey (Johnston et al. 2019). These low-redshift constraints are complemented by alignment measurements in three magnitude-limited samples at z = 1 . 0 from the Horizon AGN simulation (Dubois et al. 2014).  \nValidations of the alignment model are presented in Ho ff -mann et al. (in prep.). A comparison of the resulting signal against theory predictions are studied in Paviot et al. (in prep.), whereas predictions on the IA parameters for Euclid -like samples are investigated in Tutusaus et al. (in prep.).", '5.7. Physical parameters': 'The Euclid galaxy clustering main probe will select emission line galaxies as the tracers to sample the large-scale structure of the Universe and conduct cosmological inference. In order to compute the emission line fluxes, we first need to estimate other physical properties of the galaxies that the emission lines depend on.  \nWe start by computing the galaxy ultraviolet (UV) photon flux density from the galaxy SED. We integrate the galaxy SED, without the extinction component, derived in Sect. 5.4:  \nf no -ext UV = R f no -ext SED R UV( λ ) λ 2 d λ c R R UV( λ ) d λ , (38)  \nwhere f no -ext SED is the SED flux density, c is the speed of light and R UV is a top-hat filter response in the UV with full transmission in the wavelength range from 1500 to 2300 Å and zero transmission otherwise.  \nThe star-formation rate (SFR) is computed from the UV flux, f no -ext UV in Eq. (38) following the relation of Kennicutt (1998) for a Chabrier initial mass function (IMF, Chabrier 2003),  \nSFR M ⊙ yr -1 = 0 . 8 × 10 -28 L UV = 9 . 57 × 10 21 f no -ext UV d 2 L , (39)  \nwhere the constants are for L UV given in erg s -1 Hz -1 units; f no -ext UV , in erg cm -2 s -1 Hz -1 ; and the luminosity distance, d L, in h -1 Mpc. Equivalently,  \nlog10 SFR M ⊙ yr -1 ! = -0 . 4 GLYPH<16> M no -ext UV + 18 . 65 GLYPH<17> . (40)  \nWe derive the galaxy stellar mass, M , from the galaxy SED and its luminosity. First, we estimate the stellar mass-to-light ratio from the galaxy colours and SED. As our reference band to compute luminosities is the r 01 filter (see Sect. 5.1), we estimate the stellar mass-to-light ratio in this band. The stellar mass is then obtained simply by multiplying its r 01-band luminosity by the stellar mass-to-light ratio. The galaxy metallicity is computed from the galaxy stellar mass following the relation of Curti et al. (2020), which we slightly modify and extrapolate to cover the whole stellar mass range covered by our sample. We assume a relation for the mean value of the stellar mass-metallicity given by  \n12 + log10 O H ! ( M , SFR , z ) = 6 . 7 + 2 . 22 1 + exp [ -0 . 85 ( M-M ref)] , (41)  \nwhere M is the stellar mass and M ref is a reference value given by  \nM ref(SFR , z ) = 7 . 5 + 0 . 5 log10 SFR M ⊙ yr -1 ! + 0 . 1( z -1 . 0) , (42)  \nwhere SFR is the star formation rate and z , the galaxy redshift. Figure 28 shows the stellar mass-metallicity relation (MZR) as a function of redshift and SFR in the top panel. In the bottom panel, we show the contours of a galaxy sample selected at I E < 24 . 5 and redshift z = 1 in stellar mass-metallicity space together with three MZR lines for di ff erent SFR values.  \nFor each galaxy, we draw a realisation assuming a Gaussian distribution of the metallicity values around its mean, with a standard deviation (e.g., Tremonti et al. 2004; Curti et al. 2020)  \nσ = 0 . 01 + 0 . 03 { 8 . 92 -[12 + log10(O / H)]mean } , (43)  \nwhere [12 + log10(O / H)]mean is the mean value of the metallicity given by Eq. (41).', '5.8. Galaxy emission lines': 'The determination of the galaxy redshift with the Euclid NISP instrument (Euclid Collaboration: Jahnke et al. 2024) in slitless mode rely mostly on the identification of emission lines. Therefore, the computation of emission line fluxes is important to understand the completeness and purity of the emission lineselected galaxies that will be used for the clustering cosmological probe.  \nFirst, we compute the H α line flux, f H α , from the starformation rate using the Kennicutt (1998) relation adapted to the Chabrier IMF, using the un-extincted UV absolute AB magnitude  \nlog10 f H α erg cm -2 s -1 ! = -0 . 4 M no -ext UV -2 log10 d L h -1 Mpc ! -16 . 44 , (44)  \nwhere f H α is the flux in the H α line, M no -ext UV is the unextinguished UVabsolute AB magnitude, and dL , the luminosity distance. We add a random scatter following a Gaussian distribution of standard deviation of 0.05 in the logarithm of the flux to the computed logarithmic value of the flux.  \nOur H α fluxes are estimated from the UV fluxes that rely on the SED assignment procedure, which is based on observed optical colours (see Sect. 5.4). We then compute the dustextinguished H α flux. We get the stellar continuum extinction  \nFig. 28. Top : Mean values of the metallicity as a function of the stellar mass for di ff erent values of redshift (colour code in the figure legend) and SFR (dashed, solid, and dotted lines correspond to log10[SFR / ( M ⊙ yr -1 )] = 1 , 0 , -1, respectively). Bottom : The contours show the distribution of metallicity values as a function of stellar mass for a sample at redshift z = 1 . 0 with magnitude limit i < 24 . 5. The red lines indicate the mean value of the metallicity-stellar mass relation, Eq. (41), for the values of SFR indicated in the legend (that are the same as those in the top panel).  \n<!-- image -->  \nFig. 29. Redshift distribution (galaxies per deg 2 per unit redshift) of samples selected with a flux limit cut of f H α = 2 × 10 -16 erg s -1 cm -2 . The sold lines show the redshift distributions from models 1 (light blue) and 3 (light green) of Pozzetti et al. (2016). The dashed lines are the distributions in the Flagship catalogue for the same flux cut for the model 1 (orange) and 3 (red) simulated fluxes.  \n<!-- image -->  \nat the H α wavelength using the extinction law and colour excess value E ( B -V ) of the galaxy SED (Sect. 5.4). We then convert the stellar continuum extinction to a nebular emission line extinction using the redshift-dependent factor 0 . 44 + 0 . 2 z (Calzetti et al. 2000; Saito et al. 2020), which we clip at 0.44 and 1.0 as minimum and maximum values, respectively, and apply this extinction factor to the H α flux. We provide both the extinguished  \nFig. 30. BPT diagram at z = 1 . 20. The red contours show the density of Flagship galaxies in this diagram. The black solid line is the relation of Kewley et al. (2013a).  \n<!-- image -->  \nand un-extinguished fluxes in our catalogue. Finally, we calibrate the resulting f H α distribution to the Pozzetti et al. (2016) models applying a correction to the computed H α fluxes with abundance matching to the models. Figure 29 shows the redshift distributions for the simulated model 1 and model 3 fluxes.  \nIn the SED assignment procedure, we compute each SED as a linear combination of SEDs at fixed redshift values. In particular, we assign the extinction with distributions computed at particular redshifts given in footnote 12. The H α flux is computed from the unextinguished UV flux. In the process of subtracting the extinction from the SEDs, we introduce discontinuities in the H α flux-redshift distribution at those particular redshifts. As a temporary fix, we compute a correction of the H α fluxes to smooth the distributions at those redshifts. In Fig. 29, we can see the performance of this smoothing procedure that was applied at redshifts z > 0 . 5, but not at z = 0 . 5 as there is not enough volume in the simulation lightcone to do it properly.  \nThe other hydrogen lines of the Balmer and Paschen series are computed from the H α flux assuming case B recombination from Osterbrock & Ferland (2006). The [O ii ] and [O iii ] fluxes are computed from the H β fluxes and the metallicity of the galaxy taking into account the relations of the [O ii ] / H β and [O iii ] / H β flux ratios as a function of metallicity of Curti et al. (2020). The [N ii ] fluxes are calculated from the [O iii ] / H β flux ratio and the H α flux following the Baldwin-Phillips-Terlevich (BPT) diagram (Baldwin et al. 1981) using the relations of Kewley et al. (2013a) that include redshift evolution. 19 The [S ii ] flux is computed from the H α flux and the metallicity following the relation of Curti et al. (2020). We obtain the [S iii ] flux from the [S ii ] flux and a fit to the relation between the [S iii ] / [S ii ] and [O iii ] / H β flux ratios shown by Mannucci et al. (2021) in their Fig. 2. We generate a realisation of all fluxes adding scatter to the value obtained with the relations used. As an example, Fig. 30 shows the position of the Flagship galaxies at z = 1 . 2 and with H α flux f H α > 2 × 10 -16 erg s -1 cm -2 in this diagram, together with the Kewley et al. (2013a) relation at this redshift.  \nFig. 31. Flagship galaxy catalogue: stripes show sections of the lightcone for z < 3 for dark matter haloes (top panel), the photometric sample for I E < 24 . 5 (middle), and an H α spectroscopic galaxy sample for a flux cut f H α > 2 × 10 -16 erg cm -2 s -1 (bottom), which selects galaxies only within the range 0 . 9 < z < 1 . 8 given the NISP red grism wavelength coverage. The stripe abscissae share the same scaling in comoving distance. The colour scaling in each panel depends on the density of the sample.  \n<!-- image -->', '5.9. Galaxy lensing': "The implementation of lensing properties on mock galaxies uses the all-sky dark matter lensing maps discussed in Sect. 3. Given this set of 2D maps covering the entire redshift range of the lightcone output of the simulation, our procedure is the same as the one originally implemented in the MICE simulations (Fosalba et al. 2015). This is a simple 3-step algorithm that we describe below.  \n- 1. for a given galaxy at the 3D position in the lightcone (ˆ n , z ), where ˆ n gives its angular position in the sky and z its redshift, find the corresponding 3D pixel in the discretised lightcone, with pixel center coordinates, (ˆ n i , zj ), where the galaxy sits in (i.e., the 3D pixel in the suite of 'onion slices' or all-sky lensing maps in HEALPix tessellation described in Sect. 3),\n- 2. get the lensing values for this 3D pixel using the dark matter all-sky lensing maps, L i , j ≡ L ( ˆ n i , zj ), where the components of the lensing vector are L = ( κ, γ 1 , γ 2), which is convergence and shear, and\n- 3. assign these pixelised dark matter lensing values, L i , j , to the mock galaxy.  \nThe resulting implementation of galaxy lensing using the above method is thus limited by the pixel resolution used, HEALPix N side = 8192, which corresponds to a pixel scale of 0 . ' 43. Therefore, we only expect to accurately model lensing observables down to ∼ 0 . ' 5 scales, as we will discuss in detail when we validate the mock properties below (see Sect. 6.3). A further limitation intrinsic to this method is that di ff erent galaxies falling within a given 3D pixel in the 'onion universe' grid pattern of the lightcone will have identical lensing properties. These two limitations can be overcome using the same approach but using a finer pixel scale (i.e., higher N side) and / or using interpolation schemes.", '6. Mock validation': 'Fig 31 shows a small rectangular portion of the Flagship catalogue in the lightcone. 20 The figure illustrates the breath of the simulation and covers a region of 4368 h -1 Mpc (horizontal direction) × 300 h -1 Mpc (vertical direction) with a depth of 80, 160 and 200 h -1 Mpc for the top, middle and bottom panels respectively, with the observer located in the left-bottom corner. In the top panels we show the dark matter haloes where one can see the growth of structure progressing from right, z = 3, to left, z = 0. The middle and lower panels present likely galaxy samples to be used in the Euclid cosmological analysis. The middle panel shows the galaxies selected with an I E < 24 . 5 magnitude cut, which we expect to be representative of the weak lensing sample. In the lower panel, we show a sample selected with an H α flux cut f H α > 10 -16 erg cm -2 s -1 , trying to mimic the galaxy clustering sample.', '6.1. General galaxy properties': "In order to assess the performance and possible limitations of the Flagship galaxy catalogue, we compare the galaxy properties, distributions, and relations to observations. In particular, we compare the mock galaxies to observed quantities relevant for the Euclid cosmological probes, including observed galaxy number densities, colours, emission lines fluxes and ratios, as well as sizes at di ff erent redshifts.  \nIn Fig. 32, we compare the number counts of the Flagship simulation to literature data and COSMOS2020 (Weaver et al. 2022) counts of magnitudes corrected for Galactic extinction. We show the number density per unit area and magnitude in the I E and H E bands, as well as in similar filters observed in COSMOS, without applying further corrections to account for the di ff erences in the filter transmissions. Also, fluxes and magni-  \n<!-- image -->  \nFig. 32. Galaxy number densities in the I E and H E bands compared to literature data: COSMOS2020 (Weaver et al. 2022) counts in H VISTA and i HSC are derived from the 'farmer' version of the photometric catalogue, after removing masked regions and objects classified as stars; Driver et al. (2016); Durham Cosmology Group's counts are from the compilation available at the dedicated web page ( http://star-www. dur.ac.uk/~nm/pubhtml/counts/counts.html ).  \n<!-- image -->  \ntudes from Flagship are intrinsic, while the observed ones are a ff ected by the photometric noise, and their distribution could, therefore, be broader.  \nIn general, we find excellent agreement with observations in the optical band from the brightest to the faintest objects in the mock. We find, however, some di ff erences in the H E band, where the galaxy number density in the mock is about 40% higher than in the literature. At the EWS limit of H E < 24 . 0 (detection limit at 5 σ for point sources) the integrated number density in the mock is about 188 000 deg -2 compared to ∼ 131 000 deg -2 in the COSMOS field.  \nIn Fig. 33, we also explore how galaxy colours compare with the COSMOS2020 catalogue: the optical ( r -i ) vs. ( u -r ) diagram is quite well reproduced at the EWS magnitude limit, with a small o ff set. The evolution with redshift of the i HSC -H VISTA galaxy colour (a proxy for I E -H E), not used in the input calibration, also reproduces well the trend observed in COSMOS2020 (for which we are using the photometric redshift obtained with LePhare; Ilbert et al. 2006). The lack of galaxies with the bluest colours in Flagship compared to the COSMOS2020 distribution is likely to be the cause of the overestimation in the H -band counts.  \nFig. 33. Galaxy colour-colour and colour-redshift diagrams compared to COSMOS2020 data.  \n<!-- image -->  \nIn Fig. 34, we show the galaxy stellar mass function (GSMF) derived for the full octant of our Flagship lightcone in various redshift bins. The high resolution of the simulation allows us to derive the GSMF down to very low mass limits ( ∼ 10 8 M ⊙ ) up to z = 3. Comparison to observations in the COSMOS field (Ilbert et al. 2013) shows a good consistency. In particular, FS2 is always consistent at high mass while at z > 1 . 5, it shows an excess below the characteristic galaxy mass of 10 10 M ⊙ .  \nIn Fig. 35, we show the SFR mass relation at di ff erent redshifts from the simulation. The comparison to observational data in the COSMOS field is fairly satisfactory, with a small remaining o ff set and broader relation in COSMOS, due to observational errors. The wide Euclid survey area will allow us to map this relation up to very high mass where previous observations are limited by area and statistics.  \nBesides photometric properties, the most important quantities to validate for the galaxy clustering Euclid cosmological probe are the emission line fluxes and number densities. Fig 36 shows the emission line galaxy densities for a flux limit of f H α + [N ii ] > 2 × 10 -16 erg s -1 cm -2 , together with the global EWS limit of I E < 24 . 5 (10 σ ) and H E < 24 (5 σ ). As expected, the H α emitters spectroscopically observed by Euclid will be only a few percent of the total EWS photometric sample. Finally, as described in Sect. 5.8, we verify that the Flagship mock has been calibrated to the empirical models of Pozzetti et al. (2016). Some di ff erences are still present, when we consider the combination of H α + [N ii ], mainly due to [N ii ] contribution to the total flux, which is assumed constant in the models and not in the Flagship catalogue. In particular, the redshift distribution of fluxes calibrated on model 3 is fairly consistent with the one from Bagley et al. (2020) when we consider a sample limited in H α + [N ii ] fluxes to 2 × 10 -16 erg s -1 cm -2 . From the Flagship catalogue,  \nFig. 34. Galaxy stellar mass function at di ff erent redshifts. In red the GSMF derived from the Flagship lightcone compared to GSMF derived at z < 0 . 1 in GAMA (Driver et al. 2022, yellow line) and at higher redshifts in COSMOS from the ULTRAVISTA sample (Ilbert et al. 2013, green squares) and the COSMOS2020 release (Weaver et al. 2023, blue lines). Note that we have transformed the Flagship stellar masses to units of the Hubble constant of H = 70 h 70 km s -1 Mpc -1 to bring the units of the Flagship catalogue to the ones used by the other surveys.  \n<!-- image -->  \ncalibrated on the models 3 and 1 of Pozzetti et al. (2016), we expect that the EWS will map a density of about 3500-6300 deg -2 H α + [N ii ] emitters with fluxes above 2 × 10 -16 erg s -1 cm -2 in the redshift range 0 . 9 < z < 1 . 8, among which about 2000-  \nFig. 35. Galaxy star-formation rate - stellar mass relation at di ff erent redshifts. The Flagship lightcone catalogue (red shaded area) is compared to COSMOS2020 (Weaver et al. 2023, blue contours).  \n<!-- image -->  \nFig. 36. Expected number densities for photometric and spectroscopic emission line flux selected samples as a function of redshift in the Euclid Wide survey. Di ff erent colours show various selection cuts for the EWS (in I E, H E and in the 2 model calibrations for the line fluxes of f H α + [N ii ] > 2 × 10 -16 erg s -1 cm -2 ). The empirical model of Pozzetti et al. (2016) and the data based on slitless HST spectroscopy (Bagley et al. 2020) are also shown.  \n<!-- image -->  \n3800 deg -2 have H α flux brighter than that same flux limit. In the Euclid Deep Survey (EDS), these numbers increase to about 38 000-55 000 deg -2 in the redshift range 0 . 4 < z < 1 . 8 with H α + [N ii ] fluxes above 5 × 10 -17 erg s -1 cm -2 .  \nWeunderline, however, that the measurements of galaxy redshift rely mainly on the identification of emission lines and on their ratios. We show in Fig. 37 how our mock catalogue populates the BPT emission line ratio diagram (Baldwin et al. 1981) and how it compares with observations. In particular, we show that we recover quite well the main locus of the BPT diagram from the SDSS galaxy main sample at low redshift ( z < 0 . 3) and also the trend to higher values shown by high-redshift observation in the Fiber Multi-Object Spectrograph (FMOS)-COSMOS survey (Kashino et al. 2019) and by the MOSFIRE Deep Evolution Field (MOSDEF) survey with Keck by Kriek et al. (2015) and Reddy et al. (2015) in the redshift range 1 . 4 < z < 1 . 7. We also show other diagrams with various emission line ratios (Fig. 38), which confirms that the assignment of emission line fluxes is consistent with those observed in SDSS galaxies at low redshift and the FMOS-COSMOS and MOSDEF surveys at higher redshifts.  \nFinally, the weak lensing cosmological probe relies on the measurement of galaxy shapes and sizes. Furthermore, the spectroscopic signal-to-noise ratio of the Euclid slitless spectra depends on galaxy size and, with it, the probability of detecting and measuring emission lines and galaxy redshifts for the clustering cosmological probe. The procedure to assign galaxy sizes distribution was already described in Sect. 5.5 and its comparison to the calibration data shown in Figs. 23 and 24. Additionally, in Figs. 39 and 40 we show the distribution of disk sizes of H α emitters in the redshift range 0 . 9 < z < 1 . 6 of a sample limited in flux to f H α + [N ii ] > 2 × 10 -16 erg s -1 cm -2 and compare it to the observed distributions from slitless spectroscopic HST Wide Field Camera 3 (WFC3) data sets by Bagley et al. (2020). We note that the distribution is realistic, with pretty similar dispersion, even if the median predicted disk size in the FS2 is around 0 . '' 48, compared to 0 . '' 35 in Bagley et al. (2020). This di ff erence can be explained by a selection e ff ect, as the detection probability is dependent on galaxy size, with more compact galaxies  \n<!-- image -->  \nFig. 37. BPT diagram. Top panel: Emission line ratios in the redshift range z < 0 . 5 for Flagship (colour map) compared to SDSS (in grey) main target sample for objects classified as galaxies. Bottom panel: Emission line ratio in the redshift range 0 . 9 < z < 1 . 8 compared to FMOS stacked data in COSMOS at z ∼ 1 . 5 from Kashino et al. (2019) and MOSFDEF data from Kriek et al. (2015) and Reddy et al. (2015). The black lines are the relations by Kewley et al. (2013a) at z ∼ 0 (dotted lines) and z ∼ 1 . 5 (solid lines).  \n<!-- image -->  \n<!-- image -->  \nFig. 38. Emission lines ratios in the Flagship simulation. Flagship galaxy fluxes (colour map) compared to SDSS fluxes (grey map) from the main target sample for objects classified as galaxies in the local Universe. In blue, we also show FMOS stacked data in COSMOS from Kashino et al. (2019). Top panel: [O iii ] / H β versus [S ii ] / H α . Bottom panel: [O iii ] / H β versus stellar mass; the thin dotted curves indicate the divisions between star-forming / composite galaxies and AGN at z ≃ 0 (Juneau et al. 2014), with AGN populating the region of larger [O iii ] / H β above the dotted lines.  \n<!-- image -->  \nFig. 39. Size distribution. Disk half-light radius distribution of an H α selected sample, at the limit of the EWS in the redshift range 0 . 9 < z < 1 . 6, compared to data from Bagley et al. (2020).  \n<!-- image -->  \nFig. 40. Galaxy half-light radius versus H α flux compared to observed data of Bagley et al. (2020). In red, the data from the Flagship catalogue. In blue, the data from Bagley et al. (2020)  \n<!-- image -->  \nhaving a higher probability of being detected. We performed a preliminary comparison of the predicted size distribution under this selection e ff ect, obtaining a consistent shift in the size distribution. The description of the techniques needed to perform this test is beyond the scope of this paper. The details of the selection procedure and its dependence on galaxy properties will be presented in Euclid Collaboration: Monaco et al. (in prep.).", '6.2. Galaxy clustering': "The clustering properties of the FS2 galaxies were tested at various levels, to check the stability of clustering for objects produced in various steps along the pipeline. To this aim we selected four classes of objects, identified so as to follow the main step of the catalogue construction (the thresholds are chosen to have a similar level of clustering in the first of the redshift bins defined below): (i) dark matter haloes on the lightcone more massive than M h = 10 12 h -1 M ⊙ (the starting point), (ii) galaxies with absolute r -band magnitude brighter than Mr = -20 (after application of HOD + AM, Sect. 5.1), (iii) galaxies with apparent  \nFig. 41. Monopoles (thick continuous lines), quadrupoles (thin continuous lines) and hexadecapoles (dotted lines) of the measured redshiftspace PK for the four samples reported in the legend. The monopole includes shot noise. Thin black lines give the best fit model (convolved with the window function) for the monopole and quadrupole of the H α ELG sample.  \n<!-- image -->  \nFig. 42. Inferred cosmological parameters (redshift-space f and AP parameters) for the four samples, together with linear bias b 1, in the redshift bin [0 . 9 , 1 . 1].  \n<!-- image -->  \nmagnitude H E < 22 (after SED assignment, Sect. 5.4), 21 (iv) emission line galaxies (ELGs) with H α line flux greater than f H α = 2 × 10 -16 erg s -1 cm -2 (after emission lines assignment, Sect. 5.8). For each selected catalogue we consider the four redshift bins of the spectroscopic sample (namely [0 . 9 , 1 . 1),  \nFig. 43. Redshift-space two-point correlation function multipoles of H α galaxies, and best-fitting prediction using the EFT model. Green, red, and orange points refer to monopole, quadrupole and hexadecapole respectively. Di ff erent panels show the results for di ff erent redshifts as labelled.  \n<!-- image -->  \n[1 . 1 , 1 . 3), [1 . 3 , 1 . 5) and [1 . 5 , 1 . 8)). To isolate the e ff ect of peculiar velocities in redshift-space distortions, we select comoving angles (instead of magnified angles, which are shifted according to the gravitational lensing displacement field), observed redshift containing the contribution from background expansion and first-order peculiar velocities and an unmagnified flux cut. This ensures that our measurements will not be a ff ected by lensing e ff ects, which have been shown to be non-negligible for high-redshift 3D galaxy clustering (Jelic-Cizmek et al. 2021; Breton et al. 2022; Euclid Collaboration: Jelic-Cizmek et al. 2023) and whose e ff ect is outside the scope of this paper. For each catalogue, we generate an associated random catalogue constructed by replicating the galaxy catalogue 50 times, keeping the redshift while assigning random sky positions within the angular footprint of the mock survey (an octant of the sky). This way, the average number density of the random reproduces exactly 50 times that of the galaxy catalogue.  \nFirst, we report our validation of the galaxy clustering probe in Fourier space using the o ffi cial code developed to compute the power spectrum ( PK ) in Euclid (Euclid Collaboration: Sefusatti et al., in prep.), based on the Yamamoto-Bianchi esti-  \nFig. 44. Linear galaxy bias, b 1, as function of redshift for H α emitters in the Flagship simulation. In blue we show the estimation from the fit to the two-point correlation function multipoles in redshift space using the EFT model. In orange, we use linear theory in real space to fit the correlation function monopole. Finally, in green we fit the real-space monopole of the power spectrum using EFT at one-loop order.  \n<!-- image -->  \nmator (Bianchi et al. 2015). The code accepts a catalogue with (cosmology-independent) angular coordinates and redshift; we computed the monopole of the real-space PK , obtained by providing the true redshift of the galaxy, and the first three odd multipoles of the redshift-space PK , obtained by providing the 'observed' redshift that includes peculiar velocities. The same PK code performs a measurement of the power spectrum of the randoms, which is used to create a model of the window function.  \nThe measured power spectra were fitted with a model based on a standard perturbation theory (SPT) at 1-loop with e ff ective field theory (EFT) counter-terms using the PBJ code presented by Moretti et al. (2023). Cosmological parameters were fixed to those used in the simulation. This model contains the following parameters for galaxy bias, shot noise and EFT that were used as nuisance parameters: b 1 (linear bias), b 2 (quadratic bias), b G 2 (second order Galileon bias), c 0 and c 2 (the EFT counterterms), α P (deviation from the Poisson shot-noise), ϵ 0 , k 2 ( k 2 dependent shot-noise) and c ∇ 4 δ (Finger-of-God counter-term). In redshift space we left the growth rate f free and added AlcockPaczynski (AP) free parameters α ∥ and α ⊥ to the fiducial power spectrum. Indeed, the aim here is not to perform accurate cosmological inference from the spectroscopic catalogue but to check that we can infer the input cosmological parameters consistently when using the di ff erent catalogues along the mock production pipeline. The model was convolved with the window function following the matrix multiplication method of d'Amico et al. (2020). For the covariance matrix we used a simple analytic Gaussian, leading-order covariance; to minimise the impact of this naive choice and the lack of a proper convolution with the window function, we binned the power spectrum in bins of 8 times the fundamental frequency ( k f ) of the box used for the measurement ( L box = 10 h -1 Gpc), so that 8 k f ∼ V 1 / 3 e ff , where V e ff is the e ff ective survey volume. The model range of validity is thus from k min = 8 k f = 0 . 005 h Mpc -1 to k max = 0 . 2 h Mpc -1 . For the likelihood, we neglected the hexadecapole, which had a negligible e ff ect on the posteriors.  \nFirst, the real-space PK was fitted by the model with fixed cosmological parameters, with the aim of obtaining the b 1 linear  \nbias parameter of the four classes of objects. In the second stage, we fitted the redshift-space PK of the four catalogues, leaving f , α ∥ and α ⊥ free. Figure 41 shows the PK measurement of the monopole, quadrupole and hexadecapole of the four samples in the first redshift bin, [0 . 9 , 1 . 1]; we also report the best-fit model PK of the H α ELGs. Figure 42 shows the resulting cosmological parameters obtained from the four samples, together with linear bias b 1; posteriors are marginalised over the other nuisance parameters. The cosmological posteriors are very stable when varying the sample and are compatible to within 1 σ with the true values, amounting to 1 for the AP parameters and to f = 0 . 9 for the growth rate (computed for the median redshift z = 1); the small bias obtained for some of the parameters raises no concern as we are using a single realisation and a simplified treatment of the covariance. The variation of parameter error bars for the four catalogues is very limited, and this is due to the fact that at k = 0 . 2 h Mpc -1 the power spectrum is still well above the shot noise level in all cases, so the di ff erent number densities of the samples do not influence much the parameter error bars in this regime. Linear bias, b 1, values are compared with those obtained from real space (dashed lines). In this case, we also obtain a very moderate bias in this measurement, which is very consistent for all the samples. In conclusion, the various steps in the definition of the galaxy sample do not bias the inferred cosmological parameters.  \nWe now focus on the multipoles of the two-point correlation function for H α galaxies in redshift space. To estimate the correlation function we use the LS estimator (Landy & Szalay 1993)  \nξ ( s , µ ) = DD( s , µ ) -2 DR( s , µ ) + RR( s , µ ) RR( s , µ ) , (45)  \nwhere DD, DR and RR respectively stand for data-data, datarandom and random-random pair counts that we estimate with Corrfunc 22 (Sinha & Garrison 2020), in bins of s the comoving pair separation (between 25 and 150 h -1 Mpc) and µ the cosine of its angle (using the mid-point definition) with respect to the line of sight. To estimate the multipoles of the correlation function we integrate the 2D correlation function as  \nξℓ ( s ) = (2 ℓ + 1) Z 1 0 ξ ( s , µ ) L ℓ ( µ ) d µ, (46)  \nwhere ℓ is the multipole and L ℓ the ℓ -th order Legendre polynomial. Note that we integrate µ from 0 to 1 due to symmetry along the line of sight when using auto-correlations and the midpoint definition, that is µ = r · s / ( rs ), where r = ( r 1 + r 2) / 2 and s = r 2 -r 1, with r 1 and r 2 the vectors of the galaxy positions.  \nWe compare our measurements to a theoretical prediction based on EFTofLSS (Ivanov et al. 2020; d'Amico et al. 2020) as implemented in COMET (Eggemeier et al. 2023). To obtain the prediction, we first produced Gaussian covariance matrices for each redshift bin following the recipe from Grieb et al. (2016). Once we compute a covariance with a model which fits the data, we run Monte Carlo Markov chains (MCMCs) using PyMultinest (Buchner et al. 2014). Consistent with the power spectrum analysis presented above, we performed a fixedtemplate fit assuming the fiducial cosmology of the simulation and varying f the growth rate, α ∥ and α ⊥ the dilation parameters, the linear and quadratic biases (higher-order biases are fixed  \nassuming the local-lagrangian approximation) and EFT counterterms. In principle, it is possible to vary σ 12 (the variance of density fluctuations in spheres of 12 Mpc, see Sánchez 2020), but we instead fix it to its fiducial value so as to avoid large degeneracies in the estimation of the linear bias.  \nThe simulation measurements and best-fit model shown in Fig. 43 display a very good agreement, especially for the highestredshift bins (with mean redshift z = 1 . 4 and z = 1 . 65). For the low-redshift bins, it seems that the quadrupole is underestimated by the model at scales above 90 h -1 Mpc. This does not look particularly worrying given that these scales are not independent and that the discrepancy stays between 1 and 2 σ . We also note that for z = 1, the monopole is close to zero at 80 h -1 Mpc, which could be attributed to sample fluctuations in the redshift bin for our simulation realisation.  \nFinally, we show a simulation-based estimation of the linear bias as a function of redshift in Fig. 44. Di ff erent models (linear theory or EFT) in di ff erent configurations (real or redshift space, Fourier or configuration space) give very similar results. Because the measurement of b 1 from a galaxy sample is not a trivial task, this result can be considered as a validation not only of the simulation itself but also of our analysis pipeline. Nevertheless, our estimate of the bias in either real or redshift space appears to be higher than others reported in the literature, such as the simulation-based forecasts for future missions in Merson et al. (2019) and Zhai et al. (2021a,b). Those results are derived from semi-analytic models of galaxy formation and evolution, which have been calibrated to match the number density of H α emitters. The methodology employed in these studies di ff ers significantly from the approach used in constructing the Flagship simulation in this work. The extrapolation to clustering measurements may introduce additional uncertainty into the estimated bias, and thus we do not expect great consistency between the analyses. However, the collection of ELG data by Euclid is expected to provide tight constraints on these parameters in the near future.", '6.3.1. Main observables: shear 2-point statistics and galaxy-galaxy lensing': "Galaxy weak lensing is one of the two main cosmological probes of Euclid . The standard summary statistics used for weak lensing include the shear 2-point correlation functions, ξ ± , which in turn are related to the tangential and cross-component of the shear,  \nξ ± ( θ ) = ⟨ γ t( ϑ ) γ t( ϑ ' ) ⟩ ± ⟨ γ × ( ϑ ) γ × ( ϑ ' ) ⟩ , (47)  \nwhere θ is the separation between ϑ and ϑ ' and γ t and γ × are the tangential and cross-component of the shear, defined by γ t( ϑ ) = -ℜ [ γ ( ϑ )e -2i ϕ ] and γ × ( ϑ ) = -ℑ [ γ ( ϑ )e -2i ϕ ], where ϕ is the polar angle of ϑ -ϑ ' .  \nIn the weak lensing limit these shear correlations are related to the gradient or E-mode component of the shear angular power spectrum (i.e, the B-mode component is typically negligible: see Bartelmann & Schneider 2001; Hilbert et al. 2009),  \nξ ± ( θ ) = 1 2 π Z d ℓ ℓ J 0 / 4( ℓθ ) C γγ ℓ , (48)  \nbeing C γγ ℓ = C κκ ℓ , J 0 / 4 are the Bessel functions of the first kind of order 0 and 4 respectively. where we have assumed no B-modes, and the validity of the Limber approximation (e.g., LoVerde &  \nAfshordi 2008). Another basic observable is given by the galaxygalaxy lensing, or the average tangential shear of a background galaxy sample produced by the foreground matter distribution, ⟨ γ t ⟩ , which is directly related to the cross power spectrum of the convergence field of the background galaxies and the foreground galaxy number counts, C κg (see Jeong et al. 2009),  \nℓ  \n⟨ γ t ⟩ ( θ ) = 1 2 π Z d ℓ ℓ J 2( ℓθ ) C κg ℓ , (49)  \nwhere we have taken the flat-sky limit, which is very accurate for practically all angular scales (e.g., θ < few degrees). The exact expression can be obtained by replacing the the integral over ℓ by a sum over discrete ℓ modes, and features Wigner d matrix elements (Chon et al. 2004; de Putter & Takada 2010).  \nIn this section we provide a validation of these basic weaklensing statistics which are, along with galaxy clustering correlations, the two main cosmology probes of Euclid . Figure 45 shows the simulation measurements of the shear correlation functions at source redshifts z = 1 , 2, using the athena code 23 (Kilbinger et al. 2014) compared to theory predictions from Halofit (Takahashi et al. 2020). The mock galaxy sample has been cut at the nominal magnitude limit of the photometric galaxy sample, I E < 24 . 5. The measured correlations are in agreement with theory predictions to within 10% down to highly nonlinear scales, i.e, 1 ' for ξ + (5 ' for ξ -). A similar level of agreement between measurements and theory is found for other redshift slices 0 . 95 < z < 1 . 05 , 1 . 95 < z < 2 . 05, and 2 . 9 < z < 3, as displayed in Fig. 46. We have also compared mock measurements of the galaxy-galaxy lensing, obtained with athena , with theory predictions that use Halofit for the nonlinear matter power spectrum re-scaled with a simple linear galaxy bias model. As shown in Fig. 47 we get good agreement within 1015% for this simple model even on small (nonlinear) scales. Figure 48 shows a comparison of the galaxy-galaxy lensing estimator for di ff erent lens-source z -bin pairs across the simulation redshift range. It shows that overall there is good agreement between the simulation and nonlinear theory predictions, except for the smallest angular scales, where the linear galaxy bias model is expected to break down.", '6.3.2. Magnification bias': "Gravitational lensing by large-scale structures in the Universe changes the number density of background sources and thus it induces a cross-correlation signal between background and foreground galaxy populations (Moessner & Jain 1998; Bartelmann & Schneider 2001). For a magnitude limited survey, the cumulative number of galaxies above a flux limit f scales as N 0( > f ) ∼ Af α , where A is the area of the survey, and α is the power-law slope of the background number counts. Lensing preserves the surface brightness of galaxies by increasing the observed survey depth (i.e, decreasing the e ff ective flux limit) and the e ff ective survey area by the same amount: f → f /µ , A → A /µ , where µ is the magnification. These two competing e ff ects induce the so-called magnification bias in the cumulative number of background sources,  \nN ( > f ) ∼ 1 µ A f µ ! -α = µ α -1 N 0( > f ) . (50)  \n<  \nzs  \n<  \nFigure 49 shows the measured magnification bias b = 5 ∗ s -2, where s is the slope of the background galaxy number counts, see Eq. (51), as a function of magnitude limit in the visible ( I E) and  \n<!-- image -->  \n1.95  \n<  \nzs  \n<  \n2.05  \nFig. 45. Shear correlation functions for source samples at z = 1 (top), z = 2 (bottom). Simulation measurements (symbols) are compared to linear and nonlinear ( Halofit ) theory predictions (lines). Errors are obtained from jackknife resampling. The lower sub-panels show fractional di ff erences between simulations and nonlinear theory.  \n<!-- image -->  \nWe define the logarithmic slope of the background number counts at redshift z , for a magnitude limit m , as  \ns = 2 5 α ≡ d log10 N ( < m , z ) d m . (51)  \nArticle number, page 30 of 45  \nIn the weak-lensing limit, µ = 1 + δµ where | δµ | ≪ 1, and we can Taylor expand, µ α -1 ≈ 1 + ( α -1) δµ , and therefore the magnified overdensity of background sources is given by,  \nδ all = N -N 0 = δ m + δ p  \nN 0 (52) = ( α -1) δµ = (2 . 5 s -1) δµ = (5 s -2) δκ , (53)  \nwhere in the last equality, we have identified δµ = 2 δκ , that is valid in the weak-lensing limit. Note that in δ all we have defined the two qualitatively di ff erent contributions:  \n- 1. magnified magnitudes , δ m = αδµ ,\n- 2. magnified or lensed positions , δ p = -δµ .  \nThese two contributions cannot be separated observationally, but we define two di ff erent galaxy samples accordingly in our simulation in order to validate the two magnification contributions separately.  \nThe net magnification from these two competing e ff ects depends on how the loss of sources due to the area dilution, δ p, is compensated by the gain of sources from the flux magnification, δ m. Number counts for source populations with flat luminosity functions, such as faint galaxies, decrease due to magnification, whereas sources with steep luminosity functions, such as quasars, increase. Note that, in the particular case when s = 0 . 4, then α = 1, and there is no net magnification e ff ect.  \nFor the implementation of the magnification in the flux (or magnitudes) and the galaxy positions, we follow Fosalba et al. (2015) and we refer to that paper for further details. Below we just provide the main definitions.  \n- 1. Magnified magnitudes: flux magnification makes the mock galaxy magnitudes, m , brighter by an amount  \n∆ m = 5 2 log10 µ = 2 . 5 log10(1 -δµ ) ≃ 2 . 174 κ , (54)  \nwhere in the last equality we have Taylor expanded log10(1 -δµ ) and used the fact that δµ ≃ 2 κ in the weak-lensing limit. Therefore, knowing the value of the convergence, κ , at a given point in the source plane, it is straightforward to compute the flux magnification induced, which in turn produces the change in the background number counts, δ m.  \n- 2. Magnified or lensed positions: the 'observed' or lensed position, β , of a light ray is shifted from the 'true' or unlensed position, θ , by an angle given by the scaled deflection vector, α , according to the lens equation on the source plane (see e.g., Bartelmann & Schneider 2001). In the single-plane (or Born) approximation, the lens equation reads  \nθ = β + α , (55)  \nwhere the deflection vector, α is a tangent vector at the unlensed position of the light ray, and the lensed position is found by moving along a geodesic on the sphere in the direction of this tangent vector and for an arc length given by the scaled deflection angle, α . If we denote the unlensed position on the sphere by ( θ, ϕ ), then the lensed position, ( θ ' , ϕ +∆ ϕ ), can be simply derived by using identities of spherical triangles (Lewis 2005) .  \nFig. 46. Shear cross-correlation functions for samples at three di ff erent source bins 0 . 95 < z < 1 . 05 , 1 . 95 < z < 2 . 05, and 2 . 9 < z < 3. Panels show all possible cross-correlation z -bin pairs. We denote a z -bin pair with e.g., (1 , 2) when correlating shear amplitudes from the z ≃ 1 and z ≃ 2 source bins.  \n<!-- image -->  \nTable 3. Fitting functions for the slope of the galaxy number counts vs. magnitude cut, s = s 0 + s 1 I E + s 2 I 2 E + s 3 I 3 E , for di ff erent source redshift bins, z s.  \nof redshift. We find that the smooth evolution of s with the magnitude limit applied can be well fitted with a third-order polynomial, s = s 0 + s 1 I E + s 2 I 2 E + s 3 I 3 E for a given source redshift bin. Table 3 gives the coe ffi cients of the polynomial fit for a set of source redshift bins with width ∆ z = 0 . 1 up to z = 1 . 8.  \nIn order to validate the magnified or deflected positions of mock galaxies, we use the estimator for 'sample variance free' cross-correlations between background and foreground mock galaxy samples (see section 5.2 of Fosalba et al. 2015). As shown in Fig. 50, the measured cross-correlations are in good agreement with theory expectations and above the estimated noise level for most of the angular scales probed.", '6.3.3. Higher-order lensing statistics': 'We also validate the agreement of the third-order aperture statistics ⟨ M 3 ap ⟩ in Flagship with theoretical predictions. The ⟨ M 3 ap ⟩ values are the third moments of the aperture mass and therefore depend on the matter bispectrum. They can be inferred in two ways (Schneider et al. 2005): either from maps of the lensing convergence or by first measuring the third-order shear correlation functions and then convolving them with a suitable kernel function. By comparing the results of these two approaches with theory, we can simultaneously test the skewness of the convergence maps and the third-order statistics of the shear catalogues.  \nFor our model, we use the approach in Heydenreich et al. (2023), based on the bihalofit bispectrum (Takahashi et al.', '0.95 < zd < 1.05, 1.95 < zs < 2.05': "Figure 53 shows the 3D number density profile of halo galaxies in a slice of redshift and halo mass. We restricted the analysis to 199 haloes further than three virial radii from the edges of the studied 49 deg 2 region of Flagship, containing a total of over 10 4 member galaxies (with I E < 24 . 5 and H E < 24) within the virial radius: 6545, 1332, and 2926, Red Sequence, Green Valley, and Blue Cloud, respectively. One clearly sees the different radial distributions for the Red-Sequence, Green-Valley, and Blue-Cloud galaxies, with the latter having a concentration 4 times lower than for the Red-Sequence galaxies, as designed and shown in Table 1. The stacked radial distribution for all galaxies is also well fit by an NFW model even if it is the composite of three NFW models with very di ff erent scale radii. The bestfit concentration for all galaxies is intermediate between the extreme red and blue colour classes. One also sees that the NFW model is very well followed out to ≃ 1 . 8 r vir, beyond which only  \n<!-- image -->  \n<  \nzd  \n<  \n<  \nzs  \n<  \nFig. 47. Galaxy-galaxy lensing: cross-correlations between lens positions and source shear for di ff erent source redshift planes ( z s ≃ 2 in the top panel, and z s ≃ 3 in the bottom panel, respectively) at fixed lens (deflector) redshift, z d ≃ 1. Simulation measurements (filled blue circles) are compared to linear and nonlinear theory predictions (lines). Lower panels show fractional di ff erences between simulation and nonlinear theory.  \n<!-- image -->  \n2020) and the aperture filter function by Crittenden et al. (2002). This model has an expected accuracy of 10% for the aperture radii we consider here (Heydenreich et al. 2023).  \nWe measure ⟨ M 3 ap ⟩ on the convergence map at source redshift z s = 0 . 996. For this, the convergence is convolved with the aperture filter for (4 ' , 8 ' , 16 ' , 32 ' ) using Fast Fourier Transforms to obtain aperture mass maps. Then, a border of width 4 × 32 ' = 128 ' is cut from each aperture mass map to avoid border e ff ects. Finally, the mean of the product of three maps gives ⟨ M 3 ap ⟩ ( θ 1 , θ 2 , θ 3).  \nWe also measure ⟨ M 3 ap ⟩ from the shear of galaxies at redshifts between 0 . 984 and 1 . 01 with a Euclid -like magnitude cut of I E ≤ 24 . 5. For this we measure the shear three-point correlation function with TreeCorr (Jarvis et al. 2004) for galaxy separations between 0 . ' 1 and 400 ' , and convert to ⟨ M 3 ap ⟩ according to Schneider et al. (2005). The correlation function is measured individually for 40 patches, allowing for a jackknife estimate of the variance.  \nFigure 51 shows the measurement results and the model. The measurements for both methods agree with the model within the 10% model accuracy. They also agree with each other within the jackknife uncertainty. This confirms that the third moment of the matter distribution corresponds to theoretical expectations and that the galaxy shear catalogues retain the correct third-order moment.", '6.4. Galaxy clusters': "In this subsection, we concentrate on the properties of galaxies in clusters. For this purpose, we consider as 'clusters' all haloes more massive than 10 14 h -1 M ⊙ . For computational ease, we concentrate on statistics extracted from a 49 deg 2 square region of the Flagship catalogue, which is statistically representative of the whole population for the comparisons we present below.  \nFig 52 displays the number of galaxies as a function of cluster halo mass, computed in three ways: for the EWS limits of I E < 24 . 5 (diamonds) and H E < 24 (circles) and for the H α flux limit of f (H α ) > 2 × 10 -16 erg s -1 cm -2 (triangles). At a given redshift, the number of galaxies per halo increases with halo mass, as expected from the HOD approach implemented in the mock. Also, at given halo mass, the number of galaxies per halo decreases with increasing redshift, because of the broadband or emission line flux limits. Also, the number of galaxies with potentially detectable H α emission lines, in the red grism redshift range 0 . 9 < z < 1 . 8 is of the order of 5 to 20 per halo, amounting to roughly 6% (0 . 9 < z < 1 . 3) and 12% (1 . 3 < z < 1 . 8) of the number of galaxies with H E < 24 at those redshifts. Finally, there are no H α flux-limited selected galaxies in haloes of mass log10[ M h / ( h -1 M ⊙ )] > 14 . 5, as there are no star-forming blue galaxies in those massive clusters (see Sect. 5.2).  \nFig. 48. Average tangential shear for source samples in redshift bins (S1 , S2 , S3) = (0 . 95 < z < 1 . 05 , 1 . 95 < z < 2 . 05 , 2 . 9 < z < 3), around the lens redshift slices (L1 , L2 , L3) = (0 . 45 < z < 0 . 55 , 0 . 95 < z < 1 . 05 , 1 . 45 < z < 1 . 55. Panels show all possible combinations of cross-correlations between lens sample positions and source galaxy shear. Lens sample redshifts increase from top to bottom panels, whereas source redshift increases from left to right. Solid lines show a fiducial theory prediction for the nonlinear dark-matter power spectrum re-scaled with a best-fit linear galaxy bias for the lens samples, b (L1) = 1 . 2 , b (L2) = 1 . 45 , b (L3) = 2 . 0  \n<!-- image -->  \na fraction of galaxies are halo members because of the triaxiality of the haloes (Sect. 5.2).  \nFigure 54 compares the Flagship line-of-sight velocity dispersion as a function of halo mass with that measured in a 7 deg 2 lightcone extracted from the GAEA semi-analytical model (Hirschmann et al. 2016). The agreement between the Flagship and the semi-analytical model is excellent. This indicates that the velocities of Flagship, which follow analytical prescriptions (Sect. 5.3) are fairly realistic, because the galaxy velocities in GAEA are those of the subhaloes as extracted from the dark matter-only simulation on which the GAEA semi-analytical model was extracted. This agreement is robust to the di ff erent virial definitions (the ROCKSTAR bound halo virial mass following the Bryan & Norman 1998 definition for Flagship vs. M 200c for GAEA ). Note that the 2 3 v vir and 2 3 v 200c lines (which by definition both have slopes of 2 / 3) only serve as a reference and are not theoretical expectations.  \nFigure 55 compares the I E-band cluster luminosity functions (LFs) of Flagship in three redshift bins with the Schechter (1976) function fits performed by Sarron et al. (2018) on clusters us-  \ning i -band galaxy luminosities from the Canada France Hawaii Legacy Survey (CFHTLS). The Flagship LF reproduces reasonably well the observed one, given the intrinsic di ff erences in the LF modelling technique. In particular, the numbers of faint galaxies (absolute magnitude between -19 and -21) in Flagship are similar to those predicted by the LF fits to the CFHTLS data. However, one sees several di ff erences. First, Flagship predicts more luminous galaxies in clusters than observed in CFHTLS, regardless of colour. Second, the faint-end slopes for the red (respectively blue) galaxies is much (respectively somewhat) shallower than observed in CFHTLS. Third, the break of the red galaxy LF occurs at somewhat lower luminosities for Flagship than for CFHTLS, while it appears similar for the blue LFs of both Flagship and CFHTLS.  \nFigure 56 shows the Flagship colour-magnitude relations for field and cluster galaxies, using proxies for rest-frame U -B colours and rest-frame B magnitudes at selected redshifts where the 4000 Å break lies in between the two considered wavebands for the colour, so as to emphasize the separation between RedSequence and Blue-Cloud galaxies. One sees that Flagship re-  \nFig. 49. Magnification bias for magnitude limited source galaxy samples with redshift bin-width d z = 0 . 1.  \n<!-- image -->  \nroduces nicely the bimodality of colours at z = 0 . 1, although the Red Sequence is barely visible for the field galaxies, because the Blue-Cloud galaxies are so dominant. This trend is also visible at higher redshifts, showing that cluster galaxies in Flagship are more likely to be luminous and red at all redshifts. Again, the Red Sequence is present by construction in all galaxies (blue contours), but is overshadowed by the dominant Blue Cloud.", '7. Summary and conclusions': "We have presented the Flagship galaxy mock, a large simulated catalogue especially designed for the Euclid mission. Euclid will use weak lensing and galaxy clustering as its main observables to infer cosmological parameters. The Flagship catalogue simulates these observables in a self-consistent manner and covers the volume and depth required for Euclid . We also compute many other properties for each galaxy to increase the science cases that can be addressed with the catalogue.  \nWe ran a 4 trillion particle dark matter N -body simulations in a box of 3600 h -1 Mpc length side at the Swiss National Supercomputing Centre. We generated a lightcone on the fly out to redshift z = 3 in one octant of the sky containing 31 trillion particle positions and velocities. We have checked that the dark matter clustering properties behave as expected from theoretical predictions.  \nWe also produce an all-sky lensing map in HEALPix format with N side = 8192 resolution, corresponding to 0 ' . 43 per pixel. From these maps, we can compute the convergence, shear and displacement at any position in the simulation lightcone. The lensing statistics are in general agreement with expectations from theoretical predictions.  \nWe produce a halo catalogue with the ROCKSTAR halo finder of 15.8 billion main haloes in one octant of the lightcone. In order to push the completeness of the catalogue to faint magnitudes as to be complete for the Euclid magnitude selection limit at low redshift, we select haloes down to a threshold of 10 particles, before discarding unbound particles. After the unbinding process, some of our haloes end up fewer than 10 particles for the M bound mass definition that we use. While many of these fewparticle structures may not be individual virialised haloes, their statistical clustering strength will still be as expected, owing to  \n<!-- image -->  \n<  \nzd  \n<  \n<  \nzs  \n<  \n<  \nzd  \n<  \n<  \nzs  \n<  \nFig. 50. Lens-source cross-correlations as a test of source galaxy position deflections. We show two di ff erent lens-source z -bin pairs for galaxy samples cut at I E < 26: a case with z d ≃ 1 and z s ≃ 2 in the top panel, z d ≃ 1 . 5 and z s ≃ 3 in the bottom panel, respectively. Flagship measurements (filled blue circles) are compared to linear and nonlinear theory predictions (lines). Lower panels show the ratio of FS2 measurements over nonlinear theory prediction. Sample variance errors around nonlinear theory are displayed as dotted lines. Cross-correlation in absence of magnification (i.e, noise estimate) is shown as open red circles. Note that FS2 measurements use an estimator that roughly cancels sample-variance (see text for details).  \n<!-- image -->  \nthe uniform clustering strength as a function of mass at this low mass range. We compute the halo mass function for the di ff erent mass estimates computed with ROCKSTAR . At low redshift, the HMF roughly agrees with the T08, D16 and C17 HMFs for the same halo mass definition. However, there is a small di ff erence in the HMF slope at higher redshift z ≃ 1 . 5. We adopt the virial bound mass definition as our fiducial mass for the catalogue. We nevertheless reassign the halo mass values to correct for completeness and discreteness at low halo masses. We compute the  \nFig. 51. Top : Third-order aperture mass statistic, measured from the convergence map at z s = 0 . 996 (blue, dashed), measured from the shear of galaxies at 0 . 984 < z s < 1 . 01 (orange points), and modelled (black, solid). Bottom : Relative di ff erence of measurements to theory.  \n<!-- image -->  \nFig. 52. Arithmetic means of the number of galaxies per halo as a function of halo mass. The colours indicate the redshifts, while the symbols indicate the selection: I E < 24 . 5 (diamonds), H E < 24 (circles), and H α flux > 2 × 10 -16 erg s -1 cm -2 (triangles, limited to the redshift range, 0 . 9 < z < 1 . 8, where the red grism can see the H α line). We requested at least 2 clusters per cell of redshift and log halo mass. The error bars show the standard deviations of the numbers of galaxies per halo. The abscissae are slightly o ff set for clarity.  \n<!-- image -->  \nclustering of the haloes in the lightcone and compare them to theoretical expectations, finding reasonable agreement.  \nThe next step in our mock production pipeline is to generate a galaxy mock catalogue from the halo catalogue using a combination of halo occupation distribution and abundance matching techniques. We complement these with observed correlations to generate recipes to assign other properties to each galaxy.  \nWe start by generating the galaxy luminosities. First, we compute how many galaxies there are as a function of halo mass  \nFig. 53. 3D galaxy number density profiles of stacks of Flagship haloes, after removing all haloes lying closer than 3 r vir from the edges of the region. The maximum likelihood concentration parameters (virial radius over NFW scale radius), obtained by fitting an NFW model to the distribution of radial distances r (restricted to the shaded region), are shown in the legend.  \n<!-- image -->  \nthreshold. We use the halo mass function and the halo occupation distribution for the calculation of what we call the cumulative galaxy function (Eq. 8). For the HOD, we assume a simple parameterisation given in Eqs. (6) and (7). We compare the CGF to the cumulative luminosity function to establish the relation between halo mass and luminosity. In our luminosity assignment procedure, we apply scatter to the luminosities resulting from the AM relation. Therefore, we compute the unscattered luminosity function, which is the one representing the luminosities that, after being scattered, produce the observed luminosity function (Eqs. 9 and 10). We generate the relation between halo mass and luminosity for values of these quantities that have the same abundance (Eq. 12). We assume that this redshift-dependent relation (Fig. 10) is applicable to central galaxies and generate their luminosities with it. To assign the satellite galaxies' luminosities, we compute the global satellite luminosity function by subtracting the central LF from the LF for all galaxies (Fig. 11). We assume that all haloes share the same functional form of the satellite LF within each halo. We fit the parameters of this universal LF relation to make the sum of the individual LF of all haloes coincide with the global satellite LF. We randomly draw the satellite luminosities from these individual halo LFs that depend on the luminosity of the central galaxy.  \nWenext assign a galaxy colour. We divide the galaxy population into three colour types: red, green and blue. We fit the local SDSS colour-magnitude diagram with three Gaussian distributions, one for each colour type. We define functions determining the fraction of red and green satellites as a function of luminosity. The fraction of blue satellites and red, green and blue centrals is then determined by the CMD and the HOD. We assign a ( g 01 -r 01)HOD colour by randomly sampling the distributions.  \nWe place central galaxies at the central position of each halo. Satellites are distributed with the same triaxial NFW profile as the dark matter. We implement colour segregation, changing the concentration index of each population colour type. Red galaxies  \nFig. 54. Halo line-of-sight velocity dispersion vs. halo mass, for lightcones from the GAEA semi-analytical model (without orphan galaxies, blue circles , 7 deg 2 ) and for Flagship ( red stars , here for the 49 deg 2 region), at three representative redshifts. The lines are simple nontheoretical representations. Note that while Flagship uses virial quantities, GAEA uses R 200c for the virial radius and therefore the x axis is log10[ M 200c / ( h -1 M ⊙ )] for GAEA . Only haloes with at least 10 member galaxies with I E < 24 . 5 and H E < 24 within the virial cylinder are considered, with at least 10 such haloes per mass bin. The velocity dispersions are measured with the gapper method (Wainer & Thissen 1976), and the error bars assume Gaussian velocity distributions, i.e., ϵ ( σ ) = σ/ √ 2 ( N -1), where N is the number of galaxies per mass bin.  \n<!-- image -->  \nare assigned NFW profiles with the same concentration index as the DM. Green and blue galaxies are distributed using concentration indices that are 1 / 2 and 1 / 4 of the concentration of the DM.  \nCentral galaxies are assumed to be at rest at the centre of the halo. The velocities of the satellite galaxies are assigned by solving the Jeans equation of local dynamical equilibrium with anisotropy parameters coming from observations of local clusters and di ff erentiated according to the colour type.  \nWe choose the spectral library used by Ilbert et al. (2009) as our basis for spectral energy distribution assignment. We construct a sample of 136 templates using di ff erent extinction laws and values from the original 31 templates in the library. We rankorder them according to their ( g 01 -r 01)HOD colour. We compute the g -r colour distributions of the SDSS at low redshift and of the COSMOS catalogue at z ≥ 0 . 5. For each galaxy, we choose the template that has the same percentile in the colour distribution in SDSS or COSMOS catalogues as the ( g 01 -r 01)HOD colour already assigned. We assign a probability to the six clos-  \nFig. 55. Cluster luminosity functions of red and blue galaxies for Flagship (points) and Schechter (1976) function fits by Sarron et al. (2018) to the CFHTLS absolute magnitudes (lines). For Flagship ( I E band), the red and blue galaxies are defined by color\\_kind , while for CFHTLS ( i band) they were defined according to g -i vs. i colour-magnitude diagrams, which indicate a colour bimodality at these three redshifts. The normalisation in Flagship is per cluster, while for CFHTLS it is per cylinder of 1 h -1 Mpc radius. The abscissas of the points are slightly o ff set for clarity. The discrepancies between the CFHTLS analysis and Flagship are discussed in the text.  \n<!-- image -->  \nt templates in the rank-ordered list to this one based on their colour di ff erence. We also assign a probability depending on the extinction value of these six templates compared to the distribution of extinction values in the COSMOS galaxies. We finally choose a template randomly sampling from the probability resulting from the product of these two colour and extinction probabilities. We do this template assignment at two redshifts close to the galaxy redshift. The final SED is an interpolation between these two templates with weights based on their redshift distance.  \nWe assign the values of the convergence and shear at the galaxy positions from the values of the HEALPix maps of these quantities. We also compute a displaced position of the galaxy based on the lensing displacement field.  \nWe assign the shapes and sizes of the galaxies using measured distributions in observed HST fields. We use the GOODS South field and the CANDELS observations as calibrators. We model galaxies either as one component (bulge) or two components (bulge and disk). The bulge component is simulated as a Sérsic profile. We use a simple cut in the CMD to divide these two options. For all the galaxies, we compute the scale height and the fraction of light coming from the bulge component. For  \nFig. 56. Contours of the colour-magnitude diagrams of Flagship field galaxies (blue) and cluster ( M vir > 10 14 h -1 M ⊙ ) galaxies (red), at selected redshifts. The figure uses proxies for the rest-frame U -B and rest-frame B magnitude. The galaxies are limited to I E < 24 . 5 and H E < 24. The redshift bins have widths of ± 0 . 05 , 0 . 06 , 0 . 07 , 0 . 1 , 0 . 2 and 0.2. The contours are spaced by 0.2 dex and are KDE-smoothed. The figure highlights how the Red Sequence is more prominent in clusters.  \n<!-- image -->  \neach of the two galaxy components (or just for the bulge), we compute the half-light radius, the Sérsic index, the inclination angle, the ellipticity and the axis ratio, mimicking the distributions of our calibration samples.  \nWe compute the SFR of the galaxies from the rest-frame UV luminosity of the SEDs. We compute the stellar mass from the galaxy luminosity and the stellar mass-to-luminosity ratio. We assign a metallicity from the stellar mass, the SFR and the redshift of the galaxy.  \nWe assign the flux of the H α line based on the value of the SFR. We recalibrate the fluxes of the H α as to agree with the Pozzetti et al. (2016) models. The other Balmer lines are computed assuming the hydrogen lines case B recombination ratios. The other usually most prominent lines are assigned based on observed correlations.  \nOverall, we compute 199 quantities for each galaxy. These properties include the galaxy positions and velocities, the galaxy fluxes in several bands with and without extinction, the lensing properties at the galaxy positions, the galaxy shapes and sizes, the SED, the SFR, the stellar mass and the metallicity, the emission lines fluxes calibrated to both model 1 and model 3 of Pozzetti et al. (2016) with and without extinction, a photo-  \nmetric redshift estimate, the intrinsic alignment properties and some of the halo properties they belong to.  \nWe validate the catalogue by comparing it to observations. We compute the catalogue number counts. In the optical bands, the Flagship number counts coincide within the envelope spanned by observations. In the near-infrared, however, our number counts are slightly larger than observations. We check the colour-colour and colour redshift distributions against the COSMOS2020 catalogue. Overall, there is good agreement between the two, with the exception of a small deficit of blue galaxies in colours including a near-infrared band. We also compute the stellar mass function as a function of redshift. We find relatively good agreement with the determination from the GAMA survey and the COSMOS field. We compare the emission line fluxes in various diagnostic plots and find them to be consistent with observational data.  \nWe also validate the two main Euclid cosmology observational probes. We perform a basic validation of the clustering properties of a few samples using common selection criteria, cutting our sample in halo mass, absolute and apparent magnitude, and emission line flux. We compare our results to theoretical expectations, considering both the power spectrum and multipoles of the two-point correlation function in configuration space. We compute the expected linear bias for an emission line flux-selected sample similar to the one expected in Euclid. We also check that we can recover the input cosmological parameter from this clustering analysis. Overall, all the computed clustering statistics behave as we would expect from theoretical expectations, validating the use of the catalogue for Euclid analysis.  \nWe also validate the lensing properties of the catalogue. We compute the two-point shear correlation function and the average tangential shear and compare them to theoretical predictions. Both show good agreement. We also check the lensing magnification properties of the sample against models and found a good agreement. Finally, we also compare the higher-order lensing statistics to models corroborating that the catalogue shows the correct higher-order behaviour.  \nWecheck how some galaxy cluster properties are reproduced in the Flagship galaxy catalogue. We show the number of galaxies in clusters as a function of halo mass and redshift for magnitude and emission line flux-limited samples. We check the radial profile of galaxies within clusters as a function of their colour type, showing an NFW distribution for a stack of clusters. We present the relation between halo mass and the line of sight velocity dispersion at di ff erent redshifts and compare them to the GAEA semi-analytical model, finding good agreement. We show the LF of galaxies inside clusters for red and blue galaxies at different redshifts and compare them to data from the CFHTLS. We also show the colour-magnitude diagram for cluster galaxies and compare it to that of the overall population at several redshifts. The cluster red sequence is apparent with respect to the rest of the galaxy population.", '7.1. Future developments': "Modelling the properties of a massive and complex galaxy survey like Euclid necessarily entails making certain assumptions and approximations with respect to the true galaxy samples. Below we discuss some of the possible shortcomings in the approach we have taken in the modelling, how they can impact the di ff erent observables and, in some cases, how these could be improved in future releases of the Flagship galaxy mock.  \nIn this paper we have implemented a simple HOD approach to assign galaxies to dark-matter haloes, in which the number  \nof satellite galaxies in each halo depend only on the halo mass. Several studies have shown that other halo properties are needed to accurately simulate the galaxy population within this type of framework. This is commonly referred to as halo assembly bias (e.g., Contreras et al. 2023). While our HOD prescription is simple, the way we assign galaxy properties is more complex than the approach of typical HOD and AM models, and in a certain sense our method resembles models that introduce assembly bias as our galaxy properties assignment depends on other parameters. For example, the way we distribute galaxies within haloes depends on galaxy colour which can be related to the assembly history (e.g., Hearin & Watson 2013). Moreover, normally observations of a particular galaxy sample, selected with a flux limit or colour criteria as for example luminous red galaxies, are interpreted in the HOD framework fitting parameters of the occupancy distribution. In the Flagship catalogue, the HOD is assigned up front and corresponds to the whole population. Any galaxy subsample selected with colour and / or flux cuts will have a di ff erent HOD than the original assignment. Similarly, we have not included the e ff ects of gas physics, which are expected to significantly alter the distribution of dark matter and thus of galaxies within galaxy cluster scales (e.g., Schneider & Teyssier 2015; Gebhardt et al. 2024). Again, our galaxy properties assignment depends on luminosity and colour in order to fit observations, and therefore it must incorporate the e ff ect that baryons have on the galaxy distributions when compared to dark-matter only simulations. Nevertheless, we are exploring HOD prescriptions that depend on other parameters for implementation in future versions of the catalogue.  \nAnother improvement in the catalogue that we are considering is to increase the resolution of the HEALPix lensing maps. Currently, we use a pixel scale of 0 . ' 43, i.e., HEALPix N side = 8192, without interpolation between pixels. This resolution is enough to measure adequately the lensing correlation functions down to 1 arcmin (e.g., for ξ + and the tangential shear), which coincide with the smallest scales probed by Euclid . Therefore, the resolution of the lensing properties is good enough for the main purpose of the simulation. Moreover, in the MICE-Grand Challenge simulation (Fosalba et al. 2015), which used the same assignment scheme, we tested that interpolating linearly between the HEALPix pixels only improved marginally the e ff ective resolution on sub-arcmin scales, and at a high computational cost.  \nWe are also exploring how to improve our SED assignment. Currently, we use linear interpolation between SEDs selected from the COSMOS SED template set (see Sect. 5.4). We would like to use a wider set of templates that can cover better the range of observed SEDs as a function of redshift. Another shortcoming of our simulation is the lack of AGN. We are working on methods to include them in future versions (e.g., Allevato et al. 2021).  \nThe calibration of our recipes to assign galaxy properties are normally restricted to low redshift observed samples, z ≲ 1, which we extrapolate to higher redshifts. We are working on extending the calibration to higher redshifts using more observational data.  \nThe assignment of galaxy positions and velocities is simplistic, relying on analytic profiles without substructure. This is mainly driven by our push to simulate faint galaxies while covering a large volume. While for the relevant scales to be explored with the Euclid main cosmological probes this level of detail is enough, we are exploring keeping more information from the original N -body simulation to be able to produce more realistic galaxy profiles within haloes and therefore support more extensive cluster science analyses.  \nWe have assigned shapes and sizes to galaxies based on distributions measured with HST. The distributions of each parameter are well reproduced by our catalogue. However, we have not enforced the correlations between di ff erent morphological parameters in part due to the usage of two di ff erent samples for our calibration. We are working on a multi-parameter calibration for future implementations, that will be based on the Euclid data themselves which are going to provide a large and homogeneously measured morphological sample.", '7.2. Data availability': 'Given the comprehensive set of validation tests that we have performed on the Flagship galaxy mock catalogue, we believe it is a valuable resource to perform a wide suite of astrophysical analyses beyond its original design goals of supporting the Euclid weak lensing and galaxy clustering analyses. We make the catalogue publicly available so that everybody can potentially benefit from its usage. We distribute this catalogue from the CosmoHub platform.  \nAcknowledgements. FJC and PF acknowledge support form the Spanish Ministerio de Ciencia, Innovación y Universidades, projects PID2019-11317GB, PID2022-141079NB, PID2022-138896NB; the European Research Executive Agency HORIZON-MSCA-2021-SE-01 Research and Innovation programme under the Marie Skłodowska-Curie grant agreement number 101086388 (LACEGAL) and the programme Unidad de Excelencia María de Maeztu, project CEX2020-001058-M. We acknowledge the support of the PRACE project (Call 17) "Simulating the Euclid Universe" for the computer time on the Piz Daint supercomputer at CSCS, Lugano, Switzerland, which made possible the worldleading Flagship 2 N-body simulation. The Port d\'Informació Científica (PIC) is maintained through a collaboration between CIEMAT and IFAE. This work has been supported by MCIN / AEI grant EQC2021-007479-P and European Union NextGenerationEU grant PRTR-C17.I1 and by Generalitat de Catalunya. This work has been partially funded by Premiale MITIC 2015. The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid , in particular the Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft, funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d\'Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the State Secretariat for Education, Research, and Innovation (SERI) at the Swiss Space O ffi ce (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site ( http://www.euclid-ec.org ).', 'References': "- Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2023, Phys. Rev. D, 107, 083504\n- Aihara, H., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S8\n- Akeson, R., Armus, L., Bachelet, E., et al. 2019, arXiv e-prints, arXiv:1902.05569\n- Albrecht, A., Bernstein, G., Cahn, R., et al. 2006, arXiv / astro-ph, 0609591\n- Allevato, V., Shankar, F., Marsden, C., et al. 2021, ApJ, 916, 34\n- Amendola, L., Appleby, S., Avgoustidis, A., et al. 2018, Living Reviews in Relativity, 21, 2\n- Angulo, R. E. & Hahn, O. 2022, Living Reviews in Computational Astrophysics, 8, 1\n- Angulo, R. E., Zennaro, M., Contreras, S., et al. 2021, MNRAS, 507, 5869 Bagley, M. B., Scarlata, C., Mehta, V., et al. 2020, ApJ, 897, 98\n- Balaguera-Antolínez, A., Kitaura, F.-S., Alam, S., et al. 2023, A&A, 673, A130 Baldry, I. K., Balogh, M. L., Bower, R., Glazebrook, K., & Nichol, R. C. 2004, in American Institute of Physics Conference Series, Vol. 743, The New Cosmology: Conference on Strings and Cosmology, ed. R. E. Allen, D. V. Nanopoulos, & C. N. Pope, 106-119\n- Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5  \nBartelmann, M. & Schneider, P. 2001, Phys. Rep., 340, 291  \n- Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143\n- Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717, 379\n- Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109\n- Benson, A. J., Cole, S., Frenk, C. S., Baugh, C. M., & Lacey, C. G. 2000, MN-  \nRAS, 311, 793  \nBerlind, A. A. & Weinberg, D. H. 2002, ApJ, 575, 587  \n- Bertschinger, E. & Gelb, J. M. 1991, Computers in Physics, 5, 164\n- Bianchi, D., Gil-Marín, H., Ruggeri, R., & Percival, W. J. 2015, MNRAS, 453, L11\n- Blanton, M. R., Hogg, D. W., Bahcall, N. A., et al. 2003, ApJ, 592, 819\n- Blanton, M. R., Lupton, R. H., Schlegel, D. J., et al. 2005a, ApJ, 631, 208\n- Blanton, M. R., Schlegel, D. J., Strauss, M. A., et al. 2005b, AJ, 129, 2562\n- Breton, M.-A., de la Torre, S., & Piat, J. 2022, A&A, 661, A154\n- Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000\n- Bryan, G. L. & Norman, M. L. 1998, ApJ, 495, 80\n- Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125  \nBullock, J. S., Wechsler, R. H., & Somerville, R. S. 2002, MNRAS, 329, 246  \n- Butcher, H. & Oemler, A., J. 1978, ApJ, 219, 18\n- Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682\n- Carretero, J., Castander, F. J., Gaztañaga, E., Crocce, M., & Fosalba, P. 2015, MNRAS, 447, 646\n- Carretero, J., Tallada, P., Casals, J., et al. 2017, in Proceedings of the European Physical Society Conference on High Energy Physics. 5-12 July, 488 Chabrier, G. 2003, PASP, 115, 763\n- Chon, G., Challinor, A., Prunet, S., Hivon, E., & Szapudi, I. 2004, MNRAS, 350, 914\n- Cole, S., Lacey, C. G., Baugh, C. M., & Frenk, C. S. 2000, MNRAS, 319, 168 Collister, A. A. & Lahav, O. 2005, MNRAS, 361, 415\n- Comparat, J., Prada, F., Yepes, G., & Klypin, A. 2017, MNRAS, 469, 4157 Contreras, S., Angulo, R. E., Chaves-Montero, J., White, S. D. M., & Aricò, G. 2023, MNRAS, 520, 489\n- Cooray, A. & Sheth, R. 2002, Phys. Rep., 372, 1  \nCouchman, H. M. P., Thomas, P. A., & Pearce, F. R. 1995, ApJ, 452, 797 Crittenden, R. G., Natarajan, P., Pen, U.-L., & Theuns, T. 2002, ApJ, 568, 20 Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2020, MNRAS, 491, 944 Dahlen, T., Mobasher, B., Somerville, R. S., et al. 2005, ApJ, 631, 126 Dakin, J., Hannestad, S., & Tram, T. 2022, MNRAS, 513, 991  \n- d'Amico, G., Gleyzes, J., Kokron, N., et al. 2020, J. Cosmology Astropart. Phys., 5, 005\n- De Lucia, G., Fontanot, F., Xie, L., & Hirschmann, M. 2024, arXiv e-prints, arXiv:2401.06211\n- de Putter, R. & Takada, M. 2010, Phys. Rev. D, 82, 103522\n- de Salas, P. F., Forero, D. V., Ternes, C. A., Tórtola, M., & Valle, J. W. F. 2018, Physics Letters B, 782, 633\n- Despali, G., Giocoli, C., Angulo, R. E., et al. 2016, MNRAS, 456, 2486\n- Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168\n- Diemer, B. 2018, ApJS, 239, 35\n- Diemer, B. & Joyce, M. 2019, ApJ, 871, 168\n- Dimauro, P., Huertas-Company, M., Daddi, E., et al. 2018, MNRAS, 478, 5410 Dolag, K., Komatsu, E., & Sunyaev, R. 2016, MNRAS, 463, 1797  \nDong-Páez, C. A., Smith, A., Szewciw, A. O., et al. 2024, MNRAS, 528, 7236  \nDriver, S. P., Andrews, S. K., Davies, L. J., et al. 2016, ApJ, 827, 108  \nDriver, S. P., Bellstedt, S., Robotham, A. S. G., et al. 2022, MNRAS, 513, 439  \nDubois, Y., Pichon, C., Welker, C., et al. 2014, MNRAS, 444, 1453  \n- Eggemeier, A., Camacho-Quevedo, B., Pezzotta, A., et al. 2023, MNRAS, 519, 2962\n- Ereza, J., Prada, F., Klypin, A., et al. 2023, arXiv e-prints, arXiv:2311.14456 Euclid Collaboration: Blanchard, A., Camera, S., Carbone, C., et al. 2020, A&A, 642, A191\n- Euclid Collaboration: Cropper, M., Al Bahlawan, A., Amiaux, J., et al. 2024, A&A, this issue\n- Euclid Collaboration: Jahnke, K., Gillard, W., Schirmer, M., et al. 2024, A&A, this issue\n- Euclid Collaboration: Jelic-Cizmek, G., Sorrenti, F., Lepori, F., et al. 2023, arXiv e-prints, arXiv:2311.03168\n- Euclid Collaboration: Knabenhans, M., Stadel, J., Potter, D., et al. 2021, MNRAS, 505, 2840\n- Euclid Collaboration: Mellier, Y., Abdurro'uf, Acevedo Barroso, J., Achúcarro, A., et al. 2024, A&A, this issue\n- Euclid Collaboration: Scaramella, R., Amiaux, J., Mellier, Y., et al. 2022, A&A, 662, A112\n- Euclid Collaboration: Serrano, S., Hudelot, P., Seidel, G., et al. 2024, arXiv eprints, arXiv:2401.01452\n- Fidler, C., Kleinjohann, A., Tram, T., Rampf, C., & Koyama, K. 2019, J. Cosmology Astropart. Phys., 01, 025\n- Fosalba, P., Gaztañaga, E., Castander, F. J., & Crocce, M. 2015, MNRAS, 447, 1319\n- Fosalba, P., Gaztañaga, E., Castander, F. J., & Manera, M. 2008, MNRAS, 391, 435\n- Fosalba, P. & Szapudi, I. 2004, ApJ, 617, L95\n- Garrison, L. H., Eisenstein, D. J., Ferrer, D., Maksimova, N. A., & Pinto, P. A. 2021, MNRAS, 508, 575  \nGarrison, L. H., Eisenstein, D. J., & Pinto, P. A. 2019, MNRAS, 485, 3370 Gebhardt, M., Anglés-Alcázar, D., Borrow, J., et al. 2024, MNRAS, 529, 4896 Giavalisco, M., Ferguson, H. C., Koekemoer, A. M., et al. 2004, ApJ, 600, L93 Ginzburg, D., Desjacques, V., & Chan, K. C. 2017, Phys. Rev. D, 96, 083528 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759  \n- Grieb, J. N., Sánchez, A. G., Salazar-Albornoz, S., & Dalla Vecchia, C. 2016, MNRAS, 457, 1577  \nGrogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35 Gu, Y., Yang, X., Han, J., et al. 2024, MNRAS, 529, 4015 Habib, S., Pope, A., Finkel, H., et al. 2016, New A, 42, 49 Harnois-Déraps, J., Pen, U.-L., Iliev, I. T., et al. 2013, MNRAS, 436, 540 Hatton, S., Devriendt, J. E. G., Ninin, S., et al. 2003, MNRAS, 343, 75 Hearin, A. P. & Watson, D. F. 2013, MNRAS, 435, 1313 Heydenreich, S., Linke, L., Burger, P., & Schneider, P. 2023, A&A, 672, A44 Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140 Hilbert, S., Barreira, A., Fabbian, G., et al. 2020, MNRAS, 493, 305 Hilbert, S., Hartlap, J., White, S. D. M., & Schneider, P. 2009, A&A, 499, 31 Hirschmann, M., De Lucia, G., & Fontanot, F. 2016, MNRAS, 461, 1760 Hu, W. 2000, Phys. Rev. D, 62, 043007 Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841 Ilbert, O., Capak, P., Salvato, M., et al. 2009, ApJ, 690, 1236 Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55 Ishiyama, T., Prada, F., Klypin, A. A., et al. 2021, MNRAS, 506, 4210 Ivanov, M. M., Simonovi'c, M., & Zaldarriaga, M. 2020, J. Cosmology Astropart. Phys., 05, 042  \nIvezi'c, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111  \nJarvis, M., Bernstein, G., & Jain, B. 2004, MNRAS, 352, 338  \n- Jelic-Cizmek, G., Lepori, F., Bonvin, C., & Durrer, R. 2021, J. Cosmology Astropart. Phys., 4, 055  \nJeong, D., Komatsu, E., & Jain, B. 2009, Phys. Rev. D, 80, 123527 Jing, Y. P., Mo, H. J., & Boerner, G. 1998, ApJ, 494, 1 Johnston, H., Georgiou, C., Joachimi, B., et al. 2019, A&A, 624, A30 Juneau, S., Bournaud, F., Charlot, S., et al. 2014, ApJ, 788, 88 Kashino, D., Silverman, J. D., Sanders, D., et al. 2019, ApJS, 241, 10  \n- Kau ff mann, G., Colberg, J. M., Diaferio, A., & White, S. D. M. 1999, MNRAS, 303, 188\n- Kau ff mann, G., White, S. D. M., & Guiderdoni, B. 1993, MNRAS, 264, 201 Kennicutt, Jr., R. C. 1998, ARA&A, 36, 189\n- Kewley, L. J., Dopita, M. A., Leitherer, C., et al. 2013a, ApJ, 774, 100\n- Kewley, L. J., Maier, C., Yabe, K., et al. 2013b, ApJ, 774, L10\n- Kilbinger, M., Bonnett, C., & Coupon, J. 2014, athena: Tree code for secondorder correlation functions, Astrophysics Source Code Library, record ascl:1402.026\n- Klypin, A., Gottlöber, S., Kravtsov, A. V., & Khokhlov, A. M. 1999, ApJ, 516, 530\n- Klypin, A. & Prada, F. 2018, MNRAS, 478, 4602\n- Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36 Kovacs, E., Mao, Y.-Y., Aguena, M., et al. 2022, The Open Journal of Astrophysics, 5, 1  \nKravtsov, A. V., Berlind, A. A., Wechsler, R. H., et al. 2004, ApJ, 609, 35 Kriek, M., Shapley, A. E., Reddy, N. A., et al. 2015, ApJS, 218, 15 Lagos, C. d. P., Tobar, R. J., Robotham, A. S. G., et al. 2018, MNRAS, 481, 3573 Landy, S. D. & Szalay, A. S. 1993, ApJ, 412, 64 Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193 Lesgourgues, J. 2011, arXiv e-prints, arXiv:1104.2932 Lewis, A. 2005, Phys. Rev. D, 71, 083008 Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473 LoVerde, M. & Afshordi, N. 2008, Phys. Rev. D, 78, 123506 LSST Dark Energy Science Collaboration (LSST DESC), Abolfathi, B., Alonso,", 'D., et al. 2021, ApJS, 253, 31': '- Mamon, G. A., Biviano, A., & Boué, G. 2013, MNRAS, 429, 3079 Mamon, G. A., Cava, A., Biviano, A., et al. 2019, A&A, 631, A131 Mannucci, F., Belfiore, F., Curti, M., et al. 2021, MNRAS, 508, 1582 McDonough, B. & Brainerd, T. G. 2022, ApJ, 933, 161 Mead, A. J., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401 Menon, H., Wesolowski, L., Zheng, G., et al. 2015, Computational Astrophysics and Cosmology, 2, 1\n- Merson, A., Smith, A., Benson, A., Wang, Y., & Baugh, C. 2019, MNRAS, 486, 5737\n- Miller, L., Heymans, C., Kitching, T. D., et al. 2013, MNRAS, 429, 2858 Moessner, R. & Jain, B. 1998, MNRAS, 294, L18\n- More, S., van den Bosch, F. C., & Cacciato, M. 2009, MNRAS, 392, 917 Moretti, C., Tsedrik, M., Carrilho, P., & Pourtsidou, A. 2023, J. Cosmology Astropart. Phys., 12, 025\n- Murray, S. G., Power, C., & Robotham, A. S. G. 2013, Astronomy and Computing, 3, 23\n- Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493  \nOsterbrock, D. & Ferland, G. 2006, Astrophysics Of Gas Nebulae and Active Galactic Nuclei (University Science Books)  \nPeacock, J. A. & Smith, R. E. 2000, MNRAS, 318, 1144  \nPillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077  \n- Planck Collaboration: Aghanim, N., Akrami, Y., Ashdown, M., et al. 2020,  \nA&A, 641, A6  \n- Polletta, M., Tajer, M., Maraschi, L., et al. 2007, ApJ, 663, 81\n- Potter, D. & Stadel, J. 2016, PKDGRAV3: Parallel gravity code, Astrophysics Source Code Library\n- Potter, D., Stadel, J., & Teyssier, R. 2017, Computational Astrophysics and Cosmology, 4, 2\n- Pozzetti, L., Hirata, C. M., Geach, J. E., et al. 2016, A&A, 590, A3\n- Prévot, M. L., Lequeux, J., Maurice, E., Prévot, L., & Rocca-Volmerange, B. 1984, A&A, 132, 389\n- Reddick, R. M., Wechsler, R. H., Tinker, J. L., & Behroozi, P. S. 2013, ApJ, 771, 30\n- Reddy, N. A., Kriek, M., Shapley, A. E., et al. 2015, ApJ, 806, 259\n- Robotham, A. S. G. & Howlett, C. 2018, Research Notes of the American Astronomical Society, 2, 55\n- Rowe, B. T. P., Jarvis, M., Mandelbaum, R., et al. 2015, Astronomy and Computing, 10, 121  \nSaito, S., de la Torre, S., Ilbert, O., et al. 2020, MNRAS, 494, 199  \nSánchez, A. G. 2020, Phys. Rev. D, 102, 123511  \nSandage, A. & Visvanathan, N. 1978, ApJ, 223, 707  \n- Sarron, F., Martinet, N., Durret, F., & Adami, C. 2018, A&A, 613, A67  \nSchaye, J., Kugel, R., Schaller, M., et al. 2023, MNRAS, 526, 4978  \nSchechter, P. 1976, ApJ, 203, 297  \n- Schneider, A. & Teyssier, R. 2015, J. Cosmology Astropart. Phys., 2015, 049  \nSchneider, A., Teyssier, R., Potter, D., et al. 2016, J. Cosmology Astropart. Phys.,  \n4, 047  \nSchneider, P., Kilbinger, M., & Lombardi, M. 2005, A&A, 431, 9  \nScoccimarro, R., Sheth, R. K., Hui, L., & Jain, B. 2001, ApJ, 546, 20  \nSeljak, U. 2000, MNRAS, 318, 203  \n- Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41  \nSingh, S. & Mandelbaum, R. 2016, MNRAS, 457, 2301  \nSinha, M. & Garrison, L. H. 2020, MNRAS, 491, 3022  \n- Skibba, R. A. & Sheth, R. K. 2009, MNRAS, 392, 1080\n- Smith, A., Cole, S., Grove, C., Norberg, P., & Zarrouk, P. 2022, MNRAS, 516, 4529\n- Smith, A., Grove, C., Cole, S., et al. 2023, arXiv e-prints, arXiv:2312.08792\n- Somerville, R. S. & Primack, J. R. 1999, MNRAS, 310, 1087\n- Springel, V. 2005, MNRAS, 364, 1105\n- Springel, V., Pakmor, R., Zier, O., & Reinecke, M. 2021, MNRAS, 506, 2871\n- Springel, V., Wang, J., Vogelsberger, M., et al. 2008, MNRAS, 391, 1685\n- Stadel, J. G. 2001, PhD thesis, University of Washington, Seattle\n- Szapudi, I., Prunet, S., Pogosyan, D., Szalay, A. S., & Bond, J. R. 2001, ApJ, 548, L115\n- Takahashi, R., Nishimichi, T., Namikawa, T., et al. 2020, ApJ, 895, 113\n- Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152\n- Tallada, P., Carretero, J., Casals, J., et al. 2020, Astronomy and Computing, 32, 100391\n- Tasitsiomi, A., Kravtsov, A. V., Wechsler, R. H., & Primack, J. R. 2004, ApJ, 614, 533\n- Tatum, J. B. 1985, JRASC, 79, 302\n- Teyssier, R. 2002, A&A, 385, 337\n- Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709\n- Tiret, O., Combes, F., Angus, G. W., Famaey, B., & Zhao, H. S. 2007, A&A, 476, L1\n- To, C.-H., DeRose, J., Wechsler, R. H., et al. 2024, ApJ, 961, 59  \nTremonti, C. A., Heckman, T. M., Kau  \nff  \nmann, G., et al. 2004, ApJ, 613, 898  \nvan den Bosch, F. C., Aquino, D., Yang, X., et al. 2008, MNRAS, 387, 79  \n- Villaescusa-Navarro, F., Genel, S., Anglés-Alcázar, D., et al. 2023, ApJS, 265, 54\n- Vogelsberger, M., Marinacci, F., Torrey, P., & Puchwein, E. 2020, Nature Reviews Physics, 2, 42\n- Wainer, H. & Thissen, D. 1976, Psychometrica, 41, 9\n- Weaver, J. R., Davidzon, I., Toft, S., et al. 2023, A&A, 677, A184\n- Weaver, J. R., Kau ff mann, O. B., Ilbert, O., et al. 2022, ApJS, 258, 11\n- White, S. D. M. & Frenk, C. S. 1991, ApJ, 379, 52\n- White, S. D. M. & Rees, M. J. 1978, MNRAS, 183, 341\n- Zehavi, I., Zheng, Z., Weinberg, D. H., et al. 2011, ApJ, 736, 59\n- Zhai, Z., Chuang, C.-H., Wang, Y., Benson, A., & Yepes, G. 2021a, MNRAS, 501, 3490\n- Zhai, Z., Wang, Y., Benson, A., Chuang, C.-H., & Yepes, G. 2021b, MNRAS, 505, 2784\n- 1 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s / n, 08193 Barcelona, Spain\n- 2 Institut d\'Estudis Espacials de Catalunya (IEEC), Edifici RDIT, Campus UPC, 08860 Castelldefels, Barcelona, Spain\n- 3 Department of Astrophysics, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland\n- 4 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain\n- 5 Port d\'Informació Científica, Campus UAB, C. Albareda s / n, 08193 Bellaterra (Barcelona), Spain\n- 6 INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93 / 3, 40129 Bologna, Italy\n- 7 Institut d\'Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis boulevard Arago, 75014 Paris, France\n- 8 Institut d\'Astrophysique de Paris, 98bis Boulevard Arago, 75014, Paris, France\n- 9 Kavli Institute for the Physics and Mathematics of the Universe (WPI), University of Tokyo, Kashiwa, Chiba 277-8583, Japan\n- 10 Laboratoire Univers et Théorie, Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, 92190 Meudon, France\n- 11 Instituto de Astrofísica de Canarias, Calle Vía Láctea s / n, 38204, San Cristóbal de La Laguna, Tenerife, Spain\n- 12 Instituto de Astrofísica de Canarias (IAC); Departamento de Astrofísica, Universidad de La Laguna (ULL), 38200, La Laguna, Tenerife, Spain\n- 13 Université PSL, Observatoire de Paris, Sorbonne Université, CNRS, LERMA, 75014, Paris, France\n- 14 Université Paris-Cité, 5 Rue Thomas Mann, 75013, Paris, France\n- 15 Dipartimento di Fisica - Sezione di Astronomia, Università di Trieste, Via Tiepolo 11, 34131 Trieste, Italy\n- 16 INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34143 Trieste, Italy\n- 17 INFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste TS, Italy\n- 18 IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2, 34151 Trieste, Italy\n- 19 Institut de Física d\'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain\n- 20 Instituto de Astronomia Teorica y Experimental (IATECONICET), Laprida 854, X5000BGR, Córdoba, Argentina\n- 21 Minnesota Institute for Astrophysics, University of Minnesota, 116 Church St SE, Minneapolis, MN 55455, USA\n- 22 Institut de Ciencies de l\'Espai (IEEC-CSIC), Campus UAB, Carrer de Can Magrans, s / n Cerdanyola del Vallés, 08193 Barcelona, Spain\n- 23 Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstr. 25 / 8, 6020 Innsbruck, Austria\n- 24 Leiden Observatory, Leiden University, Einsteinweg 55, 2333 CC Leiden, The Netherlands\n- 25 Department of Astronomy, School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China\n- 26 Shanghai Key Laboratory for Particle Physics and Cosmology, Shanghai 200240, China\n- 27 Waterloo Centre for Astrophysics, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada\n- 28 Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada\n- 29 INFN Gruppo Collegato di Parma, Viale delle Scienze 7 / A 43124 Parma, Italy\n- 30 SISSA, International School for Advanced Studies, Via Bonomea 265, 34136 Trieste TS, Italy\n- 31 International Centre for Theoretical Physics (ICTP), Strada Costiera 11, 34151 Trieste, Italy\n- 32 Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544, USA\n- 33 Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, 31400 Toulouse, France  \nArticle number, page 40 of 45  \n- 34 Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK\n- 35 Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK\n- 36 Institute of Physics, Laboratory for Galaxy Evolution, Ecole Polytechnique Fédérale de Lausanne, Observatoire de Sauverny, CH1290 Versoix, Switzerland\n- 37 Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1, 85748 Garching, Germany\n- 38 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France\n- 39 Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK\n- 40 Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada\n- 41 Dipartimento di Fisica "Aldo Pontremoli", Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy\n- 42 INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy\n- 43 Department of Astronomy, University of Geneva, ch. d\'Ecogia 16, 1290 Versoix, Switzerland\n- 44 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark\n- 45 Cosmic Dawn Center (DAWN)\n- 46 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK\n- 47 INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy\n- 48 Departamento de Física Teórica, Facultad de Ciencias, Universidad Autónoma de Madrid, 28049 Cantoblanco, Madrid, Spain\n- 49 Serra Húnter Fellow, Departament de Física, Universitat Autònoma de Barcelona, E-08193 Bellaterra, Spain\n- 50 Konkoly Observatory, HUN-REN CSFK, MTA Centre of Excellence, Budapest, Konkoly Thege Miklós út 15-17. H-1121, Hungary\n- 51 MTA-CSFK Lendület Large-Scale Structure Research Group, Konkoly-Thege Miklós út 15-17, H-1121 Budapest, Hungary\n- 52 Université Paris-Saclay, CNRS, Institut d\'astrophysique spatiale, 91405, Orsay, France\n- 53 ESAC / ESA, Camino Bajo del Castillo, s / n., Urb. Villafranca del Castillo, 28692 Villanueva de la Cañada, Madrid, Spain\n- 54 School of Mathematics and Physics, University of Surrey, Guildford, Surrey, GU2 7XH, UK\n- 55 Institut für Theoretische Physik, University of Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany\n- 56 Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93 / 2, 40129 Bologna, Italy\n- 57 INFN-Sezione di Bologna, Viale Berti Pichat 6 / 2, 40127 Bologna, Italy\n- 58 INAF-Osservatorio Astronomico di Padova, Via dell\'Osservatorio 5, 35122 Padova, Italy\n- 59 Universitäts-Sternwarte München, Fakultät für Physik, LudwigMaximilians-Universität München, Scheinerstrasse 1, 81679 München, Germany\n- 60 INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese (TO), Italy\n- 61 Dipartimento di Fisica, Università di Genova, Via Dodecaneso 33, 16146, Genova, Italy\n- 62 INFN-Sezione di Genova, Via Dodecaneso 33, 16146, Genova, Italy\n- 63 Department of Physics "E. Pancini", University Federico II, Via Cinthia 6, 80126, Napoli, Italy\n- 64 INFN section of Naples, Via Cinthia 6, 80126, Napoli, Italy\n- 65 Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, PT4150-762 Porto, Portugal\n- 66 Faculdade de Ciências da Universidade do Porto, Rua do Campo de Alegre, 4150-007 Porto, Portugal\n- 67 Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, 10125 Torino, Italy\n- 68 INFN-Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy\n- 69 INAF-IASF Milano, Via Alfonso Corti 12, 20133 Milano, Italy\n- 70 Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, 52056 Aachen, Germany\n- 71 INAF-Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone, Italy\n- 72 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, Viale Berti Pichat 6 / 2, 40127 Bologna, Italy\n- 73 European Space Agency / ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma, Italy\n- 74 Université Claude Bernard Lyon 1, CNRS / IN2P3, IP2I Lyon, UMR 5822, Villeurbanne, F-69100, France\n- 75 Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland\n- 76 UCB Lyon 1, CNRS / IN2P3, IUF, IP2I Lyon, 4 rue Enrico Fermi, 69622 Villeurbanne, France\n- 77 Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, Edifício C8, Campo Grande, PT1749-016 Lisboa, Portugal\n- 78 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, 1749-016 Lisboa, Portugal\n- 79 INAF-Istituto di Astrofisica e Planetologia Spaziali, via del Fosso del Cavaliere, 100, 00100 Roma, Italy\n- 80 INFN-Padova, Via Marzolo 8, 35131 Padova, Italy\n- 81 School of Physics, HH Wills Physics Laboratory, University of Bristol, Tyndall Avenue, Bristol, BS8 1TL, UK\n- 82 Aix-Marseille Université, CNRS / IN2P3, CPPM, Marseille, France\n- 83 Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Via Irnerio 46, 40126 Bologna, Italy\n- 84 FRACTAL S.L.N.E., calle Tulipán 2, Portal 13 1A, 28231, Las Rozas de Madrid, Spain\n- 85 Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, 0315 Oslo, Norway\n- 86 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, 91109, USA\n- 87 Department of Physics, Lancaster University, Lancaster, LA1 4YB, UK\n- 88 Felix Hormuth Engineering, Goethestr. 17, 69181 Leimen, Germany\n- 89 Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark\n- 90 Cosmic Dawn Center (DAWN), Denmark\n- 91 Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany\n- 92 NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA\n- 93 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK\n- 94 Department of Physics and Helsinki Institute of Physics, Gustaf Hällströmin katu 2, 00014 University of Helsinki, Finland\n- 95 Université de Genève, Département de Physique Théorique and Centre for Astroparticle Physics, 24 quai Ernest-Ansermet, CH1211 Genève 4, Switzerland\n- 96 Department of Physics, P.O. Box 64, 00014 University of Helsinki, Finland\n- 97 Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland\n- 98 European Space Agency / ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands  \n99  \nAix-Marseille Université, CNRS, CNES, LAM, Marseille, France  \n- 100 NOVA optical infrared instrumentation group at ASTRON, Oude Hoogeveensedijk 4, 7991PD, Dwingeloo, The Netherlands\n- 101 INFN-Sezione di Milano, Via Celoria 16, 20133 Milano, Italy\n- 102 Universität Bonn, Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany\n- 103 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, via Piero Gobetti 93 / 2, 40129 Bologna, Italy  \n- 104 Department of Physics, Institute for Computational Cosmology, Durham University, South Road, DH1 3LE, UK\n- 105 Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA\n- 106 Université Côte d\'Azur, Observatoire de la Côte d\'Azur, CNRS, Laboratoire Lagrange, Bd de l\'Observatoire, CS 34229, 06304 Nice cedex 4, France\n- 107 Université Paris Cité, CNRS, Astroparticule et Cosmologie, 75013 Paris, France\n- 108 University of Applied Sciences and Arts of Northwestern Switzerland, School of Engineering, 5210 Windisch, Switzerland\n- 109 Department of Physics and Astronomy, University of Aarhus, Ny Munkegade 120, DK-8000 Aarhus C, Denmark\n- 110 Space Science Data Center, Italian Space Agency, via del Politecnico snc, 00133 Roma, Italy\n- 111 Centre National d\'Etudes Spatiales - Centre spatial de Toulouse, 18 avenue Edouard Belin, 31401 Toulouse Cedex 9, France\n- 112 Institute of Space Science, Str. Atomistilor, nr. 409 M˘agurele, Ilfov, 077125, Romania\n- 113 Departamento de Astrofísica, Universidad de La Laguna, 38206, La Laguna, Tenerife, Spain\n- 114 Dipartimento di Fisica e Astronomia "G. Galilei", Università di Padova, Via Marzolo 8, 35131 Padova, Italy\n- 115 Departamento de Física, FCFM, Universidad de Chile, Blanco Encalada 2008, Santiago, Chile\n- 116 Satlantis, University Science Park, Sede Bld 48940, Leioa-Bilbao, Spain\n- 117 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal\n- 118 Universidad Politécnica de Cartagena, Departamento de Electrónica y Tecnología de Computadoras, Plaza del Hospital 1, 30202 Cartagena, Spain\n- 119 Centre for Information Technology, University of Groningen, P.O. Box 11044, 9700 CA Groningen, The Netherlands\n- 120 INFN-Bologna, Via Irnerio 46, 40126 Bologna, Italy\n- 121 Dipartimento di Fisica, Università degli studi di Genova, and INFN-Sezione di Genova, via Dodecaneso 33, 16146, Genova, Italy\n- 122 Institut de Physique Théorique, CEA, CNRS, Université ParisSaclay 91191 Gif-sur-Yvette Cedex, France\n- 123 INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy\n- 124 Astronomical Observatory of the Autonomous Region of the Aosta Valley (OAVdA), Loc. Lignan 39, I-11020, Nus (Aosta Valley), Italy\n- 125 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK\n- 126 School of Physics and Astronomy, Cardi ff University, The Parade, Cardi ff , CF24 3AA, UK\n- 127 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, 10010, New York, NY, USA\n- 128 Université Paris-Saclay, CNRS / IN2P3, IJCLab, 91405 Orsay, France\n- 129 INFN-Sezione di Roma, Piazzale Aldo Moro, 2 - c / o Dipartimento di Fisica, Edificio G. Marconi, 00185 Roma, Italy\n- 130 Université St Joseph; Faculty of Sciences, Beirut, Lebanon\n- 131 Junia, EPA department, 41 Bd Vauban, 59800 Lille, France\n- 132 ICSC - Centro Nazionale di Ricerca in High Performance Computing, Big Data e Quantum Computing, Via Magnanelli 2, Bologna, Italy\n- 133 Instituto de Física Teórica UAM-CSIC, Campus de Cantoblanco, 28049 Madrid, Spain\n- 134 CERCA / ISO, Department of Physics, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH 44106, USA\n- 135 Departamento de Física Fundamental. Universidad de Salamanca. Plaza de la Merced s / n. 37008 Salamanca, Spain\n- 136 IRFU, CEA, Université Paris-Saclay 91191 Gif-sur-Yvette Cedex, France  \n- 137 Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy\n- 138 Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy\n- 139 Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France\n- 140 Ludwig-Maximilians-University, Schellingstrasse 4, 80799 Munich, Germany\n- 141 Max-Planck-Institut für Physik, Boltzmannstr. 8, 85748 Garching, Germany\n- 142 Thales Services S.A.S., 290 Allée du Lac, 31670 Labège, France\n- 143 Institute Lorentz, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands\n- 144 Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA\n- 145 Department of Physics & Astronomy, University of California Irvine, Irvine CA 92697, USA\n- 146 Department of Astronomy & Physics and Institute for Computational Astrophysics, Saint Mary\'s University, 923 Robie Street, Halifax, Nova Scotia, B3H 3C3, Canada\n- 147 Departamento Física Aplicada, Universidad Politécnica de Cartagena, Campus Muralla del Mar, 30202 Cartagena, Murcia, Spain\n- 148 Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK\n- 149 Department of Astronomy, University of Florida, Bryant Space Science Center, Gainesville, FL 32611, USA\n- 150 Department of Computer Science, Aalto University, PO Box 15400, Espoo, FI-00 076, Finland\n- 151 Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing (GCCL), 44780 Bochum, Germany\n- 152 DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 155, 2200 Copenhagen, Denmark\n- 153 Centre for Electronic Imaging, Open University, Walton Hall, Milton Keynes, MK7 6AA, UK\n- 154 XCAMLimited, 2 Stone Circle Road, Northampton, NN3 8RF, UK\n- 155 Department of Physics and Astronomy, Vesilinnantie 5, 20014 University of Turku, Finland\n- 156 Serco for European Space Agency (ESA), Camino bajo del Castillo, s / n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain\n- 157 ARC Centre of Excellence for Dark Matter Particle Physics, Melbourne, Australia\n- 158 Centre for Astrophysics & Supercomputing, Swinburne University of Technology, Victoria 3122, Australia\n- 159 School of Physics and Astronomy, Queen Mary University of London, Mile End Road, London E1 4NS, UK\n- 160 Department of Physics and Astronomy, University of the Western Cape, Bellville, Cape Town, 7535, South Africa\n- 161 ICTP South American Institute for Fundamental Research, Instituto de Física Teórica, Universidade Estadual Paulista, São Paulo, Brazil\n- 162 Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, Stockholm, SE-106 91, Sweden\n- 163 Astrophysics Group, Blackett Laboratory, Imperial College London, London SW7 2AZ, UK\n- 164 Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, Avenue des Martyrs, 38000, Grenoble, France\n- 165 INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125, Firenze, Italy\n- 166 Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 2, 00185 Roma, Italy\n- 167 Centro de Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto, Portugal\n- 168 Dipartimento di Fisica, Università di Roma Tor Vergata, Via della Ricerca Scientifica 1, Roma, Italy\n- 169 INFN, Sezione di Roma 2, Via della Ricerca Scientifica 1, Roma, Italy  \n- 170 Aurora Technology for European Space Agency (ESA), Camino bajo del Castillo, s / n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain\n- 171 Center for Frontier Science, Chiba University, 1-33 Yayoi-cho, Inage-ku, Chiba 263-8522, Japan\n- 172 Department of Physics, Graduate School of Science, Chiba University, 1-33 Yayoi-Cho, Inage-Ku, Chiba 263-8522, Japan\n- 173 Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The University of Edinburgh, Edinburgh EH9 3FD, UK\n- 174 Caltech / IPAC, 1200 E. California Blvd., Pasadena, CA 91125, USA\n- 175 Theoretical astrophysics, Department of Physics and Astronomy, Uppsala University, Box 515, 751 20 Uppsala, Sweden\n- 176 School of Physics & Astronomy, University of Southampton, Highfield Campus, Southampton SO17 1BJ, UK\n- 177 ASTRON, the Netherlands Institute for Radio Astronomy, Postbus 2, 7990 AA, Dwingeloo, The Netherlands\n- 178 Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands\n- 179 Anton Pannekoek Institute for Astronomy, University of Amsterdam, Postbus 94249, 1090 GE Amsterdam, The Netherlands\n- 180 Department of Physics, Royal Holloway, University of London, TW20 0EX, UK\n- 181 Department of Physics and Astronomy, University of California, Davis, CA 95616, USA\n- 182 Center for Cosmology and Particle Physics, Department of Physics, New York University, New York, NY 10003, USA\n- 183 Department of Astronomy, University of Massachusetts, Amherst, MA01003, USA\n- 184 Observatorio Nacional, Rua General Jose Cristino, 77-Bairro Imperial de Sao Cristovao, Rio de Janeiro, 20921-400, Brazil\n- 185 Instituto de Astrofísica de Canarias, c / Via Lactea s / n, La Laguna E-38200, Spain. Departamento de Astrofísica de la Universidad de La Laguna, Avda. Francisco Sanchez, La Laguna, E-38200, Spain\n- 186 Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281 S9, 9000 Gent, Belgium\n- 187 Space Telescope Science Institute, 3700 San Martin Dr, Baltimore, MD21218, USA', 'Appendix A: Halo mass function with different mass estimates': 'The Flagship halo catalogue was produced with the ROCKSTAR code adapted to the lightcone (see Sect. 4). ROCKSTAR computes di ff erent estimates of the halo mass. The mass values that we calculated are the mass of the particles linked together with a friends-of-friends algorithm of linking length b = 0 . 2, M fof; the mass contained within the virial radius, M vir; the sum of the mass of the bound particles within the virial radius, M bound; the mass of the particles within an overdensity of 200 relative to the background density, M 200b; and the mass of the particles within an overdensity of 200 relative to the critical density, M 200c.  \nIn this Appendix, we provide a comparison of the cumulative halo mass function resulting from the di ff erent mass estimates and also to the T08 M vir halo mass function as a function of redshift.  \nFigure A.1 shows the ratio of the cumulative halo mass function for all the mass estimates to the cumulative M bound HMF. It also contains the ratio of the T08 M bound HMF to the M bound HMF for comparison. The M vir and M bound HMFs coincide except for small di ff erences at the low and high mass end due to the unbound particle rejection process. The M 200b HMF is higher than the M bound HMF at low redshift z < 1, similar at 1 . 0 < z < 1 . 5, and lower at z > 1 . 5. The M 200c HMF is always lower than the M bound HMF. The M 200b and M 200c HMFs vary as expected due to the di ff erent redshift dependence of the density threshold used to define the mass. The ratio of the background density to the critical density gets lower at low redshift when the e ff ect of ΩΛ starts to be important, therefore increasing the ratio of M 200b to the M 200c values. The M fof HMF is quite di ff erent from the M bound HMF. This is expected as the friends-of-friends halo-finding technique and its mass definition is di ff erent from the rest. The T08 M vir HMF shows the same trends as described in Sect. 4.1.', 'Appendix B: The SciPIC algorithm': 'SciPIC , first described in Carretero et al. (2017) and named after Scientific Pipeline at PIC, is a suite of algorithms integrated into a powerful pipeline dedicated to the generation of massive synthetic galaxy catalogues based on halo catalogues coming from N -body dark matter cosmological simulations. SciPIC is the algorithm used to generate the Euclid Flagship galaxies. Originally, most of the algorithms were developed to the production of the MICE galaxy catalogues. In the case of MICE, the code was written in C and executed on di ff erent desktops. In particular, a significant amount of time was dedicated to input-output tasks, given that several steps were executed sequentially. However, with the Flagship parent halo catalogue containing a much larger volume of data, the methodology employed for running the MICE catalogue is unfeasible. The code has been refactored, optimised and ported to Python. Additionally, it runs on top of Apache Spark, an engine for scaleable computing. SciPIC is executed in the PIC Big Data platform, based on Hadoop, which comprises 20 nodes for a total processing power of 960 CPUs. The code runs e ffi ciently with a fast interface with the CosmoHub portal (Tallada et al. 2020), where the input halo catalogue is ingested, and the output galaxy catalogues are stored and distributed. The current implementation is able to generate a 15TB catalogue of 5B galaxies in 3 hours. The fact that it takes so little time to generate a mock allows for multiple iterations of progressive refinement.  \nRatio of Cumulative HMFs to Mbound HMF  \nFig. A.1. Natural logarithm of the ratio of the cumulative halo mass function for all mass estimates computed with ROCKSTAR and also the T08 M vir HMF to the M bound cumulative HMF as a function of redshift.  \n<!-- image -->  \nFig. C.1. Adopted model for the radial velocity dispersion profiles in virial units, for the case of dark matter concentration of 4, for the three galaxy classes. The symbols indicate the exact values, while the curves indicate the 2D-fifth-order polynomial approximations used in Flagship. One clearly sees the adopted constant values beyond the fit limits of r / r vir.  \n<!-- image -->', 'Appendix C: Details of the measurement of radial velocity dispersions of galaxies in their haloes': 'We provide here details on the mock velocities. The Jeans equation (14) is solved for the dynamical pressure, ρσ 2 r , which for the Tiret et al. (2007) velocity anisotropy profile of Eq. (15) is (Mamon et al. 2013)  \nρ ( r ) σ 2 r ( r ) = G ( r + r β ) 2 β ∞ Z r max r GLYPH<16> s + r β GLYPH<17> 2 β ∞ s 2 ρ ( s ) M ( s ) d s , (C.1)  \ntrivially yielding the radial velocity dispersion σ r for our chosen NFW model for ρ ( r ). For each of the three galaxy colour classes, we pre-computed σ r using Eq. (C.1) on a 25 × 16 grid of [log10( r / r vir) , log10( c dark)], with -2 ≤ log10( r / r vir) ≤ 0 . 4 and -1 ≤ log10( c dark) ≤ 1 . 5, in steps of 0.1 dex. We then fit a twodimensional 5th-order polynomial to allow Flagship to rapidly determine the radial velocity dispersion at given radial distances. The rms errors on log10 σ r are 0.0019, 0.0022, and 0.0026 respectively for red, green, and blue galaxies, i.e., better than 0.6% rms precision on σ r for all colour classes.'}