11institutetext: Universitäts-Sternwarte, Fakultät für Physik, Ludwig-
Maximilians-Universität München, Scheinerstr. 1, 81679 München, Germany
11email<EMAIL_ADDRESS>22institutetext: Kavli Institute for
Particle Astrophysics and Cosmology, Department of Physics, Stanford
University, Stanford, CA, USA 33institutetext: SLAC National Accelerator
Laboratory, Menlo Park, CA, USA 44institutetext: Excellence Cluster ORIGINS,
Boltzmannstr. 2, 85748 Garching, Germany

# The impact of stellar population synthesis choices on forward-modelling-
based redshift distribution estimates

Luca Tortorelli 11 Jamie McCullough 112233 Daniel Gruen 1144

(Received XXXX; accepted XXXX)

###### Abstract

Context. The forward-modelling of galaxy surveys has recently gathered
interest as one of the primary methods to achieve the required precision on
the estimate of the redshift distributions that will allow stage IV surveys to
perform cosmological tests with unprecedented accuracy. One of the key aspects
of forward-modelling a galaxy survey is the connection between the physical
properties drawn from a galaxy population model and the intrinsic galaxy
spectral energy distributions (SEDs), achieved through stellar population
synthesis (SPS) codes, e.g. FSPS (Conroy et al., 2009). However, SPS requires
many detailed assumptions about the constituents of galaxies, for which the
model choice or parameter values are currently uncertain.

Aims. In this work, we perform a sensitivity study of the impact that the
variations of the SED modelling choices have on the mean and scatter of the
tomographic galaxy redshift distributions.

Methods. We assume the Prospector-$\beta$ (Wang et al., 2023b) model as the
fiducial input galaxy population model and use its SPS parameters to build
observed magnitudes of a fiducial sample of galaxies. We then build samples of
galaxy magnitudes by varying one SED modelling choice at a time. We model the
colour-redshift relation of these galaxy samples using the KiDS-VIKING
remapped version (McCullough et al., 2023) of the Masters et al. (2015) self-
organising map (SOM). We place galaxies in the SOM cells according to the
simulated galaxy colours. We then build color-selected tomographic bins
following the prescription of McCullough et al. (2023) and compare each
variant’s binned redshift distributions against the estimates obtained for the
original Prospector-$\beta$.

Results. We find that the SED components related to the initial mass function,
to the Active Galactic Nuclei, to the gas physics, and to the attenuation law
substantially bias the mean and the scatter of the tomographic redshift
distributions with respect to those estimated with the fiducial model.

Conclusions. For the uncertainty of these choices currently present in the
literature, and regardless of any stellar mass function based reweighting
strategy applied, the bias in the mean and the scatter of the tomographic
redshift distributions is larger than the precision requirements set by next
generation Stage IV galaxy surveys, e.g. the Vera C. Rubin Observatory’s
Legacy Survey of Space and Time (LSST) and Euclid.

###### Key Words.:

(Cosmology:) large-scale structure of Universe — Galaxies: statistics —
Galaxies: stellar content

##  1 Introduction

In recent years, results from galaxy surveys have started to enter the realm
of precision cosmology, a prerogative that had previously been unique to the
studies based on cosmic microwave background data. Next-generation
cosmological surveys, like the Vera C. Rubin Observatory’s Legacy Survey of
Space and Time (hereafter, Rubin-LSST, Ivezić et al. 2019), Euclid (Racca et
al., 2016), and the Roman space telescope (Dore et al., 2019), aim at
providing the best tests yet of models for the dark energy, including from
constraints on $\sigma_{8}$ and $\Omega_{M}$ that either reconcile or further
strengthen the existing tension between the cosmological parameters estimated
from early Universe probes and those from the late Universe (e.g. Planck
Collaboration et al. 2020; Heymans et al. 2021; Abbott et al. 2022).

In order to perform the cosmological measurements with the large-scale
structure distribution, we need to measure the statistics of the density field
evolving over time. Two of the main probes used for this analysis are the weak
gravitational lensing shear field and the galaxy clustering. Weak lensing (see
Mandelbaum 2018 for a comprehensive review) is widely considered as one of the
most promising probes of the growth of dark matter structure and the dark
energy equation of state because light from galaxies gets continuously lensed
by the intervening matter and this gets imprinted on the measured shape of
galaxies. By measuring the latter for millions of objects, we can reconstruct
the intervening mass distribution and the statistics of the density field.
Galaxy clustering relies on using galaxies as a biased tracer (Desjacques et
al., 2018) of the underlying matter field at all physical scales. Galaxy
clustering provides information about the large scale structure distribution
of matter in the Universe, which in turns allows us to obtain information
about the structure growth and the expansion history of the Universe. For both
probes, the cosmological information is usually extracted using the angular
two-point correlation function that gives a glimpse of the spatial
distribution of galaxies in tomographic photometric redshift bins (e.g.
Rodríguez-Monroy et al. 2022).

Both weak lensing and galaxy clustering measurements are a challenging effort
because they require a strong control of the observational systematic effects
(Massey et al., 2013; Shapiro et al., 2013; Mandelbaum, 2015). This is
especially true for the systematics related to the accurate measurement of
galaxy shapes (Kacprzak et al., 2012, 2014), to the spatially varying survey
properties (Rodríguez-Monroy et al., 2022) and to those related to the robust
estimates of galaxy redshifts (Crocce et al., 2016; Myles et al., 2021).
Indeed, accurate knowledge of the redshift bins of the galaxy samples observed
by cosmological galaxy surveys is required to obtain precise cosmological
parameter estimates.

This is the reason why stage IV experiments (Albrecht et al., 2006) have
stringent requirements on the tomographic redshift distribution uncertainties.
The latter is generally modelled in terms of uncertainty on the mean redshift
of each tomographic bin, since cosmological parameters are most sensitive to
shifts in the redshift mean $\bar{z}$ (Amon et al., 2022; van den Busch et
al., 2022; Li et al., 2023; Dalal et al., 2023), and uncertainty in redshift
scatter $\sigma_{z}$ within a bin, which has a second-order effect on the
amplitude of lensing signals, but a first-order effect on the strength of
angular galaxy clustering. In the case of Rubin-LSST, the requirements set
down by the LSST Dark Energy Science Collaboration (The LSST Dark Energy
Science Collaboration et al., 2018) foresee that, for the weak lensing
measurements, the systematic uncertainty in the mean redshift of each source
tomographic bin should not exceed $0.002(1+z)$ in the year 1 (Y1) analysis,
which becomes $0.001(1+z)$ in the year 10 (Y10) analysis. For galaxy
clustering, the systematic uncertainty on the mean should not exceed
$0.005(1+z)$ in the Y1 and $0.003(1+z)$ in the Y10 analysis. In addition, the
systematic uncertainty in the source redshift scatter $\sigma_{z}$ should not
exceed $0.006(1+z)$ in the Y1 and $0.003(1+z)$ in the Y10 weak lensing
analysis, while for galaxy clustering these numbers lower to $0.1(1+z)$ for Y1
and $0.03(1+z)$ for Y10.

Rubin-LSST is going to measure broad-band fluxes for billions of objects up to
$i$-band depth of 26.3 for the Y10 catalogue (Ivezić et al., 2019). Obtaining
a redshift for each object via spectroscopy is therefore not a viable solution
and we need to resort to photometric redshift (photo-z) estimates.
Unfortunately, the photo-z calibration accuracies of existing methods are
still not able to achieve the Rubin-LSST requirements for the weak lensing and
galaxy clustering analysis. Although the most recent methods (Wright et al.,
2020; Myles et al., 2021) when applied to Stage III dark energy experiments
under idealized conditions yield uncertainties on the mean redshift that are
of the order of $0.001-000.6$, there are a number of effects that
significantly degrade the redshift characterisation. Indeed, effects such as
spectroscocpic incompleteness, outliers in the calibration redshift, sample
variance, blending and astrophysical systematics may cause systematic errors
in the mean redshift and scatter that are $\sim 7-10\times$ larger than the
Stage IV requirements (Newman & Gruen, 2022).

Two main factors concur with the lack of the required accuracy on the redshift
distribution estimates (see Newman & Gruen 2022 for a comprehensive review).
The first is the inability to estimate the redshift of an individual galaxy
precisely since estimates are generally based on measurements in only a few
broad noisy photometric bands where very few spectral features can be used to
constrain the galaxy redshift, leading to degeneracies between multiple fits
to a galaxy spectral energy distribution (SED, Buchs et al. 2019; Wright et
al. 2019; Wang et al. 2023b). The second factor is having a good knowledge of
the galaxy population, which, given the intrinsic uncertainty described above,
sets an informative prior on the redshift of a photometrically observed
galaxy. This is a truly fundamental problem because often the spectroscopic
redshift samples on which photo-z methods are trained, or physical models of
the galaxy population are built, systematically miss populations of galaxies,
giving rise to biased colour-redshift relations.

The forward-modelling of galaxy surveys offers an alternative approach to the
problem of accurate galaxy redshift distribution estimates. The forward
process involves the detailed modelling of all the observational and
instrumental effects of a galaxy survey. This includes the modelling of the
galaxy population, i.e. the intrinsic distribution of redshift-evolving
physical properties of galaxies, the modelling of the galaxy stellar
populations, i.e. the connection between the physical properties and the
galaxy SEDs, the mapping of intrinsic properties and fluxes to the observed
ones by means of image (see Plazas Malagón 2020 for a review) and spectra
simulators (Fagioli et al., 2018, 2020) and the characterisation of the
selection function of galaxies observed in a survey, conditioned on these
physical properties and observational and instrumental effects. This method
has already shown promising results in characterising survey redshift
distributions and galaxy population properties when applied to photometric and
spectroscopic galaxy surveys (Fagioli et al., 2018; Tortorelli et al., 2018;
Fagioli et al., 2020; Tortorelli et al., 2021; Kacprzak et al., 2020; Herbel
et al., 2017; Bruderer et al., 2016; Alsing et al., 2023; Leistedt et al.,
2023; Alsing et al., 2024; Moser et al., 2024).

The ongoing forward-modelling efforts model the galaxy population physical
properties, such as redshifts, stellar masses, star-formation histories and
metallicities, via either parametric relations (e.g. Tortorelli et al. 2020;
Wang et al. 2023b; Alsing et al. 2023; Moser et al. 2024) or machine learning
models (Alsing et al., 2024). The galaxy stellar population models, which in
this work are referred to the set of prescriptions that are used to generate a
galaxy spectrum based on its physical properties (e.g. the stellar templates,
the prescription for the dust attenuation, active galactic nuclei emission,
initial mass function, velocity dispersion), are instead created using three
main methods: empirically determined templates (e.g. Blanton & Roweis 2007;
Brown et al. 2014 templates), stellar population synthesis (SPS) codes (e.g.
Flexible Stellar Population Synthesis FSPS Conroy et al. 2009; Conroy & Gunn
2010) and, equivalently but computationally much more feasible, emulators of
SPS-generated spectra, such as SPECULATOR (Alsing et al., 2020), SPENDER
(Melchior et al., 2023) and DSPS (Hearin et al., 2023). The empirical
templates are generally controlled by a limited number of parameters, e.g. the
coefficients of the linear combination of templates (Tortorelli et al., 2021),
while the population of SPS-based SEDs, especially those that employ FSPS, are
constructed either by drawing from the stellar population model parameter
distributions employed in the Prospector-$\alpha$ (Leja et al., 2017, 2018,
2019b) and Prospector-$\beta$ models (Wang et al., 2023b, a), or by using the
default values implemented in a code like FSPS. However, it has been shown in
Conroy et al. (2009) that the choices of the parameters of the stellar
population components affect the galaxy colours up to the magnitude level.
Since the galaxy colours are used to assign a galaxy to a redshift bin for
weak lensing and galaxy clustering studies, and to estimate the redshift
distribution of a bin via the colour-redshift relation, studying the impact of
the choices of the stellar population parameters and components on the SPS-
based forward-modelling estimates of the galaxy redshift distributions is of
crucial importance.

In this work, we evaluate the impact of stellar population parameter choices
on the mean and scatter of the redshift distributions of tomographic redshift
bins to evaluate whether the variations induced by those choices are within
stage IV cosmological survey requirements. The primary purpose of this is to
identify choices that are significantly affecting redshift distribution
estimates and, hence, need to be better informed by data in order not to limit
cosmological studies from these surveys.

In order to do this, we construct galaxy observed-frame SEDs using the FSPS
code and generate model apparent AB magnitudes integrating the SEDs in the
$ugriZYJHK_{s}$ bands. Those bands are the ones used to train the Masters et
al. (2015, 2017) self-organising map (SOM, Kohonen 2001) that robustly maps
the empirical distribution of galaxies in the COSMOS field to the
multidimensional colour space spanned by the filters mentioned above. We use
the remapped version of the SOM presented in McCullough et al. (2023),
therefore the $ugriZYJHK_{s}$ bands refer to those of the KiDS-Viking-450
survey (Wright et al., 2019).

SOMs have been used to calibrate the colour-redshift relation with
spectroscopic quality redshift estimates (Masters et al., 2015, 2017, 2019;
Stanford et al., 2021; McCullough et al., 2023) and for the redshift
calibration of the weak lensing source galaxies (Hildebrandt et al., 2020;
Myles et al., 2021). In this work, we use the model apparent AB magnitudes to
place galaxies in the remapped SOM and use the redshift bin definition in
McCullough et al. (2023) to construct the redshift distributions. We then
compute the mean and scatter of redshift in each of those tomographic bins for
different stellar population parameter choices to evaluate whether their
impact is within stage IV cosmological surveys requirements.

FSPS requires physical properties of galaxies as input, e.g. star-formation
histories (SFHs), metallicities and redshifts. Therefore, we need a model that
serves both as a galaxy population model, from which these physical properties
should be drawn, and as a stellar population model, from which the stellar
population parameters used to build the SED can be drawn instead. We use the
Prospector-$\beta$ model (Wang et al., 2023b) as our input galaxy and stellar
population model, which has been proven to provide a realistic representation
of the galaxy population as shown by its use to robustly measure the physical
properties of galaxies observed with both HST and JWST (Wang et al., 2024). We
use the colours and the mean and scatter of the tomographic redshift
distributions obtained from this model as our fiducial set of properties. We
then build new samples of galaxies varying one stellar population component at
a time. Additionally, we also consider the case where we vary all of them at
the same time for each galaxy within the observationally meaningful ranges set
for every component. We vary all the relevant stellar population components
implemented in FSPS, from those connected to stellar physics to those related
to the Active Galactic Nuclei (AGN) model and the gas emission. We then
evaluate the differences in colours, means and scatters of the tomographic
bins against the fiducial ones. We also apply reweighting strategies based on
fitting the stellar mass function to test whether any SED mismodelling might
be compensated by matching the observed abundance of galaxies in a survey.

