© The Authors 2025

1. Introduction

We are entering a transformative era for cosmology and galaxy evolution, driven by current and upcoming wide-field photometric and spectroscopic surveys. The recent baryon acoustic oscillation (BAO, Blake & Glazebrook 2003; Seo & Eisenstein 2003) measurements from three years of operation of the Dark Energy Spectroscopic instrument (DESI) challenged the standard ΛCDM model by pointing towards a Universe with a time-evolving dark energy equation of state (DESI Collaboration 2025). Meanwhile, cosmic shear measurements from the complete Kilo-Degree Survey (KiDS-Legacy, Wright et al. 2024) have found the $S_8\equiv {\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$ parameter to be consistent with Planck Legacy constraints (Wright et al. 2025a; Stölzner et al. 2025), which potentially resolves the longstanding ‘S₈ tension’ (Abdalla et al. 2022; Amon et al. 2022; Dalal et al. 2023; Dark Energy Survey and Kilo-Degree Survey Collaboration 2023; Di Valentino et al. 2025) between early- and late-time cosmological probes. Euclid has recently released the Quick data release 1 (Euclid Collaboration: Aussel et al. 2025), based on 63.1 deg² of the Euclid Deep Fields (Euclid Collaboration: Mellier et al. 2025) to nominal wide-survey depth. These data span a wide range of science cases, including studies of the star-forming (Euclid Collaboration: Enia et al. 2025) and passive (Euclid Collaboration: Cleland et al. 2025) galaxy populations, as well as the identification of strong gravitational lensing systems from individual galaxies (Euclid Collaboration: Rojas et al. 2025) to massive clusters (Euclid Collaboration: Bergamini et al. 2025). Looking ahead, upcoming surveys such as the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (hereafter, Rubin-LSST, Ivezić et al. 2019) and spectroscopic programmes like 4MOST (de Jong et al. 2019) are expected to deliver the most precise constraints to date on the large-scale structure (LSS) of the Universe. These will significantly improve measurements of key cosmological parameters while also enabling detailed studies of structure growth and galaxy evolution at low (Driver et al. 2019) and intermediate redshifts (Iovino et al. 2023).

Next-generation surveys will probe a regime where the statistical uncertainties of measurements are subdominant compared to the systematic ones. Maximising the scientific return therefore requires a precise characterisation of measurement biases and incompleteness. An exemplary case is the use of weak gravitational lensing (see Mandelbaum 2018 for a comprehensive review) as a powerful probe of the LSS. Its ability to constrain key cosmological parameters, such as the clustering amplitude of the matter power spectrum, σ₈, and the matter density parameter, Ω_m, is highly sensitive to uncertainties in the redshift distribution of source galaxies, especially the mean redshift (Amon et al. 2022; van den Busch et al. 2022; Li et al. 2023; Dalal et al. 2023). Three notable recent studies highlight this challenge. The first is an improved redshift calibration method and sample that significantly enhanced the cosmic shear analysis of the KiDS-Legacy dataset (Wright et al. 2025b). The second involves the Hyper-Suprime Cam (HSC) Year 3 cosmic shear analysis, (Dalal et al. 2023) which struggled to tighten constraints over the Year 1 analysis despite a three-fold increase in sky area due to the wide prior placed on their higher redshift tomographic bins. The third example is that of the Dark Energy Survey (DES), which discarded the two highest redshift bins in the MagLim sample due to the redshift distribution estimates being the dominant systematic error in DES weak lensing measurements (Abbott et al. 2022). Given these limitations, Stage IV experiments (Albrecht et al. 2006), such as Euclid and Rubin-LSST, have set stringent requirements on the systematic uncertainty in the mean redshift of each source tomographic bin. For the Rubin-LSST weak lensing analysis, the systematic uncertainty should not exceed mean redshift errors of |δz|< 0.002 × (1 + z) in the year 1 analysis, and 0.001 × (1 + z) in the year 10 analysis (The LSST Dark Energy Science Collaboration 2018). However, current photometric redshift (photo-z) calibration methods from Stage III dark energy experiments (Myles et al. 2021; Wright et al. 2020; van den Busch et al. 2022) still fall short of these goals. They yield uncertainties on the mean redshift that are degraded by spectroscopic incompleteness and outliers in the calibration redshifts, sample variance, blending, and astrophysical systematics that cause systematic errors in the mean and standard deviation of redshift that are still one order of magnitude greater than what is required for Stage IV experiments (see Newman & Gruen 2022 for a comprehensive review). Efforts are on-going to calibrate the colour-redshift relation in order to explicitly target the colour manifold spanned by Stage IV surveys (Masters et al. 2015, 2017, 2019; Euclid Collaboration: Guglielmo et al. 2020; Stanford et al. 2021; Euclid Collaboration: Saglia et al. 2022; Gruen et al. 2023; McCullough et al. 2024). However, these observational efforts need to be complemented with an accurate assessment of their incompleteness in the targeted galaxy population, especially at faint magnitudes, where a significant gap remains in the spectroscopic coverage. Similar issues arise when using survey data to study the evolution of the galaxy population. In photometric surveys, redshift uncertainties hinder the accurate reconstruction of galaxy evolution across cosmic time, particularly at high redshift. In spectroscopic surveys, the complicated spectroscopic selections often lead to a difficult assessment of the stellar mass completeness or the contamination of a different population of objects with respect to that selected via photometry, which potentially biases the conclusions one draws from the analysis of these samples.

The forward-modelling of photometric and spectroscopic galaxy surveys is a powerful approach with which to tackle both the problem of accurate galaxy redshift distribution estimates and the characterisation of the galaxy population that underlies the observations of a survey. This methodology involves the construction of a detailed, end-to-end model of a galaxy survey. Key components are a realistic galaxy population model (i.e. the intrinsic distribution of redshift-evolving physical properties of galaxies), a model for the galaxy spectral energy distributions (SEDs) that connects the galaxy physical properties to their emitted light, a mapping of the intrinsic properties and fluxes onto the observed ones by means of image and spectra simulators, and informative observed data against which to compare and constrain the model. By forward-modelling all relevant observational and instrumental effects and calibrating model parameters using simulation-based inference (SBI, Cranmer et al. 2020) methods such as approximate Bayesian computation (ABC, Akeret et al. 2015), the galaxy redshift distributions constructed from observed data can be replaced by those that have been accurately determined via simulations in a cosmological analysis. Moreover, forward-modelling enables the simultaneous inference of the underlying relations between a galaxy physical properties, which offers a consistent framework for studying galaxy evolution across cosmic time. Compared to other modelling approaches, such as semi-analytic models or hydrodynamical simulations, forward-modelling offers two main advantages: it allows for fast realisations of mock galaxy catalogues, which enables efficient exploration of high-dimensional parameter spaces and of the detailed simulation of the observational and measurement process, and it provides property distributions that match, at all redshifts, those from the observations it aims to reproduce.

The concept of forward-modelling galaxy surveys for cosmological applications was first introduced in Refregier & Amara (2014) and further developed in a series of subsequent studies (Herbel et al. 2017; Tortorelli et al. 2018, 2020; Fagioli et al. 2018, 2020; Tortorelli et al. 2021; Bruderer et al. 2016; Moser et al. 2024; Fischbacher et al. 2025a). These works adopt a common phenomenological galaxy population model, which has been progressively refined by incorporating deeper photometric data across a broader range of wavebands, including narrow-band filters (Tortorelli et al. 2018, 2021), and tested against spectroscopic data from the Sload Digital Sky Survey (SDSS, Fagioli et al. 2018, 2020). A major strength of this approach is its ability to accurately simulate both photometric and spectroscopic datasets by means of fast and realistic image (UFIG, Bergé et al. 2013; Bruderer et al. 2016; Fischbacher et al. 2025b) and spectra simulators (USPEC, Fagioli et al. 2018, 2020; Tortorelli et al., in prep.). These tools make it possible to consistently apply the same measurement pipelines and selection cuts on data and simulations. This ensures a robust and consistent comparison of observables and properly accounts for systematics introduced by the measurement process itself, including, for example, the impact of rejecting blended sources, which are a source of bias expected to significantly affect cosmological constraints from Stage IV surveys (see e.g. Ramel et al. 2024 and references therein). Using this methodology, Kacprzak et al. (2020) performed a cosmic shear analysis with DES Y1 data, while Tortorelli et al. (2020, 2021) measured the B-band galaxy luminosity function from the Canada-France-Hawaii Telescope Legacy Survey (Cuillandre et al. 2012), and inferred galaxy physical properties from the narrow-band SEDs in the Physics of the Accelerating Universe survey (Navarro-Gironés et al. 2024).

The latest iteration of this phenomenological model, named GALSBI, has been made publicly available in Fischbacher et al. (2025a, 2024). This model builds upon the foundational works of Herbel et al. (2017), Tortorelli et al. (2018, 2020, 2021), Kacprzak et al. (2020), and Moser et al. (2024) and incorporates analytical parametrisations of galaxy luminosity functions, morphologies, and SEDs. A key innovation introduced in Fischbacher et al. (2025a) is the development of an emulator that directly generates catalogues with observed-like galaxy properties from the true input ones. This emulator, composed of a neural network classifier and a normalising flow, bypasses the computationally expensive need to simulate survey images using UFIG and significantly accelerates the forward-modelling process. This approach enabled the authors to explore 40 different configurations of the model and analysis choices. Combined with a highly informative dataset, such as the HSC deep fields (Aihara et al. 2018), the authors placed very tight constraints on GALSBI parameters, with this phenomenological model being a realistic representation of the galaxy population visible in a Stage IV precursor dataset like HSC. This is testified to by the strong agreement between in the magnitude, colour, and size distributions down to i ≤ 25 and by the agreement of the simulated tomographic and non-tomographic redshift distributions with those derived from COSMOS photo-z catalogues (Weaver et al. 2022), which are computed using both the EAZY (Brammer et al. 2008) and LEPHARE (Arnouts et al. 1999; Ilbert et al. 2006) codes, in a Stage-III setup. Some tension is still present in a Stage-IV setup, which is a gap that we aim to fill with a well-constrained physical model of the galaxy population.

Alternative forward-modelling approaches to the one adopted by GALSBI have recently emerged in the literature and offer complementary strategies that differ in several key aspects. Some of these aspects are the treatment of the noise model for the data, the modelling of galaxy SEDs, and the way galaxy properties are sampled from the underlying population model. For example, Sudek et al. (2022) used the SKYPY framework (Amara et al. 2021) to forward-model photometric data and investigate the sensitivity of redshift distributions to assumptions about galaxy demographics. Their approach is based on the works of Herbel et al. (2017), Fagioli et al. (2018, 2020), and Tortorelli et al. (2018, 2020, 2021), and relies on empirically derived SED templates, which are typically governed by a small set of parameters. A more physically motivated alternative is the use of evolutionary stellar population synthesis (hereafter, SPS, Tinsley 1980; Pickles 1985; Bruzual & Charlot 2003; Maraston 2005; Conroy et al. 2009; Vazdekis et al. 2016; Conroy 2013; Iyer et al. 2025), a technique that exploits the knowledge of stellar evolution to model the starlight component of galaxy SEDs by combining the emission due to stars using a stellar initial mass function (IMF, Salpeter 1955; Chabrier 2003; Kroupa 2002) and a star-formation history (SFH) according to given isochrones and atmosphere libraries. SPS, when combined with models for the chemical evolution, dust attenuation and emission, and gas and active galactic nuclei (AGN) emission, enables the physical properties of galaxies to be directly connected to their emitted light by accounting for all these major SED modelling components. Though computationally more expensive, SPS offer a richer and more flexible description of galaxy spectra. Well-established SPS-based codes, such as FSPS (Conroy et al. 2009; Conroy & Gunn 2010), have been integrated into galaxy population models like PROSPECTOR-α (Leja et al. 2017) and PROSPECTOR-β (Wang et al. 2023) to generate galaxy SEDs from the physical properties drawn from their model priors. PROSPECTOR-β, in particular, was used in Tortorelli et al. (2024) to forward model the nine band u, g, r, i, Z, Y, J, H, Ks photometry from the combined Kilo-Degree Survey (KiDS) and VISTA Kilo degree Infrared Galaxy (VIKING) survey (Wright et al. 2019) and evaluate the impact of SED modelling choices on forward-modelling-based redshift distribution estimates. Another notable example is the forward-modelling framework developed by Alsing et al. (2023) and Leistedt et al. (2023), which also uses SPS but accelerates magnitude predictions via an emulator of FSPS, named SPECULATOR (Alsing et al. 2020). This framework replaces full image or spectral simulations with a data-driven noise model, which enables the efficient estimation of redshift distributions for photometric surveys. The initial implementation of their model was based on analytical parametrisations of the galaxy population (Alsing et al. 2023). In follow-up works (Alsing et al. 2024; Thorp et al. 2024, 2025), the authors replaced this analytical model with a diffusion model (POP-COSMOS) that was introduced to learn the joint distribution of galaxy physical properties by matching observed magnitudes in 26 bands from the COSMOS 2020 catalogue (Weaver et al. 2022). The pre-trained POP-COSMOS model can also be used as an efficient Bayesian method for estimating individual photo-zs and galaxy properties, since it leverages a GPU-accelerated affine-invariant ensemble sampler to achieve fast posterior sampling (Thorp et al. 2024). Hahn et al. (2023a, 2024) introduced the Bayesian SED modelling framework PROVABGS to forward-model the DESI Bright Galaxy Sample (BGS, Hahn et al. 2023b), which enabled both individual- and population-level SPS parameters to be inferred using the POPSED framework (Li et al. 2024). In a parallel development, Alarcon et al. (2023) proposed a parametric model for galaxy assembly histories, which, in conjunction with the differentiable SPS code DSPS (Hearin et al. 2023), was used to construct the extra-galactic forward model DIFFSKY and to populate the OpenUniverse2024 simulation suite (OpenUniverse et al. 2025) with synthetic galaxies.

In this paper, we introduce a new SPS-based galaxy population model called GALSBI-SPS. Developing a realistic galaxy population model based on physical properties and SPS is crucial for the success of forward-modelling Stage IV cosmological experiments. Such a model allows information to be leveraged from a more constrained dataset (low-redshift samples, deep fields, or spectroscopic datasets pre-selected through photometry) and enables that information to be transferred to the fainter, magnitude-limited samples targeted by upcoming surveys. In a purely data-driven framework, such extrapolation beyond the training set is inherently uncertain. In a phenomenological model, this extrapolation is somewhat more viable, provided that the SED templates calibrated at low redshift remain valid at higher redshifts. In contrast, a physically motivated model like GALSBI-SPS enables a more principled and flexible extrapolation. By sampling galaxy physical properties from redshift-evolving galaxy stellar mass functions (GSMFs) and galaxy SFH priors, one can construct higher-redshift analogues of galaxies observed at low redshift, i.e. fainter and younger versions that retain physical consistency. This capability is a direct consequence of using SPS coupled with models for the gas, dust, and AGN emission and absorption, which allows for the generation of realistic galaxy SEDs across cosmic time based on galaxy physical properties. As a result, SPS enables the detailed simulation of the faint and high redshift galaxy population that Stage IV surveys are going to target. Moreover, the integration of SPS-based forward-modelling with image and spectra simulators such as UFIG and USPEC places us in a unique position of being able to characterise the full selection function of a survey as a function of galaxy physical properties, instrumental configuration, target selection criteria, and observing conditions.

The prescriptions implemented in GALSBI-SPS are motivated by the findings in Tortorelli et al. (2024), where we conducted a sensitivity analysis on how different SED modelling choices affect the mean and scatter of tomographic galaxy redshift distributions in SPS-based forward-modelling. In Tortorelli et al. (2024), the joint distribution of physical galaxy properties were sampled from the PROSPECTOR-β (Wang et al. 2023) model, which uses the FSPS code to compute intrinsic galaxy spectra and magnitudes. In contrast, the present work introduces a completely new set of analytical parametrisations from which we sampled galaxy redshifts and physical properties, such as galaxy stellar masses, SFHs, gas-phase metallicities, sizes, internal dust extinction, and emission from the central AGN. This is motivated by the fact that, in PROSPECTOR-β, galaxy SFHs, gas, and dust properties are sampled from either uniform or data-driven priors that are independent of other galaxy physical properties. The physical properties we sampled from GALSBI-SPS are then passed to generative codes that make use of SPS to generate galaxy SEDs and magnitudes. The generative SED code we use is the second noteworthy change with respect to Tortorelli et al. (2024). We replace FSPS with the generative SED package PROSsc pect (Robotham et al. 2022) to generate galaxy SEDs. This is motivated by PROSsc pect execution speed, high-resolution capabilities, and ability to input user-defined simple stellar populations (SSPs). The resulting intrinsic properties and magnitudes are then provided as input to the image simulator UFIG, which we employ to simulate the entire COSMOS field footprint covered by HSC Deep data. We then process the observed and the simulated images using SOURCE EXTRACTOR (Bertin & Arnouts 1996) and apply the same selection criteria as in Fischbacher et al. (2025a). This pipeline allows us not only to test the GALSBI-SPS model against observed galaxies, detected by running SOURCE EXTRACTOR on real HSC images, but also to carry out a robust cross-validation with GALSBI. Both the SPS and the phenomenological version of the model are included in the same publicly released PYTHON package GALSBI¹ (Fischbacher et al. 2024).

In this work, we adopt parameter values for our analytical prescriptions that are either modelled on data from the Galaxy and Mass Assembly survey (GAMA, Driver et al. 2011; Liske et al. 2015) and Deep Extragalactic VIsible Legacy survey (DEVILS, Davies et al. 2018) or drawn from existing literature values. We do not perform SBI to constrain the model parameters against observational data, as the primary goal of this study is to introduce and motivate the set of prescriptions used in GALSBI-SPS to model the galaxy population. SBI is able to leverage the enormous constraining power that results from jointly analysing photometric and spectroscopic data. Our forward-modelling framework uniquely enables such a joint analysis thanks to the physically motivated galaxy population model and the ability to simulate both images and spectra from a single parent galaxy catalogue. In principle, each dataset is able to provide constraining power on specific aspects of the galaxy population model; for example, photometric data are highly sensitive to the GSMF, while spectroscopic data constrain aspects of the gas and stellar population physics. Moreover, the use of emulators, such as the ones introduced in Tortorelli et al. (2025) and Fischbacher et al. (2025a), enables the rapid testing of various analytical prescriptions for galaxy properties and different analysis choices, which circumvents the computational cost of simulating images and spectra across large survey footprints.

The paper is structured as follows. In Sect. 2, we describe the datasets used to model the analytical prescriptions and whose catalogues are later used as a comparison for our simulated properties. Section 3 presents the galaxy population model and the SPS-based SED modelling code employed to generate galaxy SEDs. In Sect. 4, we detail the forward-modelling of the HSC survey and the generation of realistic simulated images covering the COSMOS field. Section 5 shows the comparison between observed and simulated photometric and physical galaxy properties. In Sect. 6, we discuss the results of the comparison and the strategies with which one can constrain the physical model. Section 7 summarises the main findings of this study. Throughout this work, we use H₀ = 67.8 km s⁻¹ Mpc⁻¹ in a flat ΛCDM cosmology with Ω_m = 0.308, which corresponds to the Planck 15 cosmology (Planck Collaboration XIII 2016) used in Bellstedt et al. (2020, 2021, 2024).

2. Data

In this section, we present the observational datasets used in this work. The data are used in two steps. We first model the analytical prescriptions between galaxy physical properties using catalogues with spectroscopic redshifts. Then, we test whether the modelled relations allow us to reproduce magnitudes, colours, and sizes from deeper photometric observed data. The data with spectroscopic redshifts are used to model SFHs, dust, AGN and morphological parameters, while the other physical properties are drawn from analytical prescriptions taken from the literature.

2.1. Galaxy And Mass Assembly Survey

GAMA (Driver et al. 2011; Liske et al. 2015) is a spectroscopic survey conducted at the 3.9m Anglo-Australian Telescope using the AAOmega fibre-fed spectrograph. Completed in 2015, GAMA represents one of the largest spectroscopic surveys of the low redshift Universe (z < 0.25), comprising over ∼230 000 spectra with reliable redshifts. GAMA covers a total area of ∼286 deg² down to an r-band Petrosian magnitude of r < 19.8 mag, distributed across three equatorial fields (G09, G12, G15) and two southern fields (G02, G23). GAMA is complemented by multi-wavelength imaging from the UV to the far-IR, making it a legacy dataset that has enabled a wide range of galaxy evolution studies. These include measurements of merger rates (Robotham et al. 2014; Davies et al. 2015), group catalogues (Robotham et al. 2011; Taylor et al. 2020), low-z galaxy population properties (Driver et al. 2012; Loveday et al. 2012; Bellstedt et al. 2020, 2021; Thorne et al. 2022a; Bellstedt et al. 2024), and constraints on the low-z GSMF down to ${10}^{6.75}{M}_{\odot }{\mathrm{h}}_{70}^{-2}$ (Baldry et al. 2012; Wright et al. 2018; Thorne et al. 2021; Driver et al. 2022; Sbaffoni et al. 2025).

In this work, we use the PROSsc pect value added catalogue² used in Bellstedt et al. (2020, 2021, 2024), which is part of the GAMA fourth data release (DR4, Driver et al. 2022). This catalogue provides stellar masses, star-formation rates (SFRs), gas-phase metallicities, and dust parameters derived with PROSsc pect in SED-fitting mode for 233 833 galaxies in the GAMA equatorial and G23 fields. The sample spans the redshift range 0.002 < z < 4.7, with median redshift Me(z) = 0.225 and 95% of galaxies below z < 0.5. The stellar population properties were obtained with PROSsc pect using SSPs produced with the Bruzual & Charlot (2003) evolutionary SPS code, the Chabrier (2003) IMF, the skewed Normal SFH parametrisation (Robotham et al. 2022), a linearly evolving gas-phase metallicity (Robotham et al. 2022), the Charlot & Fall (2000) dust attenuation model and dust emission modelled following Dale et al. (2014). To ensure our calibration sample is representative of the underlying population at fixed stellar mass, we apply the spectroscopic completeness limit from Robotham et al. (2014), which avoids biases towards bright, star-forming systems. Details of this selection are provided in Appendix A. The resulting stellar-mass-complete sample used in this work contains 101 784 galaxies.

In Sect. 3.10, we additionally use the morphological measurements from the Kelvin et al. (2012) value added catalogue³, which provides single-component Sérsic profile fits for 221 373 galaxies in the GAMA II equatorial fields. The fit was performed with SIGMA, an R wrapper that interfaces SOURCE EXTRACTOR, PSF EXTRACTOR (Bertin 2011), and GALFIT (Peng et al. 2010a). In this work, we use the Sérsic fits from the SDSS r-band imaging.

2.2. Deep extragalactic VIsible legacy survey

DEVILS (Davies et al. 2018) is a spectroscopic survey conducted with the same instrument as GAMA, the AAOmega fibre-fed spectrograph on the 3.9m Anglo-Australian Telescope, but designed to target intermediate redshift galaxies (0.3 < z < 1.0) with high spectroscopic completeness (> 85%) down to Y < 21.2 mag. The DEVILS spectroscopic sample comprises ∼60 000 objects over a combined area of ∼6 deg² across three deep extra-galactic fields: COSMOS, ECDFS, and XMM-LSS. DEVILS can be considered the higher redshift counterpart to GAMA, sharing key science goals such as characterising group and pair environments and investigating how these environments influence galaxy evolution over the last 8 Gyr.

This work makes use of the galaxy physical properties measured with PROSsc pect using the spectroscopic and photometric data from the D10-COSMOS field of DEVILS (Thorne et al. 2021, 2022b). This catalogue has been used to measure the evolution of the GSMF and the stellar mass-SFR relation out to z ∼ 4.25 (Thorne et al. 2021), and to constrain the evolution of the bolometric AGN luminosity function (Thorne et al. 2022b). The catalogue includes stellar masses, SFRs, gas-phase metallicities, dust parameters and AGN luminosities for 494 084 sources over the redshift range 0.0002 < z < 9, with a median redshift of Me(z) = 1.24 and 95% of the sample below z < 4. Roughly ∼24 k out of 494 084 sources have spectroscopic redshifts from DEVILS, while for the remaining ones the redshifts have been estimated via photometry. We apply the stellar mass completeness criterion defined in Thorne et al. (2022b), yielding a final sample of 38 066 galaxies. The selection procedure is detailed in Appendix B.

2.3. Hyper-Suprime Cam Subaru Strategic Program

The HSC Subaru Strategic Program (Aihara et al. 2018) is a deep multi-band imaging survey conducted with the 8.2-m Subaru telescope using the HSC camera, covering the g, r, i, z, y filters with exquisite image quality and depth. The survey is structured into three sub-surveys of increasing photometric depth. The Wide survey, the shallowest one (r ∼ 26 mag), covers a large equatorial area of 1470 deg², while the Deep (r ∼ 27 mag) and the Ultra Deep (r ∼ 28 mag) surveys cover a smaller area of just ∼36 deg² in four different equatorial and northern pointings. HSC was designed to address a wide range of science goals, including weak-lensing tests of cosmology, galaxy formation and evolution over cosmic time, and the detection of high redshift galaxies. Its depth, photometric precision, and high galaxy surface density make it an ideal precursor dataset for upcoming Stage IV surveys such as Rubin-LSST (Ivezić et al. 2019).

In this work, we use the publicly available data from the HSC third data release (Aihara et al. 2022) in the COSMOS field. This data are the same ones used in Fischbacher et al. (2025a) to constrain the phenomenological version of GALSBI, leveraging HSC constraining power for the higher redshift population. The HSC Deep and Ultra-deep (DUD) data overlaps with the COSMOS field and therefore with the COSMOS2020 (Weaver et al. 2022) photo-z catalogue in what is labelled as the COMBINED footprint. The latter refers to the part of the COSMOS field where there is full coverage of optical/near-infrared data. HSC DUD in the COSMOS area contains 63 patches, out of which 56 are almost fully covered. We use these 56 patches as they are the same ones used in Moser et al. (2024) and Fischbacher et al. (2025a). Each patch corresponds to a 12 × 12 arcmin² region, equivalent to a 4200 × 4200 pixel image given the HSC pixel scale of 0.168″ per pixel.

2.4. COSMOS2020

COSMOS2020 (Weaver et al. 2022) is the latest release of the COSMOS catalogue at the time of writing and comprises over 1.7 million sources with multi-band photometry from far-UV to mid-infrared and stellar mass completeness down to 10⁹ $M_{\odot }{\mathrm{h}}_{70}^{-2}$ at z ∼ 3. The COSMOS2020 catalogue represents the current state-of-the-art in photo-z catalogues thanks to its depth, wavelength coverage, and photometric accuracy. There are a total of four publicly available COSMOS2020 catalogues that differ in the code used to extract the photometry. Following Moser et al. (2024) and Fischbacher et al. (2025a), we use the CLASSIC catalogue, where the optical/near-infrared photometry is measured with SOURCE EXTRACTOR, while the photo-zs are estimated via the template fitting codes LEPHARE (Arnouts et al. 1999; Ilbert et al. 2006) and EAZY (Brammer et al. 2008). We apply several quality cuts to ensure reliable photometric and redshift estimates. We select sources lying in areas where the photometry is reliable and there is full coverage of optical/near-infrared data by means of the FLAG_COMBINED parameter. Furthermore, we select objects for which MAG_AUTO< 99 in all band (ignoring in this study possible UV dropouts), LEPHARE or EAZY photo-zs are between 0 and 8, LEPHARE object type flag set to galaxy (LP_TYPE = 0) and SOURCE EXTRACTOR FLAGS< 4.

3. Galaxy population model

In this section, we present the SPS-based galaxy population model GALSBI-SPS. The prescriptions implemented in this new model are motivated by the insights from Tortorelli et al. (2024). In that work we identified the SED modelling components that most significantly impact galaxy magnitudes and, consequently, redshift distribution estimates of colour-magnitude selected samples for cosmological applications. Unlike that previous work, we replace FSPS (Conroy et al. 2009; Conroy & Gunn 2010) with the generative galaxy SED package PROSsc pect⁴ (Robotham et al. 2022). PROSsc pect simulates far-ultraviolet-to-radio galaxy SEDs by modelling the emission from stars, gas, dust, and AGN. The SED modelling is self-consistent: the radiation produced by episodes of star formation heats up the dust grains that then cause re-emission of the absorbed energy at longer wavelengths. The attenuation of starlight by the dust is simulated as a two component process. Old stellar populations (> 10 Myr) experience attenuation from a screen-like interstellar medium (ISM), while younger stellar populations (< 10 Myr) are embedded in birth clouds, experiencing additional attenuation. A simple energy balance scheme is then used to produce nebular emission lines, where the starlight emitted by young stars (< 10 Myr) at wavelengths shorter of the Lyman limit (< 911. Å) determines the ionising flux for nebular emission. PROSsc pect is highly flexible, allowing users to input both parametric and non-parametric galaxy SFHs and to select among different stellar population libraries (SPL), like those produced with the Bruzual & Charlot (2003) and EMILES (Vazdekis et al. 2016) evolutionary SPS codes, or custom generated with the newly released PROGENY SPL software package (Robotham & Bellstedt 2025; Bellstedt & Robotham 2025)⁵. A notable feature of PROSsc pect is the ability, when used in conjunction with PROGENY, to generate single and composite stellar populations (CSPs) that make use of stellar age- or metallicity-dependent IMFs. PROGENY also enables the generation of high-resolution SSPs for PROSsc pect, a crucial feature that will enable the use of PROSsc pect in spectrophotometric mode to jointly analyse photometric and spectroscopic data from the WAVES (Driver et al. 2019) survey conducted with the 4MOST spectrograph (de Jong et al. 2019). The code’s flexibility, ongoing development, high-resolution capabilities, and fast execution speed (∼20–30 ms per galaxy SED on a Mac M1 CPU) strongly motivate its adoption as the generative SED engine for GALSBI-SPS. Importantly, the consistent use of PROSsc pect to analyse observed GAMA and DEVILS data, and to synthesise new galaxy SEDs in generative mode, ensures internal consistency in the physical assumptions underlying our galaxy population model and in the relation of physical to observable properties.

As in the phenomenological version of the model (Fischbacher et al. 2025a), galaxies in GALSBI-SPS are drawn from two overlapping populations of blue and red galaxies that share similar parametrisations, but differ in parameter values. Adding scatter to the physical relations allows the two populations to overlap, reflecting the observed continuum of galaxy properties, both in terms of SEDs and global physical characteristics. Figure 1 presents a flowchart illustrating the main components of the model. Briefly, GALSBI-SPS samples galaxy formed stellar masses and redshifts from two redshift-evolving GSMFs for blue and red galaxies. It construct galaxy SFHs using parametrisations based on GAMA and DEVILS data and samples the remaining physical properties, necessary for realistic SED generation with PROSsc pect, conditioned on their respective relations to galaxy stellar mass, redshift, and SFR. Apparent magnitudes are obtained by redshifting, dimming, and integrating the generated SEDs over the relevant filter bands. Finally, each galaxy is assigned morphological properties, including single Sérsic indices, sizes, ellipticities, and random spatial positions within the simulated field⁶, providing all the information required for realistic image simulations of un-clustered galaxies.

Fig. 1. Flowchart describing the main components of the GALSBI-SPS galaxy population model. Galaxy physical properties are sampled hierarchically from analytical distributions and conditional relations, as detailed in Sect. 3. The box colours reflect the hierarchy of dependencies. Blue boxes represent properties that are independently sampled, with the GSMF being the starting point of the sampling. Yellow boxes represent properties conditioned only on galaxy stellar mass and redshift (namely, the galaxy SFH, velocity dispersion, morphology, and AGN). Green boxes represent properties conditioned simultaneously on galaxy stellar mass, redshift, and SFH (galaxy metallicity history, dust attenuation, and gas ionisation). All sampled properties are then passed to PROSsc pect (red box) to generate galaxy SEDs and magnitudes. Together with the morphological properties they form the input catalogue for the forward-modelling. Each model component is described in its corresponding subsection of Sect. 3.

3.1. Galaxy stellar mass function

GALSBI-SPS samples formed stellar masses of galaxies logℳ (in solar mass units M_⊙) and their redshifts z from redshift-evolving GSMFs, modelled as the sum of two Schechter functions (Schechter 1976). These two components, a low mass (‘l’ subscript) and a high mass (‘h’ subscript) ones, are assigned separate parameter values for blue (‘b’ subscript) and red (‘r’ subscript) galaxy populations. This double-Schechter parametrisation provides an accurate description of the statistical distribution of galaxy stellar masses in the local Universe (e.g. Peng et al. 2010b) and up to z ≲ 3 (Leja et al. 2020). At higher redshifts, the advantage of using a double function over a single Schechter function becomes less apparent (Weaver et al. 2023). However, we use a well motivated redshift evolution of the GSMF parameters that naturally accommodate the transition from a double to a single-Schechter form at high redshifts. The GSMF defines the number density of galaxies per comoving volume V_C per logarithmic stellar mass bin at redshift z, and is expressed in terms of logarithmic stellar mass logℳ following the formalism of Weigel et al. (2016) and Weaver et al. (2023):

$$\begin{array}{cc}\hfill \mathrm{\Phi }(log\mathcal{M},z)=& \frac{\mathrm{d}N}{\mathrm{d}(log\mathcal{M})\mathrm{d}{V}_{\mathrm{C}}}=\hfill \\ \hfill & \ln (10)e^{-{10}^{log\mathcal{M}-log{\mathcal{M}}^{\ast }(z)}}\hfill \\ \hfill & \times \left[{\varphi }_l^{\ast }(z){\left({10}^{log\mathcal{M}-log{\mathcal{M}}^{\ast }(z)}\right)}^{{\alpha }_l+1}\right]\hfill \\ \hfill & \times \left[{\varphi }_h^{\ast }(z){\left({10}^{log\mathcal{M}-log{\mathcal{M}}^{\ast }(z)}\right)}^{{\alpha }_h+1}\right].\hfill \end{array}$$(1)

The characteristic stellar mass logℳ^* is assumed to be the same for both Schechter components, as it is standard practice when fitting double Schechter functions (Baldry et al. 2012; Muzzin et al. 2013; Tomczak et al. 2014). ϕ_l^*, ϕ_h^* represents the low mass and high mass number density of galaxies at the characteristic stellar mass logℳ^*, while α_l, α_h represent the faint-end slopes of the two functions. We choose to allow for a redshift evolution of logℳ^*, ϕ_l^*, ϕ_h^*, while we keep α_l, α_h fixed to the mean values reported for blue and red galaxy populations in Tables C.2 and C.3 of Weaver et al. (2023).

We parametrise the redshift evolution of the characteristic stellar mass as

$$\begin{array}{cc}\hfill log{\mathcal{M}}^{\ast }(z)=& log{{\mathcal{M}}^{\ast }}_0+log{{\mathcal{M}}^{\ast }}_1\times log(1+z)+\hfill \\ \hfill & log{{\mathcal{M}}^{\ast }}_2\times log{(1+z)}^2.\hfill \end{array}$$(2)

The logℳ^*₀, logℳ^*₁, logℳ^*₂ values for both blue and red galaxy populations are obtained by fitting the redshift evolution of logℳ^* reported in Tables C.2 (blue galaxies) and C.3 (red galaxies) of Weaver et al. (2023). This functional form is motivated by the fact that it provides a good fit to the measurements in Weaver et al. (2023). We limit our fit to z ≤ 3 galaxies, where a second Schechter component is necessary for describing the GSMF. Beyond this redshift, we extrapolate the best-fitting relations. Although extrapolation at z > 3 is sub-optimal, the vast majority of galaxies for a cosmological analysis in Stage IV surveys will lie below this redshift. Figure 2 shows the evolution of logℳ^* with redshift, based on the best-fitting parameters listed in the first three rows of Table 1, separately for blue and red galaxies.

Fig. 2. Redshift evolution of the characteristic stellar mass logℳ* for blue and red galaxies. Blue and red points represent the measurements from Weaver et al. (2023). The blue and red dashed lines are built by fitting the measurements with Equation 2. The blue and red bands represent the 1σ uncertainty band on the fit.

The redshift evolution of the normalisation parameters ϕ_l^*, ϕ_h^* follows the same prescription adopted for luminosity functions in Moser et al. (2024),

$$\begin{array}{c}\hfill {\varphi }_{l,h}^{\ast }(z)={\varphi }_{\mathrm{amp},l,h}^{\ast }\times {(1+z)}^{{\varphi }_{\exp ,l,h}^{\ast }},\end{array}$$(3)

where the ϕ_amp^*, ϕ_exp^* differs for blue and red galaxies, and for the low-mass and high-mass components of the GSMF. The parameter values used in this work are obtained by fitting Equation (3) to the values in Tables C.2 (blue galaxies) and C.3 (red galaxies) of Weaver et al. (2023). Although not all fitted parameters are constrained over the full redshift range up to z = 3 in Weaver et al. (2023), we extrapolate the best-fitting relations when necessary. Similar considerations to the case of the high-redshift extension of logℳ^* apply. Figure 3 shows the redshift evolution of ϕ_l^* and ϕ_h^* obtained with the best-fitting parameters listed in the fourth to seventh rows of Table 1, again separately for blue and red populations. The parametric form provides a good fit to the observed evolution for low- and high-mass blue galaxies and for low-mass red galaxies. A tension is observed for high-mass red galaxies at z ≳ 2, where the model under predicts the decline in number density. This could be due to either an intrinsic evolution that would require a more flexible parametrisation or the limited statistics for that population of galaxies in the COSMOS field.

Fig. 3. Redshift evolution of the characteristic density ϕ* for blue and red galaxies. Blue, deep sky blue, magenta, and red points represent the measurements from Weaver et al. (2023). The blue, deep sky blue, magenta, and red dashed lines are built by fitting the measurements with Equation 3. The bands represent the 1σ uncertainty band on the fit.

GALSBI-SPS computes the redshift distribution of galaxies, n(z), by integrating the redshift-stellar mass density φ over a specified stellar mass range. The redshift-stellar mass density, φ, is defined as

$$\begin{array}{c}\hfill \phi =\frac{\mathrm{d}N}{\mathrm{d}log\mathcal{M}\mathrm{d}z\mathrm{d}\mathrm{\Omega }},\end{array}$$(4)

where dΩ is the solid angle. Substituting the expression for the comoving volume element (Hogg 1999),

$$\begin{array}{c}\hfill \mathrm{d}{V}_{\mathrm{C}}=\frac{{\mathrm{D}}_{\mathrm{H}}{(1+z)}^2{{\mathrm{D}}_{\mathrm{A}}}^2}{\mathrm{E}(z)}\mathrm{d}z\mathrm{d}\mathrm{\Omega }\end{array}$$(5)

where D_H, D_A are the Hubble distance and the angular diameter distance, respectively, and $\mathrm{E}(z)=\sqrt{{\mathrm{\Omega }}_{\mathrm{M}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{k}}{(1+z)}^2+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}$, and if we use the definition of the GSMF Φ(logℳ, z), the redshift distribution for each Schechter component becomes

$$\begin{array}{cc}\hfill n(z)& ={\displaystyle {\int }_{log{\mathcal{M}}_{\min }}^{log{\mathcal{M}}_{\max }}}\phi (log\mathcal{M},z)\mathrm{d}log\mathcal{M}=\hfill \\ \hfill & {\displaystyle {\int }_{log{\mathcal{M}}_{\min }}^{log{\mathcal{M}}_{\max }}}\frac{{\mathrm{D}}_{\mathrm{H}}{(1+z)}^2{\mathrm{D}}_{\mathrm{A}}^2}{\mathrm{E}(z)}\mathrm{\Phi }(log\mathcal{M},z)\mathrm{d}log\mathcal{M}=\hfill \\ \hfill & \frac{{\mathrm{D}}_{\mathrm{H}}{(1+z)}^2{\mathrm{D}}_{\mathrm{A}}^2}{\mathrm{E}(z)}{\varphi }^{\ast }\mathrm{\Gamma }(\alpha +1,{10}^{log{\mathcal{M}}_{\min }-log{\mathcal{M}}^{\ast }(z)}),\hfill \end{array}$$(6)

where Γ denotes the incomplete gamma function. Stellar mass sampling proceeds in two steps: first we draw redshifts from the n(z), then we sampled stellar masses from the conditional GSMF at the selected redshift. The total number of galaxies N expected within a sky area A (in steradian) over a redshift range [z_min, z_max] is given by

$$\begin{array}{c}\hfill N={\displaystyle {\int }_{z_{\min }}^{z_{\max }}}n(z)A\frac{{\mathrm{D}}_{\mathrm{H}}{(1+z)}^2{\mathrm{D}}_{\mathrm{A}}^2}{\mathrm{E}(z)}\mathrm{d}z.\end{array}$$(7)

The total number of parameters for the low-mass and high-mass blue and red GSMFs is 18. The GSMF parameter values used in this work are listed in Table 1. It is important to note that the logarithmic stellar mass logℳ sampled from the GSMF in GALSBI-SPS represents the total stellar mass formed by the galaxy and all its progenitors throughout their lifetime. This will be used to obtain the correct galaxy SFH to pass to PROSsc pect. Care should therefore be taken when directly comparing this GSMF with those quoted in the literature, as the GSMFs derived from observations commonly refer to the galaxy surviving stellar masses. We did not sample the latter quantity as it depends on the specific SPL being used and is generally produced as an output quantity in an SPS-based code. PROSsc pect allows for the computation of the surviving stellar mass; however, since this will increase the runtime of the galaxy population sampling, we develop a machine-learning-based emulator that predicts the galaxy surviving stellar mass from its formed stellar mass, SFH parameters and gas-phase metallicity with per-mille level accuracy. We describe the emulator in Appendix C.

3.2. Stellar population model

Evolutionary SPS (Tinsley 1980; Pickles 1985; Bruzual & Charlot 2003; Maraston 2005; Conroy et al. 2009; Vazdekis et al. 2016; Conroy 2013) models the spectral emission from a group of unresolved stars, like those in a galaxy, by making use of the knowledge of stellar evolution. The simplest flavour of an evolutionary SPS model is an SSP, where it is assumed that all stars are coeval and share the same chemical composition. Following the nomenclature of Conroy (2013), when combining SSPs wth a SFH, chemical evolution and a model for dust attenuation, we can build CSPs that are complex representations of the integrated starlight emitted from galaxies.

The basic ingredients to build an SSP are the stellar evolutionary tracks (isochrones), stellar model atmospheres, and stellar IMF. Isochrones trace the evolution of stars of given mass, chemical composition and α/Fe abundance, through the various evolutionary phases that characterise the life of a star. The key predictions that isochrones must convey are basic stellar parameters, such as bolometric luminosities, surface gravities, effective temperatures (and potentially the effective physical size), as functions of evolutionary timescales. Stellar model atmospheres, which can be either empirical or theoretical, are instead the means by which we translate these basic stellar parameters into the emergent flux of stars. Commonly the predictions of stellar properties generated by isochrones have a much finer resolution than the coarser grids of available stellar spectra. Therefore, some degree of interpolation is required to generate flux predictions (see discussion in Robotham & Bellstedt 2025). Predictions from isochrones and stellar atmospheres are then complemented by the choice of a stellar IMF, that is a weighting term that encodes the probability to form a star of given mass, thereby returning the distribution of stellar masses that were generated by a single burst of star formation.

Different SPS codes exist in the literature with which one can generate SSPs, some notable ones being the evolutionary SPS codes of Bruzual & Charlot (2003) and Maraston (2005), EMILES (Vazdekis et al. 2016), and FSPS (Conroy et al. 2009; Conroy & Gunn 2010). PROSsc pect, on the contrary, does not produce SSPs, but rather combines them with complex models of SFHs (and other SED emission components) to generate spectral outputs for simulations and fit observed galaxy SEDs. Similar functionality is implemented by other SED fitting codes, such as MAGPHYS (da Cunha et al. 2008) and CIGALE (Boquien et al. 2019), each with different emphasis on specific SED modelling choices. The original implementation of PROSsc pect allows the user to choose between the SSPs generated with the evolutionary SPS codes of Bruzual & Charlot (2003) and (Vazdekis et al. 2016) (EMILES). However, Robotham & Bellstedt (2025) recently released PROGENY, an SPL software package that allows for the creation of SSPs and their associated spectra. PROGENY directly interfaces with PROSsc pect by creating custom, user-defined SSPs that are used to create CSP to model the starlight emission from galaxies. The notable addition, in comparison with other SPS-based codes, is the possibility to produce SSPs with evolving IMFs, in terms of stellar age and/or metallicity, and the ability to output spectra with higher wavelength resolution than traditional SSPs, such as those of Bruzual & Charlot (2003).

In Robotham & Bellstedt (2025) and Bellstedt & Robotham (2025), the authors compared the PROGENY-generated SSPs against those from other evolutionary SPS codes, as well as the impact of the isochrones, stellar atmospheres, and IMF choices on the recovery of galaxy properties using simulations and GAMA observations. They found that the predicted spectra are more strongly impacted by the choice of the isochrones rather than that of the stellar atmospheres or the IMF, which is also the major factor influencing the differences in the recovered SFHs for GAMA observed galaxies. They also tested the impact of Chabrier (2003) and Kroupa (2002) IMFs, either keeping them fixed or evolving them as function of metallicity (Martín-Navarro et al. 2015), finding a very small impact on recovered SFHs and a somewhat larger impact on galaxy stellar mass. This is similar to what has been found in Tortorelli et al. (2024), where changing the IMF was leading to differences in the i-band fluxes that can be interpreted as differences in galaxy stellar mass when using the i-band flux as its proxy.

In our work, we use the PROGENY web tool to create our custom SSP by incorporating the fiducial choices of isochrones, stellar spectral atmospheres, and IMFs highlighted in Robotham & Bellstedt (2025) and Bellstedt & Robotham (2025). We select the MIST isochrones (Dotter 2016) as they cover all phases of stellar evolution, including remnant stars, important for the UV upturn of older stellar populations. They have high-resolution stellar evolution age cadences of 0.05 dex between 10⁵ and 10^10.3 yr, with a 0.5 dex metallicity gridding between 10⁻⁴ and 10^0.5 Z_⊙ and a mass range of 0.1 to 3 × 10² M_⊙. The stellar atmosphere consists of a combination of ‘base’ and ‘extended’ templates to ensure a good all-around coverage in temperature, surface gravity and metallicity space. The ‘base’ templates come from the atmospheric spectra included with FSPS (Conroy et al. 2018) (sometimes labelled as C3K). They cover temperatures from 2 × 10³ to 5 × 10⁴ K, surface gravities from 10⁻¹ to 10^5.5 cm s⁻² on a logarithmic grid and logarithmic metallicities in solar units from −2.1 to 0.5. The ‘extended’ template spectra come from the Allard et al. (2012) version of PHOENIX spectra Hauschildt et al. (1999) and are meant to extend the atmospheres coverage to very low mass stars and up to temperatures of 7 × 10⁴ K. To extend the coverage of the atmospheres to the asymptotic giant branch (AGB), white dwarf and Wolf Rayet phases of stellar evolution, we further select the POWR Wolf-Rayet library (Todt et al. 2015), the Lançon & Mouhcine (2002) AGB library (also used in Maraston 2005 and Conroy & Gunn 2010) and the Werner et al. (2003) white dwarf library. We then choose the Chabrier (2003) IMF as the last critical component for our custom SSP. The latter covers a wavelength range from 5 to 3.6 × 10⁸ Å, re-binned to match the exact wavelength sampling of the 2019 version (CB19) of Bruzual & Charlot (2003) and Plat et al. (2019). More details on the invididual components are provided in Robotham & Bellstedt (2025).

3.3. Galaxy star formation history

The galaxy SFH is a key ingredient for building a physically motivated model for the galaxy population that is based on SPS, as it represents, in the most general form, the number of stars produced in a galaxy per unit time from the conversion of gas into stars. Coupled with the IMF, the galaxy SFH represents the amount of stellar mass produced per unit time. In PROSsc pect, using the Chabrier IMF and the SSP produced with PROGENY, the SFH represents the weights of the weighted sum of SSPs with which CSPs are produced. When designing a SFH parametrisation, we need to find the right trade-off between flexibility, needed to capture the diversity of galaxy SFHs, and simplicity, to avoid difficulty in constraining with observations and introducing unnecessary parameter degeneracies. The literature offers a wide range of SFH models. The simplest is an instantaneous burst model, where all the stellar mass is produced in a single burst of star formation and the galaxy stellar population coincides with an SSP, where only the evolution of stars with a single metallicity is required. The most commonly used time-evolving parametrisation is instead the τ-model, which assumes an exponentially declining SFR. However, numerous studies have shown that this form produces biased estimates of galaxy properties when applied to mock galaxies (Simha et al. 2014; Pacifici et al. 2015; Ciesla et al. 2017; Carnall et al. 2018). To overcome these limitations, researchers have developed more flexible parametrisations. These include the delayed-τ model (Simha et al. 2014), the double power-law SFH (Alsing et al. 2023), the step-wise SFHs (Leja et al. 2019; Wang et al. 2023), and the skewed Normal SFH (Robotham et al. 2022). Non-parametric SFH models are also widely used, for example basis SFHs (Hahn et al. 2023a) constructed by applying non-negative matrix factorisation to Illustris SFHs (Genel et al. 2014) or SFHs from semi-analytic models (Lagos et al. 2019).

GALSBI-SPS models individual galaxy SFHs in units of M_⊙/yr using a truncated skewed normal distribution (Robotham et al. 2022). This form ensures that SFHs are anchored at zero star formation at the beginning of the Universe, naturally rising during the first billion years with a smooth, physically motivated truncation at early times. Unlike other parametric forms that assume an unphysical finite SFR at formation, this model forces the SFR to start from zero, aligning with expectations from galaxy formation theory. This approach has proven highly successful in fitting observed galaxy populations (Bellstedt et al. 2020, 2021; Thorne et al. 2021, 2022a,b; Bellstedt et al. 2024). The functional form of the SFH reads

$$\begin{array}{c}\hfill \mathrm{SFR}(t)=\mathrm{SFR}{(t)}_{\mathrm{norm}}\times [1-\frac{1}{2}[1+\mathrm{erf}\left(\frac{t-\mu }{\sigma \sqrt{2}}\right)]],\end{array}$$(8)

where

$$\begin{array}{cc}\hfill \mu & =\mathrm{mpeak}\times \frac{|(\mathrm{magemax}-\mathrm{mpeak})|}{\mathrm{mtrunc}},\hfill \\ \hfill \sigma & =\frac{|(\mathrm{magemax}-\mathrm{mpeak})|}{2\times \mathrm{mtrunc}}.\hfill \end{array}$$(9)

The parameter mtrunc = 2 Gyr sets the sharpness of the early-time truncation (with 0 being no truncation), while magemax = 13.4 Gyr represents the maximum age of star formation, chosen to correspond to the existence of some of the earliest observed galaxies (z ∼ 11, Oesch et al. 2016) within our adopted cosmology. The naming convention for the SFH parameters follows that introduced in Robotham et al. (2022). The second term of Equation (8) represents the truncation factor, while the first term is a skewed normal distribution expressed as

$$\begin{array}{c}\hfill \mathrm{SFR}{(t)}_{\mathrm{norm}}=\mathrm{mSFR}\times e^{\frac{X{(t)}^2}{2}},\end{array}$$(10)

where

$$\begin{array}{c}\hfill X(t)=\left(\frac{t-\mathrm{mpeak}}{\mathrm{mperiod}}\right){\left(e^{\mathrm{mskew}}\right)}^{\mathrm{asinh}\left(\frac{t-\mathrm{mpeak}}{\mathrm{mperiod}}\right)}.\end{array}$$(11)

The skewed normal distribution is governed by four parameters: mSFR represents the peak SFR in units of M_⊙/yr, mpeak is the age in Gyr at which the peak SFR occurs, mperiod is the width of the star formation period in Gyr, and mskew is the skewness controlling the asymmetry of the SFH. The chosen parametrisation achieves the best results when galaxies have experienced a single dominant epoch of star formation, though it might be less accurate if the galaxy has undergone a very recent small burst of star formation, which is however a notoriously hard-to-model SFH component. As shown in Robotham et al. 2022, this parametrisation works at the galaxy population level as it is able to correctly reproduce the diversity of SFHs seen in the SHARK semi-analytic model (Lagos et al. 2019). Furthermore, in Bravo et al. (2022), the authors found that the vast majority of galaxies have a single dominant epoch of star formation, thereby making this parametrisation reliable even in the case of burstiness in the SFH.

To meaningfully sample the four free SFH parameters (mSFR, mpeak, mperiod, and mskew) and generate realistic galaxy SFHs for blue and red galaxies, we modelled the distributions of best-fitting SFH parameter values obtained from GAMA (Bellstedt et al. 2020, 2021, 2024) and DEVILS (Thorne et al. 2021, 2022b,a) data. We restrict the modelling to stellar mass-complete samples and remove objects whose best-fitting SFH parameters lie at the edges of their priors, as defined in Table 2 of Bellstedt et al. (2020). We verify that this selection does not bias the overall parameter distributions of the mass-complete samples. Because we aim at using the SFH parametrisation in a generative mode, we cannot directly model the observed distribution of best-fitting parameters. If we were to do so, we would sample an exact replica of the galaxy population observed by GAMA and DEVILS. In particular, the sampled parameters would lead to a distribution of SFHs whose formed stellar masses do not coincide with the distribution of stellar masses drawn from the GSMF. In our model, we condition the sampling of SFH parameters such that the integrated stellar mass over each individual galaxy SFH equals the formed stellar mass sampled from the GSMF for that galaxy (see Sect. 3.1 for a discussion on the stellar masses sampled from the GSMF).

To ensure that the latter is true, we need to remodel the SFHs from GAMA and DEVILS and separate the sampling process of the SFH parameters between those that define the SFH shape (mpeak, mperiod and mskew) and the SFH normalisation (mSFR). First, we define a common time grid in observer-frame look-back time, spanning the range from t = 10⁵ yr to t = 13.7 × 10⁹ yr, corresponding to the age of the Universe at z = 0 in our adopted cosmology. Since the original SFHs in PROSsc pect are defined in galaxy-frame look-back time, in order to correctly compute the stellar mass formed, we set the SFH to 0 beyond magemax′ = 13.4 − t_lb, z Gyr, i.e. beyond the age of the Universe at the galaxy redshift. We then redefine the age of the peak star formation, mpeak. Instead of using an absolute value in Gyr and sampling it from GAMA and DEVILS data, we define mpeak′ to be the fraction of the galaxy age at which the SFH peaks,

$$\begin{array}{c}\hfill \mathrm{mpeak}\prime =\frac{\mathrm{magemax}\prime -\mathrm{mpeak}}{\mathrm{magemax}\prime }.\end{array}$$(12)

Values of mpeak′> 1 correspond to rising SFHs. Instead, mperiod and mskew are unchanged from their original definition.

Observationally, it has been shown that more massive galaxies assemble the majority of their stellar mass at earlier cosmic epochs and over a shorter period than less massive galaxies that have more prolonged periods of star formation and assemble most of their mass at later times (see De Lucia et al. 2025 and reference therein for an up-to-date review). This leads us to expect a dependence of the SFH shape parameters on galaxy redshift and stellar mass. Furthermore, we expect this trend to have different strength for blue and red galaxy populations. To model this trend, we parametrise log(mpeak′) for blue and red galaxies as function of stellar mass as

$$\begin{array}{cc}\hfill log(\mathrm{mpeak}\prime )=& log{(\mathrm{mpeak}\prime )}_{\mathrm{mass},\mathrm{evo},\mathrm{intcpt}}+\hfill \\ \hfill & log{(\mathrm{mpeak}\prime )}_{\mathrm{mass},\mathrm{evo},\mathrm{slope}}\times log\mathcal{M},\hfill \end{array}$$(13)

where the subscript ‘mass, evo’ stands for mass evolution, and with slope and intercept linearly evolving with redshift as

$$\begin{array}{cc}\hfill log{(\mathrm{mpeak}\prime )}_{\mathrm{mass},\mathrm{evo},\mathrm{intcpt}}=& {\mathrm{a}}_{\mathrm{mass},\mathrm{evo},\mathrm{intcpt}}+\hfill \\ \hfill & {\mathrm{b}}_{\mathrm{mass},\mathrm{evo},\mathrm{intcpt}}\times log(1+z),\hfill \\ \hfill log{(\mathrm{mpeak}\prime )}_{\mathrm{mass},\mathrm{evo},\mathrm{slope}}=& {\mathrm{a}}_{\mathrm{mass},\mathrm{evo},\mathrm{slope}}+\hfill \\ \hfill & {\mathrm{b}}_{\mathrm{mass},\mathrm{evo},\mathrm{slope}}\times log(1+z).\hfill \end{array}$$(14)

The upper panels of Fig. 4 show the evolution of the slope log(mpeak′)_{mass, evo, slope} and intercept log(mpeak′)_{mass, evo, intcpt} as function of redshift for the combination of mass-complete GAMA and DEVILS quiescent and star-forming galaxies, selected using the specific star-formation rate (sSFR) criterion of log(sSFR) ≥ −10.5, for star-forming objects, and log(sSFR) ≤ −10.5, for quiescent objects (Leja et al. 2017). The lower panels of Fig. 4 show instead the dependence of the median of log(mperiod) and mskew in bins of 0.1 dex on stellar mass. To keep the model simple yet realistic, we also model this dependence as a linear relation:

Fig. 4. Dependence of the SFH shape parameters on galaxy redshifts and stellar mass. The fraction of the galaxy age at which the SFH peaks, log(mpeak′), linearly depends on the galaxy stellar mass and its coefficients in turn depend on the galaxy redshift. The upper panels show the dependence of the slope, log(mpeak′)mass, evo, slope, and the intercept, log(mpeak′)mass, evo, intcpt, on redshift. Red error bars refer to quiescent galaxies from the mass-complete GAMA and DEVILS sample, while blue error bars refer to star-forming galaxies. The dashed lines represent the best-fitting linear relations for the redshift evolution of these parameters, while the red and blue bands represent the 1σ uncertainty on the fit. The lower panels instead show the dependence of log(mperiod) and mskew on the galaxy stellar mass. The error bars refer to the median and standard deviation of these quantities computed in bins of 0.1 dex in stellar mass, overlayed on top of the overall distribution. Dashed lines and colour bands represent best-fitting linear relations and 1σ uncertainty on the fit for the populations of red and blue galaxies.

$$\begin{array}{cc}\hfill log(\mathrm{mperiod})=& log{(\mathrm{mperiod})}_{\mathrm{mass},\mathrm{evo},\mathrm{intcpt}}+\hfill \\ \hfill & log{(\mathrm{mperiod})}_{\mathrm{mass},\mathrm{evo},\mathrm{slope}}\times log\mathcal{M},\hfill \\ \hfill \mathrm{mskew}=& {\mathrm{mskew}}_{\mathrm{mass},\mathrm{evo},\mathrm{intcpt}}+\hfill \\ \hfill & {\mathrm{mskew}}_{\mathrm{mass},\mathrm{evo},\mathrm{slope}}\times log\mathcal{M}.\hfill \end{array}$$(15)

Based on the data at our disposal, we do not find a strong trend of the slope and intercept of these linear relations with redshift. Best-fitting parameters for these relations obtained from GAMA and DEVILS data are reported in Table 2. We sampled the SFH shape parameters for each individual galaxy from a Multivariate Truncated Normal distribution. The mean of this multivariate distribution is different for each galaxy as it is given by Equations (13), (14) and (15) conditioned on the galaxy stellar mass and redshift. The covariance is instead fixed and estimated for blue and red galaxies from the residuals between the GAMA and DEVILS original distribution of SFH shape parameters and that obtained from sampling with the best-fitting parameters in Table 2. The sampling is truncated to the observational limits defined in Table 2 of Bellstedt et al. (2020). We then rescale log(mpeak′) back to the absolute value in Gyr.

The sampled parameters define the SFH shape for each galaxy. To obtain the value of the normalisation mSFR, we impose for each galaxy that the integral of the SFH over the galaxy lifetime returns the formed stellar ℳ drawn from the GSMF:

$$\begin{array}{cc}\hfill \mathrm{mSFR}\times {\displaystyle {\int }_0^{\mathrm{magemax}}}& sf{h}_{\mathrm{shape}}(\mathrm{mpeak},\mathrm{mperiod},\mathrm{mskew})\times \hfill \\ \hfill & [1-\frac{1}{2}[1+\mathrm{erf}\left(\frac{t-\mu }{\sigma \sqrt{2}}\right)]]=\mathcal{M},\hfill \end{array}$$(16)

which yields the final correct SFH for each individual galaxy. The four sampled parameters are passed to the MASSFUNC_SNORM_TRUNC function in PROSsc pect when generating galaxy SEDs. Additionally, we computed the average SFR over the last 100 Myr for each galaxy. Redshifts, stellar masses and SFRs represent the most important physical properties for the modelling of the galaxy population. Most of the other physical properties that are described in the following sections are conditioned on these three fundamental quantities.

3.4. Metallicity history

Most studies in the literature, which are aimed at either modelling or fitting the galaxy population, treat either the gas or the stellar metallicity (or both) as fixed quantities throughout the galaxy lifetime, often assuming the two having the same value, or as a variable but constant value (see Curti 2025 for description of chemical evolution models). However, it is well known that the gas, and by reflection the stellar metallicity, evolves with the galaxy lifetime. Every new generation of stars forms from an interstellar medium (ISM) that gets enriched in metals through supernovae and stellar winds. The subsequent generation of stars will thus form from an ISM that contains a larger fraction of metals with respect to the one from which the previous generation of stars formed. This ‘closed-box’ model is complicated, however, by the possible infall of pristine, lower metallicity gas from the inter-galactic medium and by the loss of metal-rich gas into said medium. The metallicity assumption has a great impact on galaxy SEDs and consequently on the measured galaxy properties. Indeed the well-known age-metallicity degeneracy leads higher metallicity younger SSPs to appear similar to lower metallicity but older SSPs (Worthey 1994). Bellstedt & Robotham (2025) also found that an evolving metallicity, despite having little impact on the estimation of stellar mass and SFR, allows for the removal of systematic biases in the recovered cosmic star formation history (CSFH) and produces more realistic uncertainties in the recovered SFHs. Therefore, it is of great importance to introduce a physically motivated prescription for the gas metallicity evolution and how this reflects on the stellar metallicity, since new generation of stars form from the metal-enriched star-forming gas.

PROSsc pect models the star-forming gas metallicity evolution of an individual galaxy by mixing discrete SSPs with varying ages and metallicities to achieve a target stellar population age and metallicity, potentially using all metallicities available in the provided SSPs. This represents an improvement over approaches that fix the metallicity to a fiducial value, since it allows for the metallicity to be free to vary at either discrete or in between library interpolated values⁷. The functional form of the metallicity evolution we adopt is the same one used in Bellstedt et al. (2021, 2024) and Alsing et al. (2023). The basic idea is to linearly map the stellar mass evolution of each galaxy on to the shape of the gas-phase metallicity evolution (Driver et al. 2013). The metal enrichment is thus following a 1-to-1 relation with the mass-build, such that for example when half of the stellar mass of a galaxy has been assembled, half of the metal enrichment has occurred. In this prescription, the yield, defined as the fraction in mass of a given element produced per unit mass of stars formed and returned to the ISM, is not fixed, but rather evolving, declining with time. This prescription represents what we would expect in a closed-box mode of star formation, but with some pristine inflow and a constant ejecta metallicity regardless of the metallicity of the gas that formed the stars (Robotham et al. 2022).

This prescription has proven to reproduce the CSFH (Bellstedt et al. 2021) and provide robust measurements of galaxy properties when applied to GAMA and DEVILS data (Bellstedt et al. 2021, 2024; Thorne et al. 2021, 2022b). As pointed out in Robotham et al. (2022), this type of metal evolution is appropriate to describe realistic galaxy formation as it looks physically more like the galaxy metallicity histories in semi-analytic models (e.g. SHARK Lagos et al. 2019) than an analytic closed box. The linear mapping breaks down only in case of extreme inflow of pristine low metallicity gas, which is likely to be rare (Übler et al. 2014).

The analytical prescription of the gas-phase metallicity evolution is as follows:

$$\begin{array}{c}\hfill Z_{\mathrm{gas}}(t)=Z_{\mathrm{gas},\mathrm{init}}+(Z_{\mathrm{gas},\mathrm{final}}-Z_{\mathrm{gas},\mathrm{init}})\frac{1}{\mathcal{M}}{\displaystyle {\int }_0^{t^{\prime }}}\mathrm{SFR}(t)\mathrm{d}t,\end{array}$$(17)

where Z_gas, init is the initial metallicity for the earliest phases of star formation, ℳ is stellar mass in units of M_⊙ drawn from the GSMF, SFR(t) is the galaxy SFH, and Z_gas, final is the present-day gas-phase metallicity of the object. This prescription links the chemical enrichment of the galaxy to the SFR, such that increased star formation is associated with an increased rate of metal production. The initial metallicity, Z_gas, init, is taken to be equal to the lowest metallicity of the MIST (Dotter 2016; Choi et al. 2016) isochrones used in PROSsc pect, while the maximum value that Z_gas, final can assume depends on the highest metallicity of the MIST isochrones. The gas-phase metallicity is important not only to establish the metallicity of the next generation of stars, but also to model the gas emission contribution to the galaxy SED, which depends on the gas-phase metallicity itself and on the gas ionisation (see Sect. 3.5).

The relation between the gas-phase metallicity, often expressed in terms of the gas-phase oxygen abundance log(O/H), and the stellar mass is known as the mass-metallicity relation (MZR, Tremonti et al. 2004; Kewley & Ellison 2008; Foster et al. 2012; Maiolino & Mannucci 2019 and references therein). This relation points to a positive correlation between the two quantities in such a way that the more massive (or luminous) the galaxy is, the higher the fraction of metals is going to be in comparison to its lower-mass counterparts. This is due to the increased gravitational potential well depth that allows massive galaxies to retain a large fraction of their metals within their halo. The relation flattens at high stellar masses where the mass growth mechanism is not related to star formation anymore, but rather to gas-poor mergers. Red galaxies constitute the majority of objects at the high mass, high metallicity end of the relation, whereas the low-metallicity part is dominated by blue galaxies (Bellstedt et al. 2021). This distribution points to the fact that in reality the MZR represents only a projection of a more fundamental relation, known as Fundamental metallicity relation (FMR, Mannucci et al. 2010; Lara-López et al. 2013), which links together the stellar mass, the gas-phase metallicity and the SFR of a galaxy.

In GALSBI-SPS, the present-day gas-phase metallicity Z_gas, final of an individual galaxy is obtained from the FMR which is a redshift-dependent relation that involves the galaxy stellar mass and the SFR. The present-day gas-phase metallicity of a galaxy, expressed as gas-phase oxygen abundance log(O/H), is drawn from a normal distribution whose mean is described by the parametrisation introduced in Bellstedt et al. (2021):

$$\begin{array}{cc}\hfill [12+log(\mathrm{O}/\mathrm{H})]=& \alpha (t_{\mathrm{lb}})log\left(\frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}\right)+\hfill \\ \hfill & \beta (t_{\mathrm{lb}})log\left(\frac{\mathcal{M}}{M_{\odot }}\right)+\gamma (t_{\mathrm{lb}}),\hfill \end{array}$$(18)

