Horizon-AGN virtual observatory - 1. SED-fitting performance and forecasts for future imaging surveys

TL;DR

This study evaluates the performance of SED-fitting techniques on simulated galaxy data from Horizon-AGN, providing insights into redshift, stellar mass, and SFR estimation accuracy for upcoming large-scale surveys like LSST and Euclid.

Contribution

It offers a comprehensive analysis of photometric redshift and stellar mass recovery using realistic mock catalogues, and forecasts improvements from survey synergy.

Findings

Photometric redshifts are accurately recovered with dependencies similar to COSMOS2015.

Stellar masses are well estimated with typical dispersion below 0.1 dex.

Adding Euclid NIR data significantly enhances redshift and stellar mass accuracy.

Abstract

Using the ligthcone from the cosmological hydrodynamical simulation Horizon-AGN, we produced a photometric catalogue over $0<z<4$ with apparent magnitudes in COSMOS, DES, LSST-like, and Euclid-like filters at depths comparable to these surveys. The virtual photometry accounts for the complex star formation history and metal enrichment of Horizon-AGN galaxies, and consistently includes magnitude errors, dust attenuation and absorption by inter-galactic medium. The COSMOS-like photometry is fitted in the same configuration as the COSMOS2015 catalogue. We then quantify random and systematic errors of photometric redshifts, stellar masses, and star-formation rates (SFR). Photometric redshifts and redshift errors capture the same dependencies on magnitude and redshift as found in COSMOS2015, excluding the impact of source extraction. COSMOS-like stellar masses are well recovered with a dispersion typically lower than 0.1 dex. The simple star formation histories and metallicities of the templates induce a systematic underestimation of stellar masses at $z<1.5$ by at most 0.12 dex. SFR estimates exhibit a dust-induced bimodality combined with a larger scatter (typically between 0.2 and 0.6 dex). We also use our mock catalogue to predict photometric redshifts and stellar masses in future imaging surveys. We stress that adding Euclid near-infrared photometry to the LSST-like baseline improves redshift accuracy especially at the faint end and decreases the outlier fraction by a factor $\sim$2. It also considerably improves stellar masses, reducing the scatter up to a factor 3. It would therefore be mutually beneficial for LSST and Euclid to work in synergy.

Figures

Tables

Table 1. Table 1: A summary of the configurations envisaged in this study. Depths are quoted in AB magnitudes. The depths in all bands are summarised in Table  5 . A complete list of the COSMOS bands is provided in Table 1 of L16. A Horizon-AGN photometric catalogue is built for each configuration.

Name	Bands	i-band depth	NIR depth	References
COSMOS-like	26 bands from $u$ to $4.5\mu m$	26.2$\pm$ 0.1 (3$\sigma$)	24.7 $\pm$ 0.1 ($K_s$, 3$\sigma$)	Laigle et al. (2016)
LSST-like	$u$,$g$,$r$,$i$,$z$,$y$	27.0 (5$\sigma$)	NA	LSST Science Collaboration et al. (2009)
Euclid+DES	$g$,$r$,$i$,$z$,$riz$,$Y$,$J$,$H$	24.5 ($riz$,10 $\sigma$) and 24.3 ($i$, 10$\sigma$)	24.0 ($H$ band, 5$\sigma$)	Abbott et al. (2018); Laureijs et al. (2011)
Euclid+LSST	$u$,$g$,$r$,$i$,$z$,$y$,$riz$,$Y$,$J$,$H$	27.0 (5$\sigma$)	24.0 ($H$ band, 5$\sigma$)	Rhodes et al. (2017)

Table 2. Table 2: Statistical errors (NMAD) and percentage of catastrophic errors ( η 𝜂 \eta ) in different i + superscript 𝑖 i^{+} magnitude bins. Results for COSMOS2015 galaxies (L16) are compared to the outcome of our simulation.

$i^{+}$	COSMOS2015		Hz-AGN		Hz-AGN
			with IB		without IB
mag	NMAD	$\eta$ (%)	NMAD	$\eta$ (%)	NMAD	$\eta$ (%)
$($22,23$]$	0.010	1.7	0.008	0.0	0.023	0.0
$($23,24$]$	0.022	6.7	0.014	0.0	0.028	0.1
$($24,25$]$	0.034	10.2	0.026	0.5	0.037	0.8
$($25,26$]$	0.057	22.0	0.052	9.2	0.065	11.7

Table 3. Table 3: Mass completeness limits (90 and 50 per cent threshold) in the Euclid ( H < 24 𝐻 24 H<24 or r ​ i ​ z < 24.5 𝑟 𝑖 𝑧 24.5 riz<24.5 ) and LSST ( i < 25.3 𝑖 25.3 i<25.3 ) configurations as a function of redshift, as estimated from intrinsic quantities.

$z_{\mathrm{sim}}$ range	$z_{\mathrm{sim}}$ median	90% mass completeness [$\log M_{\odot }$]			50% mass completeness [$\log M_{\odot }$]
(0.00,0.25]	0.22	9.00	9.00	9.00	9.00	9.00	9.00
(0.25,0.50]	0.41	9.00	9.02	9.00	9.00	9.00	9.00
(0.50,0.75]	0.65	9.08	9.43	9.00	9.00	9.05	9.00
(0.75,1.00]	0.89	9.31	9.81	9.27	9.05	9.42	9.00
(1.00,1.25]	1.14	9.48	10.12	9.67	9.26	9.70	9.30
(1.25,1.50]	1.38	9.64	10.52	9.77	9.38	9.92	9.44
(1.50,1.75]	1.63	9.76	10.65	10.01	9.51	10.22	9.60
(1.75,2.00]	1.87	9.83	10.89	10.13	9.60	10.31	9.65
(2.00,2.25]	2.11	9.88	10.78	10.23	9.66	10.32	9.70
(2.25,2.50]	2.37	9.94	10.76	10.14	9.70	10.25	9.72
(2.50,2.75]	2.62	10.01	10.72	10.06	9.75	10.26	9.66
(2.75,3.00]	2.87	10.04	10.68	9.98	9.82	10.22	9.58
(3.00,3.25]	3.12	10.17	10.53	9.93	9.87	10.22	9.54
(3.25,3.50]	3.36	10.21	10.34	9.81	9.94	10.15	9.48
(3.50,3.75]	3.61	10.27	10.33	9.82	9.96	9.98	9.37
(3.75,4.00]	3.87	10.28	10.26	9.71	9.99	10.07	9.32

Table 4. Table 4: Statistical errors (defined as normalized median absolute deviation, NMAD) and percentage of catastrophic errors ( η 𝜂 \eta ) in different i 𝑖 i and H 𝐻 H magnitude bins, in the different configurations studied here.

$i$ band	Euclid+DES		Euclid+LSST		LSST only
mag	NMAD	$\eta$ (%)	NMAD	$\eta$ (%)	NMAD	$\eta$ (%)
$($22,23$]$	0.044	3.7	0.017	0.1	0.017	0.1
$($23,24$]$	0.081	15.7	0.018	0.2	0.020	0.7
$($24,25$]$	–	$>$40.0	0.031	2.6	0.037	4.6

Table 5. Table 5: A summary of the adopted depths in all bands. The depths of intermediate bands (IB) are detailed in Laigle et al. ( 2016 ) .

Survey	band	depth
	$u$	26.6
	$B$	27.0
	$V$	26.2
	$r$	26.5
	$i^{+}$	26.2
COSMOS-like	$z^{++}$	25.9
$3\sigma$ depth,	$Y$	25.3
extended sources	$J$	24.9
	$H$	24.6
	$K_s$	24.7
	IB	25-26
	$riz$	24.5
Euclid-like	$Y$	24.0
$5\sigma$ depth,	$J$	24.0
extended sources	$H$	24.0
	$g$	24.6
DES-like	$r$	24.1
$5\sigma$ depth,	$i$	24.0
extended sources	$z$	23.9
	$u$	26.3
	$g$	27.5
LSST-like	$r$	27.7
$5\sigma$ depth,	$i$	27.0
extended sources	$z$	26.2
	$y$	24.9

Table 6. Table 6: Statistical errors (NMAD) and percentage of catastrophic errors ( η 𝜂 \eta ) in different i 𝑖 i magnitude bins, without absorption, with IGM absorption and with both dust and IGM absorption.

$i$ band	No absorption		IGM		IGM+dust
mag	NMAD	$\eta$ (%)	NMAD	$\eta$ (%)	NMAD	$\eta$ (%)
$($22,23$]$	0.009	0.0	0.009	0.0	0.008	0.0
$($23,24$]$	0.013	0.0	0.013	0.0	0.014	0.0
$($24,25$]$	0.022	0.6	0.021	0.4	0.026	0.5
$($25,26$]$	0.049	8.7	0.045	7.6	0.052	9.2

Equations

$$F({\nu }_0)=e^{-{\tau }_{\alpha }({\nu }_0)},$$

$${\tau }_{\alpha }({\nu }_0)={\int }_0^{x_s}\mathrm{d}x\frac{{\sigma }_{\alpha }{n}_{\mathrm{HI}}(x,z)}{1+z},$$

$${\sigma }_{\alpha }=\frac{{\sigma }_{\alpha ,0}c}{b(x,z)\sqrt{\pi }}{e}^{-\frac{{\left(v(x,z)(1+z){\nu }_0-c{\nu }_{\alpha }+c(1+z){\nu }_0\right)}^2}{{\nu }_{\alpha }^2{b}^2(x,z)}},$$

$$\alpha (T)n_{\mathrm{e}}(1-x_{\mathrm{HI}})=\gamma (T)n_{\mathrm{e}}{x}_{\mathrm{HI}}+\Gamma x_{\mathrm{HI}},$$

Full text

††pubyear: 2016††pagerange: Horizon-AGN virtual observatory – 1. SED-fitting performance and forecasts for future imaging surveys–D.2

1 Introduction

Our understanding of galaxy formation, evolution, and of their distribution in the large-scale structure has taken a giant step forward in the last decade, owing to large multi-wavelength datasets. Properties of different galaxy populations, and their evolution across cosmic time, can be constrained by measuring one-point statistics, such that the luminosity and stellar mass functions (e.g. Ilbert et al., 2006; Ilbert et al., 2013; Davidzon et al., 2017; Bundy et al., 2017). Two-point statistics, i.e. measuring the spatial correlation of galaxies, make it possible to investigate the role of the local environment (e.g. Abbas & Sheth, 2006; de la Torre et al., 2010; Hatfield & Jarvis, 2017) and to infer halo properties, via simplifying assumptions such as the so-called halo model (e.g. McCracken et al., 2015; Coupon et al., 2015; Legrand et al., 2018). More generally, higher order statistics (see Moresco et al., 2017) as well as topological tools such as filament tracers, can help disentangle complex environmental effects, distinct from the isotropic influence of local density peaks (e.g. Malavasi et al., 2017; Laigle et al., 2018; Kraljic et al., 2018). When implementing such statistics, one must asses the impact of observational biases on inferring the underlying properties of the population.

In particular, investigations focusing on galaxy stellar mass assembly rely on three fundamental quantities: photometric redshifts, stellar masses and star-formation rates. Large-area surveys can significantly reduce statistical errors in these kinds of measurements, and probe a wide variety of galaxy types and environments. Therefore, the dominant source of uncertainties in state-of-the-art studies became the selection biases of the surveys, the source extraction techniques, and the physical models assumed in the analysis (when needed). Even when a high-resolution galaxy spectral energy distribution (SED) is available, inferring physical properties from it is an ill-conditioned problem (Moultaka & Pelat, 2000; Moultaka et al., 2004), which prevents a complete inversion approach to be successful (see e.g. Ocvirk et al., 2006). Difficulties are even more severe when only apparent magnitudes in broad-band filters are available. In that case, SED-fitting codes are routinely used because of their versatility. These codes fit pre-computed libraries of galaxy templates to the photometry of observed objects (see a review in Walcher et al., 2011; Conroy, 2013). Some very promising alternative techniques are also being developed (see Salvato et al., 2018, for a review), including “clustering redshift" (Newman, 2008; Ménard et al., 2013), “photo-web” (Aragon-Calvo et al., 2015) and more recently machine-learning (see e.g. Masters et al., 2015; Beck et al., 2017; Pasquet et al., 2018; Gomes et al., 2018; Hemmati et al., 2018). However these alternative techniques generally require large and representative spectroscopic samples, which is not the case of SED-fitting algorithms. In order to build a template library for the SED-fitting procedure, one relies nonetheless on several assumptions, mainly concerning star formation histories (SFHs), metal enrichment, and dust extinction and spatial distribution (Conroy, 2013). These priors inevitably introduce systematics in the recovered physical quantities, which in turn may impair the statistical measurements and bias conclusions on galaxy mass assembly scenarii. For instance Bundy et al. (2017) find that depending on the assumed SFH in the SED-fitting estimates, massive ($>3\times 10^{11}\,M_{\odot}$) galaxies between $z=0$ and 0.8 may show either a mild stellar mass growth or a lack of evolution; this systematic uncertainty is dominant, considering that their extremely large sample of galaxies ($>\!41,000$), collected across $\sim\!140$ deg2, makes shot noise and cosmic variance almost negligible.

Furthermore, in order to understand the physical processes regulating galaxy mass assembly, it is important to compare observational measurements to semi-analytical and hydrodynamical simulations, where different theoretical models of galaxy evolution have been implemented (e.g. De Lucia & Blaizot, 2007; Vogelsberger et al., 2013; Dubois et al., 2014; Schaye et al., 2015). At present, such a task is not straightforward: a fair comparison should take into account biases and uncertainties affecting the observational analysis before comparing to simulated galaxies.

Therefore it is of pivotal importance to assess the performances of photometric extraction and SED-fitting codes when recovering redshift and stellar mass in order to understand their impact on the statistical analyses of the galaxy population. Broadly speaking, observational biases can occur because of image confusion (i.e., blending between two nearby galaxies), the choice of algorithm used to extract galaxy flux, and the assumptions made in the SED-fitting procedure. Previous works have already explored some of these effects. As an example, Mobasher et al. (2015) have quantified the global performances of an exhaustive list of existing SED-fitting codes, while relying on a large observed and semi-analytical mock catalogue (see also Hildebrandt et al., 2010). Focusing on mass and age estimates, Pforr et al. (2012) and Pacifici et al. (2012) investigated the impact of the chosen template SFH, while the effect of dust and metallicity has been studied in Mitchell et al. (2013) and Hayward & Smith (2015).