In section 2 we describe the Prospector-$\beta$ galaxy population model from
which we sample the physical properties of galaxies and the stellar population
parameters that are used by FSPS to generate the galaxy SEDs. Section 3
describes the stellar population components and their parameters that we vary
in this study. The impact of each component on the rest-frame and the
observed-frame galaxy colours is presented in section 4, while in section 5 we
evaluate the impact that the same components have on the SOM-based colour-
redshift relation, and in particular on the SOM cell assignments. Section 6
describes the SED modelling impact on the tomographic redshift distribution
means and scatters, while section 7 contains the discussion on the
implications of the bias induced by the stellar population components. Section
8 summarises the main findings of this study and provides future directions on
SPS-based forward-modelling studies.

Figure 1: Distribution of physical properties of galaxies drawn from the
Prospector-$\beta$ model. The redshifts and stellar masses are jointly drawn
from the sum of two Schechter functions (Leja et al., 2020), the star-
formation rates (averaged over the last $100$ Myr) are computed from the non-
parametric star-formation histories, while the stellar metallicities are drawn
from a clipped Gaussian approximating the stellar mass-stellar metallicity
relationship of Gallazzi et al. (2005).

##  2 The Galaxy Population Model

Generating galaxy SEDs with SPS codes requires as input a number of physical
properties of galaxies, most notably the galaxy star-formation and metallicity
histories (MZHs). These histories are used by e.g. FSPS, together with a set
of stellar population prescriptions, to build-up the composite galaxy input
spectra.

To generate a realistic representation of the galaxy population properties, we
use the Prospector-$\beta$ galaxy population model (Wang et al., 2023b, a).
When used as a prior, this model has been proven to be successful in
recovering stellar ages, star-formation rates (SFR) and rest-frame colours in
the redshift range $0.2\lesssim z\lesssim 15$ of galaxies from the UNCOVER
survey (Bezanson et al., 2022) that uses Hubble Space Telescope (HST) and
James Webb Telescope (JWST) data in the Abell 2744 cluster field (Wang et al.,
2024). The Prospector-$\beta$ model is tightly linked to FSPS to generate
galaxy SEDs as both galaxy and stellar population model, with most of the
stellar population components and parameter choices related to those
implemented in FSPS itself.

The Prospector-$\beta$ model is a refinement of the Prospector-$\alpha$ model,
presented in Leja et al. (2017, 2018, 2019b) and used to measure the stellar
mass function (Leja et al., 2020) and the star-forming sequence (Leja et al.,
2022) of $0.2<z<3.0$ galaxies in the 3D-HST (Skelton et al., 2014) and
COSMOS-2015 survey (Laigle et al., 2016). The Prospector-$\alpha$ model has
also been used in Alsing et al. (2020) in conjunction with FSPS to generate
the rest-frame SEDs on which the SPECULATOR emulator has been trained.

The Prospector-$\beta$ model introduces two main modifications that makes it a
more realistic representation of the galaxy population in the Universe with
respect to the Prospector-$\alpha$ model. The first improvement is the
addition of an observationally informative prior on the galaxy stellar mass
function. Prospector-$\alpha$ employs a uniform prior from which one draws
galaxy stellar masses, in contrast to the actual, Schechter-like behaviour of
the distribution of galaxy stellar masses (Marchesini et al., 2009; Baldry et
al., 2012; Muzzin et al., 2013; Moustakas et al., 2013; Tomczak et al., 2014;
Grazian et al., 2015; Song et al., 2016; Weigel et al., 2016; Davidzon et al.,
2017; Wright et al., 2018; Weaver et al., 2023). In the Prospector-$\beta$
model, the evolution of the mass function $\Phi(\mathcal{M},z)$ is modelled as
the sum of two Schechter functions (Leja et al., 2020), where the logarithmic
form of a single Schechter function $\Phi(\mathcal{M})$ at fixed redshift $z$
is

$\Phi(\mathcal{M})=\ln{10}\phi_{*}10^{(\mathcal{M}-\mathcal{M^{*}})(\alpha+1)}\exp{(-10^{\mathcal{M}-\mathcal{M^{*}}})}\
,$ (1)

with $\mathcal{M}=\log{M_{*}}$ and $\phi_{*}$, $\mathcal{M^{*}}$ and $\alpha$
evolving with redshift. In our work, we set the minimum and maximum stellar
masses for the mass function to $M_{*,\rm min}=10^{8}M_{\odot}$ and $M_{*,\rm
max}=10^{12.5}M_{\odot}$. The stellar mass function in Leja et al. (2020) is
continuous in redshift and defined only between $0.2\leq z\leq 3.0$, while for
$z<0.2$ and $z>3.0$, the model adopts a stellar mass function fixed at $z=0.2$
and $z=3.0$ parameter values, respectively. The public release of the
Prospector-$\beta$ model allows for the use of higher redshift mass functions,
e.g. Tacchella et al. (2018). However, in our work we generate samples of
galaxies in the redshift range $0.2\leq z\leq 3.0$ that is similarly used by
other galaxy population models, e.g. CosmoDC2 (Korytov et al., 2019).

The stellar masses drawn from $\Phi(\mathcal{M},z)$ are then used to put
constraints on the SFH via a joint prior. This is the second main modification
that differentiates Prospector-$\beta$ from its predecessor. In
Prospector-$\alpha$, the SFH is non-parametric and is described by the mass
formed in seven logarithmically spaced time bins with a continuity prior that
ensures a smooth transition of the SFRs from one bin to another (Leja et al.,
2019a). The ratio between SFRs in adjacent bins
$\log{(\mathrm{SFR}_{n})/\mathrm{SFR}_{n+1}}$ is drawn from a Student-t
distribution. Prospector-$\beta$ modifies this null expectation of a flat SFH
assumption in two ways, since it is inconsistent with a global peak in SFR at
cosmic noon and a systematic trend with mass. On one side, Prospector-$\beta$
matches the expectation value in each bin in lookback time to the cosmic star-
formation rate densities (SFRD) in Behroozi et al. (2019). On the other side,
it allows for a mass dependence on the start of the age bins (Wang et al.,
2023b)

$\log{(t_{\mathrm{start}}/\mathrm{Gyr})}=\log{(t_{\mathrm{univ}}(z)/\mathrm{Gyr})}+\delta_{m}$
(2)

where $\delta_{m}$ value dependes on the galaxy stellar mass. This new prior
tries to match the expectation that high-mass (low-mass) galaxies form earlier
(later). However, it does not account for the quiescent/star-forming
bimodality, approximating the double-peaked distribution of SFRs at a given
mass and redshift as a wide single-peaked distribution. The justification for
this approximation provided in Wang et al. (2023b) is that the quiescent
fraction, computed using the specific star-formation rate (sSFR), roughly
matches with the observed trend at $z<3$ (Leja et al., 2022), which also
provides a further reason for generating mock galaxies in this redshift range.

The galaxy stellar metallicity prescription is the same as the one implemented
in the Prospector-$\alpha$ model. The prior $P(Z_{*}|M_{*})$ allows to sample
galaxy stellar metallicities $\log(Z_{*}/Z_{\odot})$ conditioned on their
stellar masses $\log(M_{*}/M_{\odot})$. Stellar metallicity affects the
optical-to-near-infrared flux ratios and helps to set the normalisation and
shape of the SED blueward of the Lyman limit, therefore properly modelling it
is of great importance. The stellar metallicity is modelled as a clipped
Gaussian distribution where the mean and the standard deviation at fixed
stellar mass are taken from the $z=0$ SDSS stellar mass-stellar metallicity
relationship of Gallazzi et al. (2005). The clipped range of possible galaxy
stellar metallicities is set to match that of the MIST isochrones (Dotter,
2016; Choi et al., 2016; Paxton et al., 2011, 2013, 2015) used in FSPS, namely
$-1.98<\log{(Z_{*}/Z_{\odot})}<0.19$. The standard deviation as a function of
stellar mass is defined as the $84\mathrm{th}-16\mathrm{th}$ percentile range
from the Gallazzi et al. (2005) relationship. This definition of the standard
deviation is roughly twice the observed one from the $z=0$ relationship, and
is adopted in order to account for potential unknown systematics or redshift
evolution. The distribution of physical properties of galaxies drawn from the
Prospector-$\beta$ model is shown in Figure 1.

The Prospector-$\beta$ model also provides the prescriptions for the stellar
population components and their parameters that are required to generate SEDs
with FSPS. For some prescriptions, Prospector-$\beta$ uses the default
implementation in FSPS, while for others, namely the emission from the gas
component, the dust attenuation and emission, and the AGN component, it uses
specific prescriptions, some of them already implemented in the
Prospector-$\alpha$ model.

The galaxy gas component is modelled as a two-component emission, the nebular
continuum and the nebular emission lines, self-consistently using the stellar
ionizing continuum to power nebular lines and continuum emission, following
the CLOUDY implementation within FSPS (Byler et al., 2017). The two parameters
that govern the gas emission are the gas-phase metallicity
$\log(Z_{\mathrm{gas}}/Z_{\odot})$ and the gas ionization parameter $\log U$,
which is the ratio of ionizing photons to the total hydrogen density. In the
Prospector-$\alpha$ model, the gas-phase metallicity is set equal to the
stellar metallicity at the time of observation, under the assumption that the
stars have the same metallicity as the gas cloud from which they had formed.
In Prospector-$\beta$, instead, the gas-phase metallicity is decoupled from
the stellar one (Shapley et al., 2015; Steidel et al., 2016), and it is drawn
from a flat prior in the range $-2.0<\log(Z_{\mathrm{gas}}/Z_{\odot})<0.5$,
following the observations that the gas in the star-forming galaxies at higher
redshifts has metallicity abundances that may differ significantly from their
stellar abundances (Shapley et al., 2015; Steidel et al., 2016). The
ionization parameter is instead kept fixed at $\log U=-2$, although, as
mentioned in Leja et al. (2019b), gas in star-forming galaxies at higher
redshift might experience a stronger ionizing radiation field that requires
raising the ionization parameter to $\log U=-1$ (Cohn et al., 2018).

The dust attenuation is modelled with two components (Charlot & Fall, 2000), a
birth-cloud and a diffuse dust screen. The birth cloud component adds extra
attenuation towards young stars, simulating their embedding in molecular
clouds and HII regions. It affects nebular emission as well. The birth cloud
attenuation component scales as

$\hat{\tau}_{\lambda,1}=\hat{\tau}_{1}(\lambda/5500\ \AA)^{-1.0}$ (3)

and attenuates only stars formed in the last $10\ \rm Myr$, the typical
timescale for the disruption of a molecular cloud (Blitz & Shu, 1980). The
diffuse component, instead, attenuates all the light from the stars and the
nebular emission of a galaxy. Its wavelength dependence is set by the Noll et
al. (2009) prescription,

$\hat{\tau}_{\lambda,2}=\frac{\hat{\tau}_{2}}{4.05}[k^{{}^{\prime}}(\lambda)+D(\lambda)]\left(\frac{\lambda}{\lambda_{V}}\right)^{n}\
,$ (4)

which combines the fixed Calzetti et al. (2000) dust attenuation curve
$k^{{}^{\prime}}(\lambda)$ with a Lorentzian-like Drude profile $D(\lambda)$
describing the UV dust bump, whose strength is linked to the diffuse dust
attenuation index $n$ (Noll et al., 2009; Kriek & Conroy, 2013). The
prescription in Noll et al. (2009) is meant to provide a robust fit to the
individual galaxy SEDs and is motivated by the observations that show that a
single attenuation curve, such as the Calzetti et al. (2000) curve, is
insufficient to describe the diversity of attenuations found in star-forming
galaxies, particularly at high redshift. Therefore, a more complex law is
required to describe the dust attenuation from a generic galaxy, with the
Calzetti et al. (2000) law providing a reasonable basis as an average for
large samples. The diffuse dust attenuation index $n$ expresses a power-law
modifier to the slope of the Calzetti et al. (2000) dust attenuation curve,
allowing for dust attenuations with a flatter slope than the Calzetti et al.
(2000) law.

The optical depth of the diffuse dust $\hat{\tau}_{2}$ is drawn from a
truncated ($\hat{\tau}_{2,\mathrm{min}}=0$, $\hat{\tau}_{2,\mathrm{max}}=4$)
normal distribution with mean $\mu_{\hat{\tau}_{2}}=0.3$ and standard
deviation $\sigma_{\hat{\tau}_{2}}=1$, while the optical depth of the birth-
cloud component $\hat{\tau}_{1}$ is drawn from a truncated normal distribution
of its ratio with $\hat{\tau}_{2}$, having mean
$\mu_{\hat{\tau}_{1}/\hat{\tau}_{2}}=1$, standard deviation
$\sigma_{\hat{\tau}_{1}/\hat{\tau}_{2}}=0.3$ and range
$0<\hat{\tau}_{1}/\hat{\tau}_{2}<2$. The diffuse dust attenuation index $n$ is
uniformly distributed in the range $-1.0<n<0.4$.

Figure 2: Distribution of model colours obtained from the sample of $10^{6}$
galaxies drawn from the Prospector-$\beta$ model. The colours are obtained
using model apparent magnitudes in the KiDS-VIKING $ugriZYJHK_{s}$ bands. We
display the 16th, 50th and 84th quantile values for each colour at the top of
the corresponding 1D histogram. The dotted lines in each 1D histogram are the
16th and 84th quantile values. The 2D contours enclose the $68\%$, $95\%$ and
$99\%$ of the values.

Prospector-$\beta$ includes also a prescription to model the emission of
infrared radiation coming the dust that has been heated by the starlight. This
makes use of the Draine & Li (2007) dust emission templates that are based on
the silicate-graphite-PAH model of interstellar dust (Mathis et al., 1977;
Draine & Lee, 1984). This model has three free parameters: $U_{\rm min}$,
representing the minimum starlight intensity to which the dust mass is
exposed, $\gamma_{e}$, representing the fraction of dust mass which is exposed
to this minimum starlight intensity and $Q_{\mathrm{PAH}}$, describing the the
fraction of total dust mass that is in PAHs. $U_{\rm min}$ and $\gamma_{e}$
have uniform priors in the ranges $0.1<U_{min}<15$ and $0<\gamma_{e}<0.15$,
respectively, while $Q_{\mathrm{PAH}}$ is drawn from a truncated normal
distribution having mean $\mu_{Q_{\mathrm{PAH}}}=2$ and
$\sigma_{Q_{\mathrm{PAH}}}=2$.

Accurate forward-modelling of a galaxy SED requires also the addition of an
AGN component. In the Prospector-$\beta$ framework this is achieved using AGN
templates from the Nenkova et al. (2008a, b) CLUMPY models. These templates do
not reproduce the broad and narrow-line region emission lines commonly found
in quasar spectra, but they are created using radiative approximations by
shining an AGN broken power-law spectrum through a clumpy dust torus medium.
This model is characterised by two variables: $f_{\mathrm{AGN}}$, describing
total luminosity of the AGN expressed as a fraction of the bolometric stellar
luminosity, and $\tau_{\mathrm{AGN}}$, which is the optical depth of an
individual dust clump at $5500\ \AA$. In the Prospector-$\beta$ model, a log-
uniform, redshift independent prior is adopted for both parameters,
specifically in the ranges $-5<\log{(f_{\mathrm{AGN}})}<\log{3}$ and
$\log{5}<\log{(\tau_{\mathrm{AGN}})}<\log{150}$.