where t_lb is the observer look-back time to the galaxy redshift in Gyr and α(t_lb),β(t_lb),γ(t_lb) evolve quadratically with look-back time:

$$\begin{array}{cc}\hfill \alpha (t_{\mathrm{lb}})=& {\alpha }_0+{\alpha }_1{t}_{\mathrm{lb}}+{\alpha }_2{t}_{\mathrm{lb}}^2\hfill \\ \hfill \beta (t_{\mathrm{lb}})=& {\beta }_0+{\beta }_1{t}_{\mathrm{lb}}+{\beta }_2{t}_{\mathrm{lb}}^2\hfill \\ \hfill \gamma (t_{\mathrm{lb}})=& {\gamma }_0+{\gamma }_1{t}_{\mathrm{lb}}+{\gamma }_2{t}_{\mathrm{lb}}^2.\hfill \end{array}$$(19)

The scatter around this mean relation also evolves quadratically with look-back time:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{FMR}}(t_{\mathrm{lb}})={\sigma }_{\mathrm{FMR},0}+{\sigma }_{\mathrm{FMR},1}{t}_{\mathrm{lb}}+{\sigma }_{\mathrm{FMR},2}{t}_{\mathrm{lb}}^2.\end{array}$$(20)

We use for these parameters the values quoted in Bellstedt et al. (2021) and report them in Table 3. These values have been obtained by fitting a sample of star-forming galaxies selected via sSFR cut, log(sSFR/yr⁻¹) > −11.5. As a first approximation for this work, we extend the validity of the relation also to red galaxies, using the same parameter values for both population of objects. The present-day gas-phase metallicity in solar units can be obtained from the one expressed in terms of oxygen abundance using the relation