Beyond the impact of simplistic SFHs (like the $\tau$-model defined in Bruzual A., 1983), metallicity or dust distribution (see also Guidi et al., 2016), the performance of SED fitting is extremely sensitive to the choice of photometric filters, the depth of the survey, and flux measurements (see Bernardi et al., 2013). Hydrodynamical simulations have already been widely used to test the impact of the photometry extraction, as they allow to work on -often high resolution- mock images of realistic galaxies. Amongst the tested effects, the choice of the apertures (Price et al., 2017) and the lack of resolution (Sorba & Sawicki, 2015, 2018, integrated photometry versus pixel-by-pixel fitting) have been found to systematically underestimate stellar masses (see also Sanderson et al., 2017) or to impair morphological estimators (Bottrell et al., 2017). All these past investigations underline the importance of understanding and quantifying biases when recovering physical parameters from surveys, which could be as large as two orders of magnitude in some particular mass and redshift ranges (see e.g. the effect of dust on stellar mass computation Mitchell et al., 2013). However, most of the literature is based either on simple phenomenological prescriptions or semi-analytical models, or when the sample is based on hydrodynamical simulation, it consists in no more than a handful of galaxies (e.g. Guidi et al., 2016). Hence we still lack a study relying on a sample that combines highly realistic baryon physics with a large cosmological volume (in order to minimize statistical uncertainties), capturing both galaxies’ internal properties and environment. Moreover, this study must be an end-to-end analysis, i.e. including the same limitations introduced by the observational strategy and data reduction pipeline in current or future surveys. The present work aims to remedy this gap.

To this end, we exploit the ligthcone from the Horizon-AGN cosmological hydrodynamical simulation (Dubois et al., 2014). From this simulation a mock catalogue of about 750,000 galaxies was extracted between $z=0$ and $4$, down to $M_{*}=10^{9}{\rm M}_{\odot}$. Hence our sample combines large number statistics over a wide redshift range with a wealth of information on galaxy properties. Our aim is to carefully understand possible systematics arising when fitting the complex photometry of the galaxies with simplified templates. For this purpose galaxy photometry has to be as realistic as possible. One advantage of using hydrodynamical simulations over SAMs for this work is to better resolve galaxies in space and time. Fluxes spatially vary across the simulated galaxies (in the limit of the resolution of the simulation) depending on metallicity enrichment and dust attenuation, and therefore the integrated photometry will present a complexity similar to the real galaxies. In addition, SFHs in the hydrodynamical simulation vary on a fine time grid and depend not only on the merger history of their host halo but also on stellar and AGN feedback, and on the detail of the gas accretion history. As emphasized in e.g. Mitchell et al. (2018), several quantities (e.g., the gas return time-scale) are naturally constrained by gravitational forces and hydrodynamics, whereas they would conversely need to be globally tuned in SAMs. Finally, the lightcone geometry mimics that of observed surveys, and allows us for instance to implement the attenuation by the intergalactic medium (IGM) for each galaxy by drawing individually lines of sight through the foreground gas distribution.

The goal of this first study is to assess the photometric redshift ($z_{\mathrm{phot}}$), stellar mass ($M_{\ast}$) and star-formation rate (SFR) uncertainties caused by the choice of the filters, the signal-to-noise ratio ($S/N$) of the photometry and the SED fitting recipe used for analyzing the real galaxies.

For this purpose, observed-frame photometry is post-processed with COSMOS-like $S/N$ for each galaxy of the Horizon-AGN ligthcone (as described in Section 2). Then photometric redshifts and physical properties ($M_{\ast}$ and SFR) of mock galaxies are measured by applying the same pipeline used in the COSMOS field (Laigle et al., 2016, hereafter L16). This procedure allows us to identify which source of uncertainty dominate the error budget (Section 3). After validation on the COSMOS2015 data, we mimic (in Section 4) the expected photometry for the Euclid mission (Laureijs et al., 2011), along with the Dark Energy Survey (DES, Abbott et al., 2018) and the Large Synoptic Survey Telescope (LSST, LSST Science Collaboration et al., 2009), to predict the expected $z_{\mathrm{phot}}$ and $M_{\ast}$ accuracy they should provide at completion. The possible synergy between these surveys is also explored. We then summarise our analysis and draw conclusions in Section 5. Additional material can be found in the Appendices, where we provides more details about how the virtual photometry has been computed (Appendix A); we further discuss dust and IGM absorption (Appendix B), zero-point magnitude offsets (Appendix C), and redshift errors (Appendix D). These virtual catalogs are going to be made publicly available at https://www.horizon-simulation.org/data.html.

Throughout this study, we use a flat $\Lambda$CDM cosmology with $H_{0}=70.4$ km s*-1* Mpc*-1*, $\Omega_{m}=0.272$, $\Omega_{\Lambda}=0.728$, and $n_{s}=0.967$ (Komatsu et al., 2011, WMAP-7). All magnitudes are in the AB (Oke, 1974) system. The initial mass function (IMF) follows Chabrier (2003). Quantities are said “observed" when they include observational noise (for magnitudes) or when they are measured through SED fitting (redshift, stellar mass and SFR). If directly derived from the simulation, they are defined as “intrinsic”.

2 Data and methods

2.1 Description of the observational surveys

The virtual photometric catalogue from the Horizon-AGN simulation is built to mimic the COSMOS2015 catalogue. It also includes the photometry expected from the Euclid space-based telescope111 https://www.euclid-ec.org (Abbott et al., 2018), the Large Synoptic Survey Telescope (LSST, LSST Science Collaboration et al., 2009) and the Dark Energy Survey (DES, Laureijs et al., 2011) in terms of filter passbands and depths. We briefly describe hereafter the different configurations investigated in this work to quantify the performances of galaxy redshift and physical property computation. Table 1 provides a summary of these surveys.

2.1.1 The COSMOS field

The COSMOS deep optical and near-infrared catalogue (COSMOS2015) described in Laigle et al. (2016, hereafter L16) is used as a reference to test the performances of our estimation of galaxy properties. The catalogue includes more than 1 million objects detected within the 2 deg2 of the COSMOS field, observed in 30 bands from UV to IR ($0.25\!-\!8\,\mu$m). Here the analysis is restricted to the “ultra-deep” stripes, i.e. four rectangular regions that in COSMOS2015 have been covered with higher near-Infrared (NIR) sensitivity (in the UltraVISTA-DR2 survey, $K_{\mathrm{s}}<24.7$, $3\sigma$) than the rest of the area.

COSMOS2015 contains far- and near-UV photometry (FUV and NUV, respectively) from GALEX (Zamojski et al., 2007), but only NUV was used for the estimation of photometric redshifts and masses. In the optical, it includes the same $u,B,V,r,i,z$ data as previous releases from the Canada-Hawaii-France and Subaru telescopes (Capak et al., 2007; Ilbert et al., 2009). This baseline is complemented with Subaru medium- and narrow-band images between $4,000$ and $8,500$ Å. In the NIR, $Y,J,H,K_{\rm s}$ images come from the second data release (DR2) of the UltraVISTA survey (McCracken et al., 2012), and the $Y$-band image from Subaru/Hyper-Suprime-Cam (HSC, Miyazaki et al., 2012). The Spitzer Large Area Survey with HSC (SPLASH, Capak et al. in prep.) provides mid-IR (MIR) coverage with the four IRAC channels centred at $3.6$, $4.5$, $5.8$, and $8.0\,\mu$m.

In order to derive photometry coherently across different bands, the point-spread function (PSF) in each filter has been rescaled using a Moffat profile modelling. After the PSF homogenisation, fluxes are extracted within fixed apertures of $3\arcsec$ using SExtractor (Bertin & Arnouts, 1996) in dual image mode. Following a reduction procedure similar to that described in McCracken et al. (2012), the detection image is a $\chi$-squared sum of the four NIR images of UltraVISTA DR2 and the $z^{++}$ band. Spitzer sources are extracted by means of the code IRACLEAN (Hsieh et al., 2012).

Estimates obtained through SED fitting (photometric redshift, stellar mass, and other physical quantities) are also provided for each entry of the catalogue. The method adopted to obtain these estimates is described in Section 2.4, where the same technique is applied to simulated galaxies. Further details about the COSMOS2015 catalogue can be found in L16.

2.1.2 Future surveys: Euclid and LSST

To compute galaxy properties from SED-fitting in comparable conditions to Euclid and LSST, Horizon-AGN galaxies are also post-processed to get the photometry in Euclid, LSST and DES filters with depths similar to the ones expected for these surveys. We note that sometimes the expected depths from the literature are provided for point sources, and as a consequence give generally too optimistic estimators of the limiting magnitudes of extended sources. Our adopted limiting magnitudes are therefore probably not exactly the ones which will be obtained in the future, but they nonetheless reflect the relative depths of these upcoming surveys. The photometric baselines is detailed below, and the adopted limiting magnitudes in all bands are summarised in Table 5.

Euclid+DES configuration

Euclid will provide photometry in one broad-band optical ($riz$ filter) and 3 NIR filters ($Y$, $J$, $H$) with expected depths at completion of 24.5 (10$\sigma$, extended sources) in the optical and 24.0 (5$\sigma$, point sources) in the NIR bands. This broadband baseline alone is not sufficient to compute photometric redshifts with a high enough accuracy (especially to constrain the Balmer break in the optical), therefore it has to be complemented with ground-based optical photometry (see e.g. Sorba & Sawicki, 2011). In particular, DES provides photometry over 5000 deg2 in the Southern sky in $g$, $r$, $i$ and $z$ with depth of 24.33, 24.08, 23.44, 22.69 (10 $\sigma$, extended sources, Abbott et al., 2018), which matches the Euclid requirements (Laureijs et al., 2011). DES photometry provides a finer sampling of the optical range than the Euclid $riz$ filter alone. Collaboration between Euclid and DES is planned. Therefore in the current work we explore what would be the expected performance of such a configuration.

LSST configuration

The survey conducted on LSST will provide photometry in the optical over $30000$ deg2. LSST single visit depth should reach 24.5 in $r$ ($5\sigma$, point sources), and the co-added survey depth should reach 26.3, 27.5, 27.7, 27.0, 26.2 and 24.9 ($5\sigma$, point sources) in $u$, $g$, $r$, $i$, $z$ and $y$ bands respectively (LSST Science Collaboration et al., 2009). For weak-lensing studies, the “gold" sample of LSST galaxies with a high $S/N$ is defined with a magnitude cut $i<25.3$. We use this cut in the present work when studying the LSST-like configuration.

In the Southern sky, Euclid and LSST will overlap over at least 7000 deg2. It is therefore natural to explore the possible gain to combine both datasets. To this end we also analyse Horizon-AGN galaxies in the Euclid+LSST-like configuration.

2.2 The Horizon-AGN simulation

This study relies on Horizon-AGN222http://www.horizon-simulation.org/ (Dubois et al., 2014), a cosmological hydrodynamical simulation in overall fairly good agreement with observations, in the redshift and mass regime of the present analysis (see Kaviraj et al., 2017).

The simulation box, run with the ramses code (Teyssier, 2002), is $L_{\rm box}=100\,h^{-1}\,\mathrm{Mpc}$ on a side, and the volume contains $1024^{3}$ dark matter (DM) particles, corresponding to a DM mass resolution of $8\times 10^{7}\,M_{\odot}$. The initially coarse $1024^{3}$ grid is adaptively refined down to $1$ physical kpc. The refinement procedure leads to a typical number of $6.5\times 10^{9}$ gas resolution elements (leaf cells) in the Horizon-AGN simulation at $z=1$.

Heating of the gas from a uniform UV background takes place after redshift $z_{\rm reion}=10$, following Haardt & Madau (1996). Gas can cool down to $10^{4}\,\rm K$ through H and He collision and with a contribution from metals that follows the rates tabulated in Sutherland & Dopita (1993). Star formation occurs in regions where gas number density is above $n_{0}=0.1\,\rm H\,cm^{-3}$, following a Schmidt law: $\dot{\rho}_{*}=\epsilon_{*}{\rho_{\rm g}/t_{\rm ff}}$, where $\dot{\rho}_{*}$ is the star formation rate mass density, $\rho_{\rm g}$ the gas mass density, $\epsilon_{*}=0.02$ the constant star formation efficiency and $t_{\rm ff}$ the gas local free-fall time. Feedback from stellar winds and supernova (both type Ia and II) are included into the simulation with mass, energy, and metal releases. Galactic black hole formation is also implemented in Horizon-AGN, with accretion efficiency tuned to match the black hole-galaxy scaling relations at $z=0$. Black hole energy is released in either quasar or radio mode depending on the accretion rate (see Dubois et al., 2012, for more details).

The lightcone has been extracted on-the-fly as described in Pichon et al. (2010). For the lightcone extraction, gas leaf-cells were replaced by gas particles, and treated as the stars and dark matter particles. All particles were extracted at each coarse time step according to their proper distance to the observer at the origin. In total the lightcone contains about 22000 portions of concentric shells. The lightcone projected area is 5 deg2 below $z=1$, and 1 deg2 above. However, we restrict ourselves to 1 deg2 over the whole redshift range considered in this study. The full lightcone up to $z=4$ contains about 19 replica of the Horizon-AGN box.

2.3 Generating a mock photometric catalogue

2.3.1 Galaxy extraction

The AdaptaHOP halo finder (Aubert et al., 2004) is run on the lightcone over $0<z<4$ to identify galaxies from the stellar particles distribution. Local stellar particle density is computed from the 20 nearest neighbours, and structures are selected with a density threshold equal to 178 times the average matter density at that redshift. Galaxies resulting in less than 50 particles ($\simeq 10^{8}\,{\rm M}_{\odot}$) are not included in the catalogue. Since the identification technique is redshift dependent, AdaptaHOP is run iteratively on thin lightcone slices (about 4000 slices up to $z=4$) of few comoving Mpc (cMpc). Slices are overlapping to avoid edge effects (i.e. cutting galaxies in the extraction) and duplicates are removed.