For our work, we sample $10^{6}$ sets of galaxy physical and stellar
population properties from the Prospector-$\beta$ model prior described above
and we use those to generate galaxy observed-frame SEDs with FSPS. We then
integrate the SEDs in the KiDS-VIKING $ugriZYJHK_{s}$ bands to obtain a sample
of model apparent AB magnitudes and colours. The colours and the redshift
distributions of samples selected by colour-magnitude from the
Prospector-$\beta$ model constitute our fiducial estimates against which the
impact of the modifications of the stellar population components will be
evaluated. Figure 2 shows the distribution of colours resulting from the
Prospector-$\beta$ model.

##  3 Description of the stellar population components

The evaluation of the impact of the stellar population modelling components on
the galaxy tomographic redshift distributions is motivated by the fact that
there are still numerous phases of stellar evolution that are subject of
intense research, being not well understood on both theoretical and
observational grounds (see Conroy et al. 2009 for an in-depth discussion).
This evaluation is carried out by modifying the stellar population components
and their parameters to generate new sets of colours, and therefore redshift
distributions through the colour-redshift relation, that are compared against
the fiducial estimates.

We generate galaxy magnitude samples in the $ugriZYJHK_{s}$ bands, separately
varying all the stellar population modelling components implemented in FSPS at
fixed physical properties of galaxies, for a total of 21 sets of apparent
magnitudes per galaxy. We also consider the case where all the modelling
components are varied at the same time within the ranges set for each
individual component. We label this augmented version ‘P-$\beta$+’. The
components we vary are the gas emission, the dust attenuation and emission,
the AGN component, the initial mass function, the velocity dispersion, the
specific frequency of blue straggler stars, the fraction of blue horizontal
branch stars, the post-asymptotic giant branch phase, the circumstellar
asymptotic giant branch dust emission, the thermally pulsating asymptotic
giant branch normalisation scheme, the presence of Wolf-Rayet stars and the
inter-galactic medium absorption. The new sets of magnitudes are then used to
generate the galaxy colours, necessary to model the colour-redshift relation
with the SOM (see section 5).

The following sections provide a description of the physical details and the
prescriptions we use to model the stellar population components we vary in our
work. We also report the observationally motivated ranges over which the
stellar population component parameters are varied.

###  3.1 Gas emission

The radiation emitted by young stellar populations ionize the surrounding gas
in galaxies, producing nebular emission. The latter has two components: a
nebular continuum, generated by Bremsstrahlung, recombination and two-photon
emission processes, and a nebular line emission, generated by recombination
and electronic transitions from forbidden and fine-structure lines. The amount
of emission from the two components depends mostly on the strength of the
ionizing field and on the metallicity of the gas. This is why many codes rely
solely on these two parameters to predict the continuum and the line emission
fluxes. FSPS includes the nebular contribution on top of the stellar
population one by means of the CLOUDY (Ferland et al., 2017) spectral
synthesis code for astrophysical plasmas. CLOUDY measures the number of
ionizing photons from single stellar populations (SSPs) younger than $10\
\mathrm{Myr}$, removes them from the final SED, and uses their energy to
compute line luminosities and nebular continua. These are slow computations
due to the complex set of equations they have to solve, therefore FSPS
implements a grid of pre-computed line luminosities and nebular continua from
Byler et al. (2017) that depends on SSP age, SSP and gas-phase metallicity,
and ionization parameter $U$.

The nebular line emission can contribute between $20\%$ and $60\%$ to the UV-
optical rest-frame broadband fluxes and is responsible for nearly all the
optical emission lines present in the spectra of star-forming galaxies (Anders
& Fritze-v. Alvensleben, 2003; Reines et al., 2010). The higher the sSFR of
the galaxy, the higher is the contribution of the emission lines to the broad-
band photometry. The gas metallicity impacts the strength of the emission
lines, while the ionization parameter has an effect on the relative intensity
of emission lines (Leja et al., 2017), contributing less to Balmer lines (e.g.
$\mathrm{H\alpha}$, $\mathrm{H\beta}$), and more to forbidden lines (e.g.
$\mathrm{[OII]}$, $\mathrm{[OIII]}$). The nebular continuum emission
contributes instead more significantly to the near-infrared fluxes (Byler et
al., 2017).

The observations that the gas-phase metallicity depends on the galaxy stellar
mass dates back to the early results of SDSS (Tremonti et al., 2004). Recent,
higher redshift studies, have shown that the gas-phase metallicity evolves
with lookback time (Bellstedt et al., 2021) and that there is a more complex
dependence between gas-phase metallicity, stellar mass and SFR (Magrini et
al., 2012; Curti et al., 2020; Bellstedt et al., 2021). In
Prospector-$\alpha$, Prospector-$\beta$ and in general SPS-based forward-
modelling, various prescriptions are employed to assign physically motivated
values to the gas-phase metallicity and to the ionization parameter. In Alsing
et al. (2023), the metallicity of a galaxy is built up with the stellar mass
production, such that the present-day value of the stellar metallicity is
equal to the gas-phase metallicity. In Leja et al. (2017), the stellar and the
gas-phase metallicity are also set to be equal, without allowing for the
chemical evolution of the galaxy. Additionally, there are other studies (e.g.
Leja et al. 2019b; Wang et al. 2023b) that decouple the two metallicities,
with the stellar one being drawn from the Gallazzi et al. (2005) mass-
metallicity relation and the gas-phase metallicity drawn from a uniform
distribution of sub-solar and super-solar abundances, with ranges taken from
Byler et al. (2017). The ionization parameter is instead generally kept fixed
at $\log{U}=-2,-1$ (Leja et al., 2019b; Wang et al., 2023b) or related to the
SFR as in Alsing et al. (2023), where the authors adopt the Kaasinen et al.
(2018) relation between the gas ionization and the SFR. Decoupling the stellar
and the gas-phase metallicity and setting the ionization parameters to a
higher value are choices that are motivated by observations suggesting that
the gas in high-redshift star-forming galaxies experiences stronger ionizing
fields (Tripodi et al., 2024), e.g. due to the higher cosmic SFRD (Madau &
Dickinson, 2014), and gas-phase metallicity abundances that differ
significantly from stellar ones (Shapley et al., 2015; Steidel et al., 2016).

In order to evaluate the impact of the different choices in the gas modelling
on the tomographic redshift distributions, we generate three different samples
of galaxies. In the first sample, which we label as ‘No gas emission’, we
evaluate the extreme case where we turn the nebular continuum and the nebular
emission lines off, such that the emission from the galaxy contains only the
stellar flux, the AGN dust torus and the dust emission. The second sample,
labelled as ‘$Z_{*}=Z_{\mathrm{gas}}$’, tests the prescription of equal values
of the stellar and the gas-phase metallicity against decoupling the two,
without allowing for a chemical evolution in the galaxy. The third sample
tests the effect of varying the ionization parameter with respect to the fixed
$\log{U}=-2$ of the Prospector-$\beta$ model. We label this sample as
‘Variable $\log{U}$’. We draw the ionization parameter for each galaxy from a
uniform distribution in the range $-4\leq\log{U}\leq-1$, a range that is used
in Byler et al. (2017) and is consistent with ionization parameters observed
in local starburst galaxies (Rigby & Rieke, 2004).

###  3.2 Dust attenuation

Dust is a key component of the interstellar medium (ISM) of galaxies. It forms
in the atmospheres of evolved stars or remnants of supernovae and gets
released into the ISM (Draine, 2003), where it presents itself primarily in
two forms, carbonaceus dust grains, like PAHs, and silicate dust grains.
Observationally, dust has two main effects: it modifies the stellar continuum,
particularly in the ultraviolet to near-infrared regime, by absorbing and
scattering starlight and re-emitting the absorbed energy in the infrared, with
continuum and various emission and absorption features.

Observations (see Salim & Narayanan 2020 for a recent review) show that
attenuation laws exhibit a wide range of behaviours and we might expect them
to vary even among similar mass galaxies at different redshifts (Salim &
Narayanan, 2020). Therefore, the modelling of the internal dust attenuation
may dramatically impact the physical properties we derive for individual
galaxies (e.g. in SED fitting, Kriek & Conroy 2013; Salim et al. 2016) since
it involves not only absorption and scattering out of the line of sight, but
also the complex star-dust geometry in a galaxy and the scattering back into
the line of sight. Additionally, we need to provide a prescriptions for the
galaxy population modelling that encompass the wide range of behaviours seen
in the observations.

In section 2, we introduce the two-component dust attenuation model used by
Prospector-$\beta$ in which the radiation from all stars is subject to the
attenuation by a diffuse dust component, with stars below the dispersal time
of the birth cloud ($10\ \mathrm{Myr}$, Blitz & Shu 1980) seeing an additional
source of wavelength-dependent attenuation. To evaluate the impact of
different dust attenuation prescription on the tomographic redshift
distribution means and scatters, we generate samples of galaxies for two
additional widely-used dust attenuation laws: the Calzetti et al. (2000) and
the Milky Way (MW) attenuation laws. We label these samples as ‘Calzetti
$k(\lambda)$’ and ‘MW $k(\lambda)$’, respectively. Both of these attenuation
laws are applied to all starlight equally in FSPS.

The Calzetti et al. (2000) attenuation law is controlled by a single parameter
$\hat{\tau_{2}}$ that is related to the ratio of total to selective extinction
$R_{V}=A_{V}/E(B-V)$. This parameter sets the overall normalisation of the
attenuation curve which has a $\lambda^{-1}$ dependence in the $6300\ \AA\
\leq\lambda\leq 22000\ \AA$ range and a third degree polynomial wavelength-
dependence in the $1200\ \AA\ \leq\lambda<6300\ \AA$ range (stronger effect
towards bluer wavelengths). We sample $\hat{\tau_{2}}$ from the same
distribution as that in the Prospector-$\beta$ model, meaning a truncated
normal distribution with mean $\mu_{\hat{\tau_{2}}}=0.3$ and standard
deviation $\sigma_{\hat{\tau_{2}}}=1$.

The Milky-Way attenuation law follows the prescription in Cardelli et al.
(1989) that is parametrised in FSPS by means of two parameters: the first one
controls the ratio of total to selective extinction $R_{V}=A_{V}/E(B-V)$,
while the second parameter controls the strength of the $2175\ \AA$ UV bump, a
spectral feature on the interstellar extinction curve that is widely seen in
the MW and nearby galaxies, but whose origin remains largely unidentified
(Wang et al., 2023c). Despite being widely assumed fixed at $R_{V}=3.1$, there
have been indications in the literature that the the ratio of total to
selective extinction deviates from this value, both in the MW itself (Cardelli
et al., 1989) and in the Large Magellanic Cloud (Maíz Apellániz et al., 2017),
as well as in samples of star-forming galaxies from SDSS (Sextl et al., 2023).
In order to conservatively select a range of values that has been measured
through the observation of individual stars line of sight, we model the
distribution of $R_{V}$ values using the latest results (Zhang et al., 2023)
from the LAMOST survey (Zhao et al., 2012). The measurements in the MW show
that the $R_{V}$ distribution can be modelled with a Gaussian having mean
$\mu_{R_{V}}=3.25$ and standard deviation $\sigma_{R_{V}}=0.25$. The strength
of the UV bump is instead kept fixed at the Cardelli et al. (1989)
determination for the MW.

###  3.3 Dust emission

More than $30\%$ of the starlight energy in the Universe is re-radiated by
dust at infrared and far-infrared wavelengths (Bernstein et al., 2002). When
modelling the dust emission, its continuum spectrum is assumed to be close to
that produced by a grey body with a single temperature (see e.g. Hensley &
Draine 2021). However, the specific grain composition (typically carbonaceous
PAHs and silicate grains) gives rise to different emission and absorption
features due to e.g. vibrational modes, in this otherwise featureless
continuum (see Draine & Li 2007 for a detailed discussion).

Prospector-$\beta$, and subsequently FSPS, implements only the Draine & Li
(2007) dust emission model that is parametrised by the radiation field
strength and the fraction of grain mass in PAHs. Given that the dust emission
starts to make a significant contribution to the galaxy flux in the mid-
infrared bands, we do not expect it to have large impact on the tomographic
bins selected by optical and near-infrared colours. In order to evaluate
whether the dust emission has any impact on the tomographic bins via the
reddest wavelength bands, we create a sample of galaxies, labelled as ‘No dust
emission’, for which the extreme modelling choice is that the contribution of
dust emission has been completely neglected.

###  3.4 Active galactic nuclei

AGNs create a stronger emission than the regular nuclei of galaxies. This
extra emission is nowadays universally accepted to be generated by the
presence of a supermassive black-hole (SMBH, $\gtrsim 10^{6}M_{\odot}$)
actively accreting matter. AGN have very high luminosities
($L_{\mathrm{bol}}\lesssim 10^{48}\ \mathrm{erg/s^{-1}}$), strong evolution of
their luminosity function with redshift (Fotopoulou et al., 2016; Kulkarni et
al., 2019), and emission that covers basically the whole electromagnetic
spectrum (Padovani et al., 2017). An AGN is detectable in the optical/UV due
to emission from the accretion disk, mostly in the form of broad and narrow
emission lines strongly ionized by the central engine, in the infrared due to
thermal emission of the dust heated by the central engine, at X-ray energies
due to the corona around the accretion disk, and in gamma-ray and radio due to
non-thermal radiation related to the jet. For our study, the most interesting
region where the AGN emission is thought to make a contribution are the
optical/UV and infrared. We expect AGN to modify galaxy colours with a degree
that depends on the relative contribution of the AGN to the overall galaxy
luminosity.