$$\begin{array}{c}\hfill \frac{Z_{\mathrm{gas},\mathrm{final}}}{Z_{\odot }}={10}^{[12+log(\mathrm{O}/\mathrm{H})]-{[12+log(\mathrm{O}/\mathrm{H})]}_{\odot }},\end{array}$$(21)

where Z_⊙ is the solar metallicity in absolute units and [12+log(O/H)]_⊙ is the solar oxygen abundance. This relation allows us to sample a relative-to-solar metallicity value, such that, depending on the specific absolute value for the gas-phase solar metallicity (e.g. Z_⊙ = 0.020 and [12+log(O/H)]_⊙ = 8.69 in Asplund et al. 2009), one can accordingly scale the absolute value of Z_gas, final. The linear metallicity evolution is implemented in PROSsc pect through the ZFUNC_MASSMAP_LIN function, while the present-day gas-phase metallicity Z_gas, final is passed to the ZFINAL PROSsc pect parameter.

As for stellar metallicity Z_*, i.e. the fraction of metals lock in stars, a similar relation with the stellar mass exists (Gallazzi et al. 2005; Zahid et al. 2017; Lian et al. 2018; Looser et al. 2024). Rather than sampling from another parametric relation, we follow the prescription also adopted in other studies (e.g. Alsing et al. 2023; Leja et al. 2017) that the stellar metallicity is equal to that of the gas-phase metallicity at all times. This assumption, known as instantaneous enrichment, postulates that the stars have the same metallicity as the ISM from which they had formed. It is important to point out, as detailed for example in Fraser-McKelvie et al. (2022), that the integrated stellar populations are generally more metal-poor than the ISM, especially in low-mass galaxies, even when the same metallicity indicator is adopted. This is also found in analytical models, such as the gas regulator or ‘bathtub’ model (Lilly et al. 2013; Peng & Maiolino 2014), which predicts a difference of roughly ∼0.2 dex. We leave the implementation of a more complex prescription that takes into account these differences to future work.