2.3.2 Galaxy SED computation and dust attenuation

Although the simulation assumes a Salpeter IMF (Salpeter, 1955) to model stellar mass losses, we have decided to post-process it with a Chabrier IMF (Chabrier, 2003). The choice of the IMF is significant, as it controls both the stellar mass loss prescription and the overall mass-to-light ratio. The Chabrier IMF brings the simulated galaxy counts in much better agreement with the COSMOS2015 galaxy counts (see Appendix A. Magnitudes are $\sim 0.4$ mag fainter with a Salpeter IMF compared to Chabrier’s). For each galaxy, each stellar particle is linked to a single stellar population (SSP) obtained with the stellar population synthesis model of Bruzual & Charlot (2003, hereafter BC03). Because stellar particle ages and metallicities vary on a much finer grid than the BC03 models, an interpolation is carried out between SSPs to reproduce the desired values. In addition, each SSP is also rescaled to match the initial stellar mass of the particle, in order to follow the same mass loss fraction for the simulated galaxies as in BC03 (see Appendix A.4). This rescaling is essential to avoid discrepancies between the intrinsic and computed galaxy properties coming from the different SSP prescriptions, which is out of topic for the present work. It should be noted that the metallicity of stellar particles in Horizon-AGN has been boosted by a empirically computed factor, to match the observed mass-metallicity relation. This factor $f_{Z}$ is redshift dependent as follows: $f_{Z}=4.08430-0.213574z-0.111197z^{2}$ (see Kaviraj et al., 2017, for more details).

Dust attenuation is also modelled for each star particle using the gas metal mass distribution as a proxy for the dust distribution. Gas metal mass is evaluated in a cube of $138$ comoving kpc3 around each galaxy and a constant dust-to-metal mass ratio is adopted. Although a mesh-based code is used to run Horizon-AGN and in particular to follow the gas distribution, gas cells have been turned in particles when extracting the lightcone. In order to get a smooth metal field around the galaxy from this gas particle distribution, a Delaunay tessellation is computed on the particles to avoid cells with null values in under-dense regions, and then interpolated on a regular grid with a resolution of $\sim 1$ ckpc. The dust column density and the optical depth along the line of sight are computed for each stellar particle in the galaxy using the ${\rm R}_{V}=3.1$ Milky Way dust grain model by Weingartner & Draine (2001). Further details on the dust computation and dust-to-metal mass ratio calibration are given in Appendix A. This dust attenuation model only takes into account absorption and does not include scattering. While this is not a problem in the rest-frame optical and NIR, the impact in the rest-frame UV is non negligible (Kaviraj et al., 2017). This effect is not corrected in our catalog and as a result, our galaxies are up to 0.8 mag brighter without scattering in the UV part of the spectrum. In fact, including scattering would have a similar impact as a steepening of the dust attenuation law in the UV. A dust-free version of the catalogue is also produced in order to isolate the impact of dust attenuation on the computation of galaxy physical properties.

Flux contamination by nebular emission lines is not included in our virtual photometry. Consistently, emission line parametrization is also turned off in the SED-fitting computation. In real surveys, emission lines can help determining the photometric redshifts. On the other hand, the dispersion in the emission line ratios is poorly modelled by the SED-fitting code and can bias the redshift estimation.

Finally, we stress that we do not model extinction by the Milky way. Photometry from observed survey is generally corrected from galactic extinction using galactic reddening maps (e.g. Schlegel et al., 1998). However, the amount of absorption in a given band will depend on the source SED. As shown in Galametz et al. (2017), the band-pass extinction can vary by up to 20% from the average correction depending on the SED. The discrepancy between the effective extinction and the average correction of the photometry can potentially lead to additional systematics which are therefore not accounted for here.

Note that Horizon-AGN only reproduces global mass assembly up to some point. Although the simulation broadly matches the mass function evolution with redshift and the SFR main sequence (Kaviraj et al., 2017), at low-mass ($\log M_{*}/{\rm M}_{\odot}>9.5$) the mass function is systematically overestimated given the present sub-grid stellar feedback and star-formation recipes. Conversely, at high redshift ($z>4$) it is also underestimated because of limited spatial resolution. To be conservative, forecasts were therefore limited within the redshift and mass range where the simulation is reliable. Appendix A.5 lists the limitations of our modelling.

2.3.3 IGM absorption

Prior to convolving the galaxy spectrum with photometric filters, attenuation by the IGM must be implemented. The knowledge of the gas distribution in the ligthcone allows us to consistently implement IGM absorption, while accounting for variation from one sight-line to the other. Conversely, it is globally accounted for with an analytical prescription at the SED-fitting stage. Therefore our photometric catalogue allows to test the effect of the inhomogeneous IGM attenuation on photometric redshifts and galaxy properties. The focus is on HI absorption in the Lyman-series (hence neglecting metal lines, such as e.g. the CIV forest). Details about the implementation of the attenuation on the line-of-sight of each galaxy is given in Appendix A. An IGM-free version of the catalogue is also produced in order to isolate the impact of IGM absorption.

2.3.4 Photometry extraction and error implementation

Eventually the integrated galaxy spectra are computed between 91 and $16\times 10^{5}$ Å (with 1221 wavelength points) by adding all the SSPs of a given galaxy together. Note that IR dust emission is not computed because not used here. Apparent total magnitudes are obtained by convolving the spectrum (redshifted to the intrinsic redshifts of the galaxies) with the same filter set as used in COSMOS2015 ($u$, $B$, $V$, $r$, $i^{+}$, $z^{++}$, $Y$, $J$, $H$, $K_{\rm s}$, [$3.6\mu$m], [$4.5\mu$m], and the 14 intermediate and narrow bands333At lower redshift, GALEX FUV and NUV filters are relevant, as they bring information about young stellar populations. However, it is more difficult to reproduce their flux extraction and the related uncertainties. Therefore, NUV and FUV are excluded from the mock catalogue.), Euclid, DES and LSST (see Section 2.1).

Photometric errors are added in each band to reproduce the $S/N$ distribution and sensitivity limit of COSMOS2015 and DES, and the ones expected for Euclid and LSST (see Table 1 and Appendix A). It should be emphasized here that the photometry has been derived from the entire distribution of star particles and not through a realistic flux extraction from images (as it is done in real observations via tools like SExtractor, Bertin & Arnouts, 1996). Consequently, our mock catalogue does not include all the associated photometric issues including potential systematic effects like blending, object fragmentation, imperfect background sky subtraction and PSF homogenisation and offsets due to the rescaling of fixed aperture to total fluxes, which will be the topic of a future work. Photometric errors as implemented in the current catalogue simply correspond to gaussian noise with standard deviation depending on the galaxy flux and the depth of the surveys.

2.4 SED-fitting: method

2.4.1 Photometric redshifts

Photometric redshifts ($z_{\mathrm{phot}}$) are computed using the code LePhare (Arnouts et al., 2002; Ilbert et al., 2006) with a configuration similar to Ilbert et al. (2013). The SED library includes spiral and elliptical galaxies from Polletta et al. (2007), along with bluer templates of young star-forming galaxies built by means of the BC03 model. Dust extinction is added to the templates according to one of the following attenuation curves: Prevot et al. (1984), Calzetti et al. (2000), or a modified version of Calzetti et al. (2000) with the addition of the “graphite bump” at $\sim 2\,175$ $\mathrm{\SIUnitSymbolAngstrom}$ (e.g. Fischera & Dopita, 2011; Ilbert et al., 2009). The $E(B-V)$ values range from 0 to 0.5. IGM absorption is implemented following the analytical correction of Madau (1995).

As strong nebular emission (such as [OII] or H$\alpha$ lines) can significantly increase the flux measured in a photometric filter, nebular emission lines were considered in the $z_{\mathrm{phot}}$ computation of real data (L16). On the other hand, such options are disabled when LePhare is run on Horizon-AGN galaxies, since nebular emission lines are not modelled for them (Section 2.3).

Each template in the LePhare library is fit to the virtual photometry of the Horizon-AGN galaxies. The code computes the goodness of fit ($\chi^{2}$) for each redshift solution and their likelihood ($\mathscr{L}$). The $z_{\mathrm{phot}}$ estimate for a given galaxy is defined as the median of the $\mathscr{L}(z)$ distribution, and the $1\sigma$ uncertainty ($\sigma_{\mathrm{z}}$) is the interval enclosing 68 percent of its area (see Section 3.1.2). In order to improve the SED-fitting performance, a mild luminosity prior is also applied to reduce the fraction of catastrophic outliers. This prior is based on the observed luminosity function in the rest-frame $B$ band (e.g. Ilbert et al., 2006; Zucca et al., 2009; López-Sanjuan et al., 2017). Given the extremely low number density expected in the (extrapolated) bright end of the luminosity function, our codes excludes solutions with absolute magnitude $M_{B}>-24.2$. No other prior on the redshift distribution is applied in LePhare.

As done in L16, fluxes rather than magnitudes are used when running LePhare. This allows to deal robustly with faint or non-detected objects.

Systematic offsets

With real datasets, an important aspect of the $z_{\mathrm{phot}}$ computation is the derivation of systematic offsets which are applied to match the predicted magnitudes and the observed ones (Ilbert et al., 2006) based on the spectroscopic sub-sample. This calibration is designed to empirically correct both for the incomplete template library and possible systematics in the galaxy magnitude extraction. However this calibration might be biased because it relies on a spectroscopic subsample. As described in Appendix C, we test the computation of systematic offsets in the COSMOS-like catalogue by using a sub-sample of galaxies matching the spectroscopic catalog on COSMOS. We find that there is no need for this calibration in the simulated catalogue. The offsets introduced by an imperfect knowledge of the templates are therefore negligible in Horizon-AGN444In fact, the virtual photometry of the simulated galaxies presents less diversity than in real datasets, because it is built by the mean of BC03 SSP models for all the galaxies, which might be a reason why these offsets are negligible..

2.4.2 Stellar mass and star formation rate

Stellar mass and star formation rate (SFR) are then derived using another template library built by using the BC03 model, in the same way as for the COSMOS2015 catalogue. In this second run, similarly to what was done in COSMOS2015 for computational reasons, only two extinction laws are used (Arnouts et al., 2013 and Calzetti et al., 2000). The initial mass function (IMF) is assumed to be Chabrier (2003)’s, and the stellar metallicity of each template to be either solar ($Z_{\odot}$) or subsolar ($0.4\,Z_{\odot}$).

In our library, the SFHs used to build the SEDs are parametrised with an analytic equation. It can be exponentially declining, i.e. ${\rm SFR}(t)\propto e^{-t/\tau}\,.$ An alternative definition is the “delayed" star formation history, more suitable to model a galaxy with gas infall: ${\rm SFR}(t)\propto e^{-t/\tau}{t}/{\tau}\,.$ In both cases $\tau$ is the $e$-folding time, varying between one tenth and several Gyr555The following timesteps are used: $\tau=0.1,0.3,1,3,4,30$ Gyr for the exponentially declining star-formation histories, and $\tau=1,3$ Gyr for the delayed star formation histories.. In the latter case, $\tau$ also represents the galaxy age (since its formation) at which the SFR peaks. In particular, such SFHs cannot reproduce multiple bursts of star formation. The SFH models start forming stars from $t=0$. When building the library of templates, $t$ is sampled from a few hundreds Myr to the age of the Universe at a given redshift, with up to 44 steps at $z=0$.

3 SED-fitting performance: present surveys

3.1 Photometric redshifts

Let us first investigate our ability to recover photometric redshifts from the photometry, using the simulated COSMOS-like catalogue.

3.1.1 Comparison between $z_{\mathrm{sim}}$ and $z_{\mathrm{phot}}$

The accuracy of our SED fitting method can be first tested by comparing galaxy redshifts in the lightcone ($z_{\mathrm{sim}}$) to those obtained by LePhare666 $z_{\rm sim}$ includes galaxy peculiar velocity.. Fig. 1 (left panel) presents such a comparison for galaxies with $K_{\mathrm{s}}<24.7$. In addition to the magnitude cut, pathological cases are also excluded with reduced chi-square values $\chi^{2}_{\mathrm{red}}>10$, which represent $<0.1$ per cent of the whole sample (namely, 387 out of 541,555 $K_{\mathrm{s}}$-selected galaxies). This kind of $\chi^{2}$ selection is also applied in real surveys and removes similar fractions of problematic objects (Davidzon et al., 2017; Caputi et al., 2015).

Overall, we do not find significant systematics affecting LePhare redshift estimates. Despite the simplistic implementation of dust extinction, the limited number of templates and the fact that they have been calibrated to represent the observed universe, our recipe captures the main features of the simulated galaxy SEDs, and recovers their redshifts with a precision comparable to what is achieved with the COSMOS2015 catalogue. In our simulation, the normalized median absolute deviation (NMAD, Hoaglin et al., 1983) for the entire sample is $1.48\times\mathrm{median}(|\Delta z|)/(1+z_{\mathrm{sim}})=0.031$. The fraction of outliers, defined as objects with $|\Delta z|>0.15\,(1+z_{\mathrm{sim}})$, is $\eta=4.6$ per cent (Fig. 1, left panel). In the real survey, with the same cut at $K_{\mathrm{s}}<24.7$, a comparison between photometric and spectroscopic redshifts results yields $\mathrm{NMAD}=0.013$ and $\eta=3.1$ per cent respectively. Such difference between simulation and observation is explained because the spectroscopic samples are not representative (as detailed below).