The AGN component implemented in FSPS comes from the Clumpy AGN dust torus
model (Nenkova et al., 2008a, b), which has proven to be successful in fitting
the mid-infrared observations of AGN in the nearby universe (Mor et al.,
2009). These templates model only the emission from the dust torus surrounding
the accretion disk. However, AGN spectra in the optical/UV are also
characterised by the power-law emission from the accretion disk and by the
presence of broad and/or narrow emission lines coming from the strongly
ionized material in the accretion disk. These lines are not modelled inside
the Nenkova et al. (2008a, b) templates and thus neither are they in FSPS.
Furthermore, the prior values assigned to the $f_{\mathrm{AGN}}$ parameter
(see section 2) in Prospector-$\beta$ lead to a fraction of galaxies hosting
AGNs (according to the criterion $f_{\mathrm{AGN}}>0.1$, Leja et al. 2018)
that is roughly $50\%$ smaller than what the observations predict (Thorne et
al., 2022b). In Leja et al. (2018), the authors point out that the
$f_{\mathrm{AGN}}$ parameter prior do not aim to model the fraction of
galaxies hosting an AGN, but rather to fit the SEDs of galaxies hosting AGNs.
Therefore, the values are motivated by the range of observed values for the
power-law distribution of black hole accretion rates (Georgakakis et al.,
2017). Furthermore, as pointed out in Leja et al. (2018), the range is quite
uncertain due to the nonlinear correlation between the AGN bolometric
luminosity and the accretion rate (Shakura & Sunyaev, 1973), and the fact that
the fraction of AGN bolometric luminosity re-processed into MIR emission by
the surrounding environment is likely highly variable (Urry & Padovani, 1995).
A more informative prior range based on the comparison between the AGN and the
galaxy luminosity function is therefore desirable, possibly based on what
observations on large samples of galaxies predict (Thorne et al., 2022b).

To evaluate the impact of AGN models on the color-redshift relation, we
generate four different samples of galaxies. In one, we modify the
Prospector-$\beta$ prior by setting the AGN luminosity fraction to zero
$f_{\mathrm{AGN}}=0$, such that no AGN emission is present. We label this
sample as ‘No AGN’. In the other three samples, we substitute the Nenkova et
al. (2008a, b) templates with the parametric quasar SED model presented in
Temple et al. (2021). This model has been calibrated using the median observed
$ugrizYJHKW12$ colours of the SDSS DR16 quasar population (Lyke et al., 2020)
and aims to cover the full AGN contribution in the $912\ \AA\ \leq\lambda\leq
30000\ \AA$ rest-frame wavelength range. It models the blue continuum in the
$900\ \AA\ \leq\lambda\leq 10000\ \AA$ range due to the low-frequency tail of
the direct emission from the accretion disc, the emission from the hot dust
torus in the $10000\ \AA\ \leq\lambda\leq 30000\ \AA$ range, the Balmer
continuum at $\sim 3000\ \AA$, the broad and narrow emission lines in the
optical/UV region, and the Lyman-absorption suppression, which has the effect
of setting the model flux to zero at all wavelengths below the Lyman limit at
$\lambda_{\mathrm{LLS}}=912\ \AA$.

In swapping the Nenkova et al. (2008a, b) templates with the Temple et al.
(2021) ones, the only parameter of the Prospector-$\beta$ model we need is
$f_{\mathrm{AGN}}$, the total luminosity of the AGN expressed as a fraction of
the bolometric stellar luminosity. We compute the bolometric luminosity for
the AGN component as
$L_{\mathrm{bol,AGN}}=f_{\mathrm{AGN}}L_{\mathrm{bol,galaxy}}$ and we follow
the prescription in the Temple et al. (2021) templates code
repository111https://github.com/MJTemple/qsogen/ to relate the AGN bolometric
luminosity to the monochromatic $3000\ \AA$ continuum luminosity and to the
absolute $i$-band magnitude at $z=2$, as defined by Richards et al. (2006),
two quantities that are necessary to the code to produce the AGN spectrum with
the right units and normalisation. We generate three samples of galaxies with
the Temple et al. (2021) templates. The first one, labelled as ‘Temple+21
QSO’, uses the Baldwin effect (Baldwin, 1977) to control the balance between
strong, peaky, systemic emission and weak, highly skewed emission of the high-
ionisation UV lines as function of redshift. In the second sample, labelled as
‘Temple+21 QSO RELB’, the balance is randomly drawn to allow for more
diversity in the population of galaxies222We follow the author suggestion in a
private communication.. In the third sample, we only use the narrow emission
line templates to reproduce the population of galaxies hosting AGNs where the
dust torus blocks the view to the broad-line region. We label this latest
sample as ‘Temple+21 Nlr’. The motivation in using the Temple et al. (2021)
templates resides in a better agreement of the latter with observed AGN
colours with respect to the use of the Nenkova et al. (2008a, b) templates
(see Appendix A), which is due to the latter containing only the contribution
from the dust torus emission.

Figure 3: Exemplary illustration of the directions in which a galaxy can move
in colour-magnitude space when we vary a single stellar population component
at a time with respect to the fiducial case (black point). We only show the
variants that induce the strongest change in the $g-r$ vs. $r$ colour-
magnitude plane.

###  3.5 Initial mass function

The stellar initial mass function (IMF) describes the distribution of birth
masses of stars and it influences most observable properties of stellar
populations and thus galaxies. The first observational characterisation of the
IMF for stars in the Milky-Way is that of Salpeter (1955), who find that the
IMF can be described as a power-law of the stellar mass with slope
$\Gamma=-2.35$ for $M_{*}\geq 0.5$. More modern IMFs have been described in
the works of Kroupa (2001) and Chabrier (2003), where the IMF power-law has
the same slope than the Salpeter (1955) one at high stellar masses, but fewer
low-mass stars. These three IMFs are the most commonly used ones in extra-
galactic studies and are generally referred to as MW-like IMFs. One of the key
assumptions in modelling them has been their universality, such that galaxies
in the close and far Universe were assumed to have the same IMF also found in
the Milky Way, meaning that they do not dependend on redshift, morphological
type, local metallicity, star-formation rate, or any other parameter. However,
observations conducted in the last decade using spectroscopic, dynamical, and
gravitational lensing measurements have instead shown the existence of a
variety of IMFs. Massive ETGs, for instance, show deviations from the Milky
Way in the form of more bottom-heavy IMFs (Smith, 2020; Parikh et al., 2024),
i.e. excess of low-mass stars. Other recent studies support instead the
existence of more top-heavy (more high mass stars produced) IMFs in nebular
dominated galaxies in the early Universe (Cameron et al., 2023), as predicted
by theoretical works (Chon et al., 2021; Sneppen et al., 2022).

In SPS, the IMF determines the relative weights that are assigned to the
different parts of the isochrones. This implies that if a broad-band colour is
dominated by a single phase of stellar evolution, then the IMF has little to
not impact to it. However, if a broad-band colour is determined by a mixture
of different evolutionary phases, then the IMF will down-weight or up-weight a
specific phase, leading to IMF-sensitive changes in the broad-band colours.
Typical colour changes induced by the IMF are of the order of 0.1 mag in the
near-IR at intermediate ages (see figure 7 in Conroy et al. 2009). In addition
to colours, the rate of luminosity evolution of a galaxy is also sensitive to
IMF, being strongly related to the logarithmic slope of the IMF (see figure 8
in Conroy et al. 2009). Therefore, different IMFs can induce different
predicted colours and mass-to-light ratios, which, in turn, lead to different
predictions for physical properties of galaxies if a wrong IMF is assumed in
the analysis.

Prospector-$\beta$ models the stellar populations with a Chabrier IMF
(Chabrier, 2003). We therefore evaluate the impact of the IMF on the mean
redshift of the tomographic bins by generating samples of galaxy properties
and apparent magnitudes with two other kinds of IMF, a Salpeter (Salpeter,
1955) and a Kroupa IMF (Kroupa, 2001), both implemented in FSPS. We label
these two samples as ‘Salpeter IMF’ and ‘Kroupa IMF’, respectively.

###  3.6 Velocity dispersion

The velocity dispersion of a galaxy, which quantifies the depth of its
potential well, is not expected to make a significant contribution to broad-
band colours. The observed velocity dispersion arises from the superposition
of many individual stellar spectra that are Doppler shifted because of the
star’s motion within the galaxy. In an integrated spectrum, such as that of a
galaxy, the observed velocity dispersion is going to be similar to that of the
stars dominating the light from the galaxy. This implies that for integrated
spectra dominated by multiple stellar populations and kinematics (e.g. bulge
and disk component), the determination of the velocity dispersion is extremely
complex. Hence, many studies determine the velocity dispersion only for
spheroidal systems whose spectra are dominated by the light of red giant
stars.

The velocity dispersion effect is modeled in FSPS as a Gaussian broadening of
the spectrum in velocity space in units of $\mathrm{km\ s^{-1}}$. The user
provides as input the widths in terms of $\sigma$ of the smoothing Gaussian.
In the fiducial case, the spectra are not smoothed by this Gaussian kernel and
therefore they are at the dispersion set by the theoretical or empirical stars
with which the stellar templates have been produced. To allow for a
variability of the stellar velocity dispersion $\sigma_{*}$ values that is
physically motivated, we use the redshift-dependent stellar mass-stellar
velocity dispersion relation introduced in Zahid et al. (2016):

$\begin{split}\sigma_{\mathrm{disp}}(M_{*})&=\sigma_{b}\left(\frac{M_{*}}{M_{b}}\right)^{\alpha_{1}}\
\mathrm{for}\ M_{*}\leq M_{b}\\\
\sigma_{\mathrm{disp}}(M_{*})&=\sigma_{b}\left(\frac{M_{*}}{M_{b}}\right)^{\alpha_{2}}\
\mathrm{for}\ M_{*}>M_{b}\end{split}$ (5)

where the redshift-dependent values of
$M_{b},\sigma_{b},\alpha_{1},\alpha_{2}$ are taken from Table 1 in Zahid et
al. (2016), obtained using Sloan Digital Sky Survey and Smithsonian Hectospec
Lensing Survey data at $z<0.7$. The sample of galaxies generated with this
prescription is labelled as ‘Variable $\sigma_{\mathrm{disp}}$’.

###  3.7 Blue horizontal branch stars

The population of horizontal branch (HB) stars consists of objects with mass
similar to that of the Sun that are burning helium in their cores via the
triple-alpha process and hydrogen in a shell surrounding the core via the
Carbon-Nitrogen-Oxygen cycle. These stars present a nearly constant bolometric
luminosity (Sweigart, 1987; Lee et al., 1990, 1994) which is largely driven by
the constant mass of the Helium core when they enter this phase. A sub-
population of the HB stars are the blue HB stars. They are a late stage in the
evolution of low-mass stars ($0.8M_{\odot}\lesssim M_{*}\lesssim
2.3M_{\odot}$) found on the blue (hot) side of RR Lyrae, at the bottom of the
instability strip (Culpan et al., 2021).

Observationally, blue HB stars are difficult to study for two main reasons:
the paucity in the solar neighbourhood (Jimenez et al., 1998) and the
characteristics of their spectra. Most blue HB stars are expected to have
spectra than resemble those of main-sequence A-type and late B-type stars (Xue
et al., 2008). The main differences are their stronger Balmer jump, lack of
metal spectral lines (indication of low-metallicity), and stronger and deeper
Balmer lines (Smith et al., 2010), which can complicate the estimation of
quantities such as stellar age and metallicity (Lee et al., 2002). These
features make blue HB stars similar in broad-band colours to A- and B-type
stars, only distinguishable with sufficient spectral resolution and signal-to-
noise.

FSPS includes the modelling of blue HB stars, however their contribution to
the galaxy SED is turned off in Prospector-$\beta$. The blue HB stars
contribution is added for stellar populations older than $5\ \mathrm{Gyr}$ and
is parametrised by a free parameter $f_{\mathrm{BHB}}$ that specifies the
fraction of HB stars that are in this blue component. Following the
prescription in Conroy et al. (2009), we draw this parameter for each galaxy
from a uniform distribution in the range $0<f_{\mathrm{BHB}}<0.5$, which we
label as ‘Variable $f_{\mathrm{BHB}}$’. The addition of blue HB stars is
expected to have a relatively stronger effect in the blue rest-frame colours
rather than in the redder ones, not significantly altering the evolution of
colors at late times, but rather giving rise to a constant blueward offset.

###  3.8 Blue stragglers

Blue stagglers are old, population II main-sequence stars extending blueward
of the main-sequence turnoff point at higher luminosity, firstly identified by
Sandage (1953). These objects have been found in all stellar systems, from
globular (Piotto et al., 2004; Salinas et al., 2012) and open clusters (Milone
& Latham, 1994; Ahumada & Lapasset, 1995, 2007; Rain et al., 2021) to the
field population of the Milky Way (Preston & Sneden, 2000; Santucci et al.,
2015). They are particularly abundant in open clusters (Mathieu & Geller,
2015) and they seem to violate standard stellar evolution theory in which a
co-eval population of stars should lie on a clearly defined curve in the
Hertzsprung-Russell diagram determined solely by the initial mass. The fact
that they lie off the main-sequence is an indication of an abnormal evolution.
The general consensus is that blue stragglers start as normal main-sequence
stars and then undergo some form of rejuvenation via acquisition of extra
mass. The probable mechanisms for the latter are mass-transfer in a close
binary (McCrea, 1964) or dynamically induced stellar collisions and mergers
(Hills & Day, 1976; Davies et al., 1994). If the first process is dominant,
then blue stragglers should have a non-negligible effect on the galaxy
integrated light since they might be common throughout the galaxy. Conversely,
if the second process is dominant, then their effect is negligible because
collisions are only important in high density environments, e.g. globular
clusters, that contain only a small fraction of the galaxy stellar mass.

Blue stragglers cover a broader range in colours than normal main sequence
stars and display A-type spectra. Li & Han (2007) found that blue stragglers
makes broad-band colours bluer than $\sim 0.05\ \mathrm{mag}$, while Xin &
Deng (2005) found that they contribute up to $\sim 0.2\ \mathrm{mag}$ in the
$B-V$ colours in old open clusters. FSPS models the presence of blue
stragglers in stellar populations older than $5\ \mathrm{Gyr}$ via a parameter
$S_{\mathrm{BS}}=N_{\mathrm{BS}}/N_{\mathrm{HB}}$, that defines the number of
blue stragglers per unit HB stars. In Prospector-$\beta$, the blue stragglers
contribution is to the galaxy SED is turned off. In this work, we let this
parameter vary uniformly in the range $0.1<S_{\mathrm{SB}}<5.0$ (Piotto et
al., 2004; Conroy et al., 2009) when generating a galaxy for our sample. This
range covers typical values for both globular clusters (Piotto et al., 2004)
and field stars (Preston & Sneden, 2000). The expectation is that they would
contribute at most at the $\sim 0.1\ \mathrm{mag}$ level (Li & Han, 2007; Xin
& Deng, 2005; Conroy et al., 2009) in the blue bands given that they are one
order of magnitude less luminous than HB stars. We label the sample of
magnitudes generated varying the fraction of blue stragglers as ‘Variable
$S_{\mathrm{BS}}$’.

###  3.9 Asymptotic giant branch dust