3.5. Gas ionisation

The origin of emission lines in optical galaxy spectra stems from nebular gas ionised either by stellar radiation or by radiation from the accretion disc of an AGN, with the latter confined to the central regions of galaxies. Nebular gas ionised by radiation from massive stars is characterised by a set of physical properties that include the electron density n_e, gas-phase metallicity Z_gas (discussed in Sect. 3.4) and ionisation parameter log U, defined as the ratio of ionising photons to the total hydrogen density. A careful modelling of the gas emission contribution to the galaxy SED therefore requires an understanding of the relation of these physical properties to the global properties of the host galaxies, such as stellar mass and SFR.

The emission line fluxes are typically computed by means of photoionisation codes such as CLOUDY (Ferland et al. 2017) or MAPPINGS (Dopita & Sutherland 1996). These codes generate grids of emission line ratios by solving radiative transfer equations for multiple ionic species. They adopt different types of input ionising radiation with a dependence on gas properties, such as ionisation parameter and gas-phase metallicity, that can be fully modelled and explored (Kewley & Dopita 2002; Dopita et al. 2013; Gutkin et al. 2016). Running photoionisation codes for large samples of galaxies is computationally expensive. Therefore, in almost all of the SPS codes available in the literature (see Fig. 1 in Thorne et al. 2021 for an overview), the emission line fluxes come from pre-computed tables that are generated by making strong simplifying assumptions about the geometry, usually a point source at the centre of a spherical shell, and the gas composition to reduce the number of free parameters (Byler et al. 2017)⁸.

PROSsc pect makes no exception, adopting photoionisation tables from Levesque et al. (2010) to produce line emission features for a range of gas-phase metallicities and ionisation parameters. The star formation nebular emission features are computed in PROSsc pect using a simple energy balance scheme. Since efficient Hydrogen ionisation is largely caused by young O and B stars continuum flux short of the Lyman limit (911.8 Å), PROSsc pect assumes that the integrated intrinsic stellar flux at λ < 911.8 Å of stars younger than 10 Myr determines the ionised flux of the nebulae (Orsi et al. 2014; Robotham et al. 2022), with no UV photon escape fraction in the default settings. This ionising flux is then redistributed across known significant emission features interpolating between the gridded values of Z_gas and log U available from the Levesque et al. (2010) tables. The latter do not include the contribution from the nebular continuum, which has been recently shown by Miranda et al. (2025) to provide a significant contribution to the galaxy broad bands only for the extreme emission line galaxies, which are increasingly important at higher redshifts.

Another strong simplifying assumption that is adopted in pre-computed tables is to assume a fixed electron density of n_e = 100 cm⁻³ (Byler et al. 2017), which is also the assumption we use in GALSBI-SPS as we are adopting the Levesque et al. (2010) fluxes implemented in PROSsc pect. This assumption is usually motivated by the interest in targeting the gas emission from H II regions due to star formation, ignoring the contribution of AGN broad line regions where the electron density is higher. However, it is well known that the typical electron density of H II regions in low redshift galaxies is of the order of n_e ∼ 30 − 120 cm⁻³ (Brinchmann et al. 2008; Davies et al. 2021), while, for high redshift galaxies, studies found values up to a few ×10³ cm⁻³ (Masters et al. 2014; Shirazi et al. 2014; Kaasinen et al. 2017; Kashino et al. 2017), with an average evolution at a rate of ∝(1 + z)^1 − 2 (Isobe et al. 2023; Abdurro’uf et al. 2024).

The ionisation parameter values change as function of redshift and galaxy physical properties. Spectroscopic surveys of objects at 1 < z < 7 have shown that there is a general increase of the ionisation parameter with redshift at fixed stellar mass (Brinchmann et al. 2008; Shirazi et al. 2014; Kashino et al. 2017; Kaasinen et al. 2018; Papovich et al. 2022). This is typically attributed to the presence of lower log(O/H), harder ionising spectra, and higher electron densities, characteristic of the high-redshift galaxy population. Indeed, there have been reports in the literature of an inverse correlation between the gas-phase metallicity and the ionisation parameter, since the increasing absorption due to a higher opacity in the stellar winds leads to a reduction at high metallicity of the stellar ionising photon flux illuminating the gas (Orsi et al. 2014; Kashino & Inoue 2019). Other studies (Kaasinen et al. 2018) instead found correlations between log U and sSFR, rather than with metallicity. Therefore, it is still unclear how exactly these parameters are connected and whether log U and Z_gas are intrinsically correlated or not.

In order to take into account the possible correlation among these parameters and let the data constrain their dependency, we modelled the gas ionisation log U following the prescription in equation 12 of Kashino & Inoue (2019), where we assume that the galaxy can be modelled as a single H II region described by a single ionisation parameter ⁹. We sampled log U for individual galaxies from a truncated Normal distribution with mean

$$\begin{array}{c}\hfill logU={\upsilon }_0+{\upsilon }_1(log(\mathrm{O}/\mathrm{H})+4)+{\upsilon }_{2}log{n}_e+{\upsilon }_3(log(\mathrm{sSFR})+9),\end{array}$$(22)

and scatter around the mean relation equal to σ_{log U} = 0.1 dex, as quoted in Kashino & Inoue (2019). The electron density is fixed to n_e = 100 cm⁻³, while log(O/H) is computed from Z_gas, final as

$$\begin{array}{c}\hfill log(\mathrm{O}/\mathrm{H})=log\left(\frac{Z_{\mathrm{gas},\mathrm{final}}}{Z_{\odot }}\right)+{[12+log(\mathrm{O}/\mathrm{H})]}_{\odot }-12.\end{array}$$(23)

We computed logsSFR from the stellar mass and the SFR drawn from the relations in Sects. 3.1 and 3.3, log(sSFR yr⁻¹) = log(SFR/ℳ). We limit the sampled parameters in the range log U = [ − 4, −1], which is the typical range used in CLOUDY and MAPPINGS. The parameters for this relation are set to the best-fitting values in Kashino & Inoue (2019) and reported in Table 4. All the parameters are well constrained in Kashino & Inoue (2019), fully rejecting the null hypothesis of no dependence on said parameters. Kashino & Inoue (2019) derived this relation by measuring the nebular conditions from the emission-line flux ratios of a sample of low redshift 0.027 ≤ z ≤ 0.25 galaxies from the Sloan Digital Sky Survey (SDSS) that have been stacked into bins of stellar mass and sSFR. Although in this sample the gas-phase metallicity has been estimated for objects having log(sSFR yr⁻¹) ≥ −10.5, which are therefore star-forming objects, we apply this relation also to the red population of galaxies. This is motivated by the fact that the impact of the gas emission on the red galaxies spectra (and therefore magnitudes) is in principle negligible.

The parameter log U is a dimensionless quantity that expresses the ratio of hydrogen ionising photons to total hydrogen density. In the literature, some studies refer to the ionisation parameter as log U (Byler et al. 2017), while others use the quantity q [cm s⁻¹] (Levesque et al. 2010). The two definitions simply differs by the speed of light in units of [cm s⁻¹], q = U × c. PROSsc pect expresses the ionisation parameter of the gas in terms of q and allows the user to relate this parameter to the provided gas-phase metallicity using equation 5 in Orsi et al. (2014). To account for the possible dependence on other physical quantities, we provide the code with our value of the ionisation parameter obtained by converting the quantity sampled from Equation 22 into the appropriate units recognised by our generative SED engine.

3.6. Dust component

In Tortorelli et al. (2024), we showed that the dust attenuation was one of the SED modelling components that had the strongest impact on the galaxy colours (up to ∼1 mag in u − g) and on the mean redshift of the tomographic bins (Δz ∼ 0.1). Therefore, a realistic modelling of dust and its relation with the galaxy physical properties is a key ingredient for a SPS-based galaxy population model. Dust forms in the atmospheres of evolved stars or in remnants of supernovae (Draine 2003) and released in the ISM where, observationally, modifies the stellar continuum through the absorption and the scattering of starlight and the re-emission of the absorbed energy in the infrared. The dust attenuation 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.

Observations from large samples of galaxies at low redshift (Salim et al. 2018) show that the dust attenuation laws vary as function of galaxy types even among similar mass galaxies at different redshifts (Salim & Narayanan 2020). They also change between lines of sight within the same galaxy (Maíz Apellániz 2024). This points to the fact that it does not exist a single dust attenuation law with fixed parameters that is generic enough that it can be applied to the whole galaxy population, but rather that the use of a dust attenuation prescription should be tightly linked to the spectral type and physical properties of the individual galaxies (Nagaraj et al. 2022).

We modelled the effect of dust in GALSBI-SPS using the Charlot & Fall (2000) two-phase dust attenuation model, which modifies the intrinsic flux f_intr into the observed attenuated one f_att at a given wavelength λ through the attenuation factor A(λ),

$$\begin{array}{c}\hfill f_{\mathrm{att}}(\lambda )=f_{\mathrm{intr}}(\lambda )\times {10}^{-A(\lambda )/2.5},\end{array}$$(24)

where A(λ) is related to the effective optical depth of attenuation, τ, by

$$\begin{array}{c}\hfill A(\lambda )=e^{-\tau {(\lambda /{\lambda }_0)}^{\nu }}.\end{array}$$(25)

The parameter τ changes the column density of dust, λ₀ = 5500 Å is the pivot wavelength and ν is the modifying power-law index, which we set to −0.7. This model splits the galaxy starlight in a component coming from stars younger than 10 Myr and another from stars older than 10 Myr, which is the typical timescale for the disruption of a molecular cloud (Blitz & Shu 1980). Older stars are solely attenuated by a dust ‘screen’, regulated by the τ_ISM parameter, which represents the dust in the diffuse ISM. Young stars are instead attenuated by the ‘dust screen’ and by an additional ‘birth cloud’ component that adds extra attenuation towards young stars, simulating their embedding in molecular clouds and H II regions. The ‘birth cloud’ component is regulated by the τ_BC parameter. From the definition, it is clear that τ_ISM is the parameter that has the strongest effect on the galaxy SED, while τ_BC mostly impacts the UV flux (see Fig. 3 in Thorne et al. 2021). An additional factor that impact the dust attenuation in the UV is the 2175 Å bump. This a prominent absorption feature that is postulated to arise from polycyclic aromatic hydrocarbon (PAH) dust grains, even though the actual carrier has not been definitively established (Battisti et al. 2025). The bump seems to be redshift-independent (Battisti et al. 2022) and in general its strength does not correlated with any galaxy property, with a median value of roughly 1/3 of the MW value (Battisti et al. 2020) among different type of galaxies. A recent work by Battisti et al. (2025) reported though an anti-correlation with surface SFR, in such a way that the bump is supposed to be applied only to old stars, since the carrier seems to be easily destroyed by UV photons. Given the size of the sample and the uncertainty concerning the dependency of the bump strength on galaxy properties, we use the default value for the bump strength in the PROSsc pect SED computation, which is set to this roughly universal value of 1/3 of the MW value.

In order to sample the dust parameters for our generative model in a physically meaningful way, we use the DEVILS data to model the dependence of τ_ISM and τ_BC on the integrated galaxy properties. We do not find any dependence of log(τ_BC) on redshift, stellar mass or sSFR. Therefore, this parameter is simply sampled from a Normal distribution that has different means μ_log(τBC) and scatters σ_log(τBC) for the population of blue and red galaxies, truncated at the limits in Table 2 of Bellstedt et al. (2020). On the contrary, we find that log(τ_ISM) does show a dependence on redshift, stellar mass and sSFR. We decide to model this dependence as

$$\begin{array}{cc}\hfill log({\tau }_{\mathrm{ISM}})=& {\tau }_{\mathrm{ISM},0}+{\tau }_{\mathrm{ISM},1}log(1+z)+\hfill \\ \hfill & {\tau }_{\mathrm{ISM},2}log\left(\frac{\mathcal{M}}{M_{\odot }}\right)+{\tau }_{\mathrm{ISM},3}log(\mathrm{sSFR}).\hfill \end{array}$$(26)

In order to obtain the best-fitting parameters of this relation, we fit an hyper-plane to the DEVILS catalogue data using the LTSFIT library (Cappellari et al. 2013). Figure 5 shows the best-fitting relation expressing the dependence of log(τ_ISM) on redshift, stellar mass and sSFR. log(τ_ISM) is then sampled from a Normal distribution with mean following 26 and scatter σ_{obs, log(τISM)} around this relation measured from the data, truncated at the limits in Table 2 of Bellstedt et al. (2020). The parameter values of the relations from which we sampled log(τ_BC) and log(τ_ISM) for the population of red and blue galaxies are reported in Table 5. The two-phase dust attenuation is directly implemented in PROSsc pect by simply providing the values of τ_BC and τ_ISM as input to the parameters TAU_BIRTH and TAU_SCREEN, respectively.

Fig. 5. Correlation between log(τISM) and redshift, stellar mass and sSFR for star-forming (left panel) and quiescent (right panel) galaxies from DEVILS data. The solid red (blue) line shows the best-fitting relation that expresses this dependence for star-forming (quiescent) galaxies. The dotted lines show the observed rms scatter of the relation. Each sample contour encloses 50%, 84%, and 99% of the values.

PROSsc pect contains also an energy balance prescription that re-emits the bolometric sum of attenuated light in the FIR using the Dale et al. (2014) templates. As discussed in Tortorelli et al. (2024) and also visible in Fig. 3 in Thorne et al. (2021), this emission impacts the galaxy SED mostly beyond the 20000 Å rest-frame, which is in the Ks band, the reddest band of interest, thereby providing an impact on galaxy colours that is negligible for our science case. Therefore, we set the ALPHA_SF_BIRTH and ALPHA_SF_SCREEN parameters in PROSsc pect to their default values (1 and 3, respectively) and we do not model the dust emission parameter distributions, leaving a more thorough modelling to a future analysis.