Specific sub-samples of galaxies may show a smaller scatter than the global one. Galaxies’ $z_{\mathrm{phot}}$ precision depends on their $S/N$, which in turn correlates with apparent magnitude, therefore it is useful to estimate the NMAD as a function of the latter. Fig. 1 (right panel) presents the $z_{\mathrm{phot}}$ vs $z_{\mathrm{spec}}$ diagram after dividing Horizon-AGN galaxies in four $i^{+}$-band magnitude bins. NMAD and $\eta$ of each sub-sample are reported in Table 2, along with the corresponding metrics obtained in L16 for COSMOS2015. The precision of Horizon-AGN $z_{\mathrm{phot}}$ is comparable to that found in the real survey. However the fraction of catastrophic failures is systematically larger in COSMOS2015. Such a discrepancy can be explained by observational uncertainties related to image source extraction (e.g., confusion noise) not considered in Horizon-AGN. In addition, it should be remembered that the observed $\eta$ and NMAD are provided only for a sub-sample of high-confidence spectroscopic galaxies, which are generally biased towards bright objects.

Comparison with a $z$COSMOS-like sub-sample

The COSMOS2015 $z_{\mathrm{phot}}$ precision was assessed by using data collected during various spectroscopic campaigns (see Table 5 of L16). One of the most important is the zCOSMOS survey, accounting for almost half of high-quality galaxy spectra at $z<3$ (Lilly et al., 2007)777The full description of the zCOSMOS-Deep sample characteristics and the evaluation of the redshift estimation performance will be published in a future paper (Lilly et al. in prep).. To reproduce a similar subset in Horizon-AGN, a sub-sample of galaxies in the lightcone is identified using selection criteria similar to zCOSMOS. We randomly extract mock galaxies at $i^{+}<22.5$ until we match both the magnitude and redshift distributions of the zCOSMOS-Bright sample (Lilly et al., 2007), and do the same at $B<25$ to mimick the zCOSMOS-Deep (Lilly et al. in prep). The zCOSMOS-Deep selection aimed at enforcing a target selection at $z>1.5$, but some faint galaxies at lower redshift were also observed. Yet, those interlopers are not included in the zCOSMOS-Deep mock sample by applying a sharp $1.5<z<2.5$ cut. It has to be noted that a large fraction of the catastrophic failures seen at $z_{\mathrm{spec}}<1$ in the real data correspond indeed to the selection that is not replicated in the simulation. The result is a set of about 5,000 galaxies for which pseudo-spectroscopic measurements are created by perturbing $z_{\mathrm{sim}}$ with a Gaussian random error having $\sigma=0.0004(1+z)$.888 The standard deviation of the Gaussian random error corresponds to the 1$\sigma$ uncertainty of zCOSMOS-Bright galaxy estimated by repeated measurements (Lilly et al., 2007). Given the order of magnitude of $z_{\mathrm{phot}}$ error, that pseudo-spectroscopic perturbation is negligible in the following analysis.

Figure 2 shows the $z_{\mathrm{phot}}$ versus $z_{\mathrm{spec}}$ comparison in Horizon-AGN (upper panel) and COSMOS2015 (lower panel). Since stars are not included in our simulation the analysis is restricted to $z>0.1$ to avoid both stellar interlopers present in the real survey and bright COSMOS2015 galaxies erroneously classified as stars. With such a zCOSMOS-like selection, the NMAD measured in Horizon-AGN is in excellent agreement with COSMOS2015 ($0.0080$ and $0.0081$ respectively). On the other hand, catastrophic errors are more numerous in the real sample, as also found in the previous test (see Table 2). Appendix D.2, focuses on the 22 outliers with $|\Delta z|>0.3\,(1+z_{\mathrm{spec}})$ to understand this difference. In most of the case, the failure arises either because of uncertain photometry (fragmented or blended objects) or spectroscopic misidentification.

Impact of IGM absorption

The Horizon-AGN lightcone allows us to quantify the importance of correctly accounting for IGM absorption, by comparing the $z_{\mathrm{phot}}$ estimate computed from the IGM and IGM-free (i.e. turning off IGM both in the photometry computation and at the SED-fitting stage) versions of the catalogue. This test is carried out with the dust-free version of the catalogue. In Horizon-AGN, IGM absorption is implemented along the line-of-sight of each galaxy knowing the foreground HI distribution, while in LePhare, absorption due to the intervening IGM between the galaxy and the observer is taken into account by applying an average correction as a function of redshift based on an analytical relation (Madau, 1995) 999As explained in Appendix A, there is a slight discrepancy between the average IGM absorption in Horizon-AGN and the correction implemented in LePhare, the latter being stronger than the former. However the present work does not aim at correcting this discrepancy as it is likely to also happen when fitting the photometry of real galaxies.. Overestimating the IGM correction (or neglecting the line-of-sight variability of IGM opacity) might impact the performance of $z_{\mathrm{phot}}$ estimate for distant galaxies (as suggested in Thomas et al., 2017). IGM absorption plays a role mostly at $z>2$, with a more dramatic attenuation of galaxy photometry from $z\sim 3$ as illustrated by Fig. 15.

The global redshift accuracy estimated by the NMAD is impacted at a below per cent level, only at the very faint end of the galaxy population. Interestingly, implementing IGM absorption slightly helps constraining the redshift of faint galaxies, when averaging it over the entire redshift range (In the bin $25<i<26$, NMAD$=0.049$ and $\eta=8.7$ % in the IGM-free version, while NMAD$=0.045$ and $\eta=7.6$ % with IGM). However, at $z>3$, the fraction of outliers populating the clump located at $\left[z_{\rm sim}>2.5\right]\cap\left[z_{\rm obs}<1.5\right]$ strongly increases when IGM is included: at $24<i<25$, only 7% of the existing outliers populate this region in the no-IGM case, while they are 58% with IGM. At $25<i<26$, in this region, they are 16% and 47% in the no-IGM and IGM cases respectively. This population of outliers also occurs when fitting observed population of galaxies (see Fig. 11 in Laigle et al., 2016) and our work suggests therefore that the way IGM is accounted for is important to mitigate it.

Impact of medium-bands

Let us now quantify the improvement of $z_{\mathrm{phot}}$ estimates due to medium-band photometry. The inclusion of these bands in a deep extragalactic survey is very useful to better constrain the redshift from spectral features occuring in the optical wavelength range (Lyman and Balmer breaks depending on the redshift, nebular emission lines in real datasets) but expensive in exposure time. It is therefore important to check whether it is worthwhile. Even though the major advantage of those filters is to find contribution from nebular emission lines (not implemented in our virtual magnitudes), they should also in principle help to constrain the galaxy continuum. Indeed, when medium-band filters are removed, the $z_{\mathrm{phot}}$ precision degrades considerably (see Table 2). For instance, the NMAD of bright objects ($22<i<23$) is degraded by almost a factor 3, going from $\mathrm{NMAD}=0.008$ (when medium-bands are included) to 0.023. At the fainter magnitudes ($24<i<26$) the difference is less remarkable, but nonetheless the absence of medium-bands results in a $z_{\mathrm{phot}}$ scatter larger by $30\!-\!40$ per cent. The outlier fraction, especially at $i>25$, also increases (see Table 1).

3.1.2 Photometric redshift errors from marginalized likelihood

When working on extragalactic surveys, an accurate knowledge of the redshift probability distribution function (PDF) is instrumental to remove observational uncertainties in galaxy statistics. For instance, the PDFs of a certain set of galaxies can be used to correct their luminosity function for the so-called Eddington bias (see Schmidt et al., 2014). Therefore, it is crucial to verify whether the SED fitting produces a reliable PDF($z$) for a given galaxy. In a Bayesian framework this PDF is the posterior probability distribution, proportional to the product of the prior distribution and the marginalized likelihood $\mathscr{L}(z)$.

From $\mathscr{L}(z)$101010 $\mathscr{L}\propto\exp(-\nicefrac{{1}}{{2}}{\chi^{2}})$,

with {\chi^{2}(z)}={\sum}_{\rm filters\,i}{(F_{{\rm obs}\,i}-F_{{\rm SED}\,i}(z,T))^{2}}\big{/}{\sigma^{2}_{i,\rm obs}}, where $F_{\rm SED\,i}(z,T)$ is the flux predicted for a template $T$ in the filter $i$ at $z$, $F_{{\rm obs}\,i}$ is the observed flux in the filter $i$ and $\sigma_{\rm obs}$ is the associated uncertainty., LePhare computes $\sigma_{\mathrm{z,fit}}$, namely the photometric redshift 1$\sigma$ error. This is defined as the redshift interval, centred at $z_{\mathrm{phot}}$, that encloses $68.27$ per cent of the $\mathscr{L}(z)$ area.

Reliability of redshift 1$\sigma$ errors from SED-fitting

The simulation provides the true $z_{\mathrm{phot}}$ uncertainties directly from the difference between SED-fitting estimates and $z_{\mathrm{sim}}$. Hence it can be checked whether the error bars provided by LePhare actually represent the redshift 1$\sigma$ uncertainty. Previous tests with spectroscopic data suggest that they could be under-estimated (see e.g. Dahlen et al., 2013, L16).

Fig. 3 shows the evolution of the median of $\sigma_{\mathrm{z,fit}}$ as a function of $z_{\mathrm{sim}}$, for galaxies in three different bins of $i^{+}$-band magnitude. From these bins galaxies with degenerate redshift solutions are excluded, i.e. whose $\mathscr{L}(z)$ function shows two peaks.111111 LePhare automatically identifies a galaxy with two acceptable redshift solutions when the secondary peak includes $>2$ per cent of the integrated $\mathscr{L}(z)$ distribution. In such a case the integrated 68.27 per cent of the $\mathscr{L}(z)$ area is strongly skewed towards the secondary solution and the $\sigma_{\mathrm{z,fit}}$ does not represent the pure statistical error but also includes systematics. Nevertheless, if those galaxies are re-introduced in the sample, results shown in Fig. 3 remain the same at $i^{+}<25$ and change less than 20 per cent in the faintest bin.

We first find that $\sigma_{\mathrm{z,fit}}$ values increases with magnitude, as shown in Fig. 3. Such behaviour is expected since fainter galaxies have a lower $S/N$ ratio and a lower constraint on the SED fit. We also find a redshift dependency of the uncertainties, mainly due to the different efficiency of optical and NIR photometry: the former being deeper, including medium-bands and can tightly constrain the Balmer break at $z<1.3$ ; whereas at $z>1.5$ the break is entirely shifted in the NIR regime, which is sampled with fewer, less sensitive bands. At $z>2.5$ the $\sigma_{\mathrm{z,fit}}$ amplitude slightly decreases as optical blue bands start to constrain the Lyman break position. Overall, LePhare predicts a symmetric scatter, with upper and lower errors such that $\sigma_{\mathrm{z,fit}}^{+}\simeq-\sigma_{\mathrm{z,fit}}^{-}$.

To establish whether the uncertainties derived by LePhare are reliable, they are compared to the 1$\sigma$ errors directly retrieved from the simulation ($\sigma_{\mathrm{z,true}}$). Using the same $i^{+}$-selected galaxies for which $\sigma_{\mathrm{z,fit}}$ was computed, we measure $\sigma_{\mathrm{z,true}}$ by means of their $\Delta z=z_{\mathrm{sim}}-z_{\mathrm{phot}}$ distribution, finding the interval that includes 68.27 per cent of it. Despite some noise, relatively good agreement is found between $\sigma_{\mathrm{z,fit}}$ and $\sigma_{\mathrm{z,true}}$, an indication that the uncertainty provided by LePhare is generally a good proxy of the actual 1$\sigma$ error dispersion (Fig. 3). However at bright magnitudes ($i^{+}<24.5$) and for $1<z_{\rm sim}<2.5$, $\sigma_{\mathrm{z,fit}}$ is generally underestimated compared to $\sigma_{\mathrm{z,true}}$. As shown in Appendix D.1, this underestimation might be due either to the underestimation of photometric errors, or to a lack of representativeness of the set of templates for this galaxy population, making $\mathscr{L}(z)$ too spiky around the median $z$. A similar trend was already discussed in L16. From a comparison with spectroscopic redshifts (their figure 13), L16 suggested indeed that the 1$\sigma$ uncertainties produced by LePhare were under-estimated. As a consequence the authors proposed a magnitude-dependent boosting factor ($f_{\sigma}$) that would enlarge $\sigma_{\mathrm{z,L16}}$ so that $\sim\!68$ percent of the COSMOS2015 $z_{\mathrm{spec}}$ would fall within $\sigma_{\mathrm{z,L16}}\times f_{\sigma}$ from $z_{\mathrm{phot}}$. This boosting factor121212 For galaxies with $i^{+}>20$, $f_{\sigma}=0.1\times i^{+}-0.8$; $f_{\sigma}=1.2$ otherwise. was however constant with redshift and increasing with magnitudes. From our analysis, it appears that this factor would generally over-correct the $z_{\mathrm{phot}}$ errors at faint magnitudes, but particular care should be given to the errors of bright galaxies at $1<z<2.5$.

Finally, note that the global behaviour of $\sigma_{\mathrm{z,fit}}$ as a function of magnitude and redshift depends on the photometric baseline available, and does not necessarily hold for different configurations (see e.g. Section 4).

zphot error comparison between COSMOS-like and COSMOS2015

Let us now compare the Horizon-AGN redshift uncertainties with those of the COSMOS2015 galaxies ($\sigma_{\mathrm{z,L16}}$) as calculated in L16 by means of LePhare. The method is the same as that applied to simulated galaxies, based on the marginalized $\mathscr{L}(z)$. Fig. 4 shows the median $\sigma_{\mathrm{z,L16}}$ as a function of redshift, for galaxies with $23.0<i^{+}<23.5$ and $25.0<i^{+}<25.5$. In the Figure $\sigma_{\mathrm{z,fit}}$ is also reported, in bins of $z_{\mathrm{phot}}$ instead of $z_{\mathrm{sim}}$ to allow the comparison with real data.