Asymptotic giant branch (AGB) stars represent the last stage in the evolution
of stars with masses in the range $0.1M_{\odot}\lesssim M_{*}\lesssim
8M_{\odot}$. They are very luminous due to helium shell flashes, contributing
up to tens of percent to the integrated light of stellar populations (Kelson &
Holden, 2010; Melbourne et al., 2012; Conroy, 2013; Melbourne & Boyer, 2013),
with photospheric spectra that reflect the mixing of inner core material with
the outer layers occurring during the dredge-up process. During this phase,
significant mass loss occurs as stars eject their envelopes and evolve towards
the white dwarf sequence. The material that surrounds the star due to the mass
loss is observed to be dust-rich (Bedijn, 1987). The degree of influence of
this AGB dust emission on galaxy spectra is still a matter of debate (Kelson &
Holden, 2010; Chisari & Kelson, 2012; Melbourne & Boyer, 2013; Silva et al.,
1998; Martini et al., 2013), but what is known is that its effect manifests
mostly in the mid-infrared. Villaume et al. (2015) found in particular that
dust shells around AGB stars affects galaxy spectra at $\lambda\gtrsim 4\
\mathrm{\mu m}$ by a factor of $\sim 10$ at intermediate ages $\sim 0.1-3\
\mathrm{Gyr}$.

FSPS models the dusty circumstellar envelopes around AGB stars using the
Villaume et al. (2015) model that makes use of the radiative transfer code
DUSTY and of empirical prescription to assign dust shells to AGB stars in
isochrones. This model is active per default in FSPS, and thus
Prospector-$\beta$, and its amplitude is controlled by a free scaling
parameter. We evaluate its effect on galaxy colours by generating a sample of
objects for which the AGB dust model has been switched off. We label this
sample as ‘No AGB dust model’.

Figure 4: Four different modelled observed-frame galaxy spectra belonging to
intrinsically the same galaxy at $z=0.04$. The black line refers to the
spectrum generated with the fiducial Prospector-$\beta$ model parameters,
while the red, green and blue curves show the spectra generated varying the
IMF (Salpeter, 1955), the attenuation law (Calzetti et al., 2000), and the AGN
prescription (Temple et al., 2021), keeping the other model parameters fixed.
The vertical bands mark the regions where some of the strongest emission lines
appear. The emission line names are reported in black. The spectra units are
_not_ arbitrarly rescaled to highlight the effect of varying the stellar
population components at fixed physical properties of galaxies on the
continuum flux.

###  3.10 Post-Asymptotic giant branch stars

Post-asymptotic giant branch (post-AGB) stars are a rapid phase in the
evolutionary path of $0.1M_{\odot}\lesssim M_{*}\lesssim 8M_{\odot}$ stars
that bridges the gap between the AGB phase and planetary nebula phase. It is
one of the least understood phases because its duration ranges between
$10^{4}-10^{5}\ \mathrm{yr}$ for low-mass stars, down to only few
decades/centuries to the most massive objects (Pereira & Miranda, 2007). After
the outer layers have been removed during the mass loss of the AGB phase,
post-AGB stars increases their surface temperature at almost constant
luminosity, moving horizontally along the Hertzsprung-Russel diagram. At some
point, the temperature becomes high enough that the post-AGB star starts
ionizing its circumstellar material and become observable as a planetary
nebula (van Winckel, 2003).

Post-AGB stars have been detected using the IRAS $[12]-[25]$ versus
$[25]-[60]$ colour–colour diagram and the $J-K$ versus $H-K$ diagram (Garcia-
Lario et al., 1997), because they were expected to be located in a colour
region between the AGB stars and the planetary nebulae. During the transition
from the tip of the AGB to the planetary nebula phase, their spectra appear as
those of M-, K-, G-, F-, A- and OB-type post-AGB supergiants for a short
period (van Winckel, 2003), surrounded by dust shells. When the temperature
increases and the dust surrounding the post-AGB stars start to diffuse, the
contribution of post-AGB stars becomes important to the UV flux from a galaxy.
Post-AGB stars are also the ionizing source of low-ionization nuclear
emission-line regions, particularly in early-type galaxies (Belfiore et al.,
2016; Byler et al., 2019).

Post-AGB stars are implemented in FSPS using the Vassiliadis & Wood (1994)
tracks. The default behaviour of the code is to include these tracks as they
are, but the user can freely change the weight given to the post-AGB phase.
Therefore, we use this flexibility to generate a sample of galaxies where we
turn off the contribution of post-AGB stars. We label this sample as ‘No Post-
AGB stars’.

###  3.11 Thermally pulsating asymptotic giant branch stars

Thermally pulsating AGB stars (TP-AGB) represent short-living phases of the
AGB which, together with AGB, HB and blue straggler stars, are important
phases for SPS modelling due to their high bolometric luminosities compared to
main sequence stars. TP-AGB arise from the instability that AGB stars undergo
during the helium shell burning. When the latter receives new helium produced
by the outer hydrogen shell fusion, the helium shell is not immediately able
to adjust to the increased energy output, leading to a thermal runaway. This
instability is typical of stellar populations older than $10^{8}\
\mathrm{yr}$, hence having initial masses $\lesssim 5\ M_{\odot}$. Therefore,
TP-AGB stars tend to dominate the energy output redward of $\sim 1\
\mathrm{\mu m}$ for stellar populations that are several Gyrs old (Maraston et
al., 2006). Their contribution has been debated for decades, but recently Lu
et al. (2024a) reported the detection of TP-AGB star signatures in the rest-
frame near-infrared spectra of three young massive quiescent galaxies at
$z=1-2$.

TP-AGB stars are included in FSPS and their contribution can be altered using
a weight factor. The default set of FSPS parameters uses the Villaume et al.
(2015) TP-AGB normalisation scheme. However, FSPS allows for other two
normalisation schemes to be used, the default PADOVA (Girardi et al., 2000;
Marigo & Girardi, 2007; Marigo et al., 2008) isochrones and the Conroy & Gunn
(2010) normalisation, which we use to generate two additional galaxy samples
labelled as ‘TP-AGB Padova 2007’ and ‘TP-AGB C&G 2010’, respectively. The
Conroy & Gunn (2010) normalisation effect can be further controlled by
shifting the bolometric luminosity $\log{L_{\mathrm{bol}}}$ and the effective
temperature $\log{T_{\mathrm{eff}}}$ of the TP-AGB isochrones with respect to
the calibrated values presented in Conroy & Gunn (2010). The hotter the TP-AGB
star, the bluer is the emission compared to its cooler sibling. Therefore,
shifting $\log{L_{\mathrm{bol}}}$ has primarily an effect on the bands
redwards of the V-band, while shifting $\log{T_{\mathrm{eff}}}$ impacts the
V-K and bluer colours. We keep these values fixed to reflect the calibration
in Conroy & Gunn (2010).

###  3.12 Wolf-Rayet stars

Wolf-Rayet (WR) stars represent a late stage in the evolution of very massive
O-type stars $\gtrsim 25\ M_{\odot}$, with surface temperatures ranging from
20,000 K to around 210,000 K, hotter than almost all other type of star. WR
occupy the same region of the Hertzsprung-Russel diagram as the hottest bright
giants and supergiants. The strong ionizing radiation from these stars,
coupled with the expanding outer envelope of material expelled by their strong
winds, leads to Doppler-broadened recombination emission lines of helium,
carbon, nitrogen, oxygen, and/or hydrogen in the WR spectrum. This outflowing
material and the explosion of WR stars in supernova Ib or Ic are important
components of the enrichment of the interstellar medium of star-forming
galaxies.

WR spectra are characterised by broad emission lines and are classified
according to the ratios of emission line strengths. Though short-lived
($\lesssim 10^{7}\ \mathrm{yr}$), the WR stars can dominate the UV emission
from young starburst galaxies, while in the UVB photometry they cannot be
distinguished from normal hot stars (Crowther, 2007). FSPS contains a spectral
library to model the WR population, which is activated by default. To test its
impact, we generate a sample of galaxies where the WR spectral library has
been turned off and the main default library is used instead. We label this
sample as ‘No WR’.

###  3.13 Intergalactic medium absorption

The bulk of baryonic cosmic matter resides in a dilute reservoir in the space
between galaxies, known as intergalactic medium (IGM). The IGM has both
cosmological and astrophysical implications: it can be used to test models of
structure formation on the smaller comoving scales (Viel et al., 2005), its
shape can bias cosmological parameter inferences from galaxy clustering
(Pritchard et al., 2007), and it is the reservoir of gas from which galaxies
feed. The higher the redshift of a galaxy, the more material its light has to
travel through. This light is absorbed by the gaseous systems that lie along
the line of sight which, at high redshift, are mostly in the form of neutral
hydrogen. This effect is known as IGM absorption and at increasing redshift it
can be so preponderant that all the light blueward of $\mathrm{Ly\alpha}$ at
$1216\ \AA$ becomes invisible to us.

The first modelling of the IGM absorption is detailed Madau (1995). The author
proposed a theoretical model which is able to reproduce the shape of the
extinction curve as function of redshift, finding that the IGM transmission
decreases with increasing redshift and that the scatter at a given redshift
ranges between $20\%$ and $70\%$. Later models by Meiksin (2006) and Inoue et
al. (2014) updated Madau (1995) prescription for the $\Lambda CDM$ model with
weaker absorption at $3.5\lesssim z\lesssim 5$ and stronger absorption at
$z\gtrsim 6$.

Figure 5: The figure panels show the rest-frame colour differences between the
fiducial galaxy sample (Prospector-$\beta$) and the ones obtained modifying
one stellar population component at the time (y-axis labels). The coloured
widths are the sum in quadrature of COSMOS median photometric errors in the
two respective bands to indicate distinguishability on a single galaxy basis.
The vertical black bars and the white dots on top of them represent the median
values of each rest-frame colour difference. The light red boxes show the
interquartile range (75th - 25th percentile), while the whiskers represent the
84th - 16th percentile ranges.

FSPS implements the IGM absorption using the Madau (1995) prescription, with a
parameter that scales the IGM optical depth. The absorption impacts only rest-
frame wavelengths below $1216\ \AA$, which becomes visible in optical broad-
bands, like the $u$-band, only from $z\sim 1.5$ onwards. The default
parameters of FSPS do not activate the Madau (1995) IGM absorption
prescription. Therefore, we generate a sample of galaxies where we include the
IGM absorption keeping the factor that scales the IGM optical depth fixed to
1. We label this sample as ‘Variable $\mathrm{IGM_{abs}}$’.

###  3.14 Combining the various prescriptions

The samples of SEDs detailed in the previous sections have been created
varying one single stellar population component at a time with respect to the
fiducial Prospector-$\beta$ model. However, the observed population of
galaxies in the Universe in reality spans all these properties ranges at the
same time. Therefore, we create an additional sample of galaxies where we aim
at evaluating the net effect on colours and redshift distributions of varying
simultaneously all the stellar population components according to the
prescriptions detailed in each subsection. We label this sample as
‘P-$\beta$+’.

In particular, we randomly generate galaxy SEDs using an IMF among the
Chabrier, Kroupa and Salpeter ones, a dust attenuation law among the Calzetti
et al. (2000), Cardelli et al. (1989) and Charlot & Fall (2000) ones and TP-
AGB stars prescription among the Villaume et al. (2015), Conroy & Gunn (2010)
and PADOVA ones. We replace the Nenkova et al. (2008a, b) templates with the
Temple et al. (2021) ones obtained by modelling SDSS DR16 quasars. The Temple
et al. (2021) templates are controlled by a single parameter
$f_{\mathrm{AGN}}$ that has the same meaning and values range as the one in
the fiducial Prospector-$\beta$ model. The gas ionization is uniformly drawn
from the Byler et al. (2017) range, while the gas metallicity is for some
objects equal to the stellar metallicity and for others decoupled and
uniformly drawn from the Prospector-$\beta$ range. The dust emission, the AGB
dust, the post-AGB stars, the WR stars and the IGM absorption are randomly
turned on and off. The velocity dispersion follows the Zahid et al. (2016)
prescription, while the fraction of blue HB stars and blue stragglers is
uniformly drawn from the observational ranges in Conroy et al. (2009).

Figure 3 shows, for an exemplary galaxy, the directions in colour-magnitude
space towards which a galaxy can move when we vary one stellar population
component at the time. The fiducial $r$-band magnitude and $g-r$ colour are
represented as a black point, while the arrows point towards the set of
magnitudes and colours generated with a different stellar population
component. The P-$\beta+$ component (magenta arrow) represents the net
movement in colour-magnitude space due to all the varied components.

##  4 The Impact of stellar population components on galaxy colours

In this section we evaluate the impact that the different stellar population
components have on both the rest-frame and the observed-frame galaxy colours.
For each of the $10^{6}$ galaxies, we compute the difference between the rest-
frame and the observed-frame colours estimated with the fiducial
Prospector-$\beta$ model and those estimated varying each stellar population
model component,
$\Delta(\mathrm{colour})=\mathrm{colour}_{\mathrm{fiducial}}-\mathrm{colour}_{\mathrm{component}}$.
We compute this difference for the $u-g$, $g-r$, $r-i$, $i-Z$, $Z-Y$, $Y-J$,
$J-H$ and $H-K_{s}$ colours. We select galaxies having $i\leq 25.3$ to
reproduce the Rubin-LSST Y10 lens gold sample cut (Mandelbaum, 2018). The use
of the rest-frame colours allows us to understand the impact of the individual
spectral features on the broad-band integration of galaxy spectra, while the
observed-frame colours have the aim of reproducing the colour variation that
we would observe in a survey that has a similar photometric precision as that
expected from Rubin-LSST. The convolution of the colour difference with the
redshift distribution of our sample of galaxies has the effect of smearing the
colour difference across the wavelength space. The cut is applied on the
$i$-band magnitude in the fiducial sample and this selection of galaxies is
then held fixed for all variants. This ensures that we evaluate the colour
differences on a per-galaxy basis.

Figure 6: Observed-frame colour differences between the fiducial galaxy sample
(Prospector-$\beta$) and the ones obtained by modifying one stellar population
component at a time. The coloured widths are the sum in quadrature of COSMOS
median photometric errors in the two respective bands to indicate
distinguishability on a single galaxy basis. The vertical black bars and white
dots on top of them represent the median values of each observed-frame colour
difference. The light red boxes show the interquartile ranges (75th - 25th
percentiles), while the whiskers represent the 84th - 16th percentile ranges.

###  4.1 The impact on an example galaxy SED

Figure 4 shows the comparison between FSPS generated observed-frame spectra
for a subset of the different choices on the stellar population components
presented in Section 3. The spectra belong to the same galaxy at $z=0.04$.
Different choices have a dramatic impact on the spectral appearance of the
galaxy which, in turn, reflect on the galaxy magnitudes and colours. In this
figure, we compare the effects of varying the IMF (red curve), the attenuation
law (green curve), and the AGN prescription (blue curve) with respect to the
fiducial Prospector-$\beta$ model (black curve), leaving all other parameters
fixed. The Salpeter IMF produces a galaxy spectrum that overall is fainter
with respect to the fiducial case that uses a Chabrier IMF. This leads to
fainter magnitudes that might impact the detection of the galaxy in e.g. a
flux-limited survey, and hence may impact the mix of physical properties of
galaxies detected at a given color. Changing the attenuation law from the
fiducial prescription in Prospector-$\beta$ to the Calzetti et al. (2000)
model instead induces a wavelength-dependent modification of the galaxy
spectrum. The latter is particularly evident towards bluer wavelengths, where
the effect of the attenuation is the strongest. Changing the AGN prescription
to the one of Temple et al. (2021), instead, gives rise to an overall brighter
spectrum due to the continuum increase, but most notably leads to the
appearance of the characteristic broad-line emission of AGN (see e.g. the
$\mathrm{[NII]}$ and $\mathrm{H{\alpha}}$ complex in the light yellow bar).
The galaxy chosen as an example has a high value of the $f_{\mathrm{AGN}}$
parameter to emphasise the dramatic effect of changing the AGN prescription.
The units are _not_ arbitrarily rescaled, in order to highlight the effect of
varying the stellar population components at fixed physical properties on the
galaxy’s flux.