3.7. Velocity dispersion

The velocity dispersion in an integrated spectrum is a light-weighted superposition of all the velocity dispersions of the individual stellar populations. Hence, the observed value is going to be similar to that of the stars dominating the light from the galaxy. From an SPS point of view, the velocity dispersion manifests itself as a broadening of the spectral lines that is usually modelled as a Gaussian smoothing in velocity space in km s⁻¹. As shown in Tortorelli et al. (2024), varying the velocity dispersion of a galaxy does not make an appreciable contribution to the broad band fluxes. However, we include a model for sampling the velocity dispersion in GALSBI-SPS as our aim is to build a model that is realistic enough to be used to also forward-model spectroscopic surveys.

As in Tortorelli et al. (2024), we modelled the velocity dispersion σ_disp for our population of galaxies following the stellar mass-velocity dispersion relation introduced in Zahid et al. (2016):

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{disp}}(\mathcal{M})& ={\sigma }_b{\left(\frac{\mathcal{M}}{M_b}\right)}^{s_1}\mathrm{for}\mathcal{M}\le M_b\hfill \\ \hfill {\sigma }_{\mathrm{disp}}(\mathcal{M})& ={\sigma }_b{\left(\frac{\mathcal{M}}{M_b}\right)}^{s_2}\mathrm{for}\mathcal{M}>M_b\hfill \end{array}$$(27)

where the values of M_b, σ_b, s₁, s₂ are reported in Table 6 and are taken from the first row of Table 1 in Zahid et al. (2016). These values are redshift independent, as Zahid et al. (2016) showed that the Sloan Digital Sky Survey and Smithsonian Hectospec Lensing Survey data at z < 0.7 are consistent with no evolution in redshift. We sampled the velocity dispersion from a Normal distribution whose mean follows the relation in Equation (27) and whose scatter is taken from the values quoted in Zahid et al. (2016). The minimum velocity dispersion is set to 10 km s⁻¹. We smooth the galaxy spectrum as a Gaussian broadening in velocity space, where the Gaussian FWHM is FWHM = 2.35482 × σ_disp. The smoothing is applied prior to adding the AGN flux contribution.

3.8. Active galactic nuclei

Active galactic nuclei are powered by supermassive black-holes (M_BH ≳ 10⁶ M_⊙) in the centre of galaxies that are actively accreting mater. They have high luminosities (L_bol, AGN ≲ 10⁴⁸ erg s⁻¹) and strong evolution of their luminosity function with redshift (Kulkarni et al. 2019; Shen et al. 2020; Thorne et al. 2022b). Their emission covers the whole electromagnetic spectrum (Padovani et al. 2017), with our wavelength range of interest, from the UV to the infrared, being powered by the accretion disc, mostly in the form of broad and narrow emission lines strongly ionised by the central engine, and by the thermal emission of the dust heated by the central engine. When actively accreting matter and contributing more than 10% to the galaxy bolometric luminosity, an AGN can dramatically modify galaxy colours, inducing strong changes in the redshift distribution of galaxies detected in a survey (Tortorelli et al. 2024).

To model AGN hosts in GALSBI-SPS, we adopt the prescription of López-López et al. (2024). We first computed the probability that a galaxy hosts an AGN as function of redshift and stellar mass by taking the ratio between the redshift-evolving GSMFs of blue and red galaxies and the AGN host galaxy mass function (HGMF) derived in Bongiorno et al. (2016). Then, we statistically assign an AGN to a host galaxy with a Bernoulli trial proportional to the computed probability. As in López-López et al. (2024), we assume fixed AGN fractions outside the redshift ranges specified by Bongiorno et al. (2016). In particular, the probability of hosting an AGN is 1% at all stellar masses below z = 0.55, while at z > 2 we maintain a constant fraction taken from the last redshift bin in Bongiorno et al. (2016). These assumptions are support by recent observational studies (Birchall et al. 2022; Williams et al. 2022; Osorio-Clavijo et al. 2023).

To assign an AGN spectral contribution, we adopt a different set of AGN templates with respect to those implemented in PROSsc pect, which are the Dale et al. (2014) model, the extended AGN template from Andrews et al. (2018) and the Fritz et al. (2006) model. These templates model the emission from the AGN as composition of power laws with different spectral indices; however, AGN spectra in the optical/UV are also characterised by the presence of broad and/or narrow emission lines coming from the strongly ionised material in the accretion disc. Similarly to what we did in Tortorelli et al. (2024), we replace the AGN model in PROSsc pect with the parametric quasar SED model presented in Temple et al. (2021), calibrated to match the observed u, g, r, i, z, Y, J, H, K, W1, W2 colours of the SDSS DR16 quasar population (Lyke et al. 2020). This model covers the full AGN contribution in the 912 Å≲λ ≲ 30 000 Å rest-frame wavelength range by modelling the low-frequency tail of the direct emission from the accretion disc (900 Å ≲ λ ≲ 10000 Å), the emission from the hot dust torus (10000 Å ≲ λ ≲ 30000 Å), the Balmer continuum (∼ 3000 Å), the broad and narrow emission lines in the optical/UV region with the Baldwin effect (Baldwin 1977) applied, and the Lyman-absorption suppression at λ < 912 Å. The Temple et al. (2021) model will be implemented directly into PROSsc pect in a future release of the code¹⁰.

The Temple et al. (2021) model requires as input the monochromatic 3000 Å continuum luminosity and the absolute i-band magnitude at z = 2, as defined by Richards et al. (2006), to produce the AGN spectrum with the right units and normalisation. We relate the bolometric luminosity of the AGN L_bol, AGN to these two quantities by following the prescription in the Temple et al. (2021) templates code repository¹¹. The bolometric AGN luminosity, in turns, is related to the galaxy bolometric luminosity L_{bol, galaxy} through the f_AGN parameter, L_bol, AGN = f_AGN × L_{bol, galaxy}, which is the only parameter we are sampling in GALSBI-SPS to add the contribution from the AGN to the galaxy SED. We modelled the f_AGN distribution from the DEVILS data used in Thorne et al. (2022b). We do not use the values of f_AGN reported in this paper though, as its definition is different from the one highlighted above. We instead recomputed f_AGN to match our definition by re-running PROSsc pect with the exact same setup as in Thorne et al. (2022b). We modelled log(f_AGN) as a two component mixture of Normal distributions. Figure 6 shows the best-fitting distributions for red and blue galaxies in DEVILS separated by sSFR superimposed on the newly computed f_AGN values. The best-fitting mean and standard deviation of the Normal distributions are reported in Table 7. We then assign AGN flux contributions to all galaxies, regardless of whether they are classified as AGN hosts. However, f_AGN is sampled from different distributions depending on the AGN classification. Galaxies hosting AGNs draw from the high-log(f_AGN) distribution, while non-AGN galaxies draw from the low-log(f_AGN) one. We adopt a threshold of log(f_AGN) =  − 1 to separate these two regimes (Thorne et al. 2022b).

Fig. 6. Histogram of the ratio between the AGN bolometric luminosity and the galaxy luminosity in log-space, log(fAGN), for star-forming (blue histogram) and quiescent (red histogram) galaxies. The dashed lines represent the samples drawn from the sum of two normal distributions with which we modelled the log(fAGN) values from DEVILS data.

3.9. Magnitude computation

The galaxy physical properties sampled in the previous sections are essential to create galaxy SEDs with PROSsc pect. The parameters we provide as input to the code are mSFR, mpeak, mperiod and mskew for the SFH, τ_BC and τ_ISM for the two component dust attenuation, Z_gas, final for the present day value of the linearly evolving gas-phase metallicity, the ionisation parameter q, and the SSP library described in Sect. 3.2. We set the redshift to z = 0 in order to create a rest-frame SED for each galaxy in units of L_⊙/Å. We then apply a Gaussian broadening to the rest-frame spectra to model the effect of the galaxy velocity dispersions (see Sect. 3.7). We computed the observed-frame spectra in units of erg s⁻¹ cm⁻² Å⁻¹ using the luminosity distance D_L for each object obtained with PYCOSMO (Refregier et al. 2018). The addition of an AGN contribution to the flux is left as a choice to the user through a specific flag in the code. The AGN contribution is computed with the Temple et al. (2021) templates as detailed in Sect. 3.8. The output AGN SED is also in the observed-frame and it can directly be added to the galaxy flux contribution. We then apply reddening due to galactic extinction to the observed-frame spectra using the Fitzpatrick (1999) extinction law. The extinction map of the Milky-Way is the one presented in Schlegel et al. (1998), where the E(B − V) values are obtained based on the galaxy RA and Dec coordinates. The final spectra are then integrated in the g, r, i, z, y HSC wavebands to obtain the observed-frame AB magnitudes. We also computed the rest-frame magnitudes by de-redshifting the galaxy spectra before applying the galactic extinction. The Python package GALSBI further allows the user to skip the rest-frame SED computation and save computational time by obtaining directly the observed and rest-frame magnitudes from the sampled physical quantities using the PROMAGE magnitude emulator presented in Tortorelli et al. (2025).

3.10. Distribution of morphological parameters

The modelling components described in the previous sections have been used to generate realistic galaxy SEDs and, consequently, magnitudes. A crucial aspect of our forward-modelling approach is the ability to simulate realistic survey images with UFIG. This requires not only magnitudes, but also assigning each galaxy a light profile, which we assume to be a Sérsic light profile (Sérsic 1963), and sampling additional morphological parameters, such as sizes, Sérsic indices, ellipticities and positions. We leave the use of more complex morphologies, such as bulge and disc, to a future work, as their relation with galaxy SFHs and physical properties require more careful modelling as shown in Robotham et al. (2022) and Bellstedt et al. (2024).

Galaxy sizes are key observational quantities to understand the formation pathways of galaxies since they imprint on the galaxy structure (e.g. Mo et al. 1998). While it is well established that galaxy sizes evolve with redshift (e.g. Daddi et al. 2005; Trujillo et al. 2007; van Dokkum et al. 2008), the degree with which they do so and whether this size growth is due to major mergers (e.g. Naab et al. 2007), minor mergers (e.g. Naab et al. 2009), progenitor bias (e.g. Fagioli et al. 2016), feedbacks (e.g. Damjanov et al. 2009), or other processes, is still under debate. Galaxy physical sizes in kpc exhibit a tight scaling relation with the stellar mass known as size-stellar mass relation (e.g. van der Wel et al. 2014). This is well characterised at low redshift and at high mass end (e.g. Shen et al. 2003), while the behaviour at higher redshift (e.g. Yang et al. 2021) and the low-mass end (e.g. Lange et al. 2015) is still subject to intense research.

In Fischbacher et al. (2025a), the authors model galaxy sizes as a log-normal distribution whose mean and scatter were related to the galaxy absolute magnitude (Shen et al. 2003). In our work, galaxy sizes, expressed in terms of the half-light radii in kpc, are also sampled from a log-normal distribution, but the mean of this distribution is given by size-stellar mass relation, while the standard deviation is taken from Shen et al. (2003). The size-mass relation has a different parametrisation for the population of blue and red galaxies. We adopt the parametrisations introduced in Nedkova et al. (2021), where the authors measured the size-stellar mass relation for a large sample of blue and red galaxies in clusters and fields with HST photometry in the redshift range 0.2 < z < 2, extending the sample down to stellar masses of ∼10^7.5 M_⊙. The dependence of the galaxy size with stellar mass is parametrised for blue galaxies as a single power law,

$$\begin{array}{c}\hfill R_e=A{\left(\frac{\mathcal{M}}{5\times {10}^{10}{M}_{\odot }}\right)}^B,\end{array}$$(28)

where log A and B are assumed to evolve linearly with redshift:

$$\begin{array}{cc}\hfill logA& =log{A}_{\mathrm{intcpt}}+log{A}_{\mathrm{slope}}\times z,\hfill \\ \hfill B& =B_{\mathrm{intcpt}}+B_{\mathrm{slope}}\times z.\hfill \end{array}$$(29)

This assumption is motivated by fitting the redshift evolution of the log A and B values from Table 3 in Nedkova et al. (2021) for the sample of blue galaxies covering the whole mass range probed (left panel of Fig. 7). We report the best-fitting parameters of the linear fit of log A and B in Table 8. The size-mass relation for red galaxies tends instead to flatten at low stellar masses; therefore, it is parametrised as a double power law,

Fig. 7. Best-fitting linear evolution in redshift of the size-stellar mass parameters for star-forming (left panel) and quiescent (right panel) galaxies. The scatter points represent the observed values of the parameters and are taken from Tables 2 and 3 in Nedkova et al. (2021). The dashed lines show the best-fitting relations, while the bands represent the 1σ uncertainty on the fit.

$$\begin{array}{c}\hfill R_e=\zeta {(\mathcal{M})}^{\eta }{(1+\frac{\mathcal{M}}{{10}^{\delta }})}^{\theta -\eta },\end{array}$$(30)

where η and θ are the slope at the low and high mass end, ζ is the normalisation, while 10^δ is the stellar mass at which the second derivative of the function is at a maximum. All these parameters also evolve linearly with redshift:

$$\begin{array}{cc}\hfill \eta & ={\eta }_{\mathrm{intcpt}}+{\eta }_{\mathrm{slope}}\times z,\hfill \\ \hfill \theta & ={\theta }_{\mathrm{intcpt}}+{\theta }_{\mathrm{slope}}\times z,\hfill \\ \hfill log\zeta & =log{\zeta }_{\mathrm{intcpt}}+log{\zeta }_{\mathrm{slope}}\times z,\hfill \\ \hfill \delta & ={\delta }_{\mathrm{intcpt}}+{\delta }_{\mathrm{slope}}\times z.\hfill \end{array}$$(31)

The linear evolution fit was obtained from the values in Table 2 of Nedkova et al. (2021) for the whole mass range of the red sample (right panel of Fig. 7), while the obtained best-fitting parameters are reported in Table 8. In this work, we assume the size in kiloparsec to be wavelength independent, although observations do show an half-light radii dependence on wavelength (e.g. Tortorelli et al. 2023) due to the different stellar populations probed. The assumption of a single Sérsic profile with no radial dependence of the stellar population justify the choice of a wavelength independent size for this work.

Fig. 8. Stellar mass dependence of the r-band Sérsic index n. The deep sky blue and orange contours are obtained by matching the star-forming and quiescent galaxies selected via sSFR of the GAMA mass complete sample with the morphological catalogue in Kelvin et al. (2012). The blue and red error bars represent the r-band Sérsic index median trend in stellar mass bins for star-forming and quiescent galaxies, respectively, with the size of the error bars representing the standard deviation in said bin. The dashed lines represent the evolution with stellar mass obtained from the best-fitting parameters of Equation 32 on the observed data. The bands represent the 1σ uncertainty on the best-fitting relations. Each sample contour encloses 50%, 84%, and 99% of the values.

In Fischbacher et al. (2025a), galaxy Sérsic indices are drawn from a redshift-independent beta prime distribution, since extending the parametrisation of the Sérsic index by allowing its evolution with redshift did not improve the distance metrics between the simulations and the observations. This is consistent with findings in the literature (Martorano et al. 2025) that show little to no evolution of the median Sérsic index in the redshift range 0.5 < z < 2.5. Martorano et al. (2025) also found that high-mass galaxies have higher Sérsic indices than lower-mass galaxies at all redshift and that star-forming galaxies have lower Sérsic indices than quiescent galaxies. Therefore, in our work, we introduce a stellar mass dependence of the Sérsic index for both blue and red galaxies. To model this dependence, we match the GAMA mass complete sample with the single-component Sérsic fit catalogue presented in Kelvin et al. (2012) and downloaded from the GAMA II database, separating blue and red galaxies via the sSFR. We use the r-band Sérsic index n to model its dependence on stellar mass and we remove from the catalogue bad fits utilising the primary χ² of both the galaxy and the PSF, as suggested by the authors. We modelled the relation as

$$\begin{array}{c}\hfill n=n_0+n_1\times {\left(\frac{log\mathcal{M}}{log({10}^{10.5}{M}_{\odot })}\right)}^{n_2},\end{array}$$(32)

where 10^10.5 M_⊙ is the median stellar mass of the sample. Figure 8 shows the best-fitting evolution of n as function of stellar mass. As in Martorano et al. (2025), we also find that high-mass galaxies tend to have higher n than lower mass ones and that star-forming galaxies have lower n than quiescent galaxies at all masses. We sampled n from a truncated normal distribution, with mean given by Equation 32 and standard deviation σ_n given by the median scatter across the stellar mass range. The sampled values are within the range 0.2 < n < 10. The best-fitting parameters of Equation (32) and the scatter across stellar mass for blue and red galaxies are reported in Table 9.

Galaxy positions are randomly assigned within the UFIG simulated image. Future works are aiming at introducing a realistic clustering prescription, as developed in Berner et al. (2022, 2024) and Fischbacher et al. (2025c). This is an important point to keep in mind for the use of forward-modelling for cosmological analysis. Indeed, Fischbacher et al. (2025a) found that by not including clustering, and therefore not being affected by cosmic variance, the phenomenological model returns tomographic redshift distributions that are much smoother than those obtained from COSMOS2020 photo-zs (Weaver et al. 2022). The authors attributed this difference to the fact that the reported uncertainty from GALSBI is underestimating the real one since it is only including the uncertainty of the galaxy population model.

Finally, we sampled both blue and red galaxies absolute ellipticities from the same modified Beta distribution as in Moser et al. (2024) and Fischbacher et al. (2025a). The α_ell and β_ell parameters of the standard Beta distribution are related to the mode e_mode and spread e_spread parameters of the modified Beta via e_mode = (α_ell − 1)/(e_spread − 2) and e_spread = α_ell + β_ell. We set e_mode = 0.2 and e_spread = 2.9. The two ellipticity components e₁ and e₂ are computed by assigning a random phase to the absolute ellipticity, i.e. no shear or intrinsic alignment prescription applied.

4. HSC image simulation

Hyper-Suprime Cam is an excellent precursor dataset for Stage IV surveys and it has been used in Moser et al. (2024) and Fischbacher et al. (2025a) to both constrain and validate the phenomenological version of GALSBI. For these reasons, we consistently use this dataset to validate the performance of GALSBI-SPS against its phenomenological version and the observed data.

The catalogue of intrinsic galaxy properties sampled from the galaxy population model described in Sect. 3 is used to create simulated HSC DUD images with the image simulator UFIG (Bergé et al. 2013; Fischbacher et al. 2025b). UFIG needs as input the catalogue of galaxy properties and the instrumental characteristics of the HSC images, namely image sizes, pixel scales, zero-points, and central pixel coordinates. We take the instrumental properties from the metadata provided by the HSC database¹². We simulated the HSC images in the g, r, i, z, y broad bands. The r and i filters have been replaced by more uniform r2 and i2 filters, whose throughputs are the ones we use in our simulations. We did not simulate single exposures, but rather the image co-added versions through the use of exposure time maps, derived from the metadata, and 2D background maps, derived from each individual observed HSC image. The background noise is added to the simulated images in the form of a Gaussian noise plus a Lanczos resampling (Duchon 1979) of order 3 to create correlated noise, mimicking the contribution from sky brightness, readout and errors during data processing. The Gaussian noise mean is read off the image header, while the standard deviation is estimated by generating a 2D background image from the observations and then applying a 3σ sigma-clipping with PHOTUTILS (Bradley et al. 2024), ensuring a different standard deviation per pixel. The 2D background image is estimated from background-subtracted observed images, so no background subtraction is performed on simulations. The magnitude zeropoint is set to 27 mag/ADU for the HSC co-adds, while the CCD gain of a single exposure multiplied by the number of exposures per pixel gives us a rough estimate of the effective gain to convert between ADUs and number of photons. Galaxies are randomly distributed on the image and rendered via photon shooting according to the galaxy Sérsic profile, naturally including Poisson noise (Bergé et al. 2013; Bruderer et al. 2016). Galaxy profiles are then convolved with a point spread function (PSF) whose model is generated with a Convolutional Neural Network, presented in Herbel et al. (2018) and updated in Kacprzak et al. (2020). The model PSF is evaluated at the position of bright unsaturated stars (18 ≤ i ≤ 22) and then interpolated across the co-add using Chebyshev polynomials as detailed in Appendix C of Kacprzak et al. (2020). The observed star positions are taken from GAIA DR3 (Gaia Collaboration 2016, 2023). Simulated bright stars (i < 18) are not randomly drawn, but placed at the positions where bright stars are in real HSC DUD images. This allows us to apply the same bright object masks of HSC to our simulated images. Fainter stars are instead drawn from the Besançon model of the Milky Way (Robin et al. 2003). Lastly, we convert photons to ADU by dividing out the effective gain and saturate pixels that are above the maximum value of the observed data. More details on the HSC DUD image simulation procedure are provided in section 3.2 of Moser et al. 2024, while Figs. 2 and 3 of the same work shows the good agreement between the simulated and observed image background pixel distribution and co-added image.