The trend of $\sigma_{\mathrm{z,L16}}$ and $\sigma_{\mathrm{z,fit}}$ are remarkably similar, both showing an increase between $z=1.5$ and 2, as discussed above. However, COSMOS2015 galaxies have a median $\sigma_{\mathrm{z}}$ that is about 50 per cent larger than Horizon-AGN. This difference is likely to be driven either by simplifications in the modelling of the photometry itself, or by failures in the photometry extraction of real data which are not modelled in the simulated catalog. On the one hand, the simulated photometry includes indeed less variety than realistic galaxies: we use a single and constant IMF, a constant dust-to-metal ratio, single SSP model for stellar mass losses, and emission lines are not modelled. These simplifications can naturally reduce the scatter of the $z_{\rm phot}$ estimate compared to the observed catalogue. On the other hand, although magnitude errors in the simulated catalogue are implemented consistently with COSMOS2015, “catastrophic detections" (such that blended or fragmented objects) and systematics (astrometry calibration issues, background removal, lack of modelling of the PSF variation within the field, mis-centering of the galaxies …etc.) are not considered in the virtual photometry. These errors will propagate in the $z_{\mathrm{phot}}$ errors. Exploring these effects is out of the scope of this paper and will be the topic of a future work.

3.2 Physical quantities: mass and star-formation rate

3.2.1 Stellar mass estimate

The overall comparison between the intrisic stellar masses ($M_{\rm sim}$) and those retrived via SED-fitting ($M_{\rm phot}$) is shown in Fig. 5 in photometric redshift bins up to $z\sim 4$. The observed stellar masses are in very good agreement with the intrinsic ones. The left and right panels of Fig. 6 present respectively the median and the dispersion around the median $\sigma_{M}$ of $\log M_{\rm phot}/M_{\rm sim}$ as a function of $\log M_{\rm phot}$ in different redshift bins. The dispersion around the median value being potentially asymetric, we measure $\sigma^{+}_{M}$ and $\sigma^{-}_{M}$ as the value which encloses 34% of the full population respectively above and below the median. The values $\sigma^{+}_{M}$ and $\sigma^{-}_{M}$ are displayed as the upper and lower lines on the right panels of Fig. 6.

Impact of $z_{\rm phot}$ uncertainties

In order to determine how much of the trend is driven by the propagation of uncertainties from the photometric redshift estimation, stellar mass computation through SED-fitting is reproduced in a second step while fixing the redshift at $z_{\rm sim}$ instead of $z_{\rm phot}$. The dashed lines in the top panels of Fig. 6 correspond to the median and dispersion of $\log M_{\rm phot}/M_{\rm sim}$ using the photometric redshifts in the computation of $M_{\rm phot}$, while the solid lines correspond to the same quantities but using the intrinsic redshifts from the simulation in the computation of $M_{\rm phot}$. Comparing the solid and dashed lines therefore allows to quantify the impact of the photometric redshift uncertainty propagation in the stellar mass computation, which is very limited. Finally, it should be noted that the dashed lines provide a direct comparison with observations, as the galaxy population is split in bins of $z_{\rm phot}$, and the computation includes both dust and redshift uncertainties. As a complement, the top panel of Fig. 19 can also be compared with Fig. 5. Overall, the propagation of $z_{\rm phot}$ uncertainties has only a small impact on retrieving stellar mass. The scatter is relatively stable over the redshift and mass ranges and is generally smaller than 0.1 dex. $M_{\rm phot}$ is preferentially underestimated up to $z\sim 2$ by at most $\sim$0.12 dex. At $z>2$ and $\log M_{\rm phot}>9.5$ the trend tends to reverse and $M_{\rm phot}$ ends up slightly overestimated.

Impact of dust attenuation

In order to isolate the role played by attenuation in driving this behaviour, the same computation is performed on the attenuation-free catalogue, while using the intrinsic redshift from the simulation. The median and dispersion of $\log M_{\rm phot}/M_{\rm sim}$ in the dust-free case are shown in the bottom panels of Fig. 6. Without attenuation, one is left with a weak systematic underestimation of the stellar mass especially at low redshift, which can be driven either by the too simplistic (single-burst) SFHs (see e.g. Leja et al., 2018) or the discretisation of the metallicity in galaxy template (Mitchell et al., 2013). At higher redshift, these assumptions are more likely to correctly represent the actual SFHs and metallicity distributions. It should be noted that the impact of attenuation (overestimation) and of these simplified assumptions (underestimation) tends to compensate each other. For example, at $1.5<z<2$, $M_{\rm phot}$ is closer to $M_{\rm sim}$ when dust is included.

3.2.2 Star formation rate estimate

It is known that SFR derived from SED-fitting has to be considered with caution, given the simplistic shape of the SFHs assumed in the templates, which for instance cannot account for recent bursts of star formation. Ilbert et al. (2015) predicted an overall offset of 0.25 dex and a scatter up to 0.35 dex, from a comparison of SFR derived on the one hand from SED-fittig and on the other hand from IR+UV flux.

The bottom panel of Fig. 5 illustrates this point, as it presents the overall comparison between the intrinsic star-formation rate (SFRsim) and that derived from SED-fitting (SFRphot) in photometric redshift bins up to $z\sim 4$. Up to $z\sim 3$, SFRphot presents a bimodal behavior, with a systematic underestimation for a large fraction of the galaxy population up to $z\sim 1.5$ and an overestimation of low-mass galaxies above. Moving towards high-$z$, the bi-modality tends to disappear but the scatter remains very large. The left and right panels of Fig. 7 present the median of $\log{\rm SFR}_{\rm phot}/{\rm SFR}_{\rm sim}$ and the dispersion around the median $\sigma_{\rm SFR}$. The dispersion around the median value being potentially asymetric, we measure $\sigma^{+}_{\rm SFR}$ and $\sigma^{-}_{\rm SFR}$ as the value which encloses 34% of the population respectively above and below the median. The values $\sigma^{+}_{\rm SFR}$ and $\sigma^{-}_{\rm SFR}$ are displayed as the upper and lower lines on the right panels of Fig. 7. The median evolves between -0.6 at $\log{M}_{*}/{\rm M_{\odot}}>10.5$ and $z<1.5$ and 0.3 at $\log{M}_{*}/{\rm M_{\odot}}<10.$, while the scatter varies from $\sim$ 0.6 at $0.8<z<1.1$ to $\sim$ 0.15 at high-$z$.

Impact of $z_{\rm phot}$ uncertainties

In order to determine the role of redshift uncertainties in driving the trend, the SFR is also computed while fixing the galaxy redshift at their intrinsic values $z_{\rm sim}$ instead of $z_{\rm phot}$. The dashed lines in the top panels of Fig. 7 correspond to the median and dispersion of $\log{\rm SFR}_{\rm phot}/{\rm SFR}_{\rm sim}$ using the photometric redshifts in the computation of ${\rm SFR}_{\rm phot}$, while the solid lines correspond to the same quantities but using the intrinsic redshifts from the simulation in the computation of ${\rm SFR}_{\rm phot}$. Comparing the solid and dashed lines therefore allows to quantify the impact of the photometric redshift uncertainty propagation in the SFR computation. As a complement, the top panel of Fig. 20 can also be compared with the bottom panel of Fig. 5. This comparison highlights that working with the simulated redshift removes the bimodality in the lowest redshift range $0.2<z<0.5$, which therefore is driven by redshift degeneracies. However the bimodality remains in all the other redshift bins.

Impact of dust attenuation

One expects the impact of dust on the precision of SFR to be much stronger than on the mass. Indeed it attenuates preferentially blue bands, which are a tracer of recently formed stars, hence directly connected to the SFR. The comparison of the SFR in the run with and without attenuation confirms this fact. The median and dispersion of $\log{\rm SFR}_{\rm phot}/{\rm SFR}_{\rm sim}$ in the dust-free case are shown in the bottom panels of Fig. 7, and the overall comparison of ${\rm SFR}_{\rm phot}$ versus ${\rm SFR}_{\rm sim}$ without dust is shown in the bottom panel of Fig. 20. At $z>2$ and in a dust-free Universe, the SFR is very well recovered without any bimodal behaviour, while it tends to be slightly underestimated at $z<2$. As for the mass, this remaining underestimation is likely to be driven by the oversimplified SFH and metallicity underlying models, which cannot fully render the complexity of low redshift galaxy SEDs. On the contrary, when dust is accounted for in the virtual Universe, a bimodal behaviour appears, due to the SFR-dust degeneracy. The cause of this trend is investigated in Appendix B.3. In particular, there is a direct correlation between the attenuation in the rest-frame $NUV$ and the SFR. Overestimating the attenuation $A_{NUV}$ at the SED-fitting stage yields an overestimation of the SFR. None of the two extinction curves used in LePhare are a good fit for the one used in Horizon-AGN, and this discrepancy is likely to be the main driver of this bimodality.

3.3 Performance of current surveys: summary

The virtual photometric catalogue, calibrated to mimic COSMOS2015, has allowed a fully consistent test of the performance of LePhare when computing galaxy redshifts, masses and SFR from broad-band photometry. We summarise below our main findings:

• •

In the same configuration as COSMOS2015, photometric redshifts are retrieved with the same overall precision (as estimated from NMAD and 1$\sigma$ uncertainties) in the virtual dataset as in the observed one. When binning the datasets in apparent magnitudes, the simulation yields as precise estimates as the observations. In particular, the 1$\sigma$ uncertainties measured from $\mathscr{L}(z)$ represent on overall a good estimate of the intrinsic errors (as measured from the difference between $z_{\rm sim}$ and $z_{\rm zphot}$), except for the bright galaxies at $1.<z<2.5$ which have in general their errors underestimated. However the averaged correcting factor for redshift errors proposed in earlier works (see e.g. Ilbert et al., 2013; Laigle et al., 2016) would generally over-correct the errors of faint galaxies. Redshift errors are generally smaller in Horizon-AGN as in COSMOS, as the mock catalog does not include systematics in the extraction of photometry from noisy images, and presents less diversity in terms of photometry. For the same reasons, although the simulated catalogue allows us to retrieve the overall redshift distribution of the catastrophic population of outliers, it systematically underestimates their fraction. Intermediate bands allows to improve redshift accuracy;

• •

Stellar masses are very well recovered, despite the use of single-burst SFH model and discrete metallicity in the SED-fitting templates, which do not a priori represent the complex SFHs of simulated galaxies. Only a small underestimation of at most $\sim 0.12$ dex persists at low redshift, and an overall scatter of the order of 0.1 dex. Conversely, dust induces a slight overestimation of the mass at high redshift. The impact of redshift uncertainties in driving the scatter is very limited;

• •

Unsurprinsgly, the SFR directly derived from the SFH are a quite poor proxy of the intrinsic SFR. The simplistic SFH and metallicity enrichment induce an underestimation, while the dust modelling (mainly the choice of the attenuation curve) induces a bimodality. As a result, the dispersion around the median values evolve between 0.2 and 0.6 dex.

4 SED-Fitting performance: Forecasts

In this Section, we use the mock catalogues reproducing Euclid-like, LSST-like, and DES-like photometry in order to predict the performance of these surveys in the three configurations presented in Table 1.

4.1 Redshift and mass completeness in Euclid and LSST

Let us first present the expected completeness of each survey, estimated from virtual catalogues using the intrinsic redshift, intrinsic stellar masses and total unperturbed magnitudes. The fraction of “detected” galaxies (i.e., those brighter than the magnitude limit of the survey) is measured as a function of redshift and stellar mass. The magnitude cuts correspond to those used for weak lensing galaxy selection in Euclid ($riz<24.5$) and LSST ($i<25.3$); the completeness at $H<24$ is also computed, namely the 5$\sigma$ detection limit expected at the completion of the Euclid mission. It is argued in the following that the latter threshold is the most suited for galaxy evolution science cases.

Fig. 8 compares the intrinsic redshift distribution of all the objects in the lightcone with $M_{\mathrm{sim}}>10^{9}{\rm M}_{\odot}$ to the sub-sample of galaxies detected in the Euclid-like catalogue, in the case of either $H$-band or optical selection. The $i$-band sample expected to be detected in LSST is also included, showing a redshift distribution similar to the Euclid $H$-band selected sample. The cut applied in the $riz$ band results in a lower completeness, below 50 per cent already at $z_{\mathrm{sim}}=1$. Both the $H<24$ and $i<25.3$ completeness drop below 50 per cent at $z_{\mathrm{sim}}=1.5$. The fact that Fig. 8 does not include galaxies below $10^{9}M_{\odot}$ has a negligible impact at $z>0.5$ because such a low-mass population is generally fainter than the magnitude limits considered here (see L16). Even at $z<0.5$, where all selections are complete, one does not expect the addition of $M_{\mathrm{sim}}<10^{9}M_{\odot}$ galaxies to significantly impact our results.

The stellar mass completeness ($M_{\rm lim}$) is shown in Fig. 9. This is a lower limit, as a function of redshift, above which $>$90 per cent of galaxies are detected in the selection band ($i$, $riz$, or $H$). The results are summarized in Table 3, along with a less conservative $M_{\rm lim}(z)$ threshold (50 per cent completeness).

As shown in Fig. 9, the 90 per cent stellar mass completeness of LSST galaxies is below $10^{9}\,{\rm M}_{\odot}$ at $z<0.5$. It increases at $z>0.5$, because of dimming, reaching $M_{\rm lim}=1.5\times 10^{10}{\rm M}_{\odot}$ at $z\sim 2$. This threshold decreases at higher redshift as galaxies within our mass range in the early universe have higher SFRs (see e.g., Speagle et al., 2014) so they become brighter in the rest-frame UV probed by the $i$ band. Conversely, at $z\lesssim 2$ the SFR starts to decline while more stellar mass is assembled, allowing an easier detection in the $H$ band (see a similar discussion in Davidzon et al., 2017, for a Spitzer/IRAC selected sample).