###  4.2 The impact on rest-frame galaxy colours

Figure 5 shows the distribution of the rest-frame colour differences for all
the SED modelling components. We report in the figure the median photometric
errors of the COSMOS2020 sample (Weaver et al., 2022) that has a similar
photometric precision to that expected by Rubin-LSST (see Appendix C). The
distributions are plotted using a boxplot representation to capture the
relevant quantile ranges of the distribution of colour differences. The wider
the distribution and the more offset its median with respect to the fiducial
value, the larger is the effect on colours of the SED model component we vary.
For the rest-frame case, the SED components that seem to have an appreciable
impact on the colour distributions are those related to the dust attenuation,
the AGN prescription, the gas physics and, consequently, the combined
P-$\beta$\+ sample.

Changing the dust attenuation prescription from the Prospector-$\beta$
fiducial two components Charlot & Fall (2000) to the Calzetti et al. (2000)
(‘Calzetti $k(\lambda)$’ label) or Cardelli et al. (1989) (‘MW $k(\lambda)$’
label) ones induces a wavelength-dependent change in the rest-frame colours,
especially in the bluer bands, where the effect of the dust attenuation is the
strongest. The distributions of colour differences for the Calzetti et al.
(2000) attenuation law are moved towards bluer colours than the fiducial case.
The median value of the distributions is beyond the median COSMOS2020 errors
for all the optical colours, while it is within the errors for the near-
infrared ones. The 84th - 16th percentile values also follow the same trend.
For the Cardelli et al. (1989) MW dust attenuation prescription, we have a
similar, but mitigated behaviour. The median value of the distributions are
outside the median COSMOS errors for the $u-g$ and the $g-r$ colours, while
they are within the errors for the other redder colours. Interestingly, the
$i-z$ rest-frame colour distribution is biased towards redder colours than the
fiducial model. The 84th - 16th percentile values are instead wider than the
errors only for the $u-g$ and the $g-r$ colours, and only marginally for the
$i-z$ rest-frame colour.

The AGN prescription effect differs depending on whether we completely remove
the AGN component or on whether we use the Temple et al. (2021) model.
Removing the AGN component (‘No AGN’ label) has an appreciable effect only in
the rest-frame $J-H$ and $H-K_{s}$ bands since the AGN prescription in FSPS
aims at reproducing the dust torus emission that is stronger in the infrared
bands. Replacing the FSPS AGN model with the Temple et al. (2021) model
(‘Temple+21 QSO’ and ‘Temple+21 QSO RELB’ labels) leads to noteworthy changes
in the rest-frame near-UV/optical colours and in the $J-H$ and $H-K_{s}$
colours. The colour difference in the near-UV/optical is mainly due to the
fact that the Temple et al. (2021) prescription aims at modelling also the
emission coming from the gas in the broad-line region that is subject to the
intense ionizing radiation field of the accretion disk. The difference in the
$J-H$ and $H-K_{s}$ colours is instead due to a different treatment of the
dust torus emission that leads to the Temple et al. (2021) model producing
slightly redder colours than the fiducial Prospector-$\beta$ model. By
modelling the emission from the ionized gas in the narrow-line region only
(‘Temple+21 Nlr’ label), the rest-frame optical colour differences are mostly
within the median COSMOS2020 photometric errors, while the different
prescription becomes detectable in the $J-H$ and $H-K_{s}$ colours. This is
due to the lack of dust torus emission when modelling the narrow-line region
only, which predicts bluer colours with respect to the FSPS AGN prescription.

SED component |  Fraction of cell change for the full sample |  Fraction of cell change for the Rubin-LSSY Y10 sample
---|---|---
Kroupa IMF |  1.4% |  2.3%
Salpeter IMF |  14.9% |  23.4%
No AGN |  6.6% |  6.5%
Temple+21 QSO |  57.4% |  43.6%
Temple+21 Nlr |  11.2% |  12.5%
Temple+21 QSO RELB |  57.3% |  43.7%
No gas emission |  43.6% |  63.1%
Variable $\log{U}$ |  17.2% |  25.0%
$Z_{*}=Z_{\mathrm{gas}}$ |  25.2% |  36.5%
No dust emission |  1.5% |  3.1%
Variable $\sigma_{\mathrm{disp}}$ |  0.1% |  0.1%
Variable $S_{\mathrm{BS}}$ |  0.8% |  1.8%
Variable $f_{\mathrm{BHB}}$ |  5.4% |  11.1%
No Post-AGB stars |  0.5% |  0.5%
No AGB dust model |  0.5% |  0.9%
Calzetti $k(\lambda)$ |  89.6% |  80.9%
MW $k(\lambda)$ |  88.6% |  76.3%
No WR |  11.3% |  12.2%
Variable $\mathrm{IGM_{abs}}$ |  16.4% |  6.0%
TP-AGB C&G 2010 |  0.0% |  0.0%
TP-AGB Padova2007 |  0.0% |  0.0%
P-$\beta$+ |  84.2% |  82.2%
Table 1: For each SED modelling component that we vary (left column), the
table shows the fraction of galaxies whose assigned cell in the SOM changes
with respect to the one assigned with the fiducial SED model. The central
column refers to the case where we consider all the $10^{6}$ galaxies in the
sample. The right-hand column refers, instead, to the case where galaxies are
selected using the Rubin-LSST Y10 cut of $i\leq 25.3$ (right-hand column). We
apply the same cut to each component based on the $i$-band magnitudes of the
fiducial sample in order to compare always the same set of galaxies.

The different gas physics prescriptions induce changes in the intensity or
presence of emission lines which in turn reflect on the galaxy rest-frame
colours. Removing completely the modelling of emission lines (‘No gas
emission’ label) results in colour differences that are larger than the median
photometric errors for the $u-g$, $g-r$ and $r-i$ optical colours. In
particular, the $u-g$ and $g-r$ distributions are biased towards bluer
colours, while the $r-i$ towards redder colours. The disagreement with the
fiducial model is particularly prominent for the $r-i$ rest-frame colour,
because this is the wavelength range where the $\mathrm{H{\alpha}}$,
$\mathrm{[NII]}$ and $\mathrm{[SII]}$ lines fall. The effects in the $u-g$ and
$g-r$ colours are instead mainly driven by the complete absence in modelling
the $\mathrm{[OII]}$, $\mathrm{H{\beta}}$ and $\mathrm{[OIII]}$ lines.
Changing the ionization state of the gas results in a variable relative
intensity of the emission lines, which is particularly evident for the
$(\mathrm{[OIII]}\ \lambda\lambda 4959,5007\ \AA\ /\ \mathrm{[OII]}\
\lambda\lambda 3726,3729\ \AA)$ and the $(\mathrm{[OIII]}\ \lambda 5007\ \AA\
/\ \mathrm{H\beta})$ emission line ratios. These ratios fall into the $u$ and
$g$ bands as testified by the 84th - 16th percentile distribution values in
the $u-g$ and $g-r$ rest-frame colours being outside of the median photometric
errors and leading to bluer and redder colours, respectively, if compared to
the fiducial model. The use of a different prescription to assign the gas-
phase metallicity to galaxies leads in a change of the overall intensity of
the nebular continuum, as well as the intensity and presence of certain
complexes of emission lines. The $u-g$ and $g-r$ colours are the ones most
affected by the change in the gas prescription, with distributions that are
biased in the opposite sense with respect to the ‘Variable $\log{U}$’ case.
Besides the lines already mentioned falling into these bands, a visual
inspection of the rest-frame spectra shows that other lines contribute to the
variation in rest-frame colours, for instance the presence or absence of the
$\mathrm{[OIII]}\ \lambda 4363\ \AA$ line.

The distribution in the rest-frame colours for the P-$\beta$\+ case can be
considered as an average of the effects of all the variable SED components
that we implement when creating this sample of galaxies. The optical colours
are dominated by the combined effect of the dust attenuation law and AGN
prescription that bias the P-$\beta$\+ colours towards bluer ones than the
fiducial model. The $Z-Y$, $Y-J$ and $J-H$ colour distributions are within the
median photometric erros, while the $H-K_{s}$ infrared colour is biased
towards bluer colours being an interplay between the Temple et al. (2021)
prescription for the dust torus emission and the lack of dust torus
prescription in the narrow-line AGN modelling when compared to the fiducial
model.

###  4.3 The impact on observed-frame galaxy colours

Figure 6 shows the same set of SED modelling components and colour
differences, but in the observed-frame, thus mimicking the effects we would
observe in a survey that has a similar photometric precision than that
expected from Rubin-LSST. The figure shows that the AGN prescription and the
dust attenuation law are the main drivers of the change in observed-frame
colours also in this case. Their impact exceeds the median photometric errors,
implying detectability on a single galaxy basis. Furthermore, the redshift
leads to dilution of the colour differences across the wavelength space. This
leads, on one hand, to the mitigation of the effects of some of the
components, for instance the prescription for the gas ionization and the gas
metallicity that are now within the photometric errors, while on the other
hand, a wider distribution of colour difference values for the AGN
prescriptions and dust attenuation laws. The P-$\beta$\+ observed-frame
colours follow the trends of the main drivers of the colour changes, namely
the dust attenuation and the AGN prescription. In particular, all the
distributions are biased towards bluer colours with respect to the fiducial
case.

Figure 7: The left-hand panel shows the distributions of the relative
difference in the cell occupation number with respect to the fiducial case.
For each pair of fiducial and variable component model, we consider only the
cells that are populated by both galaxy samples. The relative differences are
shown as
$\left(N_{j,\mathrm{comp}}-N_{j,\mathrm{fiducial}}\right)/N_{j,\mathrm{fiducial}}$
for clarity. The black vertical bars and the white dots on top of them
represent the median value of
$\left(N_{j,\mathrm{comp}}-N_{j,\mathrm{fiducial}}\right)/N_{j,\mathrm{fiducial}}$,
the boxes represent the 75th - 25th percentile ranges and the whiskers
represent the 84th - 16th percentile ranges. The light green and light blue
areas show the regions where the relative difference is within
$-0.1\leq\left(N_{j,\mathrm{comp}}-N_{j,\mathrm{fiducial}}\right)/N_{j,\mathrm{fiducial}}\leq
0.1$ and
$-0.01\leq\left(N_{j,\mathrm{comp}}-N_{j,\mathrm{fiducial}}\right)/N_{j,\mathrm{fiducial}}\leq
0.01$, respectively. The right-hand panel shows instead the distributions of
the differences in the mean redshift per cell between the fiducial case and
the SED components variation cases normalised by $1+z$. The elements on the
plot have the same meaning, but the light green and light blue areas show
instead the regions where the relative difference of cell mean redshift is
within
$-0.001\leq\left(\bar{z}_{j,\mathrm{comp}}-\bar{z}_{j,\mathrm{fiducial}}\right)/(1+\bar{z}_{j,\mathrm{fiducial}})\leq
0.001$ and
$-0.01\leq\left(\bar{z}_{j,\mathrm{comp}}-\bar{z}_{j,\mathrm{fiducial}}\right)/(1+\bar{z}_{j,\mathrm{fiducial}})\leq
0.01)$, respectively.

Concerning the dust attenuation, the distribution of colour differences for
the Calzetti et al. (2000) attenuation law is moved towards bluer colours than
the fiducial case. The median value of the distributions is beyond the median
photometric errors for all the optical colours, while it is within the errors
for the near-infrared ones. The 84th - 16th percentile values are instead all
beyond the errors except for the $H-K_{s}$ colour. For the Cardelli et al.
(1989) MW dust attenuation prescription, we have a similar, but mitigated
behaviour. The median value of the distributions are outside the median
COSMOS2020 errors for the $u-g$ and the $g-r$ colours, while they are within
the errors for the other redder colours. The 84th - 16th percentile values are
instead wider than the errors only up to the $Z-Y$ colour.

Removing the AGN modelling does not lead to changes in the observed colours
that are beyond the photometric errors, while replacing the FSPS AGN model
with the Temple et al. (2021) model leads to noteworthy changes in the
observed-frame colours from the $u-g$ down to the $Y-J$ colour. In this case,
the median of the colour difference distributions is only marginally offset
from zero, but the 84th - 16th percentile distributions are heavily biased
towards bluer colours. By modelling the emission from the ionized gas in the
narrow-line region only, the observed-frame colour differences are within the
photometric errors. Therefore they contribute only marginally to the colour
variation with respect to the fiducial model. The other SED modelling choices
contribute little to the colour variation, with the case of the modelling of
the TP-AGB phase that does not impact at all the studied bands. Therefore, we
will not further discuss these two modelling components in the rest of this
work.

Figure 6 is interesting because it provides us with a quantitative estimate of
how much the colours of a population of galaxies change when we vary the SED
modelling components, thus guiding the evaluation of which SED components can
be left at the fiducial values and which instead require more careful
modelling and constraints from the observations. In turn, this is tightly
linked to the fact that not accurately knowing the colour distributions makes
the modelling of the colour-redshift relation highly uncertain. Indeed, even
if the effect of varying the SED component cannot be detected on a colour for
a single galaxy when using median COSMOS photometric errors, it could however
be of significance on the ensemble colour distribution of a survey and thus on
its redshift distributions.

##  5 The impact of SED modelling on the SOM-based Colour-Redshift relation

In this work, we model the colour-redshift relation using the Masters et al.
(2015, 2017) SOM, which maps the empirical distribution of galaxies in the
multidimensional colour space spanned by the COSMOS $u$ band imaging from
CFHT, the $griz$ imaging from Subaru Suprime Cam, and the $YJH$ imaging from
the UltraVista survey (McCracken et al., 2012). In its original
implementation, the SOM was trained on the COSMOS 2015 (Laigle et al., 2016)
$ugrizYJHK_{s}$ photometry. We use the KiDS-VIKING $ugriZYJHK_{s}$ remapped
version (McCullough et al., 2023) of the Masters et al. (2015, 2017) SOM. In
the previous section, we have assessed the dramatic impact that the different
stellar population modelling choices have on colours. We show in this section
that this directly affects the colour-redshift relation as modelled by the
SOM, by evaluating the impact that different stellar population modelling
choices have on a galaxy SOM cell assignment and on the mean redshift
estimation per SOM cell.