One crucial aspect of forward-modelling is the ability to apply the same measurement steps and selection cuts on both data and simulations. Therefore, we do not use the photometric catalogue from COSMOS2020, but we measured galaxy photometric properties with SOURCE EXTRACTOR using the same settings for observed and simulated images. The use of the same settings ensures that all uncertainties in the extraction process are consistently modelled between data and simulations, resulting in the constraints on the galaxy population model being independent of the particular choice of the extraction software. We run SOURCE EXTRACTOR in dual-image mode, with the i-band image as detection image, using the hyper-parameters settings reported in Appendix C of Moser et al. (2024). While galaxies detected on the real images are matched to the COSMOS2020 catalogue by positions and magnitudes, galaxies detected on the simulated images are matched to the true injected galaxies by means of the SOURCE EXTRACTOR segmentation map. For each detection, we find the overlapping simulated object that minimises the sum of the differences between the measured magnitude (MAG_AUTO_b) and the true magnitude (mag_b) in all bands,

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{mag}={\displaystyle \sum _{b\in g,r,i,z,y}}|{\mathrm{mag}}_{\mathrm{b}}-\mathrm{MAG}\_\mathrm{AUTO}\_\mathrm{b}|.\end{array}$$(33)

This procedure improves the matching especially in the case of a crowded field. The matching procedure allows us to obtain catalogues of detected objects with measured quantities, such as magnitudes and sizes, for both data and simulations. These properties are then accompanied by high quality photo-z information, in the case of data, and true galaxy physical properties, including redshifts, in the case of simulations.

5. Results

In this section, we compared the photometric properties (magnitudes, colours and sizes), the redshift distributions and the scaling relations between physical properties of the galaxies that have been detected on the simulated and the observed images. The photometric measurements have been performed on the 56 HSC DUD patches fully overlapping with the COSMOS2020 photo-z catalogue, with the exact same SOURCE EXTRACTOR setup between observations and simulations. The simulated HSC images created in this work have the same instrumental characteristics as those in Fischbacher et al. (2025a). Therefore, we compared the detected galaxies drawn from GALSBI-SPS not only against observations, but also against those drawn from the phenomenological version of GALSBI. It is important to point out that in this section we are comparing the observed galaxy properties against those sampled from a model, GALSBI-SPS, whose parameter have not been explicitly constrained using SBI, as it is the case, instead, for the GALSBI phenomenological model.

5.1. Comparison of photometric properties

The SOURCE EXTRACTOR configuration file we use in our work is described in appendix C of Moser et al. (2024). The same setup was also used in Fischbacher et al. (2025a). We apply selection cuts on the catalogues of detected galaxies to avoid spurious detections and contaminations, separate galaxies from stars, and select objects with reliable photometry in every waveband. To ensure a meaningful comparison with GALSBI, the cuts we apply are the same ones that have been used in Fischbacher et al. (2025a) to compare the photometric properties and the redshift distributions between observations and simulations. The cuts are applied on the SOURCE EXTRACTOR measurements for every waveband with strict AND conditions:

$$\begin{array}{c}\hfill 0\le \&\mathrm{FLAGS}\le 3\\ \hfill 14<\&\mathrm{MAG}\_\mathrm{APER}3<30,\\ \hfill 16<\&\mathrm{MAG}\_\mathrm{APER}3\_\mathrm{i}<25,\\ \hfill 0<\&\mathrm{MAG}\_\mathrm{AUTO}<99,\\ \hfill 0<\&\mathrm{z}<6,\\ \hfill 0.1<\&\mathrm{FLUX}\_\mathrm{RADIUS}<30\\ \hfill -3<\&log(\mathrm{SNR})\equiv log\left(\frac{\mathrm{FLUX}\_\mathrm{AUTO}}{\mathrm{FLUXERR}\_\mathrm{AUTO}}\right)<4,\\ \hfill 0<\&\mathrm{ELL}<1,\\ \hfill 0<\&\mathrm{N}\_\mathrm{EXPOSURES}\\ \hfill 0.5<\&r_{50}/\mathrm{PSF}\_\mathrm{FWHM},\\ \hfill 0\le \&\mathrm{CLASS}\_\mathrm{STAR}\_\mathrm{i}<0.95.\end{array}$$(34)

FLAGS is an indicator for the photometric extraction quality, MAG_APER3 is the magnitude in a 3″ aperture, MAG_AUTO is the total magnitude (which is closest to the UFIG true magnitude), z is the redshift, FLUX_RADIUS is the radius containing 50% of the total light, SNR is the signal-to-noise ratio computed using the total flux FLUX_AUTO and its error FLUXERR_AUTO, ELL is the absolute ellipticity computed using the SOURCE EXTRACTOR windowed moments, N_EXPOSURES is the number of exposures in the co-added images at the position of the object, r₅₀ is the object size computed using the SOURCE EXTRACTOR windowed moments as in Kacprzak et al. (2020), PSF_FWHM is the point-spread-function full width at half-maximum, and CLASS_STAR is the SOURCE EXTRACTOR stellarity index. Since in the image simulations we assume that bright saturated stars are in the same position as those from the observations, we apply the HSC bright object mask to both observed and simulated images. The mask allows us to avoid sources whose photometry is impacted by the higher background value caused by nearby bright objects. We visually inspect galaxies close to bright objects in observed images, finding that the HSC bright object mask efficiently removes galaxies contaminated by starlight. The bright object induces in its proximity a higher background value that leads the extraction software to artificially increase the measured size of the detected galaxy.

Fig. 9. Fractional difference between the number of detected galaxies in the simulations, nsim, and in the observations, nobs, as a function of the i-band magnitude cut. Red scatter points refer to the GALSBI-SPS model, while blue box plots to the phenomenological version in Fischbacher et al. (2025a). The GALSBI-SPS error bars consider Poisson noise on the counts, while GALSBI box plots have been constructed using the posterior samples in Fischbacher et al. (2025a). The dashed black line marks the perfect agreement in the number counts with observations.

The first metric we compared is the number of detected galaxies, expressed as fractional difference between the number of detected galaxies in the simulations n_sim and in the observations n_obs:

$$\begin{array}{c}\hfill d_{\mathrm{frac}}=\frac{n_{\mathrm{sim}}-n_{\mathrm{obs}}}{n_{\mathrm{obs}}}.\end{array}$$(35)

Since we are using SOURCE EXTRACTOR in forced photometry mode and we are considering only those objects that have photometric measurements in all bands (we do not consider possible UV dropouts), the number of detected galaxies is band-independent. Figure 9 shows the fractional difference of detected galaxies as function of the i-band magnitude cut (i ≤ j, j = {21, 22, 23, 24, 25}) for both versions of GALSBI. The error bars for GALSBI-SPS consider Poisson noise on the counts, while the boxplots for GALSBI represent the distribution of counts estimated for each posterior sample. Since we determine the model parameters for GALSBI-SPS using mostly bright galaxy samples (see e.g. Fig. B.2 in Appendix B for the relations involving GAMA and DEVILS data), it is not surprising to see that the fractional difference of detected galaxies is at the percent level at the bright-end cut, progressively increasing at the faint-end cut (i ≤ 24,  i ≤ 25) where the model is mostly uninformed. By being constrained against HSC DUD deep data using SBI, the fractional difference of detected galaxies for GALSBI is always at the percent level.

The second metric we compared is the distribution of magnitudes in the five HSC filters g, r, i, z, y and sizes measured in the i-band. We focus on the i-band sizes because our current single Sérsic modelling does not account for radial gradients in the galaxy size. Galaxy magnitudes are measured on the HSC images using the SOURCE EXTRACTOR MAG_AUTO parameter, while galaxy sizes, defined as half-light radii in arcsec, are derived from the SOURCE EXTRACTOR FLUX_RADIUS parameter. Figures 10, 11, and 12 present comparisons of magnitude, size, and colour distributions among GALSBI-SPS (red contours), GALSBI (blue contours) and observed measurements (black contours) for various i-band magnitude cuts. We focus on the i ≤ 21, 22, 23 samples, as this is the magnitude range spanned by the data we use to create our forward model. For completeness, comparisons for the fainter i ≤ 24, 25 samples are included in Appendix D.

Fig. 10. Comparison of magnitude, size, and colour distributions between the observed HSC galaxies (black contours) and those from GALSBI (blue contours) and GALSBI-SPS (red contours) models for the i ≤ 21 cut. The lower left contours refer to the magnitude and size distributions, while the upper right to the colour distributions. MAG_AUTO refers to the galaxy magnitudes as measured by SOURCE EXTRACTOR, FLUX_RADIUS to the galaxy effective radii, while the colours are obtained as the difference between the SOURCE EXTRACTOR magnitudes. In the figure we report the GALSBI contours obtained from each posterior distribution sample in Fischbacher et al. (2025a).

Fig. 11. Comparison of magnitude, size, and colour distributions between the observed HSC galaxies (black contours) and those from GALSBI (blue contours) and GALSBI-SPS (red contours) models for the i ≤ 22 cut. The lower left contours refer to the magnitude and size distributions, while the upper right to the colour distributions. MAG_AUTO refers to the galaxy magnitudes as measured by SOURCE EXTRACTOR, FLUX_RADIUS to the galaxy effective radii, while the colours are obtained as the difference between the SOURCE EXTRACTOR magnitudes. In the figure we report the GALSBI contours obtained from each posterior distribution sample in Fischbacher et al. (2025a).

Fig. 12. Comparison of magnitude, size, and colour distributions between the observed HSC galaxies (black contours) and those from GALSBI (blue contours) and GALSBI-SPS (red contours) models for the i ≤ 23 cut. The lower left contours refer to the magnitude and size distributions, while the upper right to the colour distributions. MAG_AUTO refers to the galaxy magnitudes as measured by SOURCE EXTRACTOR, FLUX_RADIUS to the galaxy effective radii, while the colours are obtained as the difference between the SOURCE EXTRACTOR magnitudes. In the figure we report the GALSBI contours obtained from each posterior distribution sample in Fischbacher et al. (2025a).

Figure 10 shows that, for the brightest galaxy sample (i ≤ 21), the g, r, i, z, y magnitude distributions exhibit very good agreement among the three samples. This is quantified in Table 10, which reports the median, 50th-16th, and 84th-50th percentile values for each distributions. At this magnitude cut, the GALSBI-SPS median magnitudes differ from the observed data by at most ∼0.1 mag in the g and y-bands, and remain below this threshold in all other wavebands. GALSBI shows sub-0.1 mag median magnitude agreement across all bands, except for the g-band, although this difference is not visually noticeable in the magnitude histogram. In the magnitude-size planes (bottom-left row of Fig. 10), both GALSBI-SPS and GALSBI captures the bulk of the observed magnitude-size distribution. However, GALSBI-SPS tends to generate a population of galaxies with larger median apparent size, whereas GALSBI achieves a notably better match, likely due to being directly constrained by HSC data through SBI. This improved agreement is reflected in the median size differences, which is of only 0.03″ for GALSBI, compared to 0.24″ for GALSBI-SPS.

The i ≤ 22 galaxy sample (Fig. 11) displays similar trends. Magnitude distributions remain consistent between observations and simulations, with all median differences between GALSBI-SPS and the data remaining below 0.1 mag. The median size discrepancy for GALSBI-SPS is slightly reduced to 0.19″, but still significant. In contrast, GALSBI shows a ∼0.03 mag agreement with the observations in the median r, i, z, y bands, while the g-band median magnitude differs from the observed value by an amount (0.14 mag) that is larger than the GALSBI-SPS case. Also at this magnitude cut, GALSBI size distribution perfectly agrees with the one from the data.

For the fainter i ≤ 23 sample (Fig. 12), the g, i, z, and y-band median magnitudes from GALSBI-SPS continue to show sub-0.1 mag agreement with the data. The r-band, instead, increased its median difference, with respect to the previous i-band magnitude cuts, to roughly 0.11 mag. As for GALSBI, the median magnitudes in the r, i, z, and y bands differ less than 0.04 mag, while the g-band median magnitude differs by 0.13 mag with respect to that from the data. It is also worth noticing that while the median magnitudes in the g-band show a better agreement with data for GALSBI-SPS, the g-band histograms in the top-left box of Fig. 12 show a distribution for GALSBI-SPS galaxies that is skewed towards fainter values with respect to that of GALSBI and the data. The inclusion of fainter galaxies, which dominate the number counts and tend to have smaller apparent sizes, helps reduce the median size mismatch between GALSBI-SPS and the data down to 0.16″. GALSBI shows instead a very good agreement with data since the model has been constrained against observations by minimising distance metrics on the ≤23 galaxy sample.

We complete the photometric comparison among GALSBI-SPS, GALSBI and the observations by considering the g − r, r − i, i − z, z − y colour distributions. This comparison is very informative as colours are a proxy of galaxy SEDs and of the modelling choices that go into creating them. The comparisons for the i ≤ 21, 22, 23 galaxy samples are reported in the upper right panels of Figs. 10, 11, and 12, while the comparison for the fainter samples are reported in Appendix D. The median, 50th-16th, and 84th-50th percentile values for each colour distribution, are reported in Table 11. The three samples span a very similar broad band colour space at all i-band magnitude cuts. For the i ≤ 21 sample (Fig. 10), the difference between the GALSBI-SPS and the observed median colours is ≲0.7 mag for the r − i, i − z, z − y colours, while the GALSBI-SPS g − r median colour differs by 0.11 mag with respect to the observed value. As for GALSBI, the median values of the colours are more consistent with the observed ones than what GALSBI-SPS shows. While the 2D distributions show a similar colour space coverage among the three samples, the 1D colour distributions show a more marked visual agreement for GALSBI rather than for GALSBI-SPS with the data. Both GALSBI and the data show a more pronounced bimodal g − r colour distribution than GALSBI-SPS. While GALSBI-SPS tends to produce a larger fraction of bluer galaxies with respect to the observations, thereby leading to a lower g − r median colour, GALSBI has the opposite trend, producing an excess of redder galaxies with respect to the observations. The r − i colour distribution shows a similar trend for GALSBI-SPS, GALSBI and the data, while for the i − z colour, GALSBI-SPS tends to produce a population of bluer objects compared to GALSBI and the data. The z − y colour distribution is similar among the three samples instead.

For the i ≤ 22 sample (Fig. 11), we observe little evolution in the median colours for GALSBI-SPS with respect to the i ≤ 21 case. For GALSBI and the data, the r − i, i − z, z − y colours also show little evolution in their median values, while the g − r median colour becomes bluer than the i ≤ 21 sample. Notably the difference in the g − r median colour for GALSBI-SPS and the data is of the order of 0.02 mag, while for GALSBI this difference is 0.08 mag. The other three colours show differences in the median values with respect to the observations of ∼0.06 mag for GALSBI-SPS and 0.02 mag for GALSBI. The colour space spanned by three samples is similar albeit with some differences with respect to the i ≤ 21 case. In the g − r colour, observed data are dominated by bluer colours, as we would expect for a sample containing higher redshift objects. The bimodal distribution predicted by GALSBI shows instead a similar number of objects with bluer or redder g − r colours, while GALSBI-SPS predicts a distribution that is very similar to the previous case. Similar considerations to the i ≤ 21 case apply for the r − i and z − y colours, while the i − z colour distribution for GALSBI and the data shows an excess of redder objects with respect to GALSBI-SPS.

In the i ≤ 23 sample (Fig. 12), the observed g − r colour distribution becomes increasingly bluer, whereas GALSBI-SPS does not show the same trend, maintaining a roughly constant median of ∼0.89 mag. This leads to a discrepancy of 0.14 mag between GALSBI-SPS and the data, while for the same colour GALSBI median differs by 0.08 mag. In this magnitude cut, GALSBI has an excess of galaxies with red g − r colours with respect to the observations, while the distribution of blue g − r colours is similar. The difference of the median r − i, i − z, z − y colours between GALSBI-SPS and the data is smaller than the previous case, of the order of 0.05 mag, while for GALSBI is of the order of 0.02 mag. The colour space spanned by three samples is similar, with the r − i, i − z, z − y 1D colour distributions show a better agreement among the three samples than the previous magnitude cut.

To summarise, GALSBI-SPS produces a fairly realistic galaxy population up to i ≤ 23. It yields a percent-level agreement with the HSC data in the number of detected galaxies up to a i-band cut of i ≤ 23. Beyond this threshold, the model is uninformed and the number counts differ up to 20% at i ≤ 25. GALSBI-SPS magnitude distributions are in agreement with those from GALSBI and the observed HSC data at all magnitude cuts. The size distribution, instead, is skewed towards larger apparent size galaxies for GALSBI-SPS compared to GALSBI and the data, with the tension decreasing from 0.24 to 0.16″ when including fainter galaxies. GALSBI-SPS galaxy colours span a very similar colour space as their observed counterparts. However, the fraction of bluer to redder g − r colours differ both from the data and GALSBI. At increasingly fainter magnitude cuts, the g − r colours of HSC data tend to become bluer, while GALSBI-SPS g − r colours do not evolve appreciably. While GALSBI produces a more clear bimodal g − r colour distribution, the ratio of bluer to redder objects is different from that of the observations. Overall, the better agreement that GALSBI has with the observations compared to GALSBI-SPS is motivated by GALSBI model parameters being constrained against HSC data with SBI. GALSBI-SPS model parameters are instead a combination of prescriptions fitted to low-redshift observations and values taken from literature studies that use very differently selected samples of objects. A SBI run with appropriate distance metrics can jointly constrain the model parameters, mitigating the residual discrepancies we see in the current setup.

Fig. 13. Comparison of the redshift distributions for different magnitude cuts between GALSBI-SPS (red lines), GALSBI (blue lines), and the photo-z estimates with EAZY (green line) and LEPHARE (orange line) of the observed galaxies in the COSMOS field. In the figure we report the GALSBI redshift distributions obtained from each posterior sample in Fischbacher et al. (2025a) and the GALSBI-SPS redshift distributions obtained from each bootstrap sample.

5.2. Comparison of redshift distributions

In this section we compared the redshift distributions derived for simulations and observations. While redshifts in simulations are drawn from the input galaxy population model, and are by definition without error, apart from matching errors, the redshift estimates for observations are taken from the state-of-the-art COSMOS2020 photo-z catalogue (see Sect. 2.4). These estimates have been obtained using two different template fitting codes, EAZY and LEPHARE. The comparisons between the simulated and observed redshift distributions are reported in Fig. 13 for different i-band cuts (using SOURCE EXTRACTOR MAG_APER3 magnitudes as in Fischbacher et al. 2025a), while the mean redshifts, and their uncertainties, median, 50th-16th, and 84th-50th percentile values are reported in Table 12. Given that we generate a single realisation of the COSMOS field with GALSBI-SPS, the errors on the mean redshift estimates are obtained using bootstrap resampling. We repeatedly sampled with replacement from the original set of galaxies to create multiple bootstrap samples, we apply the magnitude selection, and then we computed the mean for each bootstrap sample and the standard deviation of these bootstrap means. We generate a total of 10³ bootstrap samples.

For all the magnitude-limited samples, we note that GALSBI-SPS redshift distributions visually peak at a lower redshift compared to the estimates from GALSBI, EAZY and LEPHARE. For the i ≤ 21 sample (left-most panel of Fig. 13), the median redshift of GALSBI-SPS differs by just 0.002 and 0.007 with respect to those of EAZY and LEPHARE, respectively. However, the presence of high-redshift bright objects at i ≤ 21 in GALSBI-SPS leads to a noticeable shift in the mean redshift, which is higher by 0.077 and 0.084 when compared to the EAZY and LEPHARE estimates, respectively. GALSBI shows instead a better agreement with the observations in the mean redshift, with a 0.016 and 0.009 difference for EAZY and LEPHARE, respectively.

As we move to the i ≤ 22 sample (central panel of Fig. 13), the redshift distributions of EAZY, LEPHARE, and GALSBI shifts towards higher redshifts, an effect that is less pronounced for GALSBI-SPS. GALSBI-SPS median redshift that is now 0.033 and 0.023 smaller than the EAZY and LEPHARE estimates, respectively, with the mean redshift still remaining higher than the other estimates, differing by 0.049 from EAZY and 0.058 from LEPHARE. As for GALSBI, the mean redshifts now differs by 0.033 from EAZY and 0.024 from LEPHARE.

In the i ≤ 23 sample (right-most panel of Fig. 13), it is visually clear how the i ≤ 23 GALSBI-SPS redshift distribution is skewed towards lower redshift objects than the observations. This is quantified by the difference in the median redshift, which is now of the order of 0.11 for EAZY and 0.093 for LEPHARE, with GALSBI-SPS median redshift being smaller than the other two estimates. As for the mean redshift, GALSBI-SPS estimate is larger by a value of 0.012 and 0.026 with respect to EAZY and LEPHARE estimates, respectively. This is a consistent trend across the magnitude cuts that is due to bright high redshift galaxies passing the i-band magnitude cut. GALSBI shows instead a mean redshift that is lower by 0.036 and 0.022 with respect to EAZY and LEPHARE estimates, respectively. Table 4 in Fischbacher et al. (2025a) reports the comparison of GALSBI mean redshifts with EAZY and LEPHARE estimates for fainter magnitude cuts down to i ≤ 25. In the table, the authors include also EAZY and LEPHARE mean redshift estimates obtained using the re-weighted mocks from Moser et al. (2024). Using those re-weighted values leads to GALSBI mean redshift being consistent within the errors with EAZY and LEPHARE estimates.