For galaxy evolution studies relying on the Euclid photometry one should prefer an $H$-band selection, instead of the nominal $riz<24.5$, if the scientific goals require a sample highly complete in stellar mass. Given an average galaxy SED, the Euclid optical selection would correspond to a cut at $H\sim 23$, while the survey will go deeper in the NIR by about 1 mag. Therefore, in the following sections, the analysis is carried on with the $H<24$ selected sample. Modulo observational uncertainties, the nominal mass completeness for the Euclid $H<24$ sample is well fit by the function $M_{\mathrm{lim}}(z)=4.5\times 10^{8}(1+z)^{2.4}\,M_{\odot}$, reaching a maximum of $2\times 10^{10}M_{\odot}$ at $z=4$ in the studied redshift range (see Fig. 9). Although the limit based on $riz$ detections is generally higher (e.g., $8\times 10^{10}M_{\odot}$ at $z\sim 2$) it is enlightening to note that it starts to decline at $z>2$, as already discussed for the optical selection in LSST131313Recall that the star forming main sequence in Horizon-AGN reproduces well the observed one at $z\sim 4$, but at lower redshifts the simulation underestimates it by $\lesssim 0.3$ dex (see figure 3 in Kaviraj et al., 2017). Consequently, observed galaxies at $2<z<4$ are expected to be brighter in the rest-frame UV, which would therefore enhance their $riz$ magnitude by at most $\sim$0.5 mag. Such a magnitude offset corresponds to a stellar mass limit approximately $\sim$0.2 dex lower than that displayed in Fig. 9. A $H<24$ selection would still provide a higher completeness. The underestimation of the SFR in the virtual catalogue being independent of mass, this remark is also valid for the LSST catalogue. .

4.2 Forecasts for Euclid and LSST photometric redshifts

Photometric redshifts are computed in the same way as in the COSMOS-like case (Section 2.4). The performance of SED-fitting in the Euclid+DES, LSST-only and Euclid+LSST configurations for $z_{\rm phot}$ estimation is presented in Fig. 10. The sample is split in different $i$-band bins ($i$ taken from either LSST or DES photometry). The results in terms of NMAD and catastrophic outliers are presented in Table 4, for both $H$ and $i$ magnitude bins.

In all configurations, the usual population of outliers at high redshift and faint magnitudes $i>24$ is found (see the discussion in Section 3.1 about the impact of IGM absorption). Euclid (combined with optical DES photometry) performs relatively well at bright magnitudes ($i<23$). However, because of the lack of blue optical bands to constrain the Balmer break, the accuracy at very low redshift ($z<0.5$) is lower than in COSMOS, even for bright galaxies. At fainter magnitudes, the main limitation of this configuration is the shallow depth of the survey.

In the LSST configuration, without NIR photometry, a large fraction of catastrophic outliers is present at $1.2<z<2.5$ at all magnitudes with a relatively symmetric patterns. The reason is the same as for the increase of $z_{\rm phot}$ uncertainties in this redshift range (see e.g. Fig. 4). At this redshift, the Balmer break is not constrained anymore by the optical bands and enters NIR. Without NIR bands to properly constrain its position, determining the redshift is challenging (see also Fotopoulou & Paltani, 2018; Gomes et al., 2018). The situation improves from $z>2.5$ when the Lyman break enters the optical bands.

The LSST-like catalogue alone performs therefore less well than the LSST+Euclid one. While the NMAD improves only by $\sim$ 1 percent in the faintest magnitude bin, adding the Euclid NIR bands allows to reduce the fraction of outliers by a factor 2. The benefit to combine Euclid and LSST is further illustrated in Fig. 11. It presents the photometric errors computed by LePhare in the LSST-only and Euclid+LSST configurations in different $H$ bins. The $z_{\rm phot}$ errors dramatically improve in the redshift range $1.5<z<2.5$, especially in the faintest magnitude bins. This is true as long as galaxies are bright enough in the NIR. At fainter magnitudes $H\gtrsim 25$, Euclid is not deep enough to properly constrain the Balmer break. From the Euclid perspective, adding the LSST optical bands to the Euclid+DES baseline considerably decreases redshift uncertainties and fraction of outliers (see Table 4) especially at faint magnitudes, at is provides deeper photometry in the $g$, $r$, $i$ and $z$ band. Furthermore the addition of the $u$-band is considerably useful from $z\sim 2.5$, when the Lyman-break enters the $u$-band.

4.3 Forecast for Euclid and LSST stellar masses

Figure 12 presents the overall comparison between intrinsic and reconstructed stellar masses in five photometric redshifts bins between 0.2 and 3, in the Euclid+DES, LSST-only and LSST+Euclid configurations. The number of objects in each redshift bin varies as the performance of $z_{\rm phot}$ estimation varies from one configuration to the other.

With the LSST-like catalogue, the performance is much poorer than with the Euclid+DES or Euclid+LSST configurations, with a very large scatter from $z>0.7$. Indeed, without NIR photometry, the stellar mass will be determined on the basis of the photometric filters which trace the young stellar populations. For example, the massively star-forming galaxies at high redshift (e.g. the massive galaxies in the bin $1.7<z_{\rm phot}<2.2$) will generally get their mass overestimated, which drives the very large scatter above the median. On the contrary, passive galaxies generally get their mass underestimated (e.g. in the bin $0.7<z_{\rm phot}<1.2$), which drives the very large scatter below the median. The resulting scatter (as defined from the RMS of $\log M_{\rm phot}/M_{\rm sim}$) can be as large as 0.5 at $z>2$ (see e.g. Figure 13).

It can also be noted that the Euclid+LSST configuration performs in general better than the Euclid+DES one at $z<2$: the additional LSST optical bands help to constrain the mass reconstruction. However more bands do not always yield a better fit. At $z>2$ and $\log M_{\rm sim}>10.5$, the scatter is larger in the Euclid+LSST configuration than in the Euclid+DES one (namely with NIR photometry only). Although conter-intuitive, this discrepancy might be a consequence of the very different depth in optical and NIR. Flux errors being much smaller in the optical, the blue part of the spectrum will provide a stronger constraint to the fit compared to the NIR. When no template can fit well both the optical and NIR photometry, the preference will be given to the optical since the errorbars are smaller. As a result, the error on the mass might be higher, because the optical part of the spectrum is a poorer proxy for stellar mass than the NIR. In the case of star-forming galaxies, it can lead to an overestimation of the mass. Removing the LSST $u$-band is in general not sufficient to bring better agreement, and the other optical bands still contribute a lot to this discrepancy.

As a summary, Fig. 13 presents the evolution of completeness with redshift and mass, along with the evolution of the median and RMS of $\log M_{\rm sim}-\log M_{\rm phot}$.

4.4 Performance of future surveys: summary

We can draw the following conclusion concerning the expected performance of future surveys.

• •

With the depth of the surveys for weak lensing galaxy selection in Euclid ($H<24$) and LSST ($i$<25.3), one can expect at $z=2$ a 90 percent completeness at $\log M_{*}/{\rm M}_{\odot}>9.9$ and $10.2$ respectively.

• •

The Euclid+DES $z_{\rm phot}$ accuracy is of the order of several percent even at low redshift and for bright objects, due to the absence of deep optical photometry to constrain the Balmer break position ($z<1.5$); the fraction of catastrophic outliers dramatically increases with fainter objects. The LSST $z_{\rm phot}$ accuracy is of the order of 2 percent at bright magnitudes ($i<24$). The absence of NIR photometry does not allow to properly constrain the Balmer break at $1.5<z<2.5$, leading to a significant fraction of catastrophic outliers in this redshift range. A dataset which would include Euclid NIR photometry in addition to the LSST optical photometry would provide a better $z_{\rm phot}$ accuracy than LSST alone or Euclid alone; it would also decrease the fraction of outliers by 2.

• •

Our LSST-like catalogue benefits a lot from NIR Euclid-like photometry for stellar mass reconstruction, reducing the scatter up to a factor 3. There is therefore a mutual benefit for LSST and Euclid to work in synergy.

5 General summary and conclusion

Using the realistic photometric catalogue extracted from the Horizon-AGN hydrodynamical simulation, we investigated the performance of SED-fitting algorithms to compute galaxy properties. Compared to previous studies, the additional value of the present modelling relies on the use of an hydrodynamical lightcone from which the photometry has been consistently post-processed. This lightcone contains a large diversity of galaxies in terms of masses and star formation activity (but also orientation with respect to the line-of-sight), over a representative cosmological volume. Galaxy photometry therefore naturally accounts for the diversity of SFHs, metallicity enrichment and dust distribution which result from the complex history of their formation, driven by the combination of pristine gas infall, stellar and AGN feedback, mergers, etc, which is consistently followed in the simulation. This lightcone also allowed for the self-consistent implementation of in-homogeneous IGM absorption within each galaxy spectrum in order to test its impact on $z_{\rm phot}$ estimation.

We used the well-calibrated COSMOS2015 dataset to assess the performance of the SED-fitting software LePhare in extracting galaxy $z_{\rm phot}$ from our mock catalogue. We also quantified our ability to estimate the corresponding $z_{\rm phot}$ uncertainties. We then estimated the biases in galaxy masses and SFR estimation through SED-fitting relying on our ability to turn on and off various physical processes in the mocks (see Table 2 and Section 3.3 for a detailed summary). Finally, we quantified the expected performance for the upcoming imaging surveys Euclid and LSST, given the available photometric baseline and expected depths (see Table 3, Table 4 and Section 4.4 for a detailed summary).

In addition to these findings specific to some survey configurations, this work has allowed to draw the following general conclusions on the process of measuring galaxy properties from their photometry:

Choice of the photometric baseline:

The added value of having medium-bands in the photometric baseline to improve redshift precision is obvious when comparing $z_{\rm phot}$ computed with and without medium-bands (Table 2). At the faint end of the galaxy population, better estimating the galactic continuum with these bands improves NMAD and $\eta$ by $\sim 50$ per cent. One can expect that the gain is even larger in the real universe, when nebular line emission can be used to constrain redshift more efficiently;

Deriving $z_{\rm phot}$ without deep optical photometry is challenging below $z<1.5$ (bottom panel in Fig. 10).

NIR photometry is mainly driving the performance of $z_{\rm phot}$ at $1.5<z<2.5$ and of stellar mass computation (e.g. compare top and bottom lines in Fig. 12). Adding optical bands to the NIR photometry helps reducing the scatter (compare middle and bottom lines in Fig. 12);

Impact of dust and IGM attenuation on the $z_{\rm phot}$ estimates

While the impact of dust is significant, it globally does not bias much the $z_{\rm phot}$ reconstruction performance with the current method to include it in SED-fitting. Overestimating the IGM absorption at the SED-fitting stage also impacts the $z_{\rm phot}$ estimate and can explain a large fraction of the population of catastrophic outliers at $z_{\rm sim}>2.5$ and $z_{\rm phot}<1.5$, as seen in real data. There are however other possible explanations for this observed population of outliers, including systematics at the stage of photometry extraction, which we do not test in the present work;

Uncertainties on $z_{\rm phot}$ estimates:

1-$\sigma$ uncertainties estimated from the SED-fitting code LePhare are a good representation of the intrinsic $z_{\rm phot}$ errors, except for the bright galaxies at $1<z<2.5$ for which the errors are generally underestimated.. The remaining discrepancies can be understood in the limits of our end-to-end pipeline (e.g. no emission lines, no modelling of the possible failures in the extraction of the photometry);

Uncertainties on $M_{\rm phot}$ and SFR estimates:

The scatter and systematics in the stellar mass and SFR computation are a combination of three effects: 1) the inherently limited SFH and metallicity enrichment pattern used to build the template library, which usually drives a global underestimation of stellar mass and SFR 2) the way dust is accounted for (in particular the choice of the dust extinction curves at the SED-fitting stage) and the degeneracy between dust and SFR in the blue bands 3) the propagation of $z_{\rm phot}$ errors in the mass and SFR estimates, which increase the scatter. The net result exhibits a complex trend, which also depends on the photometric baseline available (e.g. compare Fig. 5 and Fig. 12). As mentioned before, the impact of these effects on SFR is more dramatic, with in particular a bimodal behaviour mainly driven by dust.

Amongst possible actions to improve the redshift, mass and SFR estimates, building a template library with an additional parametrizable double burst star-formation history could temper the remaining systematic offsets at low redshift. When computing masses and SFR, one could also try to build the template library in LePhare with additional dust extinction curves, which could mitigate the bimodal behaviour in the SFR computation. Finally, it would be worth testing if allowing the mean IGM absorption to slightly vary at $z>3$ (in order to account for the line-of-sight variability of IGM opacity or the uncertainty on the model) can reduce the fraction of catastrophic outliers.

Our study does not account for systematics in the photometry extraction, and therefore our estimates of $\sigma_{z}$, $\eta_{z}$ and $\sigma_{M}$ must be understood as lower limits. However in the light of our results, we can discuss if LSST and Euclid will, at face value, fulfill their requirements.

For LSST, the redshift errors quantified from the root-mean-square scatter ($\sigma_{z}^{\rm rms}=$rms$(z_{p}-z_{s})/(1+z)$) and fraction of outliers $\eta_{3\sigma}$ (i.e the fraction of objects with $(z_{\rm p}-z_{\rm s})/(1+z_{\rm s})>3\sigma_{z}^{\rm rms}$) must be respectively smaller than 0.05 (with a goal of 0.02) and 10% at all redshifts, as specified by LSST Science Collaboration et al. (2009). For $z<3$ and $i^{+}<24.5$ (resp. $i^{+}<25.0$), we found $\sigma_{z}^{\rm rms}=0.043$ (resp. 0.060) and $\eta_{3\sigma}=1.03$% (resp. 1.50%). Note that $\sigma_{z}^{\rm rms}$ is quite sensitive to the presence of outliers, and computing it by excluding the outliers (as defined by $\eta_{3\sigma}$) yields $\sigma_{z}^{\rm rms}=0.026$ (resp. 0.032). At face value, the requirements are fulfilled, but as cautioned previously, these errors might be increased because of the systematics in the photometry extraction. When adding the Euclid photometric baseline to LSST, the errors decreases to $\sigma_{z}^{\rm rms}=0.031$ (resp. 0.044) (and $\sigma_{z}^{\rm rms}=0.025$ (resp. 0.030) when excluding the outliers).