For each galaxy belonging to the different samples that we build by varying
the stellar population components, we perform the SOM cell assignment
procedure using the observed-frame model colours, without adding photometric
noise and assuming equally sized errors in all colours, which is equivalent to
a perfect cell assignment. We show in Appendix B that the use of COSMOS2020
photometric errors has an impact on the SOM cell assignment of a galaxy that
is negligible compared to its variation in colours due to stellar population
choices.

Figure 8: The left-hand panel shows the differences
$\Delta(\bar{z})_{\mathrm{comp}}$ between the mean of the fiducial magnitude-
limited ($i\leq 25.3$) redshift distribution and those obtained from the
magnitude-limited redshift distributions for each varied SED modelling
component. The right-hand panel shows instead the differences
$\Delta(\sigma_{z})_{\mathrm{comp}}$ on the redshift distributions scatter.
The light blue stars refer to the case where the total number of galaxies for
each magnitude-selected sample is matched to the fiducial number of objects
via an overall magnitude offset (Fitted SMF), while the orange circles refer
to the unmatched case (Fixed SMF). We report the mean errors on
$\Delta(\bar{z})_{\mathrm{comp}}$ and $\Delta(\sigma_{z})_{\mathrm{comp}}$
estimated via bootstrap resampling in the black boxes in the bottom right-hand
corner of each subplot. The blue and red bars show the systematic uncertainty
requirements on the redshift distribution mean and scatter for the galaxy
clustering and weak lensing 3x2 pt analysis of Rubin-LSST Y10, as detailed in
Mandelbaum (2018).

Table 1 shows the fraction of galaxies whose assigned cell in the SOM change
with respect to the one assigned with the fiducial model. The left-hand column
shows the fraction when considering the full sample of $10^{6}$ galaxies,
while the right-hand column shows the case of galaxies belonging to the Rubin-
LSST Y10 sample, $i\leq 25.3$. We apply the Y10 cut based on the $i$-band
magnitudes of the fiducial sample in order to always compare the same set of
galaxies. The Y10 cut selects 191025 galaxies. The fractions track the
variation in the colour distribution reported in Sect. 4, with the largest
fraction of cell change in the full sample arising from the attenuation laws.
Indeed, the Calzetti et al. (2000) and Cardelli et al. (1989) attenuation laws
lead to a change of $89.6\%$ and $88.6\%$, respectively, in the assigned cells
with respect to the fiducial case for the full sample, while for the Y10 cut
the values drop to $80.9\%$ and $76.3\%$. The AGN prescription also leads to a
SOM cell assignment variation, which is less dramatic than the attenuation law
case, but still leading to a fraction of cell assignment change of roughly
$57\%$ ($\sim 43\%$) variation in the full sample (Y10 sample) when adopting
the Temple et al. (2021) AGN templates. The fractions are instead lower ($\sim
11\%$) when using only the narrow-line region template. The P-$\beta$\+ case
leads instead to a fraction of cell change of $84.2\%$ for the full sample and
$82.2\%$ for the Y10 sample. Another SED component that generates a sensible
change in the SOM cell assigment is that provided by the removal of the gas
component emission, ‘No gas emission’. Although the effect on the galaxy
colours distribution of this component is smaller than other components, such
as AGN prescription and dust attenuation laws, the colour change they lead to
causes a variation in the SOM cell assignment with respect to the fiducial
case of $43.6\%$ for the full sample and $63.1\%$ for the Y10 sample. This is
due to the unrealistic galaxy population that this case models, which is
consistent with real galaxies only in the case of quiescent objects.

Other stellar population components that give a large ($\gtrsim 10\%$) SOM
assignment change are the ‘Salpeter IMF’ ($14.9\%$ for the full sample and
$23.4\%$ for the Y10 sample), the ‘Variable $\log{U}$’ ($17.2\%$ for the full
sample and $25.0\%$ for the Y10 sample), the ‘$Z_{*}=Z_{\mathrm{gas}}$’
($17.2\%$ for the full sample and $25.0\%$ for the Y10 sample) and the ‘No WR’
($11.3\%$ for the full sample and $12.2\%$ for the Y10 sample) components,
with the latter inducing a change due to its impact on the rest-frame u-band
magnitudes. The other stellar population components lead to variations which
are $\lesssim 10\%$.

Having assigned each galaxy to a SOM cell, we aim at characterising how much
the occupation number and the mean redshift of each cell change with respect
to the fiducial case when using samples generated with different stellar
population components. The left-hand panel of Figure 7 shows the distribution
of the relative difference in the occupation numbers of the SOM cells. For
each pair of fiducial and variable component model, we consider only the cells
that are populated by both galaxy samples. We report the relative difference
as
$\left(N_{j,\mathrm{comp}}-N_{j,\mathrm{fiducial}}\right)/N_{j,\mathrm{fiducial}}$,
with the light green and light blue areas in Figure 7 showing the regions
where the relative difference is within $\pm 1\%$ and $\pm 10\%$,
respectively. From this figure, we see that the only components that do not
lead to a substantial variation in the cell occupation number are those
related to the advanced stages of stellar evolution (blue stragglers, blue HB
stars, AGB, post-AGB, TP-AGB), as well as the ‘Kroupa IMF’, ‘No dust
emission’, ‘Variable $\sigma_{\mathrm{disp}}$’ and ‘Variable
$\mathrm{IGM_{abs}}$’ components.

The right-hand panel of Figure 7 shows instead how the mean of the redshift
distribution in each cell change when one varies the stellar population
components with which galaxy spectra are built. In particular, we report the
difference in the mean redshift between the fiducial model and those obtained
by varying one stellar population component at the time normalised by $1+z$,
$\left(\bar{z}_{j,\mathrm{comp}}-\bar{z}_{j,\mathrm{fiducial}}\right)/(1+\bar{z}_{j,\mathrm{fiducial}})$.
The light green and light blue areas in Figure 7 show the regions where the
relative difference is within $\pm 0.1\%$ and $\pm 1\%$, respectively. The
same components that have percent level variation in the occupation number are
those that have permille level variation in the mean redshift of the SOM
cells, while the other components induce relative variations in the mean
redshift that are in some cases even greater than the percent level.

The results shown in Figure 7 highlight the direct impact that varying a
single stellar population component has on the colour-redshift relation. For
the same galaxy, varying the spectrum generates a variation in colours which
in turn associates the galaxy to a different SOM cell. The combined effect of
the population of galaxies results in a change of cell occupation numbers and
mean redshift per cell that at the cell level is much greater than the percent
and permille nominal precision requirements by Stage III and IV surveys,
respectively.

##  6 The SED modelling impact on the redshift distribution means and scatters

In this section we evaluate the impact of varying the SED modelling components
on the redshift distribution of our sample of simulated galaxies. First, we
evaluate how the entire redshift distribution of the sample of galaxies
changes between SED model variants when we apply a pure magnitude selection.
Then, we use these magnitude-selected samples to create tomographic bins
following the prescription in McCullough et al. (2023). This procedure allows
us to evaluate the impact on the mean and the scatter of the tomographic
redshift bins induced by the variation in the SED modelling components.

In forward-modelling the galaxy population, any misspecification of the SED
model would result in a partial compensation of this impact by the fitted
distribution of physical properties of galaxies, particularly the stellar mass
function (SMF). What effect this has in practice depends on the galaxy data
being used as a constraint and on the degrees of freedom that are given to the
galaxy population model. We introduce three limiting cases with which we try
to evaluate to what degree this may compensate the impact of SED mismodelling
on the colour-redshift relation. We assume that only the SMF that is being
fitted is used to compensate for the SED mismodelling, while other physical
properties such as metallicity and star formation history are kept as-is. The
SMF is indeed the most flexible aspect of the galaxy population model for
fitting the observed abundance of galaxies, which is the primary constraint in
the forward-modeling of photometric surveys. Degrees of freedom on other
galaxy physical property distributions may also have an interplay with SED
modelling, which however is beyond the scope of this work to investigate. The
three SMF cases we consider are:

  * •

Fixed SMF: the different variants keep the distribution of all physical
properties of galaxies fixed and we apply the $i$-band selection separately on
each variant. This leads to different galaxies, and different total numbers of
galaxies, to be selected in each variant.

  * •

Fitted SMF: we apply an overall multiplicative correction of stellar mass that
compensates the SED modelling choice dependent stellar mass-to-light ratio.
This is done in such a way that the overall abundance of our apparent $i$-band
magnitude limited sample is kept fixed.

  * •

Fitted SMF of (type,z): we assume a highly flexible fit of the SMF as a
function of galaxy type and redshift. We realise this practically such that
the exact same set of galaxies is selected into the magnitude limited sample,
regardless of which variant of the SED modelling is being considered.

###  6.1 The impact on the magnitude-limited sample

A magnitude-limited selection of target galaxies is generally performed using
a sharp cut in colours or magnitudes. Since the results presented in Section 4
highlight the dramatic effect that varying the SED components has on the
galaxy magnitudes and colours, we expect a similar effect to be reflected on
the redshift distributions resulting from a magnitude-limited selection. In
the case of a Stage IV galaxy survey like Rubin-LSST, the LSST Dark Energy
Science Collaboration requirements (Mandelbaum, 2018) for the weak lensing and
galaxy clustering analysis on the systematic uncertainty in the mean and
redshift scatter of each source tomographic bin refer to a sample of galaxies,
the Y10 lens gold sample, that has a limiting magnitude of
$i_{\mathrm{lim}}=25.3$. We use this value to select magnitude-limited samples
of objects from each of the 22 sets of galaxies generated by varying the SED
modelling components.

Figure 9: Redshift distribution $N(z)$ of two sets of five tomographic bins
built with the prescription highlighted in Sect. 6. The filled histograms
refer to the $N(z)$s built with the fiducial sample, while the histograms with
the dashed outlines represent the $N(z)$s built using the ‘Calzetti
$k(\lambda)$’ component. The solid and dashed vertical lines represent the
mean of the tomographic bins for the fiducial and the ‘Calzetti $k(\lambda)$’
components, respectively. The legend shows the tomographic bin means and
$1\sigma$ errors on the means. We use the ‘Calzetti $k(\lambda)$’ component to
demonstrate our analysis methodology as it induces the largest shift in the
$N(z)$. The shift is caused by the brighter rest-frame UV emission that the
‘Calzetti $k(\lambda)$’ component induces (see Figure 4). This leaves the low-
redshift galaxy observables unaffected, while it brightens high-redshift
galaxies in the observed optical bands, hence making them more numerous in a
magnitude-limited sample.

We compute the mean $\bar{z}_{\mathrm{comp}}$ and the standard deviation
$\sigma_{z,\mathrm{comp}}$ of each of the 22 resulting magnitude-limited
redshift distributions $N(z)_{\mathrm{comp}}$ and compare those against the
mean $\bar{z}_{\mathrm{fiducial}}$ and the standard deviation
$\sigma_{z,\mathrm{fiducial}}$ obtained for the magnitude-limited fiducial
sample $N(z)_{\mathrm{fiducial}}$,

$\begin{split}\Delta(\bar{z})_{\mathrm{comp}}&=\bar{z}_{\mathrm{fiducial}}-\bar{z}_{\mathrm{comp}},\\\
\Delta(\sigma_{z})_{\mathrm{comp}}&=\sigma_{z,\mathrm{fiducial}}-\sigma_{z,\mathrm{comp}}.\end{split}$
(6)

The errors on these estimates are obtained using bootstrap resampling. The
procedure involves repeatedly sampling with replacement from the original set
of galaxies to create multiple bootstrap samples, applying the magnitude
selection, computing the mean for each bootstrap sample and then the standard
deviation of these bootstrap means. We generate a total of $10^{3}$ bootstrap
samples.

For the Fixed SMF sample, the variation in the stellar population components
generates samples of galaxies that contain a different number of objects
brighter than the $i$-band magnitude cut. The number of galaxies at $i\leq
25.3$ ranges from 148447 galaxies in the case of the ‘Salpeter IMF’ to 243329
galaxies in the case of ‘Temple+21 QSO’. However, in modelling a galaxy
population, an important property that is fitted to the observations is the
total number of galaxies that the model produces. To account for this, the
Fitted SMF sample matches the $\mathrm{d}n/\mathrm{d}m_{i}$ of each of the 22
components to the fiducial one by applying a constant magnitude shift to the
$i$-band. This corresponds to an overall multiplicative correction of stellar
mass that compensates the SED modelling choice dependent stellar mass-to-light
ratio. We estimate the $i$-band magnitude shift for each component variant
such that the total number of galaxies having $i\leq 25.3$ is equal to that of
the fiducial sample (191025). The shifted $i$-band magnitude distributions for
the various SED modelling components are then used to select a sample and
compute the mean and the standard deviation of the magnitude-limited redshift
distributions. We compute the errors on these estimates as in the previous
case, i.e. via bootstrap resampling.

Figure 8 shows the differences $\Delta(\bar{z})_{\mathrm{comp}}$ and
$\Delta(\sigma_{z})_{\mathrm{comp}}$ between the mean (left-hand panel) and
the scatter (right-hand panel) of the magnitude-limited fiducial redshift
distributions and those obtained for each varied SED modelling component. The
light blue stars refer to the $\mathrm{d}n/\mathrm{d}m_{i}$ matched case
(Fitted SMF), while the orange circles to the unmatched case where the
$i$-band selection produces samples of different numbers of galaxies (Fixed
SMF). The light blue and orange error bars in the black boxes in the bottom
right corners show the bootstrap-estimated errors on the mean and the scatter
of each $\Delta(\bar{z})_{\mathrm{comp}}$ and
$\Delta(\sigma_{z})_{\mathrm{comp}}$. In particular, we report the error
averaged across all the evaluated components. The blue and the red bars show
the systematic uncertainty requirements on the redshift distribution mean and
scatter for the galaxy clustering and weak lensing 3x2 pt analysis of Rubin-
LSST Y10, as detailed in Mandelbaum (2018). The requirements state that the
systematic uncertainty in the mean redshift of each tomographic bin should not
exceed $0.001(1+z)$ for the weak lensing 3x2 pt analysis and $0.003(1+z)$ for
the galaxy clustering analysis, while the systematic uncertainty in the
photometric redshift scatter should not exceed $0.003(1+z)$ for the weak
lensing 3x2 pt analysis and $0.03(1+z)$ for the galaxy clustering analysis.

The Fixed SMF case (orange circles) shows the dramatic impact that a sharp
magnitude cut has on the mean and standard deviation of the magnitude-selected
galaxy redshift distributions. We find that a magnitude selection leads to
changes in the mean of the redshift distributions up to $\Delta z\sim 0.1$.
The ‘Kroupa IMF’, the ‘Salpeter IMF’ and the ‘No gas emission’ samples lead to
redshift distributions that are biased low with respect to the fiducial case,
since they tend to shift the $\mathrm{d}n/\mathrm{d}m_{i}$ towards fainter
galaxies that do not make the magnitude cut. On the other hand, the ‘Temple+21
QSO’, the ‘Temple+21 QSO RELB’, the ‘Calzetti $k(\lambda)$’, the ‘MW
$k(\lambda)$’ and the ‘P-$\beta$+’ samples tend to shift the
$\mathrm{d}n/\mathrm{d}m_{i}$ towards brighter galaxies than the fiducial
case, leading to redshift distributions that are biased high.