One notable feature is the smoothness of the simulated redshift distributions compared to real ones. As discussed in Moser et al. (2024) and Fischbacher et al. (2025c), this is due to the absence of clustering in our simulations, as well as the fact that we do not run the template fitting codes on simulated galaxies, hence leading to smoother than the COSMOS2020 distributions. This also implies that the simulated uncertainties on the mean redshift underestimate the real ones since they only include the uncertainty of the galaxy population model. Moser et al. (2024) analysed the effect of sample variance on the COSMOS field, finding a mean redshift offset between the COSMOS field and the other deep fields of roughly ∼0.015 when cut at i ∼ 23. This effect is similar in magnitude to the systematic offset between EAZY and LEPHARE.

Fig. 14. Comparison of the stellar mass-SFR and the size-stellar mass planes between the observed galaxies (green and orange contours) and GALSBI-SPS (red contours) model galaxies for the i ≤ 21 (upper panel), i ≤ 22 (central panel), and i ≤ 23 (bottom panel) cuts. The lower left contours refer to the stellar mass-SFR planes, while the upper right to size-stellar mass planes. Stellar masses and SFRs for the observations are taken from the values in the COSMOS2020 catalogue, which have been estimated using EAZY (green contours) and LEPHARE (orange contours) photo-zs. The physical size estimates are provided by the SOURCE EXTRACTOR FLUX_RADIUS parameter in the r-band converted to kiloparsec using the adopted cosmology.

5.3. Galaxy physical properties

One of the main features of GALSBI-SPS is the ability to sample galaxy physical properties from the model and provide precise estimates of the relations between these properties that can be compared against those obtained from observations. The COSMOS2020 catalogue contains stellar mass and SFR estimates obtained using EAZY and LEPHARE photo-z estimates. In particular, the LEPHARE photo-zs are used in combination with Bruzual & Charlot (2003) templates and a delayed tau model for the SFH, while EAZY photo-zs are used to fit galaxy SEDs with FSPS. In the following, we refer to the stellar masses and the SFRs obtained from these fitting procedures as EAZY and LEPHARE estimates of said properties. We can compare the distribution of stellar masses and SFRs against those from the GALSBI-SPS model for the detected population of galaxies, allowing us to evaluate how realistic is the distribution of simulated galaxies physical properties we sampled from our model. In particular, we qualitatively compared the stellar mass-SFR and the size-mass planes, two important proxies of the galaxy population in a survey. A quantitative comparison is left to a future work where the model parameters are jointly constrained against observations with appropriate distance metrics and where the stellar mass and SFR are consistently estimated across the three datasets.

Figure 14 shows the stellar mass-SFR (lower left contours) and the size-stellar mass (upper right contours) planes for the same three i-band magnitude cuts adopted in the previous sections. The figure shows how GALSBI-SPS detected galaxies are able to reproduce both the star-forming main sequence (Brinchmann et al. 2004; Noeske et al. 2007; Daddi et al. 2007; Salim et al. 2007; Whitaker et al. 2012; Popesso et al. 2023) and the passive cloud in the stellar mass-SFR plane, overlapping with the distributions from the observations. The stellar mass distribution of GALSBI-SPS galaxies seems to agree more with the EAZY estimates than with the LEPHARE ones, with the latter two being discrepant due to the different modelling choices. The figure shows also that the simulations tend to produce objects with an overall higher SFR distribution compared to observations in the i ≤ 21 cut. The trend reverses instead for the fainter samples. This trend could explain what we see in the colour and redshift distribution plots. The i ≤ 21 sample contains higher SFR simulated objects, which tends to be brighter, bluer, and visible at higher redshifts, thereby producing the high redshift tail in the GALSBI-SPS i ≤ 21 redshift distribution. At fainter i-band cuts, the observations have more objects with higher SFR, leading to bluer colours than those predicted by GALSBI-SPS.

The upper right panels of Fig. 14 compares the size-stellar mass planes for the three i-band cuts. The physical size estimates are provided by the SOURCE EXTRACTOR FLUX_RADIUS parameter in the r-band converted to kpc using the adopted cosmology. The sizes have been consistently measured on data and simulations with the same SOURCE EXTRACTOR setup. The figure shows very good qualitative agreement at all i-band cuts. We also observe a similar trend as in Figs. 10, 11, and 12, where GALSBI-SPS tends to produce a size distribution that is skewed towards larger values with respect to the observations. We also observe, as in Fig. 12, that the distribution of physical sizes tend to be more in agreement in the faint sample rather than in the bright sample.

6. Discussion

The results of this work demonstrate the ability of GALSBI-SPS to provide a realistic and robust representation of the galaxy population up to a depth of i ≤ 23. Within this range, the model successfully reproduces the observed galaxy number counts and captures the multidimensional distribution of magnitudes, sizes, colours, and physical properties of observed galaxies. This level of agreement is particularly noteworthy, considering that the model parameters were taken from different literature sources, each based on a different sample selections.

However, there are some tensions with the observations. Apparent galaxy sizes in arcsec are systematically larger than the observed counterparts, with median offsets of by 0.24″, 0.19″, 0.16″ at i ≤ 21,  i ≤ 22,  i ≤ 23, respectively. Furthermore, at i ≤ 21, GALSBI-SPS under-predicts the fraction of redder objects in the bimodal g − r colour distribution with respect to observations, whereas the phenomenological model (GALSBI) over-predicts the number of this population. At fainter cuts, i ≤ 22,  i ≤ 23, GALSBI-SPS still under-predicts the fraction of redder objects in the bimodal g − r colour, while also producing an overall distribution of g − r colours that is redder than what the observations suggest. GALSBI-SPS redshift distributions tend to have a mean redshift that is higher than what the COSMOS2020 photo-zs seem to suggest. Furthermore, the galaxy population in GALSBI-SPS has a higher average SFR at i ≤ 21 and a lower average SFR at i ≤ 23 compared to the values quoted in the COSMOS2020 catalogue.

These discrepancies could be due to the estimated values for the mean and scatter of the SFH shape parameter distributions and to the fact that we modelled these complex distributions from GAMA and DEVILS data using a single truncated multivariate normal distribution. This choice is motivated by the necessity to not dramatically increase the number of model parameters. We also assume linear redshift and stellar mass evolutions of the SFH shape parameters. While this modelling captures the bulk of the population, the parameter values used thus far are not able to reproduce bimodal or strongly non-Gaussian features in the SFH parameter space. As a result, the predicted SFR distribution at i ≤ 21 is skewed towards higher values than what is reported in the COSMOS2020 catalogue. This excess star formation makes galaxies appear brighter and therefore detectable at higher redshifts, leading to a larger redshift distribution mean compared to the observations. At fainter limits (i ≤ 22,  i ≤ 23), the redshift evolution of the SFH parameters leads to a SFR distribution that has lower values relative to the data. This under-prediction of strongly star-forming galaxies results in a lack of evolution in the g − r colour distribution, inconsistent with observational trends.

A key avenue for future improvements is to enforce consistency in galaxy formation histories, ensuring, for example, that the stellar mass assembled by galaxies at z = 1 reflects the stellar mass distribution of the z = 0 progenitors at that epoch. This is a non-trivial task, as galaxies do not evolve through smooth SFHs alone, but are affected by mergers, feedback, and environmental processes. While unimodal SFHs, like those adopted in this work, are effective at capturing the long-term star formation trends of most galaxies, many of which are dominated by a single formation epoch, they struggle to reproduce short-lived bursts of recent star formation. These bursts are particularly important for modelling galaxies that lie above the star-forming main sequence, but are notoriously hard to model and constrain. A more flexible SFH model that incorporates a burst component, such as those available in PROSsc pect, could help disentangle long-term trends from recent star formation episodes (e.g. Ciesla et al. 2015; Małek et al. 2018). Indeed, as shown in Fig. 14, GALSBI-SPS does not reproduce the population of highly star-forming galaxies above the main sequence at stellar masses of ∼10^{9 − 10.5} M_⊙, a shortcoming that could likely be addressed by including a burst component. In this work, however, we do not implement such a component, as the SFH parameters are calibrated on best-fitting values from GAMA and DEVILS data, both of which assume unimodal SFHs. Although the lack of a burst model may partially contribute to the observed discrepancies in magnitudes and colours, it is unlikely to be the sole driver. In particular, it cannot account for the systematic offset in apparent sizes, which are consistently overestimated across all magnitude cuts. This suggests that some of the discrepancy is due to the use of parameter values obtained from observations, while in reality the intrinsic distribution of the SFH or size parameters might be different than that observed.

A second key limitation to the predictive power of GALSBI-SPS arises from the use of literature-based parameter values. As the primary goal of this work is to motivate and illustrate the modelling choices for our galaxy population description, we adopt parameters from published studies, each based on galaxy samples with different selection criteria and properties, such as redshift and star formation activity. For example, we infer SFH, metallicity, and dust attenuation parameters for blue and red galaxies in GAMA and DEVILS using a sSFR cut of log(sSFR/yr⁻¹) > −10.5. In contrast, Weaver et al. (2023) adopt a colour-based selection in the NUVrJ diagram when deriving the GSMF for the two populations. Even among studies that use sSFR-based classifications, thresholds vary. Some studies use our same threshold to separate blue and red galaxies (Kashino & Inoue 2019), while others adopt a more conservative cut at log(sSFR/yr⁻¹) > −11.5 (Bellstedt et al. 2021).

Additional discrepancies may stem from the treatment of dust attenuation and the inclusion of an AGN component, both of which primarily affect rest-frame UV and blue optical magnitudes and colours (Tortorelli et al. 2024). To assess the impact of these components, we conduct two tests. First, we consider a limiting case where no dust attenuation is applied to any galaxy. If dust were the primary cause of the observed discrepancies, we would expect a notable increase in UV and blue optical flux, leading to brighter g-band magnitudes, bluer colours, and a larger number of detections at faint flux levels. However, we find that this limiting case is not alleviating the observed tensions. Similarly, we test removing the AGN contribution to the galaxy flux. While this alters the shape of the multidimensional colour distribution by removing AGN-like colours, the trend in the g − r colour evolution remains unchanged.

Beyond refining the modelling of SFH parameter distributions, the most significant improvement in the predictive power of GALSBI-SPS is set to come from constraining the model parameters with informative observational data using SBI. A core strength of our forward-modelling framework is the ability to adjust model parameters, guided by physically motivated priors and carefully chosen distance metrics, to achieve closer agreement with observations. We have not performed this optimisation in the present work to focus on motivating the modelling choices for GALSBI-SPS. Optimising model parameters is computationally intensive due to the run-time of generating large numbers of catalogue realisations from broad prior distributions for the model parameters. While PROSsc pect is considerably more time-efficient at generating SEDs than other SPS codes, its execution speed is still not sufficient to function as a scalable generative model across the wide prior space required for SBI. This calls for the inclusion in our framework of a machine-learning-based emulator to enable fast computation of synthetic magnitudes. For this reason, we developed PROMAGE a feed-forward neural network that emulates the computation of observer- and rest-frame magnitudes from the generative galaxy SED package ProSpect. This emulator will be crucial in enabling a full-scale SBI run to optimise the model parameters. Generating full-survey image simulations for every prior sample also presents a major computational bottleneck. As shown in Fischbacher et al. (2025c), this can be mitigated by using an emulator that maps intrinsic galaxy properties directly to detected ones, bypassing image-level simulations. Such emulators are also critical for testing model extensions to GALSBI-SPS. In keeping with the forward-modelling philosophy, additional complexity should only be introduced when residual tensions with the data remain after the SBI procedure. This could include a recent burst in the SFHs, the use of a single rather than a two component (blue/red) population, or the inclusion of realistic clustering of galaxies with the SHAM-OT method (Fischbacher et al. 2025c). This last aspect could further expand the scientific reach of our framework, enabling studies of the stellar-to-halo mass relation, clustering-based cosmological parameter estimation, and the environmental dependence of galaxy physical properties.

A major strength of the forward-modelling approach lies in its ability to combine heterogeneous datasets to constrain the model parameters via SBI. Each dataset carries distinct constraining power towards different components of the model. As shown in Tortorelli et al. (2024), several SED modelling ingredients, such as dust attenuation, AGN emission, gas properties, IMF, and the SFH (as shown in this work), can significantly alter broad band galaxy colours. These affect, in turn, forward-modelling-based redshift distributions, often producing shifts larger than what is required for Stage IV surveys. To mitigate this, it is crucial to design distance metrics and use the appropriate datasets that are sensitive to these SED modelling components. For instance, HSC DUD data provide strong constraints on components that affect broadband magnitudes and colours, such as GSMF, SFH, dust attenuation, galaxy sizes and the AGN contribution. This dataset, already used in Fischbacher et al. (2025a) to calibrate the phenomenological version of GALSBI, serves as an ideal precursor for Stage IV photometric surveys, thereby providing a constrained model that is able to reproduce the colour space spanned by these future data. However, several key components of galaxy SED models have little effect on broad band photometry (Tortorelli et al. 2024). These include the current gas-phase metallicity, ionisation state of the gas, along with other stellar population features. Constraining these features requires spectroscopic data with higher spectral resolution, as provided by surveys like DESI and 4MOST. These surveys enable more detailed SED modelling, making it possible to create galaxy population models applicable to both photometric and spectroscopic datasets. Examples of valuable spectroscopic constraints from DESI include the use of the quasar sample (Chaussidon et al. 2023), which would help in refining AGN modelling in the near-UV/optical, DESI BGS (Hahn et al. 2023b) high signal-to-noise spectra, which could help in linking Balmer decrement-estimated extinction with galaxy global physical properties, and the emission line galaxy samples (Raichoor et al. 2023), which can be used to trace gas ionisation and metallicity up to high redshift and link them to the galaxy global physical properties. On the 4MOST side, programmes such as 4MOST-STePS (Iovino et al. 2023), targeting long-integration spectra of intermediate-redshift galaxies, offer further opportunities to constrain stellar population models. In order to forward-model these spectroscopic surveys and constrain the relevant parameters, it is essential to develop realistic and efficient spectral simulators that map intrinsic galaxy SEDs to observed spectra under realistic observational conditions. The USPEC (Fagioli et al. 2018, 2020) and USPEC2 (Tortorelli et al., in prep.) codes are designed for this purpose, supporting the simulation of fibre-based spectra for SDSS, DESI, and 4MOST. Finally, to make SBI over large parameter spaces feasible, it is necessary to accelerate the generation of synthetic spectra from PROSsc pect. This will require the development of a dedicated spectral emulator of PROSsc pect outputs, analogous what we have already achieved with PROMAGE.

The choice of the SBI scheme is critical for obtaining robust and unbiased posterior estimates of model parameters. Several recent studies (Tortorelli et al. 2020, 2021; Kacprzak et al. 2020; Moser et al. 2024; Fischbacher et al. 2025c) have successfully employed ABC, one of the earliest and most widely used SBI techniques, to derive posterior estimates. The use of ABC, one of the earliest SBI methods developed, is justified by its inherent conservative nature. It mitigates the risk of convergence to local minima and is particularly well suited for models with a large number of parameters, such as GALSBI and GALSBI-SPS. The robustness of ABC stems from its reliance on physically motivated priors and well-chosen distance metrics tailored to be sensitive to specific aspects of the model. The construction of informative priors plays a key role in guiding the inference process, anchoring it to observationally supported values while allowing for uncertainty through appropriate prior variances. Despite its strengths, ABC is computationally intensive. It converges more slowly to the approximate posterior than more modern machine learning–based SBI methods, such as neural density estimators (Papamakarios & Murray 2016; Lueckmann et al. 2017; Papamakarios et al. 2018; Alsing et al. 2019) or full-field inference approaches (Zeghal et al. 2025). These techniques are capable of modelling complex posteriors using fewer simulations, but may yield overconfident results, particularly when sample sizes are small or priors are misspecified (Hermans et al. 2021; Tam et al. 2022). This can lead to the premature exclusion of plausible regions of parameter space, an issue ABC avoids by its incremental acceptance process. Nevertheless, ABC is not without limitations. One of the primary challenges is the sharp drop in acceptance rate as the distance threshold ϵ approaches zero, which can render the method inefficient for fine-grained posterior estimation (Alsing et al. 2018).

A well-constrained galaxy population model, coupled with a realistic forward-modelling framework for photometric and spectroscopic surveys, enables a wide range of scientific applications. From a galaxy evolution perspective, a major strength of our approach is its ability to accurately measure the luminosity and the GSMFs for large galaxy samples, as well as providing an unbiased estimate of the intrinsic scaling relations among global galaxy physical properties. This is enabled by the simulation’s perfect knowledge of intrinsic quantities such as stellar mass, luminosity, and redshift, alongside a precise characterisation of survey completeness. Using emulators to test model extensions, we can evaluate which physical prescriptions best match the observed trends, for example the redshift evolution of gas-phase metallicity or the dependence of dust attenuation on global galaxy properties. From a cosmological standpoint, our framework offers a robust avenue for meeting the stringent redshift calibration requirements of Stage IV surveys, thereby enabling precise measurements of LSS cosmological parameters. Once the model parameters are constrained via SBI, the resulting redshift distributions can be used to re-analyse existing cosmological datasets, quantify the gains in parameter precision from improved calibration, and support future Stage IV analyses. These capabilities are made possible by the forward-modelling framework’s core feature: the ability to apply identical selection functions and analysis pipelines to both simulations and observational data.

To conclude, we have presented the first implementation of GALSBI-SPS, a comprehensive SPS–based galaxy population model designed for forward-modelling applications. This model represents a major step forward in predictive galaxy population modelling, as it consistently generates both physical and morphological galaxy properties. For this first implementation, despite adopting parameter values drawn from literature sources, each based on different datasets and redshift regimes, GALSBI-SPS offers a realistic description of the galaxy population down to i ≤ 23. The model reproduces the observed distributions of magnitudes, colours, and sizes from HSC data with good fidelity, although moderate discrepancies remain, particularly in the g − r colour and in the median apparent sizes. We benchmark GALSBI-SPS against its phenomenological counterpart, GALSBI, which has undergone extensive calibration against both broad and narrow band photometric data across multiple studies. Although GALSBI-SPS is still in an early development phase and has not yet been optimised via SBI, it achieves a comparable level of realism, demonstrating the strength of its underlying physical modelling framework. The predictive power of GALSBI-SPS is expected to increase significantly once its parameters are jointly constrained using SBI and informative photometric and spectroscopic datasets. A calibrated version of the model will provide forward-modelled redshift distributions with the accuracy required for Stage IV cosmological analyses. It will also allow for robust measurements of fundamental scaling relations, such as the luminosity function, GSMF, size–mass relation, and stellar mass–SFR relation, derived from a model where the true physical parameters and their redshift evolution are precisely known. GALSBI-SPS is publicly available as part of the existing GALSBI Python package, offering the community a unified framework for generating both physical and phenomenological galaxy population models.

7. Conclusions

The forward-modelling of photometric and spectroscopic galaxy surveys represents a robust approach to characterise the galaxy population observed in a survey and in doing so obtain an accurate estimate of galaxy redshift distributions for cosmological applications. Key to its success is the development of a realistic galaxy population model and a mapping of galaxy physical properties to their SEDs. In this work, we present the first implementation of GALSBI-SPS, a new SPS-based galaxy population model that is designed for forward-modelling applications in both galaxy evolution and cosmology. This model represents a major step forward in predictive galaxy population modelling, as it consistently generates both physical and morphological galaxy properties. The prescriptions we implement are motivated by the sensitivity analysis conducted in Tortorelli et al. (2024). The model generates realistic, survey-independent galaxy catalogues, which we use to forward-model HSC DUD observations in the COSMOS field. This work aims to introduce GALSBI-SPS, to explain the motivation for the modelling choices adopted for the physical galaxy properties and for their SEDs, and to evaluate the performance of the model with literature values for the parameters against its phenomenological version, GALSBI (Fischbacher et al. 2025a), and observed data.

The model sampled galaxy physical properties from analytical parametrisations. Galaxies are drawn from two overlapping populations of blue and red objects that share similar parametrisations but different parameter values. The sampling started by jointly drawing formed stellar masses of galaxies and their redshifts from double Schechter GSMFs, with different parameter values for blue and red objects. SFHs are modelled using the truncated skewed normal distribution (Robotham et al. 2022). The shape parameters for the SFH are sampled from modelling the distribution of best-fitting values in GAMA and DEVILS, while the overall normalisation came from forcing the integrated SFH to produce the formed stellar mass sampled from the GSMF. We then linearly mapped the stellar mass evolution of each galaxy on to the shape of the gas metallicity evolution, with the final gas-phase metallicity at the time of observation being dependent on the galaxy stellar mass, redshift and SFR averaged over the last 100 Myr. We then sampled gas ionisations, dust parameters, and velocity dispersions, and we assigned each galaxy an AGN template whose flux contribution depends on a sampled f_AGN parameter. We also sampled morphological properties for our galaxy samples, whose light distributions are modelled as single Sérsic light profiles. Galaxy physical sizes and Sérsic indices are sampled conditioned on their redshifts and stellar masses, with different analytical prescriptions for red and blue galaxies, while the ellipticities are sampled for both blue and red galaxies from the same modified Beta distribution as in Moser et al. (2024) and Fischbacher et al. (2025c). The galaxy physical properties are then passed to the generative galaxy SED package PROSsc pect, which enables self-consistent modelling of stellar, nebular, dust, and AGN emission. We then computed apparent magnitudes from PROSsc pect SEDs after adding the AGN flux contribution.