As for Euclid, the expected redshift error $\sigma_{z}^{\rm rms}$ must be smaller than 0.05 (with a goal of 0.03) and the fraction of outliers $\eta_{0.15}$ (same definitition as in this paper) is required to stay below 10%, with a goal of 5% (Laureijs et al., 2011). In the Euclid+DES configuration, at $riz<23.5$ (resp. 24.5) we get $\sigma_{z}^{\rm rms}=0.09$ (resp. 0.17) and $\eta_{0.15}=3.45$% (resp. 9.67%). Excluding the outliers in the computation yields $\sigma_{z}^{\rm rms}=0.46$ (resp. 0.057). Photometry deeper than DES in bands narrower than the actual Euclid $riz$ filter (e.g. the photometric baseline provided by LSST) will be required to improve these performances and extend them at fainter magnitudes.

Although measuring stellar mass is not pivotal for weak-lensing based cosmology, it is of great interest for galaxy evolution science, i.e. to make the best of the legacy programs of Euclid and LSST, and therefore to fullfill their secondary science goals. In particular, the huge area of these surveys will allow to drastically decrease the statistical errors on mass functions, two-point correlation functions or any other environmental measurements (e.g. groups and clusters, cosmic web analysis). For these studies, the NIR coverage provided by Euclid will be of prime importance to extract accurate galaxy masses. On the other hand, without deep optical photometry in narrow optical filters, Euclid will be unable to separate galaxy populations from their colours, which is pivotal e.g. to study galaxy bimodality and the growth of the population of quiescent galaxies. To this end, combining Euclid and LSST would be a powerful configuration, which would benefit to both surveys and allows for the first time to address some of the most pressing questions in the field of galaxy formation today.

In following works we will pursue this discussion by exploring how redshift and mass errors propagate into one and two-point statistics, and we will quantify the effect of imperfect photometry extraction from mock images.

Acknowledgements