The mentioned components bias the mean of the redshift distribution with
respect to the fiducial case by an amount that exceeds the requirements for
the Y10 weak lensing and galaxy clustering analysis. The other components
instead cause biases that are within the Y10 requirements. In the case of the
redshift scatter, all biases are always within the requirements for the galaxy
clustering Y10 analysis, while for the weak lensing analysis the components
mentioned above bias the redshift scatter by an amount that is beyond the Y10
requirements.

Matching the total number of galaxies to that of the fiducial case only
partially mitigates the bias on the mean and the scatter of the redshift
distributions. In particular, the only variants that are made compatible by
the matching with both the weak lensing and the galaxy clustering Y10
requirements on the redshift distribution mean are the ‘Kroupa IMF’ and ‘No
gas emission’. Matching the total number of galaxies also leads the ‘Kroupa
IMF’ variant to become compatible with both the weak lensing and the galaxy
clustering Y10 requirements on the redshift distribution scatter.

Figure 10: The figure shows the impact that varying the choices of the SED
modelling component has on the mean of the tomographic redshift distribution
bins. Each panel shows a different tomographic bin, from the lowest (leftmost)
to the highest (rightmost) redshift one. Each row represents a different SED
modelling component we vary. The scatter points represent the difference
between the mean redshift per bin obtained with the fiducial
Prospector-$\beta$ model and that obtained by varying one SED modelling
component at the time. Orange scatter points are obtained by weighting the
$P(z|c_{i})$ for the number of simulated galaxies per SOM cell, while light
blue points are obtained by weighting the $P(z|c_{i})$ for the observed
abundances in the Masters et al. (2015); McCullough et al. (2023) SOM cells.
Circle-shaped, star-shaped and diamond-shaped scatter points refer to the
cases where we perform an $i$-band magnitude-selection of the sample of
galaxies that is different for each component, matched for the total number of
galaxies and based on the fiducial $i$-band distribution, respectively. We
report the mean errors on the difference of the tomographic bin means
estimated via bootstrap resampling in the black boxes in the bottom right-hand
corner of each subplot. The blue and red coloured bands refer to the Rubin-
LSST Y10 galaxy clustering and weak lensing requirements, respectively.

###  6.2 The impact on the tomographic redshift distributions

In order to evaluate whether varying the stellar population modelling
components introduce biases in the mean and scatter of the tomographic
redshift bins, we use the galaxies assigned to the SOM cells to infer the
redshift distributions for five tomographic bins built following the
prescription in McCullough et al. (2023).

####  6.2.1 Tomographic redshift distributions construction

After selecting the magnitude-limited sample of galaxies, we perform the SOM
cell assignment algorithm, such that the $i$th cell $c_{i}$ in the SOM is
populated by an associated number of galaxies $n_{i}$, with a redshift
probability distribution $P(z|c_{i})$. We use the histogram of redshifts of
galaxies in a cell as the model for the $P(z|c_{i})$. We then sort cells
according to their median redshift and create 5 equally populated tomographic
bins by placing galaxies in each tomographic bin according to the median
redshift of the cell they reside in. The redshift distribution of galaxies in
each tomographic bin is then given by the redshift probability distribution in
each cell weighted by the cell occupation,

$N(z|\mathrm{bin})=\sum_{i\ \in\ \mathcal{C}_{i}}n_{i}P(z|c_{i}),$ (7)

where $\mathcal{C}_{i}$ is the set of cells used in the tomographic bin. We
populate the lowest redshift bin first and then we move to the consecutive one
once the total number of galaxies in each bin $N_{\mathrm{tot}}$ is roughly
equal to one fifth of the full sample of magnitude-selected objects populating
the SOM cells. The mean of the galaxy redshift distribution in each
tomographic bin depends on the mean redshift in each colour cell
$\bar{z}_{c_{i}}$ weighted by the number of galaxies per cell $n_{i}$,

$\begin{split}\bar{z}_{\mathrm{bin}}&=\frac{\int z\ N(z|\mathrm{bin})\
\mathrm{d}z}{N_{\mathrm{tot}}}\\\
&=\frac{1}{N_{\mathrm{tot}}}\left[n_{1}\bar{z}_{c_{1}}+...+n_{\mathcal{N}}\bar{z}_{c_{\mathcal{N}}}\right],\end{split}$
(8)

where $N_{\mathrm{tot}}$ is the total number of galaxies in the tomographic
bin. The variance of each tomographic bin is given by

$\sigma^{2}_{z}=\frac{\sum_{i=1}^{\mathcal{N}}n_{i}\sigma_{z,c_{i}}^{2}}{N_{\mathrm{tot}}}+\frac{\sum_{i=1}^{\mathcal{N}}n_{i}(\bar{z}_{c_{i}}-\bar{z}_{\mathrm{bin}})^{2}}{N_{\mathrm{tot}}},$
(9)

where $\sigma_{z,c_{i}}^{2}$ is the redshift distribution variance per cell.
In Figure 9, we show the five tomographic bins built for the fiducial sample
and for the ‘Calzetti $k(\lambda)$’ sample using the prescription highlighted
above. We choose to plot this stellar population variant against the fiducial
to highlight the strong impact that stellar population components might have
on the tomographic redshift distributions. Indeed, the ‘Calzetti $k(\lambda)$’
is one of the components providing the largest impact on the tomographic
redshift distributions. This is due to the brighter rest-frame UV emission
that the ‘Calzetti $k(\lambda)$’ variant induces (see Figure 4). This leaves
the low-redshift galaxy observables unaffected, while it brightens high-
redshift galaxies in the observed optical bands, hence making them more
numerous in a magnitude-limited sample.

Figures 10 and 11 represent the main result of this work. Figure 10 shows the
bias induced on the tomographic redshift distribution bin means when modelling
the galaxy SEDs with different stellar population component choices, while
Figure 11 shows the bias induced on the tomographic bins’ scatter. The blue
and red vertical bars represent the requirements for the Rubin-LSST Y10 galaxy
clustering and weak lensing analysis, respectively, as detailed in section
6.1. The error bars in the black boxes are obtained similarly to section 6.1,
i.e. via bootstrap resampling.

The different points for a given stellar population variant represent
different ways in which the redshift probability distribution in each cell is
weighted towards the bin average. In particular, the orange scatter points are
obtained when weighting the redshift probability distributions $P(z|c_{i})$ in
each cell by the number of simulated galaxies $n_{\mathrm{sim}}$ that occupy
that cell in the same variant. The light blue scatter points, instead, are
obtained when weighting the $P(z|c_{i})$ by the observed KiDS-VIKING abundance
$n_{\mathrm{KV}}$ in the remapped SOM (Masters et al., 2015; McCullough et
al., 2023) cells, similar to what would be the case in a scheme like that of
Myles et al. (2021).

At fixed variation of stellar population variant and scatter point colour, we
generate three different samples of magnitude-limited objects with which we
construct the tomographic bins. The circles represent the case where, for each
different stellar population component, we select the magnitude-limited sample
of galaxies using the $i$-band magnitude distribution from that specific
component (Fixed SMF). This process selects different numbers of galaxies
among the various stellar population components, as pointed out in section
6.1. The star-shaped scatter points represent instead the case where we match
the $\mathrm{d}n/\mathrm{d}m_{i}$ of each of the 22 components to the fiducial
one by applying a constant magnitude shift to the $i$-band magnitude
distributions (Fitted SMF). This process aims at matching the total number of
magnitude-selected galaxies with that of the fiducial case by matching a model
to the observed abundances. The diamond-shaped scatter points instead
represent the case where, for each variation of the stellar population
components, we select the magnitude-limited sample of galaxies using the
$i$-band magnitude distribution of the fiducial case (Fitted SMF of (type,z)).
This process effectively employs always the same population of galaxies, but
with a colour-redshift relation that is established by the SOM and therefore
changes for every stellar population component.

Figure 11: The figure shows the impact that varying the choices of the SED
modelling component has on the scatter of the tomographic redshift
distribution bins. Each panel shows a different tomographic bin, from the
lowest (leftmost) to the highest (rightmost) redshift one. Each row represents
a different SED modelling component we vary. The scatter points represent the
difference between the redshift scatter per bin obtained with the fiducial
Prospector-$\beta$ model and that obtained by varying one SED modelling
component at the time. The orange scatter points are obtained by weighting the
$P(z|c_{i})$ for the number of simulated galaxies per SOM cell, while the
light blue points are obtained by weighting the $P(z|c_{i})$ for the observed
abundances in the Masters et al. (2015); McCullough et al. (2023) SOM cells.
Circle-shaped, star-shaped and diamond-shaped scatter points refer to the
cases where we perform an $i$-band magnitude-selection of the sample of
galaxies that is different for each component, matched for the total number of
galaxies and based on the fiducial $i$-band distribution, respectively. We
report the mean errors on the difference of the tomographic bin scatters
estimated via bootstrap resampling in the black boxes in the bottom right-hand
corner of each subplot. The blue and red coloured bands refer to the Rubin-
LSST Y10 galaxy clustering and weak lensing requirements, respectively.

Figures 10 and 11 show that the variation of the stellar population component
choices with respect to those adopted in the fiducial Prospector-$\beta$ model
affects the derived colour-redshift relation. This, in turn, induces biases in
both the mean and the scatter of the redshift distributions of tomographic
bins that are beyond the requirements set by Stage IV galaxy surveys for
several of the variants. As already shown throughout the work, the components
that make the greatest contribution to the bias are those related to the IMF,
AGN, gas, and dust attenuation modelling, including the P-$\beta$\+ component
that encapsulates the net effect of all the varied SED modelling components.
These variants mostly bias the redshift to higher values.

In the following sections, when presenting the biases induced by each SED
modelling component, we test the consistency with the weak lensing and
clustering requirements using errors estimated via bootstrap resampling. These
errors are somewhat different for every SED component, however in Figures 10
and 11 we show only the errors averaged across all the components for plot
readability.

####  6.2.2 The impact of the initial mass function

The first SED modelling component of which we evaluate the impact on the
tomographic redshift distributions is the IMF. In our work, we test two
different widely used IMFs, the ‘Kroupa IMF’ and ‘Salpeter IMF’ that are
compared against the Chabrier IMF adopted in the fiducial model. The ‘Kroupa
IMF’ component (first row of Figures 10 and 11) does not induce a dramatic
shift in the mean redshift and scatter of the tomographic bins given its
similarity with the Chabrier IMF. When considering the case of different
magnitude-limited samples of galaxies (Fixed SMF, circle-shaped points), the
shift in the means induced by this component is consistent with the fiducial
model within the galaxy clustering requirements in the second and third
tomographic bin for $n_{\mathrm{sim}}$ and in the first and third tomographic
bin for $n_{\mathrm{kv}}$. It is instead consistent with the weak lensing
requirements in the first, fourth and fith tomographic bin for
$n_{\mathrm{sim}}$ and in the fourth and fith tomographic bin for
$n_{\mathrm{KV}}$. If we account for SED mismodelling and for the observed
abundances (Fitted SMF, star-shaped points, and Fitted SMF of (type,z),
diamond-shaped points), the ‘Kroupa IMF’ component becomes consistent with the
fiducial model within the Stage IV weak lensing requirements on the mean
redshift for every tomographic bin. As for the requirements on the redshift
scatter, the ‘Kroupa IMF’ component is consistent with the fiducial model
within the weak lensing requirements for every tomographic bin, with the only
exception of the Fixed SMF weighted by $n_{\mathrm{sim}}$ case that is
consistent with the clustering requirements only in the first tomographic bin.

The ‘Salpeter IMF’ component (second row of Figures 10 and 11) has instead a
slightly different behaviour. Since we have seen in Figure 8 that the
‘Salpeter IMF’ has a strong impact on the $i$-band magnitude and thus
selection, this component in the Fixed SMF case (circle-shaped points) selects
a very different sample of objects with respect to the fiducial ones, thereby
biasing the mean redshift of the tomographic bins well beyond the Stage IV
requirements. Upon matching the total number of galaxies selected to the
fiducial model (Fitted SMF and Fitted SMF of (type,z) cases), the impact is
strongly mitigated in every tomographic bin in such a way that for the highest
degree of SED mismodelling compensation (Fitted SMF of (type,z)) we have
consistency with the fiducial model within the weak lensing requirements on
the redshift mean in every tomographic bin. The impact of the ‘Salpeter IMF’
component on the redshift scatter with respect to the fiducial model is always
at least within the clustering requirements. Also in this case, the highest
degree of SED mismodelling compensation (Fitted SMF of (type,z)) makes this
component consistent with the fiducial model within the weak lensing
requirements.

The results for the IMF show that fitting the SMF in order to match the
observed abundance of galaxies to compensate for any possible IMF-related SED
mismodelling is crucial to obtain redshift distribution estimates that are not
biased beyond the requirements set by Stage IV surveys.

####  6.2.3 The impact of active galactic nuclei

AGN are the SED modelling component that induces the largest bias on the mean
and scatter of the tomographic redshift bins. In our work, we test the impact
of AGN by replacing the Nenkova et al. (2008a, b) templates implemented in
FSPS with the more realistic model presented in Temple et al. (2021). The
‘Temple+21 QSO’ and the ‘Temple+21 QSO RELB’ components (fourth and sixth rows
of Figures 10 and 11) produce redshift distribution estimates that have mean
redshifts biased high with respect to the fiducial ones. In the majority of
cases they are not consistent with the Stage IV requirements. For the Fixed
SMF (circle-shaped points) and the Fitted SMF (star-shaped points) cases, the
‘Temple+21 QSO’ and the ‘Temple+21 QSO RELB’ components are consistent with
the fiducial model neither within the clustering nor within the weak lensing
requirements on the mean redshifts, independently of how we weight the
$P(z|c_{i})$ by the cell occupation. In the Fitted SMF of (type,z) case
(diamond-shaped points), weighting by $n_{\mathrm{sim}}$ (orange points) leads
the ‘Temple+21 QSO’ component to be consistent with the clustering
requirements in the fourth tomographic bin and not consistent with any of the
requirements in the other four bins, while the ‘Temple+21 QSO RELB’ component
becomes consistent with the clustering requirements in the second and fourth
bins, and not consistent in the remaining three. Weighting the Fitted SMF of
(type,z) case by $n_{\mathrm{KV}}$ (light blue diamond-shaped points) leads
the ‘Temple+21 QSO’ to be not consistent with any of the requirements in the
third bin, consistent with the clustering requirements in the fifth bin and