The catalogue of intrinsic galaxy properties, which is sampled from the galaxy population model, is then used to create simulated HSC DUD images in the COSMOS field with the UFIG image simulator. We simulated images in the five broad bands g, r, i, z, y following the prescriptions adopted in Moser et al. (2024) and Fischbacher et al. (2025c) and we ran SOURCE EXTRACTOR on the images to produce catalogues of detected galaxies. We consistently applied the same SOURCE EXTRACTOR setup on both real and simulated HSC images.

After we applied the same selection on both data and simulations to obtain a sample with reliable photometric information, we compared GALSBI-SPS model galaxies against those sampled from GALSBI and against the observed data using a number of distance metrics. We first compared the relative number of detected galaxies between the data and the simulations for five different i-band cuts. We find that both GALSBI and GALSBI-SPS yield a similar number of detected galaxies to the observations, within few percent, up to i ≤ 23, a threshold that is similar to a Stage III-like experiment (e.g. KiDS, Wright et al. 2025a). However, at fainter magnitudes, the relative difference between GALSBI-SPS and the data increases beyond 10%, while for GALSBI it stays at the few percent level. This different behaviour is due to two factors. The first is that GALSBI has been constrained against HSC DUD data, while GALSBI-SPS relies entirely on literature-based parameters and is not yet tuned via SBI. The second is that the model is uninformed beyond i ≤ 23 given the GAMA, DEVILS, and literature data being used. We then compared the magnitude, size, and colour distributions among the three samples up to the threshold where the relative number of detected objects agrees within few percent (i ≤ 23). The magnitude distributions show a very good agreement at all i-band thresholds, with median differences always below 0.14 mag. Model galaxies cover the full colour space spanned by HSC data.

However, some tensions are visible in the 1D g − r colour distributions. GALSBI-SPS does not show the same fraction of redder colour objects in the bimodal g − r distribution compared to the observations and GALSBI. Furthermore, the g − r colour distribution does not evolve with the i-band cut, which leads to an overall distribution of redder g − r colours at i ≤ 22,  i ≤ 23. Apparent sizes are systematically overestimated by ∼0.2″ in the median across all magnitude cuts. Furthermore, GALSBI-SPS predicts redshift distributions that have higher mean values than what the observations and GALSBI return. This difference goes from ∼0.08 at i ≤ 21 to 0.012 at i ≤ 23. The comparison of the galaxy physical properties drawn from GALSBI-SPS with those from the COSMOS2020 catalogue show the ability of GALSBI-SPS to qualitatively reproduce the distribution of galaxies in the stellar mass-SFR and size-stellar mass planes. GALSBI-SPS tends to produce a population of higher SFR galaxies at i ≤ 21 and lower SFR at i ≤ 23 compared to observed data.

In Sect. 6 we discuss the strengths and limitations of the current version of GALSBI-SPS. While it provides a realistic, robust, and survey-independent representation of the galaxy population up to a depth of i ≤ 23, the model predictive power can be significantly increased in the future. We mostly use values taken from the literature for the model parameters, which introduces inconsistencies across different physical components of the model, as these are often derived from samples with differing redshift coverage, selection criteria, and sSFR definitions. In this work, we choose to not constrain the model with SBI to focus on motivating our choices for the GALSBI-SPS modelling components. Future work will focus on constraining the model parameters with SBI and leveraging the joint constraining power of photometric and spectroscopic datasets such as HSC, LSST, DESI, and 4MOST. This improvement in the modelling would contribute to the generation of more realistic g − r colours and a bluer/redder fraction of objects. We will use the machine-learning galaxy magnitude emulator PROMAGE to accelerate the generation of large numbers of catalogue realisations from broad prior distributions for the model parameters. The use of emulators coupled with physically motivated priors and carefully chosen distance metrics will allow GALSBI-SPS model parameters to be constrained with informative photometric and spectroscopic data. These improvements will enhance the model’s ability to meet the stringent requirements on the redshift distributions for a Stage IV cosmological analysis. Furthermore, they will enable robust inference of the evolution of galaxy properties across cosmic time, such as the luminosity function, GSMF, size-stellar mass relation, and the stellar mass-SFR relation, by using a model where the true physical properties are precisely known. GALSBI-SPS is publicly released as part of the GALSBI Python package (Fischbacher et al. 2024) and provides the community with a unified framework for generating physically and phenomenologically motivated galaxy population models suitable for forward-modelling wide-field galaxy surveys.

Acknowledgments

LT is grateful to Andrew Battisti, Luke Davies, Claudia Lagos, Adam Marshall, Risa Wechsler, Sucheta Cooray for useful discussions. This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. This project was supported in part by grant 200021_192243 from the Swiss National Science Foundation. SB acknowledges funding by the Australian Research Council (ARC) Laureate Fellowship scheme (FL220100191). ASGR acknowledges funding by the Australian Research Council (ARC) Future Fellowship scheme (FT200100375). GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalogue is based on data taken from the Sloan Digital Sky Survey and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programmes including GALEX MIS, VST KiDS, VISTA VIKING, WISE, Herschel-ATLAS, GMRT and ASKAP providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA website is https://www.gama-survey.org/. DEVILS is an Australian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The DEVILS input catalogue is generated from data taken as part of the ESO VISTA-VIDEO (Jarvis et al. 2013) and UltraVISTA (McCracken et al. 2012) surveys. DEVILS is part funded via Discovery Programs by the Australian Research Council and the participating institutions. The DEVILS website is https://devilsurvey.org. The DEVILS data is hosted and provided by AAO Data Central (https://datacentral.org.au/).

References

Appendix A: GAMA completeness-selected sample

Fig. A.1. Stellar mass-redshift selection limit for the GAMA sample. At fixed stellar mass slice of 0.1 dex width, each galaxy is colour-coded by the g − i colour quantile in that slice. The black line represents the low redshift 95% extreme of the g − i distribution at all stellar masses.

Fig. A.2. Redshift distributions for the mass complete samples of GAMA (red) and DEVILS (blue) galaxies.

The Value Added Catalogue used in Bellstedt et al. (2020, 2021) contains gas, dust and stellar population properties measured with PROSsc pect for 233707 galaxies. We match this catalogue to the GKVSCIENCECATV02 catalogue (trimmed-down, science-ready version of table GVKINPUTCAT) in the GAMA DR4 database via the CATAID keyword to assign the corresponding u, g, r, i, Z, Y, J, H, Ks photometry to each object. Following the GAMA DR4 catalogue documentation, we select sources having SC ≥ 7, which is the GAMA III Main Survey 95% spectroscopic completeness limit. This cut reduces the number of objects to 195260 and it is meant to include sources that are not classified as ambiguous, that lie within the survey region and outside of the star-masked regions, and that have r < 19.65. In order to select a sample of objects that contains all galaxy types at a given stellar mass without being biased to bright, star-forming galaxies, we define the lower stellar mass range for GAMA using the spectroscopic completeness limit defined by Robotham et al. (2014). This completeness cut takes into account the fact that at higher redshifts we are able to see only the bluer galaxies at fixed stellar mass, thereby creating an unrepresentative sample of objects that is strongly biased towards blue galaxies. We therefore define a variable stellar mass redshift selection limit by computing, at fixed stellar mass slice of 0.1 dex width, the g − i colour quantile of each galaxy in that slice. Figure A.1 shows the stellar mass redshift plane. Each galaxy is colour coded by the g − i quantile it belongs to at fixed stellar mass slice, from 0% to 100%, depending on whether the galaxy is the bluest or the reddest in the stellar mass bin selected. The stellar mass completeness limit is the low redshift 95% extreme of the g − i distribution at all stellar masses (black line in Fig. A.1). The final mass complete sample of galaxies contains 101784 objects. We report in Figs. A.2 and B.2 the redshift and i-band magnitude distributions of this sample.

Appendix B: DEVILS completeness-selected sample

Fig. B.1. Stellar mass-redshift selection limit for the DEVILS sample. The black line represents stellar mass completeness selection as function of redshift defined in Thorne et al. (2022b).

Fig. B.2. Histogram of i-band magnitudes for the mass complete samples of GAMA and DEVILS galaxies. GAMA i-band magnitude refers to GAMA-reprocessed version of the KiDS-VIKING images, while DEVILS i-band magnitude refers to the HSC version.

The DEVILS data in the COSMOS field contains stellar masses, SFRs, gas metallicities, dust parameters and AGN luminosities for 494, 084 sources, of which ∼24k have spectroscopic redshifts from DEVILS, while for the remaining ones have photometric-estimated redshifts. In order to select an equally unbiased sample, we apply the stellar mass completeness limit reported in Fig. 4 of Thorne et al. (2022b). The sample of galaxies above the completeness limit (black line) is shown in Fig. B.1. This mass completed sample contains 38, 066 galaxies. We report in Figs. A.2 and B.2 the redshift and i-band magnitude distributions of this sample. GAMA galaxies are 13, 6, 4, 3, 3 times more than DEVILS galaxies at i ≤ 21, 22, 23, 24, 25, respectively.

Appendix C: Galaxy surviving stellar mass emulator

Fig. C.1. Histogram showing the prediction accuracy of the galaxy surviving stellar mass emulator. Red and blue bands represent the ranges containing 95% and 99% of the samples from the test set, respectively, while dashed and solid vertical lines represent the 95th and 99th percentile values. The trained emulator produces per-mille level accurate predictions for 99% of the test set.

The galaxy surviving stellar mass is a quantity that is more closely tied to what we estimate from the observed galaxy light, since it is only the existing stars that contribute to it. PROSsc pect is able to compute the galaxy surviving stellar mass; however, this increases the computational time of individual galaxy SEDs. We therefore trained a machine-learning-based emulator to accelerate the computation of galaxy surviving stellar masses. The emulator predicts galaxy surviving masses from their formed masses, redshifts, SFH parameters and gas-phase metallicities at the time of observation. We trained the emulator on a sample of 10⁶ objects, split in 80% for training, 10% for validation, and 10% for testing. The results in Fig. C.1 shows the prediction accuracy of the emulator. We are able to recover the input galaxy surviving stellar masses with a per-mille level accuracy on 99% of the test set.

Appendix D: Faint sample comparison

Fig. D.1. Comparison of magnitude, size, and colour distributions between the observed galaxies (black contours) and those from GALSBI (blue contours) and GALSBI-SPS (red contours) models for the i ≤ 24 (left panel) and i ≤ 25 (right panel) cuts. The lower left contours refer to the magnitude and size distributions, while the upper right to colour distributions. MAG_AUTO refers to the galaxy magnitudes as measured by SOURCE EXTRACTOR, FLUX_RADIUS to the galaxy effective radius in arcsec, while the colours are obtained as difference between the SOURCE EXTRACTOR magnitudes. We report in the figure the GALSBI contours obtained from each posterior distribution sample in Fischbacher et al. (2025a).

We report in Fig. D.1 the magnitude, size, and colour distributions for the fainter samples at i ≤ 24, 25. By going at fainter magnitudes, we are including in the sample galaxies that are probing a regime where GALSBI-SPS is uninformed by the data we use to create the model. At both i ≤ 24 and i ≤ 25, the i, z, y magnitude distributions are in very good agreement, while for the g and r-band, GALSBI-SPS magnitude distributions are marginally fainter than their observed counterparts. The samples span a similar colour space but the GALSBI-SPS g − r colour distributions produce a population of redder g − r colours than the observations. As we have seen in Fig. 14, the SFR distribution does not scale with fainter magnitude cuts as simulations do. The lack of strongly star-forming objects means redder colours than what the observations produce. Since in the i ≤ 25 galaxy sample we include more higher redshift objects whose bluer emission moves in redder wavebands, this lack of strongly star-forming objects is also marginally visibile in the r − i colour. GALSBI magnitudes, sizes, and colours are in a better agreement with the observations since the distance metrics used to match simulations to observations have been computed using samples cut at these i-band thresholds.

All Tables

All Figures

Fig. 1. Flowchart describing the main components of the GALSBI-SPS galaxy population model. Galaxy physical properties are sampled hierarchically from analytical distributions and conditional relations, as detailed in Sect. 3. The box colours reflect the hierarchy of dependencies. Blue boxes represent properties that are independently sampled, with the GSMF being the starting point of the sampling. Yellow boxes represent properties conditioned only on galaxy stellar mass and redshift (namely, the galaxy SFH, velocity dispersion, morphology, and AGN). Green boxes represent properties conditioned simultaneously on galaxy stellar mass, redshift, and SFH (galaxy metallicity history, dust attenuation, and gas ionisation). All sampled properties are then passed to PROSsc pect (red box) to generate galaxy SEDs and magnitudes. Together with the morphological properties they form the input catalogue for the forward-modelling. Each model component is described in its corresponding subsection of Sect. 3.

In the text

	Fig. 2. Redshift evolution of the characteristic stellar mass logℳ* for blue and red galaxies. Blue and red points represent the measurements from Weaver et al. (2023). The blue and red dashed lines are built by fitting the measurements with Equation 2. The blue and red bands represent the 1σ uncertainty band on the fit.
In the text

	Fig. 3. Redshift evolution of the characteristic density ϕ* for blue and red galaxies. Blue, deep sky blue, magenta, and red points represent the measurements from Weaver et al. (2023). The blue, deep sky blue, magenta, and red dashed lines are built by fitting the measurements with Equation 3. The bands represent the 1σ uncertainty band on the fit.
In the text

Fig. 4. Dependence of the SFH shape parameters on galaxy redshifts and stellar mass. The fraction of the galaxy age at which the SFH peaks, log(mpeak′), linearly depends on the galaxy stellar mass and its coefficients in turn depend on the galaxy redshift. The upper panels show the dependence of the slope, log(mpeak′)mass, evo, slope, and the intercept, log(mpeak′)mass, evo, intcpt, on redshift. Red error bars refer to quiescent galaxies from the mass-complete GAMA and DEVILS sample, while blue error bars refer to star-forming galaxies. The dashed lines represent the best-fitting linear relations for the redshift evolution of these parameters, while the red and blue bands represent the 1σ uncertainty on the fit. The lower panels instead show the dependence of log(mperiod) and mskew on the galaxy stellar mass. The error bars refer to the median and standard deviation of these quantities computed in bins of 0.1 dex in stellar mass, overlayed on top of the overall distribution. Dashed lines and colour bands represent best-fitting linear relations and 1σ uncertainty on the fit for the populations of red and blue galaxies.

In the text

Fig. 5. Correlation between log(τISM) and redshift, stellar mass and sSFR for star-forming (left panel) and quiescent (right panel) galaxies from DEVILS data. The solid red (blue) line shows the best-fitting relation that expresses this dependence for star-forming (quiescent) galaxies. The dotted lines show the observed rms scatter of the relation. Each sample contour encloses 50%, 84%, and 99% of the values.

In the text

	Fig. 6. Histogram of the ratio between the AGN bolometric luminosity and the galaxy luminosity in log-space, log(fAGN), for star-forming (blue histogram) and quiescent (red histogram) galaxies. The dashed lines represent the samples drawn from the sum of two normal distributions with which we modelled the log(fAGN) values from DEVILS data.
In the text

	Fig. 7. Best-fitting linear evolution in redshift of the size-stellar mass parameters for star-forming (left panel) and quiescent (right panel) galaxies. The scatter points represent the observed values of the parameters and are taken from Tables 2 and 3 in Nedkova et al. (2021). The dashed lines show the best-fitting relations, while the bands represent the 1σ uncertainty on the fit.
In the text

Fig. 8. Stellar mass dependence of the r-band Sérsic index n. The deep sky blue and orange contours are obtained by matching the star-forming and quiescent galaxies selected via sSFR of the GAMA mass complete sample with the morphological catalogue in Kelvin et al. (2012). The blue and red error bars represent the r-band Sérsic index median trend in stellar mass bins for star-forming and quiescent galaxies, respectively, with the size of the error bars representing the standard deviation in said bin. The dashed lines represent the evolution with stellar mass obtained from the best-fitting parameters of Equation 32 on the observed data. The bands represent the 1σ uncertainty on the best-fitting relations. Each sample contour encloses 50%, 84%, and 99% of the values.

In the text

Fig. 9. Fractional difference between the number of detected galaxies in the simulations, nsim, and in the observations, nobs, as a function of the i-band magnitude cut. Red scatter points refer to the GALSBI-SPS model, while blue box plots to the phenomenological version in Fischbacher et al. (2025a). The GALSBI-SPS error bars consider Poisson noise on the counts, while GALSBI box plots have been constructed using the posterior samples in Fischbacher et al. (2025a). The dashed black line marks the perfect agreement in the number counts with observations.

In the text

Fig. 10. Comparison of magnitude, size, and colour distributions between the observed HSC galaxies (black contours) and those from GALSBI (blue contours) and GALSBI-SPS (red contours) models for the i ≤ 21 cut. The lower left contours refer to the magnitude and size distributions, while the upper right to the colour distributions. MAG_AUTO refers to the galaxy magnitudes as measured by SOURCE EXTRACTOR, FLUX_RADIUS to the galaxy effective radii, while the colours are obtained as the difference between the SOURCE EXTRACTOR magnitudes. In the figure we report the GALSBI contours obtained from each posterior distribution sample in Fischbacher et al. (2025a).

In the text

Fig. 11. Comparison of magnitude, size, and colour distributions between the observed HSC galaxies (black contours) and those from GALSBI (blue contours) and GALSBI-SPS (red contours) models for the i ≤ 22 cut. The lower left contours refer to the magnitude and size distributions, while the upper right to the colour distributions. MAG_AUTO refers to the galaxy magnitudes as measured by SOURCE EXTRACTOR, FLUX_RADIUS to the galaxy effective radii, while the colours are obtained as the difference between the SOURCE EXTRACTOR magnitudes. In the figure we report the GALSBI contours obtained from each posterior distribution sample in Fischbacher et al. (2025a).

In the text

Fig. 12. Comparison of magnitude, size, and colour distributions between the observed HSC galaxies (black contours) and those from GALSBI (blue contours) and GALSBI-SPS (red contours) models for the i ≤ 23 cut. The lower left contours refer to the magnitude and size distributions, while the upper right to the colour distributions. MAG_AUTO refers to the galaxy magnitudes as measured by SOURCE EXTRACTOR, FLUX_RADIUS to the galaxy effective radii, while the colours are obtained as the difference between the SOURCE EXTRACTOR magnitudes. In the figure we report the GALSBI contours obtained from each posterior distribution sample in Fischbacher et al. (2025a).

In the text

Fig. 13. Comparison of the redshift distributions for different magnitude cuts between GALSBI-SPS (red lines), GALSBI (blue lines), and the photo-z estimates with EAZY (green line) and LEPHARE (orange line) of the observed galaxies in the COSMOS field. In the figure we report the GALSBI redshift distributions obtained from each posterior sample in Fischbacher et al. (2025a) and the GALSBI-SPS redshift distributions obtained from each bootstrap sample.

In the text

Fig. 14. Comparison of the stellar mass-SFR and the size-stellar mass planes between the observed galaxies (green and orange contours) and GALSBI-SPS (red contours) model galaxies for the i ≤ 21 (upper panel), i ≤ 22 (central panel), and i ≤ 23 (bottom panel) cuts. The lower left contours refer to the stellar mass-SFR planes, while the upper right to size-stellar mass planes. Stellar masses and SFRs for the observations are taken from the values in the COSMOS2020 catalogue, which have been estimated using EAZY (green contours) and LEPHARE (orange contours) photo-zs. The physical size estimates are provided by the SOURCE EXTRACTOR FLUX_RADIUS parameter in the r-band converted to kiloparsec using the adopted cosmology.

In the text

	Fig. A.1. Stellar mass-redshift selection limit for the GAMA sample. At fixed stellar mass slice of 0.1 dex width, each galaxy is colour-coded by the g − i colour quantile in that slice. The black line represents the low redshift 95% extreme of the g − i distribution at all stellar masses.
In the text

	Fig. A.2. Redshift distributions for the mass complete samples of GAMA (red) and DEVILS (blue) galaxies.
In the text

	Fig. B.1. Stellar mass-redshift selection limit for the DEVILS sample. The black line represents stellar mass completeness selection as function of redshift defined in Thorne et al. (2022b).
In the text

	Fig. B.2. Histogram of i-band magnitudes for the mass complete samples of GAMA and DEVILS galaxies. GAMA i-band magnitude refers to GAMA-reprocessed version of the KiDS-VIKING images, while DEVILS i-band magnitude refers to the HSC version.
In the text

	Fig. C.1. Histogram showing the prediction accuracy of the galaxy surviving stellar mass emulator. Red and blue bands represent the ranges containing 95% and 99% of the samples from the test set, respectively, while dashed and solid vertical lines represent the 95th and 99th percentile values. The trained emulator produces per-mille level accurate predictions for 99% of the test set.
In the text

Fig. D.1. Comparison of magnitude, size, and colour distributions between the observed galaxies (black contours) and those from GALSBI (blue contours) and GALSBI-SPS (red contours) models for the i ≤ 24 (left panel) and i ≤ 25 (right panel) cuts. The lower left contours refer to the magnitude and size distributions, while the upper right to colour distributions. MAG_AUTO refers to the galaxy magnitudes as measured by SOURCE EXTRACTOR, FLUX_RADIUS to the galaxy effective radius in arcsec, while the colours are obtained as difference between the SOURCE EXTRACTOR magnitudes. We report in the figure the GALSBI contours obtained from each posterior distribution sample in Fischbacher et al. (2025a).

In the text