CL is supported by a Beecroft Fellowship, and thanks the Korean Astronomy and Space science Institute for hospitality when this work was finalized. ID was supported in part by NASA ROSES grant 12-EUCLID12-0004. OI acknowledges the funding of the French Agence Nationale de la Recherche for the project “SAGACE”. JD and AS acknowledges funding support from Adrian Beecroft, the Oxford Martin School and the STFC. HJMCC acknowledges support from the Programme national cosmologie et galaxies (PNCG) and the Domaine d’intérêt majeur en astrophysique et conditions d’apparition de la vie (DIM–ACAV+). The authors thank A. L. Serra for her suggestions to improve graphic rendering of the figures. This work relied on the HPC resources of CINES (Jade) under the allocation 2013047012 and c2014047012 made by GENCI and on the Horizon and CANDIDE clusters hosted by Institut d’Astrophysique de Paris. We warmly thank S. Rouberol for maintaining these clusters on which the simulation was post-processed. This research is part of Spin(e) (ANR-13-BS05-0005, http://cosmicorigin.org), ERC grant 670193 and horizon-UK. This research is also partly supported by the Centre National dEtudes Spatiales (CNES). This work is based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO programme ID 179.A-2005 and on data products produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium.

Appendix A Mocks additional features

Let us provide more details about how the virtual photometry has been computed.

A.1 Calibration of the dust attenuation

A.1.1 Dust-to-metal mass ratio

Dust attenuation is implemented assuming that the distribution of gas metallicity is a good proxy for the dust distribution. This computation implies to choose a value for the dust-to-metal mass ratio, i.e. to define which fraction of metals are locked into dust grains. For the sake of simplicity, the dust-to-metal mass ratio is assumed constant, though some works have shown that it could vary with redshift or within a same galaxy as a function of metallicity (e.g. Galametz et al., 2011; Mattsson et al., 2012; De Cia et al., 2013; Fisher et al., 2014). Most of the time, the implementation of dust attenuation in simulations uses a dust-to-metal ratio of 0.4 (e.g. Jonsson, 2006), which is the Milky-Way value (Dwek, 1998). Nevertheless, there is no evidence that this factor, derived from high-resolution models of our Galaxy, should be used at face value in the simulation, especially because the spatial resolution ($\sim 1$ physical kpc) of the simulation implies that dust scattering and absorption occurs at the subgrid scale. One expects therefore that the emergent dust attenuation will depend on the smaller scale distribution of dust and metals, which are not resolved. Given the low resolution of the simulation, we do not implement a two-component dust attenuation models to account separately for dust obscuration in both birth clouds and diffuse interstellar medium, although this has been successfully implemented elsewhere (Trayford et al., 2015).

In addition, prior to converting the metal mass into dust mass, it is important to reproduce the correct gas phase metallicity in the simulation. As discussed in Kaviraj et al. (2017), the relatively low resolution reached in Horizon-AGN implies a delayed enrichment of star-forming clouds, which underestimates the gas phase metallicity compared to observations. To correct for this, a redshift-dependent boosting factor (varying from 4 at $z\sim 0$ to 2.4 at $z\sim 3$) has been computed in order to brings the simulated mass-gas phase metallicity relation in agreement with obsevations from Mannucci et al. (2010) at $z=0$, 0.7, 2.5 and Maiolino et al. (2008) at $z=3.5$. However, the normalization of the mass-gas metallicity relation undergoes large variations depending on the chosen observable used to measure gas metallicity, up to a factor of 5 (see e.g. Andrews & Martini, 2013). In particular, the renormalisation of Horizon-AGN gas phase metallicity tends to align onto the highest values derived from observations (namely using the R23 method, see e.g. Lian et al., 2015). As a consequence, one can expect the boosting factor derived in Kaviraj et al. (2017) to be an upper limit. In fact, while comparing the simulated galaxy counts with COSMOS2015 in various bands as a function of redshift after renormalizing the gas-phase metallicity, we find that a dust-to-metal mass ratio of 0.4 is excessive as it results in too few galaxy counts compared to observations. Therefore, we empirically choose a dust-to-metal mass ratio of 0.2, which yields an overall better agreement with the observed counts in all bands at all redshifts.

A.1.2 Attenuation curve

Beyond the dust mass distribution in the galaxy, the amount of extinction at a given wavelength will be dependent upon the chosen attenuation curve. In this work, the ${\rm R}_{V}=3.1$ Milky Way dust grain model by Weingartner & Draine (2001) is used for post-processing the simulated galaxies. This model includes in particular the prominent 2175 $\mathrm{\SIUnitSymbolAngstrom}$-graphite bump. In the simulation, the spectrum of each stellar particle, assumed to be a SSP, is attenuated with this model and knowing the specific dust column density in front of each particle. However, when summing up the contribution of all the SSPs, we find that the overall attenuation curve of the resulting galaxy spectrum is less steep and the bump tends to reduce. As discussed in Fischera & Dopita (2011), turbulence can further reduce the bump and flatten the extinction curve. Fig. 14 presents the extinction curve in Horizon- AGN as averaged over one thousand galaxies randomly selected, and the two curves used in LePhare for the computation of stellar masses and SFR. $k_{\lambda}$ is defined as $A_{\lambda}/{\rm E}(B-V)$. None of the two curves used in LePhare can correctly reproduce the one in Horizon-AGN, and this discrepancy is likely to be the reason for the bimodality in the SFR, as further discussed in Appendix B.

A.2 IGM absorption

In order to implement the IGM absorption, the Lyman-$\alpha$ forest is implemented on each galaxy line-of-sight from the gas density, velocity and temperature in the IGM in front of the galaxy.

Let us consider the line-of-sight of a background source emitting at the observed wavelength $\lambda_{0}$. Here $\lambda_{\alpha}=1215.7$ $\mathrm{\SIUnitSymbolAngstrom}$ corresponds to the transition from the ground state to the first excited state of the Hydrogen atom. The wavelength of the photons emitted by the background source spectra are redshifted by a factor of $(1+z)$. At some point the light from the source will be redshifted at $\lambda_{\alpha}=1215.7$ $\mathrm{\SIUnitSymbolAngstrom}$. At this point, it may be absorbed by HI from the IGM. The probability of transmission of the light at the observed frequency $\nu_{0}$ will be given by:

[TABLE]

where $\tau_{\alpha}(\nu_{0})$, the Ly-$\alpha$ optical depth at the observed frequency $\nu_{0}$ is given by

[TABLE]

where $x$ is the comoving coordinate of the comoving point varying along the line of sight between the observer ($x=0$) and the source ($x=x_{s}$), $z$ is the corresponding redshift, $n_{\rm HI}$ the neutral hydrogen density at point $x$ and redshift $z$, and $\sigma_{\alpha}$ the Ly-$\alpha$ cross-section. Here $\sigma_{\alpha}$ is a function of the frequency $\nu$ of the photon with respect to the rest frame of the neutral hydrogen at position $x$. Here $\nu=\nu_{0}(1+z)(1+v/c)$, where $v$ is the peculiar velocity along the line-of-sight. So $\sigma_{\alpha}$ may be then written as

[TABLE]

where $b(x,z)=\sqrt{2k_{B}T(x,z)/m_{p}}$, $\sigma_{\alpha,0}=(3\pi\sigma_{\rm T}/8)^{1/2}f\lambda_{\alpha}$, with $\sigma_{\rm T}=6.25\times 10^{-25}$cm2 the Thomson cross-section, and $f=0.4162$ the oscillator strength. As we do not save the neutral hydrogen outputs for Horizon-AGN, the neutral hydrogen density is computed in post-processing by considering that the fraction $x_{\rm HI}=n_{\rm HI}/n_{\rm H}$ is a balance between photo-ionizations, collisional ionizations and recombinations. At equilibrium with the cosmic UV background field, it yields:

[TABLE]

where $\alpha$ and $\gamma$ are the collisional recombination and ionisation rates, $\Gamma$ is the photoionisation rate, and $n_{e}$ is the free electron number density. Considering a uniform background radiation field as implemented in the simulation, the photoionisation rate is assumed to be spatially uniform. Its overall normalization is quite uncertain and is adjusted in order to match the PDF of the transmitted flux at $z=1.5$, $z=2$, $z=2.5$ and $z=3$ (see and e.g. Theuns et al., 1998; Bolton et al., 2005; Lukić et al., 2015; Bolton et al., 2017). The chemical composition of the IGM is close to primordial, so it can safely be assumed that $n_{e}$ receives contribution only from ionised Hydrogen and Helium (assumed entirely in its fully ionized form), which allows to determine its composition (Choudhury et al., 2001). Prescriptions from Black (1981) are used to determine the collisional recombination and ionisation rate as a function of gas temperature.

In order to take into account the full Lyman series absorption, we assume that the only difference between Lyman-$\alpha$ and other Lyman transition comes from different absorption cross-sections. This modelling hence neglect different widths of Lorentz profiles from one transition to an other, but this is expected to be a secondary effect (Iršič & Viel, 2014).

Eventually, galaxy spectra are then multiplied with the Lyman-series absorption lines. Fig. 15 shows the median absorption by the IGM in the $u$ and $B$ bands as a function of redshift and compares with the litterature (Meiksin, 2005; Inoue et al., 2014). The median absorption in a given band will also depends on the hardness of the galaxy UV spectrum, and on the presence or absence of Lyman-limit systems (Meiksin, 2005). On overall, our implementation of IGM absorption matches well the litterature. As noted in e.g. Inoue et al. (2014), the Madau model tends to slightly overestimate the correction at $z>3$ compared to observations, and therefore also overestimate the correction of our virtual galaxies.

A.3 Flux error implementation

Implementing realistic errors on the flux is crucial for the accuracy of our forecasts to retrieve correct redshift and masses. To this end, we compute flux errors for our mock galaxies and perturb their fluxes accordingly.

For the COSMOS-like catalogue, this implementation is done in each band while relying on the COSMOS2015 catalogue. The real dataset is split in small bins of flux141414 In COSMOS2015 the total fluxes are computed from the corrected aperture magnitudes (using the offset defined in Equation 4 of L16), and then accounting for the Galactic foreground extinction (following their Equation 10). . In each bin the flux error distribution is fitted with a Gaussian function. The Horizon-AGN catalogue is then divided in the same way. For each simulated galaxy in a given bin, and for each COSMOS filter, an error is randomly chosen according to the COSMOS2015 error distribution. Galaxy apparent fluxes, which initially corresponded “exactly” to the star particles’ content, are now perturbed according to their $1\sigma$ error. We note that for this implementation, the intrinsic fluxes from Horizon-AGN are confronted to the -already perturbed- observed ones. As a result, the horizontal width of the faint-end tail tends to be larger in the simulation as in observations.

For the Euclid-like, LSST-like and Euclid+LSST catalogues, errors are implemented according to the COSMOS2015 flux error distribution in the closest filters pass-bands, and shifted according to the expected depth of the surveys at similar number of $\sigma$. We emphasize that, although reasonable, these errors might no reflect the specific noise of the survey, and in any case it does not take into account possibly systematics in the photometry.

Fig. 17 shows the distribution of $1\sigma$ errors, as a function of total magnitudes, in several bands in Horizon-AGN (grey area) and compares them to observed data (red area) when available (for the COSMOS-like sample). The adopted depths in all bands are summarised in Table 5.

Fig. 16 presents the simulated galaxy count in the $r$ and $K_{\rm s}$ bands in bins of redshift and compares them to the counts from COSMOS2015. Several features can be noticed. First of all the Horizon-AGN simulated catalogue is mass limited while COSMOS2015 is not, which explains the drop-off of the count at faint magnitudes in the $K_{\rm s}$ band at low-redshift. Part of the discrepancy at bright magnitudes is driven by the smaller area covered by Horizon-AGN (1deg2) compared to COSMOS2015 ($\sim$1.4deg2). Finally, we find that the galaxy counts are overestimated at the faint end of the distribution. This effect has already been discussed in Kaviraj et al. (2017) and is mostly driven by stellar feedback not being strong enough in Horizon-AGN.

A.4 Stellar mass loss

Let us discuss the prescriptions used in both hydrodynamical simulations and stellar population synthesis models to take into account stellar mass losses due to galactic winds, remnants and supernovae. The impact on the comparison between simulated and observed data is not negligible if the two samples rely on different prescriptions.

In Horizon-AGN, the stellar mass loss is modelled as a function of time and metallicity assuming that stars are distributed with a Salpeter (1955) IMF and supernovae type Ia occur with the frequency described in Greggio & Renzini (1983), assuming a binary fraction of 5 percent. To compare to SED-fitting estimates, one may either rescale the $M_{*}$ values in the simulation by a factor $\sim$1/1.7 (the usual conversion from Salpeter to Chabrier IMF, e.g. Santini et al., 2011) or perform the SED fitting with BC03 templates that assume a Salpeter IMF. Neither of these solutions is sufficient to fully remove the bias because even when the IMF is the same, the resulting fraction of ejected stellar mass may significantly differ. In other words, the mass of Horizon-AGN stellar particle also account for remnant mass. Different SSP models will implement their formation differently, and therefore at a given time and metallicity the remnant mass will change from one model to the other. Fig. 18 presents the mass evolution of a SSP as a function of time in Horizon-AGN and using the Maraston (2005) and BC03 models.

To account for this discrepancy, one possible solution is to choose a SED-fitting library based on a stellar population synthesis model whose features are in better agreement with that used in the simulation. However, we cannot modify LePhare as we want to compare our findings to COSMOS2015’s and therefore be consistent with the set-up used for that catalogue. Instead, we prefer to correct the Horizon-AGN virtual photometry by matching a BC03 SSP to each hydrodynamical particle at the time of its formation. We let the SSP evolve so that at any age of the stellar particle we can compute the mass loss fraction according to the BC03 model.

We stress that these details, related a-priori to the “sub-grid physics”, should be reconciled before comparing the simulation to real data (e.g., Henriques et al., 2015). However, although most of the studies take into account the IMF conversion, it is difficult to find in the literature comparisons that correct for the different stellar mass loss parametrisation (see e.g., Davidzon et al., 2018).

A.5 Limitations of our modelling

It should be finally emphasized that our end-to-end modelling still has shortcomings, which are listed below.

No systematics in the photometry

Although statistical photometric errors are consistently implemented in the virtual dataset (see Section A.3), systematics arising when extracting the photometry from the images (blending of objects, clumpy objects possibly split at the extraction, lensing magnification and PSF, image artefacts) are not accounted for here and will be the topic of future work.

Spatially constant IMF and stellar mass loss prescriptions

Our modelling also ignores specific aspects of galaxy evolution which can modify the photometry. For example, nebular emission lines are not implemented in the photometry. Furthermore, the pipeline implicitly assumes that the IMF does not spatially vary within the galaxies, and is perfectly known at the SED-fitting stage. A Chabrier IMF (Chabrier, 2003) is de facto used both when computing the photometry from the simulated catalogue and to build the galaxy template library for SED fitting. Fitting the galaxy photometry with a different IMF from the one chosen to compute this photometry would obviously lead to new systematics in the stellar mass computation. In addition to the choice of the IMF, which controls the amount of stars formed as a function of their mass (and therefore in particular the overall mass-to-light ratio), the chosen prescription for stellar mass losses as a function of time and metallicity is important. In the work presented here, the simulated photometry is computed with BC03 SSP templates (see Section A.4), and fitted with a SED library which includes a higher diversity of stellar mass loss prescriptions. Therefore, we do not assume to know a priori the stellar mass loss prescriptions of the simulated galaxies, and effectively the best-fit is not always a BC03 template. However, in practice the simulated galaxies have all the same stellar mass loss prescriptions, and these prescriptions are spatially constant within the galaxies: this is unlikely to be the case in the real Universe. Therefore the simulated galaxy population present less diversity than the observed one, and one could expect the SED-fitting to perform much better on the Horizon-AGN Universe than on the real one.

In despite of these limitations, we found that the $z_{\rm phot}$ accuracy of the virtual catalogue is comparable to that of COSMOS2015. This suggests that varying IMF and stellar mass loss prescriptions should not dramatically impact redshift reconstruction, which indeed relies on galaxy colours (i.e. the relative values of flux in different bands, and not the absolute value of the flux). Only the redshift errors are lower in the simulated catalogue with respect to the real one. However a much larger and systematic impact on stellar mass and SFR is expected, and the reconstruction accuracy quoted in Section 3.2 should be understood as an optimistic case. On the other hand, we have shown that with a very good photometry stellar masses are well retrieved (within 0.12 dex) with the COSMOS-like configuration when the IMF is fixed. Therefore, one can expect that all additional systematics will be driven by the uncertainty on the IMF. Comparing the mass from the SED-fitting with an independent measurement (kinematics, small-scale lensing) can therefore be a way to constrain the IMF.

Appendix B Impact of absorption on SED

Dust attenuation is implemented in the virtual photometry as post-processing at the scale of $\sim$1 kpc, following the prescription detailed in Appendix A. In particular, at equal dust mass, the resulting attenuation of the total spectrum will depend on the geometry of the galaxy and the angle under which it is seen. Conversely, at the SED-fitting stage, dust attenuation is kept as simple as possible and therefore does not depend on galaxy geometry: all the SSPs are assumed to undergo the same attenuation given a specific attenuation curve. In a similar way, IGM absorption is implemented in the Horizon-AGN lightcone independently for each galaxy, knowning the foreground distribution of HI. However, at the SED-fitting stage, IGM absorption is accounted for with an average redshift-dependent relation for all galaxies. These differences can play a role in driving the scatter and the systematic trends observed in mass and star-formation rate comparison. Further details on this effect are presented now.

B.1 $z_{\rm phot}$ performance

LePhare is run with 3 different input photometric catalogues, in order to show the impact of dust and IGM on $z_{\mathrm{phot}}$. The main run (presented in Section 3.1) includes both inter-galactic and intra-galactic (dust) absorption. An additional run is performed on an input photometry that does not include absorption from either components, and the last one is performed on an input photometry that includes only inter-galactic absorption. The comparison of the performance of these runs directly tells us about the impact of inhomogenous IGM and spatially varying dust absorption on the SED-fitting performance.

Results of the $z_{\rm phot}$ performance for no absorption, IGM-absorption only and IGM+dust absorption are presented in Table 6. On overall IGM slightly helps to constrain the $z_{\rm phot}$. However, as discussed in Section 3.1, we note that overestimating the IGM absorption at the SED-fitting stage can make a large fraction of the catastrophic outliers falling in $\left[z_{\rm sim}>2.5\right]\cap\left[z_{\rm obs}<1.5\right]$. Adding dust to the photometry reduces slightly the performance. Table 6 has to be interpreted with caution however. In a given magnitude bin, the dust-free catalogue probes a galaxy population which can overlap a fainter magnitude bin in the dusty catalog. In other words, not only the presence/absence of dust drives the difference in the performance of the $z_{\rm phot}$ in a given magnitude bin, but also the possible differences in the intrinsic SED of the galaxies.

B.2 Stellar mass estimation

As the impact of IGM is limited at the very high redshift, the subsequent analysis focus on the impact of absorption as a whole on stellar mass and SFR estimates. In order to isolate the effect of absorption, the physical parameter estimation is done at the intrinsic redshift of the galaxies (i.e. the mass and SFR uncertainties do not include the $z_{\rm phot}$ uncertainties) on two input photometric catalogues (with and without absorption).

Fig. 19 presents the comparison between the intrinsic and observed stellar mass, for galaxies selected at $K_{\rm s}<24.7$. Horizontal orange arrows indicate the mass completeness in each redshift bin. Without absorption, mass completeness is naturally better as galaxies are brighter. Dust attenuation induces a negligible overestimation of the observed stellar mass especially at high redshift and for massive galaxies. This is expected, as the mass estimates is essentially provided by the NIR bands which are barely affected by the dust. However, when only optical bands are available (e.g. in the LSST-only configuration), dust is expected to have a more dramatic effect, similar to the one observed in the SFR computation (see the discussion below).

Without any kind of absorption (bottom panel), there is a persistent underestimation of stellar mass by at most 0.1 dex at low redshift, owing to the simplistic SFHs and metallicity distributions in the template library. As noted before (see Section 3.2) the effect of dust (overestimation) and SFHs (underestimation) can act in opposite direction and therefore tend to compensate each other.

B.3 SFR estimation

Fig. 20 presents the comparison between the intrinsic and observed SFR, for galaxies selected at $K_{\rm s}<24.7$ and $0.2<z<4$. The effect of dust is dramatic: it drives a large scatter and a bimodality, with a systematic overestimation of the SFR for a fraction of the population up to at least $z\sim 2.5$. This is a consequence of the degeneracy between dust and SFR. Fig. 21 isolates what in the dust modelling drives the effect. The comparison between SFR estimated from the best-fitted template and the intrinsic SFR is presented at $1.1<z<1.5$. The diagram is colour-coded by the extinction curve used in the fit (left panel) and $\Delta_{A_{NUV}}$ (right panel), where $\Delta_{{\rm A}_{NUV}}={\rm A}_{NUV}^{\rm sim}-{\rm A}_{NUV}^{\rm phot}$. Qualitatively, when the best-fit template is derived using Arnouts et al. (2013)’s extinction curve, $A_{NUV}$ is overestimated and the SFR is consequently overestimated. The reverse happens when the best-fit template is attenuated using Calzetti et al. (2000)’s extinction curve.

Appendix C Zero-point magnitude offsets

It is common in the literature to apply an offset to apparent magnitudes in order to correct for systematics due to calibration discrepancies between different filters (e.g., Ilbert et al., 2013). These offsets are computed by LePhare using a spectroscopic sub-sample: the code fits galaxy SEDs after fixing their redshifts at the spectroscopic value, then it compares the magnitude observed in each filter to the one of the best-fit model (i.e., the template minimizing the reduced $\chi^{2}$). An offset is added in each filter to reduce the difference between predicted and observed magnitudes. The code iterates the procedure until convergence, finding the final values of the zero-point offsets that will be add by default in the next LePhare runs.

Although such a procedure is quite efficient at improving SED-fitting results, they might also introduce a bias if the used spectroscopic sub-sample is not representative of the whole population. Moreover zero-point offsets correct, at least partly, for possible incompleteness in the template library. For this reason some authors prefer not apply them, e.g. when deriving stellar masses, because besides solving calibration issues they may affect the physical interpretation of the SED fitting results (see Moutard et al., 2016).

We compute the zero-point offsets in the COSMOS filters, not for a spectroscopic-subsample but for the whole galaxy catalogue. All the offsets are found to be smaller than $0.02$ mag, which is of the order of the minimal photometric errors, with the exception of 0.029 and 0.039 mag in the $z^{++}$ and $K_{\rm s}$ band respectively. In COSMOS2015 these offsets are generally much larger (see e.g. Table 4 in Laigle et al., 2016), an indication that in real datasets they are mainly coming from calibration issues. Because such issues are not present in the virtual photometry, there is no need to include zero-point offsets to our virtual magnitudes.

Appendix D Estimating the redshift errors

D.1 Robustness of $\sigma_{z}^{\rm fit}$

By measuring the cumulative distribution of $|z_{p}-z_{s}|/(1+\sigma_{z}^{\rm fit})$ in bins of magnitudes (figure 13 in L16) for the high-confidence spectroscopic sub-sample in COSMOS, L16 concluded that $\sigma_{z}^{\rm fit}$ derived by LePhare underestimates $\sigma_{z}^{\rm true}$, with a trend increasing with fainter magnitudes. In an effort to better quantify this trend, Figure 22 displays $\log\sigma_{z}^{\rm true}/\sigma_{z}^{\rm fit}$ in bins of magnitude and redshift for the COSMOS spectroscopic sample (extreme left panel) and the Horizon-AGN (middle left panel) simulated sample. This plot highlights that Horizon-AGN presents also an underestimation of $\sigma_{z}^{\rm fit}$ at bright magnitudes in the redshift range $1<z<2.5$. This underestimation might be due either to a remaining underestimation of magnitude errors, or to a lack of representativeness of the template library. The distribution of the reduced $\chi^{2}$ in this plane (middle right panel) speaks in favour of this claim, as the regions with higher $\chi^{2}$ broadly match those where $\sigma_{z}^{\rm fit}$ is underestimated the most. The underestimation is more severe in COSMOS, because the real photometry presents more diversity than the simulated one.

D.2 Catastrophic outliers

Let us now investigate what causes the higher fraction of catastrophic outliers in the observed zCOSMOS sample (Lilly et al., 2007) with respect to the virtual photometric catalogue, as displayed in Fig. 2. It is important to check if this discrepancy is driven by additional observational limitations (spectroscopic redshift misidentification151515This is a potential issue specific to the zCOSMOS-Deep sample, zCOSMOS-Bright being much more secure., crowed photometry, etc.) independent on photometry modeling, or if the virtual photometric catalogue actually misses some essential components of the real galaxy population, which would make it a poor predicator of the accuracy of SED-fitting performance.

Let us therefore focus on the failure for the objects marked as red squares on Fig. 2. After individual inspection of the spectroscopy and the photometry for each of these objects, we conclude that in the large majority of the cases, the failure arises because of one of the two reasons:

Uncertain photometry

This case happens for $\sim 35$ per cent of the outliers. The photometry extraction might be uncertain because of clumpy galaxies which might be over-split or even identified as two different objects, or on the contrary because of blended objects (which is an issue particularly severe for $NUV$ or IRAC bands given the confusion limit). In the latter case, removing the IRAC bands in the SED fitting improves the match with the spectroscopic redshifts. This process does not impact the virtual catalogue because the identification of the virtual galaxies is done directly in 3D.

Spectroscopic misidentification

For another $\sim$35 per cent of the outliers, a second object on the VIMOS slit of a length of $>10\arcsec$ is misidentified as the target in the zCOSMOS-Deep observations. Therefore, the spectroscopic redshift attributed to the original target is erroneous. In Figure 23, an example of the misidentification issue is shown. We found that these $z_{\rm spec}$ often agree well with the photometric objects of the second on-slit objects. The full description on the zCOSMOS-Deep sample and the redshift evaluation will be given in a future paper (Lilly et al., in prep).

Finally note that for a few outliers in COSMOS at low redshift, removing the $NUV$ photometry in the SED-fitting improves the match with the spectroscopic redshifts. Recall that the $NUV$ band is not used in the SED fitting of the Horizon-AGN COSMOS-like catalogue, which could also be a reason for the lower fraction of outliers.