EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 707 (March 2026) A&A, 707 (2026) A231 Full HTML
Open Access
Issue
A&A
Volume 707, March 2026
Article Number A231
Number of page(s) 23
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202556464
Published online 17 March 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1. Introduction

Reconstructing the star formation histories (SFHs) of galaxies is a critical step toward deciphering galaxy evolution across cosmic time. An accurate reconstruction of a galaxy’s SFH enables reliable estimates of its total stellar mass and star formation rate (SFR). Of course, there are multiple factors that regulate the process of star formation and influence the overall shape and amplitude of the SFH. First, environmental factors such as galaxy mergers or intense disk instabilities can induce episodes of star formation and active galactic nuclei (AGN) activity in the core (Barnes & Hernquist 1991, 1996; Springel et al. 2005; Snyder et al. 2011; Dekel & Burkert 2014; Zolotov et al. 2015; Wellons et al. 2015; Lapiner et al. 2023). Second, feedback mechanisms, such as supernova explosions and AGNs, can prevent the hot gas in a galaxy’s halo from cooling down (e.g., Crenshaw et al. 2003; Best et al. 2005; Croton et al. 2006; Fabian 2012; Cheung et al. 2016; Harrison 2017; Barišić et al. 2017; Terrazas et al. 2017; Henden et al. 2018; Semenov et al. 2021; Dome et al. 2024). Finally, dynamical features such as spiral arms and stellar bars can regulate the process of star formation locally (Elmegreen 2011; Baba et al. 2013; Dobbs & Baba 2014; Sellwood & Carlberg 2019; Pettitt et al. 2020; Shin et al. 2023; Durán-Camacho & Duarte-Cabral 2025). Therefore, the study of SFHs not only sheds light on when and how galaxies formed their stars but also helps link their past to present-day properties, such as their morphology, stellar populations, and chemical composition.

The most widely used method for inferring the SFHs of individual galaxies is modeling their observed spectral energy distribution (SED). This approach usually involves a set of multiwavelength photometric measurements that is compared to a library of theoretical and empirical models of galaxy evolution (see reviews by Walcher et al. 2011; Conroy 2013; Pacifici et al. 2023). Ideally, broadband photometry should be combined with spectroscopic information to reduce the uncertainties in the SFHs of galaxies, particularly at earlier look-back times (e.g., Worthey et al. 1994; Bruzual & Charlot 2003; Trager et al. 2000; Gallazzi & Bell 2009; Salim et al. 2016, 2018; Chauke et al. 2018; Johnson et al. 2021; Thorne et al. 2021; Tacchella et al. 2022a; Nersesian et al. 2024; Iglesias-Navarro et al. 2024; Csizi et al. 2024; Siudek et al. 2024; Kaushal et al. 2024; Nersesian et al. 2025). Nevertheless, in the local Universe, reconstructing the stellar age distribution of individual galaxies through high-spectral-resolution spectroscopy remains a significant challenge, as most galaxies are older than 5 Gyr (Gallazzi et al. 2005), resulting in stellar spectra that are often very similar in terms of absorption features.

One of the assumptions that goes into SED fitting is the choice of the prior function for the SFH. This choice can significantly influence the reliability and precision of the derived physical properties of galaxies (e.g., Carnall et al. 2019; Leja et al. 2019a; Siudek et al. 2024). Most SED fitting frameworks adopt simplified parametric SFH models to describe the evolution of SFR with time. These models typically include parameters that control normalization, amplitudes, and characteristic timescales (e.g., the e-folding time). Common examples include exponentially declining or rising SFHs, lognormal forms, and delayed-τ models (Papovich et al. 2011; Diemer et al. 2017; Boquien et al. 2019). More flexible parametric models also exist, allowing for multiple episodes of star formation and quenching (Ciesla et al. 2017; Morishita et al. 2019; Abdurro’uf et al. 2021).

While these models are computationally less intensive, their fixed functional forms often fail to capture the expected diversity of SFHs. Numerous studies have demonstrated that this limitation can introduce systematic biases in the inferred galaxy properties (e.g., Conroy 2013; Diemer et al. 2017; Carnall et al. 2019; González Delgado et al. 2021; Suess et al. 2022; Whitler et al. 2023; Pacifici et al. 2023). To overcome these challenges, nonparametric SFH models have been developed. These flexible models do not impose a specific functional form and are designed to capture the stochastic nature of SFHs (e.g., Ocvirk et al. 2006; Leja et al. 2017, 2019a; Iyer & Gawiser 2017; Iyer et al. 2019; Euclid Collaboration: Corcho-Caballero et al. 2026).

A drawback of the nonparametric SFH models is that they are more computationally demanding compared to parametric SFHs and need broader data coverage. Nevertheless, nonparametric models have the potential to accurately reconstruct the complex SFHs of galaxies, while yielding more reliable estimates of their physical properties. Both Iyer et al. (2018) and Lower et al. (2020) validated the galaxy properties inferred from nonparametric SFH models by benchmarking them against mock observations derived from semi-analytic models and cosmological simulations. Moreover, Leja et al. (2019a) assessed the effectiveness of nonparametric SFHs using the SED fitting framework Prospector (Leja et al. 2017; Johnson et al. 2021), showing that these models yield more accurate results than their parametric counterparts. In a follow-up study, Leja et al. (2019b) found that galaxies in the 3D-HST catalog (Skelton et al. 2014) appear systematically older and more massive when using nonparametric SFHs instead of parametric models, consistent with theoretical expectations (e.g., van Dokkum 2008; Davé 2008; Davé et al. 2011; Lilly et al. 2013). The results from these studies strongly suggest that nonparametric priors have a clear advantage over traditional parametric SFHs, enabling more reliable results. However, a notable challenge still remains in the dependence of SFHs on the choice of physically motivated prior assumptions, as highlighted by Suess et al. (2022), Tacchella et al. (2022b), and Whitler et al. (2023).

Another potential issue in reconstructing SFHs arises from the modeling and fitting of SEDs from integrated photometry, which is susceptible to the outshining effect (Sawicki & Yee 1998; Papovich et al. 2001; Shapley et al. 2001; Daddi et al. 2004; Fontana et al. 2004; Pozzetti et al. 2007; Trager et al. 2008; Zibetti et al. 2009; Maraston et al. 2010; Graves & Faber 2010; Pforr et al. 2012; Conroy 2013; Sorba & Sawicki 2015, 2018; Suess et al. 2022; Whitler et al. 2023; Giménez-Arteaga et al. 2023, 2024; Narayanan et al. 2024; Jain et al. 2024). Outshining can occur when a ultraviolet (UV) bright young stellar population (younger than 10 Myr) dominates the light of a galaxy and obscures the light of older stars (older than 200 Myr). This effect can significantly bias the inferred stellar masses and SFHs of galaxies (e.g., Zibetti et al. 2009; Pforr et al. 2012; Sorba & Sawicki 2015; Jain et al. 2024). While single-aperture spectroscopy can mitigate this bias to some extent, with the inclusion of spectral features, many datasets rely solely on broad- or narrowband photometry.

A more robust alternative involves reconstructing SFHs at spatially resolved scales. Unlike integrated photometry, which captures only global averages, resolved analyses can trace local variations in stellar populations and dust content, enabling more accurate estimates of stellar mass, SFR, and formation timescales, such as the time by which 50% of a galaxy’s stellar mass was formed (Mosleh et al. 2020; Suess et al. 2019; Abdurro’uf et al. 2022a; Jain et al. 2024; Mosleh et al. 2025). Studies such as Giménez-Arteaga et al. (2023, 2024), and Jain et al. (2024) demonstrate that accounting for spatial structure reduces the impact of outshining, resulting in more accurate stellar mass measurements. The results by Jain et al. (2024) also motivate the use of flexible SFH models on global scales, a result particularly important at higher redshifts, where SFHs are often more bursty (Tacchella et al. 2020; Narayanan et al. 2024) and thus more prone to outshining effects.

The advent of the Euclid space telescope (Laureijs et al. 2011; Euclid Collaboration: Mellier et al. 2025) offers an unprecedented view of the local Universe (e.g., Hunt et al. 2025; Saifollahi et al. 2025; Cuillandre et al. 2025), providing high spatial resolution and sensitivity in one optical band (Euclid VIS instrument; Euclid Collaboration: Cropper et al. 2025), and three near-infrared bands (Euclid NISP instrument; Euclid Collaboration: Jahnke et al. 2025). By combining Euclid observations with ancillary imaging data from ground-based and space observatories, and modeling them through spatially resolved SED fitting, we can thoroughly analyze local galaxies, their SFHs, and their physical properties at kiloparsec and sub-kiloparsec scales.

However, before applying any SED fitting methods to Euclid observations, we first leverage cosmological hydrodynamical simulations. These simulations serve as powerful tools for exploring the origin and evolution of galaxies, while providing crucial benchmarks for testing and refining SED fitting techniques. In a series of Euclid papers, we employed two complementary approaches to infer the resolved stellar properties of simulated galaxies at z = 0 using synthetic Euclid observations. Euclid Collaboration: Kovačić et al. (2025) applied machine learning techniques to map the distribution of physical parameters within well-resolved simulated galaxies. In addition, Euclid Collaboration: Abdurro’uf et al. (2025) tested a pipeline for spatially resolved SED fitting of local galaxies with the piXedfit package (Abdurro’uf et al. 2021), exploring different stellar population synthesis (SPS) models and incorporating physically motivated priors, such as mass–age and mass–metallicity relations.

In this study, we also employed SED fitting with the goal of reconstructing SFHs on both global and spatially resolved scales for a sample of 25 TNG50-SKIRT Atlas galaxies (Baes et al. 2024). We analyze the SFHs derived using two methods: (1) pixel-by-pixel SED fitting and (2) integrated broadband photometry, employing both parametric and nonparametric models. The goal is to assess the feasibility of reconstructing accurate SFHs from Euclid-like data. Specifically, we tested whether we can improve upon the standard approach to SFH reconstruction, which typically involves applying either resolved rather than global SED fitting, or nonparametric rather than parametric SFHs. While each of these improvements is commonly applied separately, our approach uniquely combines both, leveraging their joint strengths to provide a more comprehensive and robust analysis of galaxy SFHs. Another key objective of this paper is to demonstrate the added value of the Euclid bands in constraining spatially resolved physical properties. This discussion complements the main analysis by showcasing the direct scientific benefits of incorporating Euclid data, which are highly relevant to the broader goals of the Euclid Collaboration.

This paper is organized as follows. In Sect. 2, we provide a brief overview of the datasets and the methodology adopted to prepare the synthetic images from TNG50-SKIRT Atlas for the pixel-by-pixel SED fitting. In Sect. 3, we describe the construction of our physical model within Prospector and provide definitions of the various SFH models used in our analysis. In Sect. 4, we compare the inferred physical properties from our fiducial SED fitting-run to the ground truth properties from the simulations. The main analysis of this paper is in Sect. 5, where we discuss the implications of how the assumed SFH models and the spatial resolution affect the inferred physical properties of the galaxies. We summarize our results in Sect. 6.

2. Simulations, sample selection, and image processing

In this section, we present a short overview of the TNG50 cosmological hydrodynamical simulation (Pillepich et al. 2019; Nelson et al. 2019) and the TNG50-SKIRT Atlas (Sect. 2.1), from which we drew our working sample of galaxies (Sect. 2.2). The TNG50-SKIRT Atlas provides synthetic images of 1154 galaxies, generated using the radiative transfer code SKIRT (Camps & Baes 2015, 2020). We also give a brief description of the image processing (Sect. 2.3) performed with piXedfit (Abdurro’uf et al. 2021; Abdurro’uf et al. 2022b), to create realistic representations that mimic the observational data.

2.1. Synthetic images of TNG50 galaxies with SKIRT

To perform the pixel-by-pixel SED fitting analysis, we used galaxies from the TNG50 simulation (Pillepich et al. 2019; Nelson et al. 2019), the highest-resolution version of the IllustrisTNG cosmological magnetohydrodynamical suite (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2018; Springel et al. 2018). With a baryonic mass resolution of 8.5 × 104 M⊙ and spatial resolution of 70–140 pc in star-forming regions, TNG50 combines high resolution with a statistically meaningful galaxy population in a cosmological volume, making it well suited for spatially resolved studies.

A sample of 1154 galaxies from the z = 0 snapshot of the TNG50 simulation was defined by Baes et al. (2024), forming the TNG50-SKIRT Atlas database. The database features high-resolution, panchromatic images of these galaxies, generated using the radiative transfer code SKIRT (Camps & Baes 2015, 2020), a 3D Monte Carlo radiative transfer code designed to model the interaction of the stellar radiation field with dust, in different astrophysical systems. Every galaxy in TNG50 is represented by a collection of stellar and gas particles. The stellar particles serve as the primary source of radiation, and depending on their mass, age, and metallicity, SKIRT assigns an SED to each one of them. Stellar particles older than 10 Myr are assigned an SED from the Bruzual & Charlot (2003) SPS models, whereas stellar particles younger than 10 Myr are assigned an SED from the MAPPINGS III templates for star-forming regions (Groves et al. 2008).

To account for the effects of diffuse dust (absorption and scattering) on the stellar radiation, a constant dust-to-metal fraction (fdust = 0.2) was assumed (Trčka et al. 2022), based on observational results in nearby galaxies (De Vis et al. 2019; Galliano et al. 2021; Zabel et al. 2021). The dust properties were modeled using The Heterogeneous Evolution Model for Interstellar Solids (THEMIS; Jones et al. 2017) framework. The main goal of the TNG50-SKIRT Atlas was to produce synthetic images spanning the UV to near-infrared (NIR) wavelengths. Consequently, the modeling of dust emission was not included in the radiative transfer analysis.

The TNG50-SKIRT Atlas originally included synthetic images in 18 wavebands, covering the UV to NIR spectrum (0.15–4.6 μm): GALEX far-UV (FUV) and near-UV (NUV) bands, Johnson UBVRI, LSST ugrizy, 2MASS JHKs, and WISE W1 and W2. Euclid Collaboration: Kovačić et al. (2025) extended this dataset by adding synthetic images for the four Euclid broadband filters (IE, YE, JE, HE). Each galaxy in the TNG50-SKIRT Atlas was simulated from five different viewing angles, resulting in 22-band images for each projection, with a consistent field of view of 160 kpc × 160 kpc. These synthetic images are in units of megajanskys per steradian and have a physical pixel size of 100 pc. For the pixel-by-pixel SED fitting analysis, only a subset of these bands is utilized: GALEX FUV and NUV, LSST u, g, r, i, z, and Euclid IE, YE, JE, and HE. Although LSST observations from the Rubin Observatory (Ivezić et al. 2019) are not yet available, other large ground-based campaigns are being carried out to complement Euclid’s observations and can be used to perform spatially resolved SED fitting on nearby galaxies. In particular, the Ultraviolet Near Infrared Optical Northern Survey (UNIONS; Gwyn et al., in prep.) in the Northern Hemisphere provides observations in the u, g, r, i, and z bands, while the Dark Energy Survey (DES; Abbott et al. 2021) in the Southern Hemisphere provides observations in the g, r, i, and z bands.

For this work, we used a subset of 11 bands spanning the UV to NIR (GALEX FUV/NUV, LSST ugriz, and Euclid IE, YE, JE, and HE), matching expected multiwavelength coverage from ongoing and upcoming surveys. Alongside synthetic images, the TNG50-SKIRT Atlas provides maps of key physical properties such as stellar mass surface density (Σ★), stellar metallicity (Z★), mass-weighted stellar age (t★), dust mass surface density maps (Σdust), and SFR surface density (ΣSFR). These maps serve as ground-truth references for testing our SED fitting.

2.2. Sample of simulated galaxies

Table 1.

Free parameters and their associated prior distribution functions in the Prospector physical model.

The simulated galaxies of TNG50-SKIRT Atlas span a wide range of stellar mass (M★), from 109.8 M⊙ to 1012 M⊙. Figure 1 shows the star-forming sequence (SFS; Strateva et al. 2001; Baldry et al. 2004; Daddi et al. 2007; Noeske et al. 2007; Elbaz et al. 2007; Whitaker et al. 2012), that is the relationship between SFR and M★, as obtained from the simulations. Galaxies can be separated into two broad galaxy populations: quiescent and star-forming. The two most common criteria to separate these two classes are based on (1) their locus on color-color diagrams (e.g., Labbé et al. 2005; Wuyts et al. 2007; Williams et al. 2009; Arnouts et al. 2013; Schawinski et al. 2014; Euclid Collaboration: Enia et al. 2026), and (2) their specific SFR (sSFR = SFR/M★). In this work, we choose to define quiescence based on the sSFR as a more robust method. We define the boundary between quiescent and star-forming galaxies based on sSFR, adopting a threshold of log 10 ( sSFR / yr − 1 ) = − 11 Mathematical equation: $ {\log_{10}} \left(\mathrm{sSFR} / \mathrm{yr}^{-1} \right) = -11 $ (Brinchmann et al. 2004; Fontanot et al. 2009; Cecchi et al. 2019; Donnari et al. 2019; Paspaliaris et al. 2023; Nersesian et al. 2025), which yields 285 quiescent and 869 star-forming galaxies.

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

SFR as a function of M★ for the TNG50-SKIRT Atlas galaxies at z = 0. Galaxies are separated as star-forming (blue stars) and quiescent (red points) based on their sSFR. The dashed orange line indicates this separation between star formation and quiescence (sSFR = 10−11 yr−1). The black squares indicate the selected sample of 25 galaxies used in our analysis.

From the parent sample of TNG50-SKIRT Atlas galaxies, a subset of 25 galaxies was selected by Euclid Collaboration: Abdurro’uf et al. (2025) for benchmarking spatially resolved SED-fitting. The galaxies were selected based on stellar mass and sSFR criteria, spanning a broad range of sSFRs (10−13–10−9 yr−1) and covering the same stellar mass range as the parent sample, while also representing a diversity of morphological types. Another criterion of the selection process was to prioritize galaxies with an almost face-on orientation, minimizing systematic effects caused by increased dust attenuation in edge-on views. Focusing on this small subsample enables a proof-of-concept analysis of spatially resolved SED fitting with nonparametric SFHs, while keeping the computational demands tractable. Figure 1 illustrates the location of these 25 galaxies on the SFR–M★ plane, where 21 are classified as star-forming and 4 as quiescent.

2.3. Observational effects and image processing with piXedfit

The synthetic images from the TNG50-SKIRT Atlas were post-processed to include realistic observational effects before applying pixel-by-pixel SED fitting. Observational noise (sky background and photon shot noise) was added following prescriptions for the Euclid Wide Survey and LSST (Euclid Collaboration: Martinet et al. 2019; Euclid Collaboration: Merlin et al. 2023), as implemented by Euclid Collaboration: Abdurro’uf et al. (2025). The adopted 10σ (S/N = 10) limiting magnitudes within a 2″ aperture for the IE, YE, JE, HE, and LSST u, g, r, i, and z bands are 24.6, 23.0, 23.0, 23.0, 23.6, 24.5, 23.9, 23.6, and 23.4, respectively (see Table 1 of Euclid Collaboration: Abdurro’uf et al. 2025). The limiting magnitudes of the synthetic images quoted here closely match the magnitude limits reported in Euclid’s Quick Data Release (Q1; e.g., Euclid Collaboration: Enia et al. 2026; Euclid Collaboration: Tucci et al. 2026). For the UV bands, background noise was added based on GALEX MIS observations (Bianchi et al. 2014), which provide deeper coverage than AIS and thus stronger constraints on recent star formation.

Each galaxy was placed at z = 0.03, and images were resampled to the pixel scales of their respective instruments. Following noise addition, all bands were convolved to their instrument point spread functions (PSFs), and subsequently homogenized in spatial resolution and sampling using the piXedfit pipeline (Abdurro’uf et al. 2021; Abdurro’uf et al. 2022b). This procedure ensures that all multiwavelength images share the resolution of the lowest-quality dataset, in this case GALEX, with a full width half maximum (FWHM) of 5 . ″ 05 Mathematical equation: $ 5{{\overset{\prime\prime}{.}}}05 $, and a pixel size of 1 . ″ 5 Mathematical equation: $ 1{{\overset{\prime\prime}{.}}}5 $. Despite the lower resolution, GALEX provides critical UV coverage to constrain recent star formation.

To test the effects of wavelength coverage and variations in imaging datasets, Euclid Collaboration: Abdurro’uf et al. (2025) analyzed three sets of imaging data cubes: (1) a combination of GALEX, LSST, and Euclid (11 bands); (2) LSST and Euclid (9 bands); and (3) Euclid-only images. They showed that Euclid images alone can accurately recover stellar mass surface density. However, to constrain SFHs, shorter-wavelength information is critical. Therefore, for the purposes of our analysis we use the combined images of GALEX, LSST, and Euclid (11 bands). To improve pixel-level S/N, piXedfit applies adaptive binning by grouping neighboring pixels with similar SED shapes, while preserving spatial information. Across the 25 selected galaxies, this yields 6958 pixel bins used in our analysis (see Fig. 2 for an example). We refer to Euclid Collaboration: Abdurro’uf et al. (2025) for further details on the procedure.

Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Example of a simulated galaxy processed with piXedfit. We show the synthetic images of GALEX, LSST, and Euclid for TNG275545 with the O1 orientation index. Observational effects were introduced to the original synthetic images from the TNG50-SKIRT Atlas database in the form of simulated noise and convolution with the PSF of the corresponding cameras. The FWHM of the PSF as well as the pixel size of the images are 5 . ″ 05 Mathematical equation: $ 5{{\overset{\prime\prime}{.}}}05 $ and 1 . ″ 5 Mathematical equation: $ 1{{\overset{\prime\prime}{.}}}5 $, respectively. The pixels are binned with piXedfit to increase the S/N of the spatially resolved SEDs.

3. SED fitting with Prospector on resolved and global scales

We applied different physical models to reconstruct the SFHs of the 25 TNG50-SKIRT Atlas galaxies, on both global and resolved scales. In this section, we provide a description of those physical models generated with the Prospector inference framework (Leja et al. 2017; Johnson et al. 2021), to fit the photometric SEDs on spatial scales and global scales.

3.1. Standard physical model within Prospector

Prospector is an SED fitting code that utilizes Bayesian forward modeling and Markov chain Monte Carlo (MCMC) sampling to explore the parameter space. Prospector can generate gridless, “on-the-fly” SEDs by combining models of stellar, nebular, and dust components into composite stellar populations. Within Prospector, there are numerous methods for treating the SFH of a galaxy. One of Prospector’s key strengths is its ability to employ nonparametric SFHs with various prior distributions and parameterizations.

One of our main goals is to determine which approach returns more meaningful results for the SFH. Since we treat the SFH in different ways, establishing a standard physical model that is consistently applied across all SED fitting runs is essential. Our physical model largely follows the model presented in Leja et al. (2017, 2018, 2019b). We also used the nested sampler dynesty (Skilling 2004, 2006; Koposov et al. 2022) to dynamically sample the parameter space and maximize a chosen objective function as the fit proceeds.

Regarding the spectra of the stellar populations, Prospector utilizes the Flexible Stellar Populations Synthesis (FSPS) code (Conroy et al. 2009). We adopted the default SPS parameters in FSPS, that is the MILES stellar library (Sánchez-Blázquez et al. 2006) and the MIST isochrones (Choi et al. 2016). MILES is an empirical stellar library of high spectral resolution (FWHM = 2.5 Å), with a rest-frame wavelength range from 3525 to 7400 Å. The FSPS code also covers the UV and NIR wavelength range but with lower spectral resolution. The MIST models are based on the open-source stellar evolution package MESA (Paxton et al. 2011, 2013, 2015, 2018). Baes et al. (2024) used the Bruzual & Charlot (2003, BC03) SPS models to generate the synthetic images of the TNG50-SKIRT Atlas. Euclid Collaboration: Abdurro’uf et al. (2025) tested whether a systematic difference exists when using different SPS models (FSPS, BC031) with SED fitting and found only a negligible effect on the results. Therefore, the use of FSPS models in this work remains valid. Throughout this paper, we adopted the Chabrier (2003) initial mass function (IMF).

A common practice in SED fitting is to use a flat prior in logarithmic space when sampling stellar metallicity values (Leja et al. 2017, 2018). However, this approach tends to favor the exploration of lower metallicity values, which can influence the fitting results. When fitting only photometric data, a broader parameter space needs to be considered, requiring wider steps for stellar metallicities. To mitigate this issue, a flat prior in linear space can be used (Nersesian et al. 2025). In Appendix A, we show that a flat prior in linear space yields improved Z★ estimates in respect to the flat prior in logarithmic space. Thus, in our analysis we adopted a flat prior in linear space for the stellar metallicity. However, as Nersesian et al. (2024) demonstrated, reliable stellar metallicity constraints remain challenging, even with both broadband and narrowband photometry. Alternatively, employing more physically motivated priors such as a mass–metallicity or age–metallicity relation based on empirical studies of the local Universe (e.g., Bundy et al. 2015) could potentially improve metallicity estimates (e.g., Euclid Collaboration: Corcho-Caballero et al. 2026; Euclid Collaboration: Abdurro’uf et al. 2025), but it may also introduce biases into the inferred mass–metallicity relationship.

Since our analysis relies solely on broadband photometry without narrowband or spectroscopic data, and to optimize computational efficiency, we chose not to include nebular emission in our modeling. Finally, the impact of dust grains on the light from stellar populations at UV and optical wavelengths is modeled using a variable dust attenuation law. The strength of the UV bump is constrained according to the findings of Kriek & Conroy (2013), while the effects of dust attenuation from diffuse dust in the interstellar medium (ISM) and from birth clouds are accounted for separately. The birth-cloud component attenuates the stellar emission from stars with an age up to 10 Myr, and its optical depth is given by

τ ̂ dust , 1 = τ ̂ 1 ( λ / 5500 Å ) − 1 . Mathematical equation: $$ \begin{aligned} \hat{\tau }_{\mathrm{dust} , 1} = \hat{\tau }_{1} \left(\lambda /5500\ \AA \right)^{-1}. \end{aligned} $$(1)

The diffuse dust component has a flexible function that attenuates both stellar and nebular emission of a particular galaxy, using a power-law-modified starburst curve (Calzetti et al. 2000), extended with the Leitherer et al. (2002) curve. We use (Noll et al. 2009)

τ ̂ dust , 2 = τ ̂ 2 4.05 [ k ′ ( λ ) + D ( λ ) ] ( λ λ V ) n , Mathematical equation: $$ \begin{aligned} \hat{\tau }_{\mathrm{dust} , 2} = \frac{\hat{\tau }_{2}}{4.05} \left[k^\prime \left(\lambda \right) + D\left(\lambda \right)\right]\left(\frac{\lambda }{\lambda _V}\right)^{n}, \end{aligned} $$(2)

where k′(λ) is the original Calzetti et al. (2000) attenuation curve, and D(λ) is a Lorentzian-like Drude profile characterizing the UV bump at 2175 Å in the attenuation curve. The normalization of the optical depth of the birth-cloud and diffuse dust components, along with the power-law index (n), which describes the shape of the attenuation curve for the diffuse dust component, are set as free parameters in our model (see Leja et al. 2017, for more details).

We should note here that the SKIRT radiative transfer calculations produce effective attenuation curves that depend on the underlying dust geometry, composition, and viewing angle. In contrast, Prospector assumes a simplified, parametric attenuation law such as the one described above. This choice is intentional, as it mirrors the conditions of real observations where the true dust geometry is unknown and must be approximated through empirical or semiempirical laws. Consequently, part of the dust-age degeneracy discussed in Sect. 4.1 may be driven by this modeling simplification, reflecting an intrinsic limitation of SED fitting under uncertain dust geometries.

While the TNG50-SKIRT Atlas provides dust attenuation maps derived from the radiative transfer simulations, these maps trace only the diffuse dust component. The attenuation associated with birth clouds in star-forming regions is modeled separately in SKIRT, making it nontrivial to construct a total, self-consistent attenuation map. Therefore, a direct comparison between the dust attenuation inferred from our SED fitting and that produced by SKIRT would be inconsistent and potentially misleading.

3.2. SFH treatments

For our fiducial run, we performed a pixel-by-pixel SED fitting, using a flexible nonparametric SFH with a Student’s-t prior distribution (also known as the “continuity” prior) described thoroughly in Leja et al. (2019a). Based on the regularization schemes by Ocvirk et al. (2006) and Tojeiro et al. (2007), the “continuity” prior favors a piecewise constant SFH without sharp transitions in SFR(t). We used eight time elements to describe the nonparametric SFH, specified in look-back time. The recent variations in the SFH were captured by fixing the first two time bins at 0–30 Myr and 30–100 Myr. We also fixed the oldest time bin, placed at 0.85 tuniv − tuniv where tuniv is the age of the Universe at the observed redshift (z = 0.03). The remaining five bins were spaced equally in logarithmic time between 100 Myr and 0.85 tuniv. The stellar metallicity was assumed to be constant across all eight time bins. In other words, a single metallicity parameter describes the entire stellar population. Our fiducial model includes 12 free parameters in total. A list of the free parameters and their associated prior distribution functions are given in Table 1.

In addition to our fiducial model, we performed a pixel-by-pixel SED fitting, using a parametric SFH with a delayed exponential of the following functional form:

SFR ∝ t age e − t age / τ , Mathematical equation: $$ \begin{aligned} \mathrm{SFR} \propto t_\mathrm{age} \, \mathrm{e}^{-t_\mathrm{age} /\tau }, \end{aligned} $$(3)

where tage is the look-back age when SFH began, and τ is the e-folding timescale. In order to account for the uncertainty in the stellar age tage and in τ, we followed the methodology described in Jain et al. (2024). Specifically, we created two datasets of SFHs for each pixel of a particular galaxy: one by varying tage while keeping the stellar mass and τ fixed, and one by varying τ while keeping the stellar mass and tage fixed. To retrieve the total parametric SFH of each pixel, we took the median SFH of the respective sets containing N SFHs, then we calculated the final SFH by averaging the SFHs accounting for the uncertainties in tage and τ.

For the purposes of this analysis, we also fit the integrated photometry of our galaxy sample, assuming the same physical model and parameter space given in Table 1, but again with different treatments of the SFH. In particular, we fit (1) a flexible nonparametric SFH and (2) a delayed τ model, similar to those SFH models adopted on spatially resolved scales. In Table 2, we summarize the different SFH models and provide the nomenclature used in our analysis.

Table 2.

Definitions and nomenclature of the different SFH models in our analysis.

Finally, we sought to investigate the impact of the four Euclid bands on the estimation of the SFHs and stellar properties. For that reason, we performed an additional fitting run using our fiducial SED fitting setup (Table 1), fitting only the GALEX and LSST bands (seven bands). The results of this analysis are presented in the Appendix B of this paper. A more comprehensive analysis of how Euclid observations affect the recovery of galaxy physical properties can be found in Euclid Collaboration: Kovačić et al. (2025) and Euclid Collaboration: Abdurro’uf et al. (2025).

4. Results

4.1. Recovering the stellar properties at spatially resolved scales

We present the results from the fiducial spatially resolved SED fitting performed with Prospector for the sample of 25 TNG50 galaxies. For each of the 6958 pixel bins analyzed, we adopted the median (50th percentile) of the posterior distributions for the inferred physical properties. These values were then used to construct 2D maps of stellar properties, including ΣSFR, Σ★, Z★, and t★,  mw. Here, the mass-weighted stellar age is defined as the average age of stellar populations, weighted by their mass. In particular, t★,  mw is given by the following functional form:

t ★ , m w = ∑ i M i t i ∑ i M i , Mathematical equation: $$ \begin{aligned} t_\star ,\,mw = \frac{\sum _i M_i \, t_i}{\sum _i M_i}, \end{aligned} $$(4)

where Mi is the mass of the i-th component (e.g., a stellar population of a certain age), and ti is the age of that component. The denominator is the total stellar mass.

Figure 3 showcases a few examples of these stellar property maps for the first five galaxies in our sample, sorted by their sub-halo ID. From these maps, we observe that the internal structures of the galaxies are adequately recovered, with detailed resolution of features such as bulges and spiral arms. The results also reproduce well-established trends in the spatial distribution of stellar properties as reported by observational studies (e.g., Sánchez-Blázquez et al. 2014; González Delgado et al. 2016; Ibarra-Medel et al. 2016; Casasola et al. 2017; García-Benito et al. 2017; Zheng et al. 2017; López Fernández et al. 2018; Zhuang et al. 2019; Dale et al. 2020; Parikh et al. 2021; Smith et al. 2022; Pessa et al. 2023) and mock observations (e.g., Nanni et al. 2022, 2023, 2024). Specifically, we find that the majority of stellar mass is concentrated in the bulge regions, where stellar populations tend to be older and more metal-rich compared to those in the disk or outskirts. In contrast, ΣSFR is generally elevated in the spiral arms and outer regions of the galaxies. Inter-arm regions, on average, exhibit lower values of both Σ★ and ΣSFR compared to the spiral arms (e.g., González Delgado et al. 2014).

Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

Examples of the stellar property maps of five TNG50-SKIRT Atlas galaxies, obtained from the spatially resolved SED fitting of GALEX, LSST, and Euclid photometry with Prospector. The maps include the stellar mass surface density (Σ★), the SFR surface density (ΣSFR), stellar metallicity (Z★), and mass-weighted stellar age (t★,  mw). Next to each galaxy’s sub-halo ID, we indicate whether a galaxy is classified as star-forming (SF) or quiescent (Q), based on its global sSFR.

To assess the quality of our fits and evaluate the ability of Prospector to recover stellar properties from broadband photometry, we compared our spatially resolved SED fitting results with the 2D maps of ground truth stellar properties from TNG50, generated using SKIRT. In Fig. 4, we show the residual maps of stellar properties for two example galaxies: a quiescent and a star-forming one. In both cases, the residuals in Σ★ remain within 0.15 dex. For the star-forming galaxy, Σ★ is slightly underestimated in the spiral arms. In the quiescent galaxy, the small number of star-forming gas particles led to only a few usable pixel bins in the ΣSFR map, with an overall underestimated SFR. For the star-forming galaxy, we observe a general pattern: ΣSFR is underestimated in the outer disk and overestimated in the inner spiral arms. These residual patterns likely arise from degeneracies between dust attenuation, stellar age, SFH, and metallicity (see also Fig. 6). For Z★, we find a systematic underestimation across both galaxies. Similarly, t★,  mw is underestimated in the quiescent galaxy. In contrast, for the star-forming galaxy, the ages are well recovered in the inner disk but tend to be underestimated toward the outer edges.

Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

Examples of residual maps of the stellar properties of two TNG50-SKIRT Atlas galaxies, obtained from the nonparametric SFH configuration. The residual maps include the stellar mass surface density (Σ★), the SFR surface density (ΣSFR), stellar metallicity (Z★), and mass-weighted stellar age (t★,  mw). Each residual map is computed as the logarithmic ratio of the inferred property from Prospector to the ground truth value from TNG50, expressed as log10(Prospector/TNG50).

In Fig. 5, we present the comparison of the four stellar properties for all spatial pixel bins across the 25 galaxies in our sample. We color-coded the scatter plots by the diffuse dust optical depth ( τ ̂ dust , 2 Mathematical equation: $ \hat{\tau}_{\mathrm{dust,\, 2}} $) in the V band. It is well known that SED fitting with photometry alone can suffer from the age-dust-metallicity degeneracy, as highlighted by previous studies (e.g., Worthey et al. 1994; Silva et al. 1998; Devriendt et al. 1999; Pozzetti & Mannucci 2000; Bell & de Jong 2001; Walcher et al. 2011; Díaz-García et al. 2015). We derived the average trend of τ ̂ dust , 2 Mathematical equation: $ \hat{\tau}_{\mathrm{dust,\, 2}} $ with the stellar properties with the Locally Weighted Regression (LOESS) method (Cleveland & Devlin 1988) as implemented in the LOESS routine2 by Cappellari et al. (2013). To facilitate quantitative comparisons, we computed the mean offset (μ) and scatter (standard deviation, σ) for the logarithmic ratio between the best-fit parameters and the ground truth. We also calculated the Spearman rank-order correlation coefficient (ρ), indicative of the significance of the relationship between two datasets. The results for each stellar property are summarized at the top-left of each panel in Fig. 5.

Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

Comparisons of the spatially resolved stellar population properties for 6958 pixel bins in the 25 sample galaxies derived from spatially resolved SED fitting with a nonparametric SFH model. The depicted physical properties include the stellar mass surface density (Σ★), the SFR surface density (ΣSFR), the stellar metallicity (Z★), and the mass-weighted stellar age (t★,  mw). The colors indicate the average trend with the dust optical depth in the V band ( τ ̂ dust , 2 Mathematical equation: $ \hat{\tau}_{\mathrm{dust,\, 2}} $), calculated with the LOESS method. Each hexbin includes a minimum of five points. Contours enclose 20%, 50%, and 80% of the total data. The dashed orange line represents the one-to-one relation. The Spearman’s correlation coefficient (ρ), measured bias (average offset μ), and scatter (standard deviation σ) relative to the ground truth from TNG50 are indicated in the legend of each panel.

For Σ★, we find excellent agreement with the ground truth, as evidenced by a small offset (μ = 0.01 dex), low scatter (σ = 0.1 dex), and a high correlation (ρ = 0.98). This suggests that stellar mass can be reliably recovered on spatially resolved scales using photometry, which covers the rest-frame NIR. In particular, this demonstrates the strong potential of Euclid observations for mapping Σ★ in local galaxies.

For ΣSFR, the recovery is similarly strong (ρ = 0.88), but with an increased mean offset (μ = 0.1 dex), and scatter around the mean (σ = 0.3 dex). In particular, the offset is stronger at the lower end of the distribution, where values of ΣSFR are very low and comparisons become less meaningful. At these lower values, dust attenuation appears to increase the inferred SFR due to the dust-age degeneracy. At the high end of the distribution, we notice that the ΣSFR values derived using Prospector tend to be lower than the ground truth. These discrepancies in the ΣSFR maps may stem from the limited number of stellar particles in the latest time bin of the TNG50 ΣSFR maps, which can result in spurious SFR values.

For Z★, we measure a significant systematic offset of μ = −0.15 dex, with a low scatter (σ = 0.1 dex) and a moderate correlation (ρ = 0.59). This is an expected result, simply because there is not enough information in broadband photometry to provide meaningful constraints on stellar metallicity (Nersesian et al. 2024; Csizi et al. 2024). Euclid Collaboration: Abdurro’uf et al. (2025) demonstrated that a better recovery of Z★ can be achieved when using mass-metallicity and age-metallicity priors, highlighting the importance of these prior functions in improving SED fitting with Bayesian methods.

The recovery of the stellar ages is limited due to a narrow dynamic range, with a smaller systematic offset (μ = −0.04 dex), and scatter (σ = 0.09 dex), yet with a more moderate correlation coefficient (ρ = 0.41). Notably, the stellar age estimates show a closer alignment with the ground truth than for Z★. The largest deviations from the true values occur for older stellar populations, while younger populations are better recovered. These findings are consistent with those by Euclid Collaboration: Abdurro’uf et al. (2025), where more informed priors were used. This suggests that the use of specialized priors may not significantly improve age estimates for this dataset, which covers only the rest-frame UV-optical-NIR regime and lacks optical spectroscopy, an important constraint for age determination (Nersesian et al. 2024, 2025). The absence of spectroscopy, and particularly the lack of spectral features that are key age indicators, such as Hδ, likely contributes to this limitation. Finally, color-coding by τ ̂ dust , 2 Mathematical equation: $ \hat{\tau}_{\mathrm{dust,\, 2}} $ reveals that the secondary peak in stellar ages is linked to the dust-age degeneracy, where Prospector infers younger ages, and the red colors are primarily attributed to dust attenuation rather than the presence of older stellar populations.

To investigate these degeneracies, we examined a representative pixel bin with an overestimated SFR and underestimated stellar age and inspected the joint posterior distributions of its inferred properties (Fig. 6). The corner plot highlights several well-known correlations: a strong positive correlation between τ ̂ dust , 2 Mathematical equation: $ \hat{\tau}_{\mathrm{dust,\, 2}} $ and sSFR (ρ = 0.69), and moderate anticorrelations of τ ̂ dust , 2 Mathematical equation: $ \hat{\tau}_{\mathrm{dust,\, 2}} $ with t★,  mw (ρ = −0.49) and Z★ (ρ = −0.38). In addition, Σ★ correlates tightly with t★,  mw (ρ = 0.81), while strongly anticorrelating with both sSFR (ρ = −0.69) and dust attenuation (ρ = −0.63). These results confirm that the well-known degeneracies between dust, age, star formation, and metallicity persist in the analysis, contributing to the observed biases in the recovered stellar properties. Breaking these degeneracies remains a challenge and likely requires the inclusion of spectroscopic information as well as information in the infrared wavelength regime.

Thumbnail: Fig. 6. Refer to the following caption and surrounding text. Fig. 6.

Joint posterior distributions of the main physical properties in our analysis, the best-fit SED, and the recovered SFH for a selected pixel bin of TNG414917 with the O4 orientation index. The corner plot displays the posterior distributions of τ ̂ dust , 2 Mathematical equation: $ \hat{\tau}_{\mathrm{ dust,\, 2}} $, sSFR, t★,  mw, Z★, and Σ★. Contours enclose the 20%, 50%, and 80% of the total data, while the vertical dashed black lines mark the median values for each parameter. The median and the 16th–84th percentile range for each physical property are given at the top of the corresponding posterior distribution. The Spearman’s correlation coefficient, ρ, quantifies the strength of correlation between two posterior distributions. Top right: Comparison of the input photometric data (red points) to the best-fit SED model (blue line). Bottom right: Subpanel illustrating the absolute residuals defined as ΔFν = log10(fit/ground truth). The photometric uncertainties are not visible in the plot because they are smaller than the symbol size. Bottom: Posterior distribution of the SFH.

4.2. Recovering the shapes of galaxy SFHs

In this section, we present the SFHs of our sample of TNG50 galaxies and evaluate the different approaches for reconstructing the SFH of a galaxy. We analyze the differences between the SFHs obtained from our fiducial spatially resolved SED fitting with Prospector (using a nonparametric model) to those derived from resolved fits assuming a parametric SFH, as well as to global-scale fits using both nonparametric and parametric models. These are further compared to the ground truth SFHs from TNG50 on resolved scales. This comparison allows us to assess how modeling assumptions and spatial resolution impact the inferred SFHs. To compute the total SFH of each galaxy from the spatially resolved data, we sum the SFRs across all pixel bins associated with that galaxy, at the specific time bin. This procedure is applied consistently to the spatially resolved SED fitting results, and the TNG50 spatially resolved SFHs.

Figure 7 presents the ΣSFR maps alongside a comparison of the five global SFH models for the first ten galaxies in our sample, sorted by their sub-halo ID. For each galaxy, the left panel shows the ΣSFR map derived using the nonparametric SFH model, highlighting the spatial distribution of the present SFR. The right panel compares the four global SFH models reconstructed with Prospector to the ground truth SFH from TNG50. The definitions and nomenclature of the different SFH models are given in Table 2.

Thumbnail: Fig. 7. Refer to the following caption and surrounding text. Fig. 7.

Examples of SFR surface density (ΣSFR) maps and SFH models for ten galaxies in our TNG50-SKIRT Atlas sample. Left panel: ΣSFR map for each galaxy, as estimated with Prospector. Right panel: Comparison of various global SFH models to the ground truth SFH from TNG50. Specifically, we show the global nonparametric SFH inferred from our fiducial spatially resolved run (SFHres,  np; orange), the global parametric SFH inferred from the spatially resolved map (SFHres,  τ; pink), the nonparametric SFH derived from fitting the integrated photometry (SFHglob,  np; yellow), the global parametric SFH based on a simple τ-model (SFHglob,  τ; blue), and the true global SFH from the TNG50 resolved map (SFHres,  TNG50; green).

A first qualitative inspection of the SFHs shown in Fig. 7 reveals that the SFHres,  np and SFHglob,  np align more closely with the ground truth SFHres,  TNG50 than the parametric SFHglob τ and SFHres τ. While all four SFH models reproduce in general the shape of SFHres,  TNG50, the nonparametric models seem to describe better the extended period of star formation. On the other hand, both parametric SFH models fail to match the SFR amplitude of SFHres,  TNG50, while the peak of SFR is skewed at a later cosmic time. This mismatch can lead to underestimated stellar masses and younger inferred stellar ages, highlighting the limitations of relying on parametric SFH assumptions.

A more quantitative way to compare SFH models with different parameterization schemes is through an analysis of their cumulative stellar mass growth and formation times. This approach is particularly effective for evaluating how the various SFH models assemble the stellar mass of galaxies over time. Figure 8 presents the recovered cumulative stellar mass growth curves, normalized to the total stellar mass estimated by each method,3 as a function of look-back time for the full sample of 25 TNG50-SKIRT Atlas galaxies. It is evident that SFHres,  np aligns with the ground truth SFHres,  TNG50, capturing both the earlier formation time and the extended duration of star formation. Similarly, SFHglob,  np is in very good agreement with both the ground truth SFHres,  TNG50 and the resolved SFHres,  np. In contrast, the parametric model, in both resolved (SFHres,  τ) and integrated scales (SFHglob,  τ), shows quite large differences in respect to SFHres,  TNG50, with delayed and less extended SFHs.

Thumbnail: Fig. 8. Refer to the following caption and surrounding text. Fig. 8.

Comparison of the normalized cumulative stellar mass growth curves of all 25 galaxies in our sample. We plot the cumulative stellar mass growth curves from the global nonparametric SFH inferred from our fiducial spatially resolved run (SFHres,  np; orange), the global parametric SFH inferred from the spatially resolved map (SFHres,  τ; pink), the nonparametric SFH derived from fitting the integrated photometry (SFHglob,  np; yellow), the global parametric SFH based on a simple τ-model (SFHglob,  τ; blue), and the true global SFH from the TNG50 resolved map (SFHres,  TNG50; green). The dotted lines indicate the different formation times by which 20%, 50%, and 80% of a galaxy’s stellar mass was formed. Both nonparametric SFHs (SFHres,  np and SFHglob,  np) seem to be in better agreement with SFHres,  TNG50, while most of the parametric SFHs (both SFHres,  τ and SFHglob,  τ) are skewed toward a later cosmic time, thus lacking the contribution from the older stellar populations.

One key observation is that, for most galaxies, the cumulative stellar mass growth predicted by the parametric SFHs is systematically biased toward later formation times. To further quantify these discrepancies, Fig. 9 compares the formation times at which 20% (t20), 50% (t50), and 80% (t80) of a galaxy’s stellar mass was assembled across the different SFH models. We quantify the differences by calculating the mean offset (μ) and scatter (standard deviation, σ) relative to the ground truth SFHs from TNG50. The values of μ and σ for each formation time are given in Table 3.

Table 3.

Mean offset and standard deviation of formation times and global properties from the TNG50 reference for each SFH model.

Thumbnail: Fig. 9. Refer to the following caption and surrounding text. Fig. 9.

Comparison of the formation timescales for the four different SFH models with respect to the ground truth values from the TNG50 simulation. Each marker represents the formation time of a galaxy in our sample, color-coded according to the corresponding SFH model. Star-forming galaxies are indicated with stars, while quiescent galaxies with points. From left to right: Time for 20% (t20), 50% (t50), and 80% (t80) of a galaxy’s stellar mass formation. The four models include the spatially resolved model from pixel-by-pixel SED fitting with a nonparametric (orange) and parametric (pink) SFH, the integrated SED fitting with a nonparametric SFH (yellow), and a simple τ model within Prospector (blue). The measured bias (average offset μ) and scatter (standard deviation σ) relative to the ground truth from TNG50 are indicated in the legends of the panels. The formation times inferred from the spatially resolved and integrated nonparametric SFHs are in best agreement with the ground truth values.

The nonparametric SFHs, both spatially resolved and global, exhibit strong agreement with the ground truth across all timescales. Our fiducial model reproduces t20 within 0.1 dex, while the global nonparametric model (SFHglob,  np) shows an even smaller offset of 0.07 dex. Interestingly, both models tend to recover slightly earlier t20 values than the reference SFHres,  TNG50. In contrast, the parametric SFH models significantly delay the early phases of star formation, and thus underestimate the contribution of the older stellar population.

For t50, both nonparametric SFH models achieve an excellent agreement with the ground truth (μ = 0.01 dex), and even for t80, the offsets remain small (within 0.05 dex). These results highlight that spatially resolved fitting with a nonparametric SFH prior can yield accurate reconstructions of galaxy SFHs. Furthermore, we report that nonparametric SFHs at global scales still perform remarkably well. This result suggests that the use of flexible SFHs at global scales can return meaningful formation histories and are less susceptible to the outshining effect compared to parametric SFHs. This provides a significant advantage, especially when spatially resolved data are not available. Crucially, this remains true even in the absence of spectroscopic information, underscoring the power of spatially resolved, and most importantly nonparametric approaches in SFH reconstruction.

4.3. Measuring the bias in global properties

The results presented in this paper suggest that nonparametric SFHs, either at spatially resolved or integrated scales, describe the variation of star formation across cosmic time more accurately. On the other hand, the use of a delayed τ SFH leads to biased fits, often skewing the SFH peak to later cosmic times and thus underestimating the contribution from the older stellar populations.

To understand better the impact of the systematic biases introduced by the different SFH models, we compare three key global stellar properties (M★, sSFR, and t★,  mw) to the ground truth values from the TNG50 simulation. For the spatially resolved fits, we derived the global properties of a particular galaxy by summing or averaging across its pixel bins. We obtained M★ and SFR of a particular galaxy, by summing the values from each pixel bin, while t★,  mw is computed as the stellar-mass-weighted average of the pixel ages. In Fig. 10, we show the absolute differences in these properties relative to the TNG50 ground truth (top row), as well as the distribution of the absolute differences (bottom row) along with the mean offset (μ), and scatter (standard deviation, σ). The values of μ and σ for each global physical property are given in Table 3.

Thumbnail: Fig. 10. Refer to the following caption and surrounding text. Fig. 10.

Differences in the global stellar properties for the four different SFH models with respect to the ground truth global values from the TNG50 simulation. From left to right: M★, sSFR, and t★,  mw. The properties were derived through SED fitting under four different configurations: spatially resolved fitting with a nonparametric SFHres,  np (orange), spatially resolved SED fitting with a simple τ-model SFHres,  τ (pink), global fitting with a nonparametric SFHglob,  np (yellow), and a global fitting with a simple τ-model SFHglob,  τ (blue). Typical uncertainties for each configuration are shown in the bottom-left corner of each panel. Star-forming galaxies are indicated with stars, while quiescent galaxies with points. The bottom panels show the residual distributions of each stellar property. The measured bias (average offset) and scatter (standard deviation) are indicated in the legends of the panels. The properties inferred from the spatially resolved and integrated nonparametric SFHs are in best agreement with the ground truth values.

Focusing on the left column of Fig. 10, we report that the global stellar masses inferred from the nonparametric model, at both resolved (SFHres,  np) and integrated (SFHglob,  np) scales, are closest to the ground truth, with a mean offset near zero. Conversely, the stellar masses inferred from the resolved parametric SFHres,  τ and global parametric SFHglob,  τ are systematically underestimated by 0.12 and 0.06 dex, respectively. The underestimation of the stellar mass at global scales is expected (e.g., Jain et al. 2024; Mosleh et al. 2025), largely due to the outshining effect, where the younger stellar populations (< 10 Myr) dominate the observed SED and outshine the older stars (> 200 Myr). The underestimation of M★ is more prominent in the case of the delayed τ SFH model. However, a more flexible parametric model, such as a double power-law SFH, could be a reliable option at integrated scales (e.g., Harvey et al. 2025).

Interestingly, we measure the strongest deviation in the resolved parametric model despite the expectation that spatial resolution should mitigate outshining. Iglesias-Navarro et al. (2025) similarly found that applying a delayed-τ model in a pixel-by-pixel analysis of JWST data led to significant deviations in inferred galaxy properties, highlighting the limitations of simple parametric SFHs. Other studies have demonstrated that resolved parametric SFHs can yield accurate global M★ when more flexible functional forms are used. For instance, Euclid Collaboration: Abdurro’uf et al. (2025) employed a double power-law model at resolved scales along with physically motivated priors, resulting in very accurate global M★ for the same sample of 25 TNG50 galaxies. Similarly, Mosleh et al. (2025) reported that the double power-law SFH at resolved scales consistently outperformed other parametric SFH models, concluding that such parameterization is essential for accurate resolved SED fitting. In contrast, our analysis uses a simpler delayed-τ model, which limits accuracy and likely explains the larger deviations from the ground truth.

The middle column of Fig. 10 shows the comparison between the fitted sSFR and the ground truth TNG50 values. The present SFR is used here to calculate the sSFR. The sSFR values inferred from SFHres,  np are in excellent agreement with the TNG50 ground truth, with a mean offset of μ = 0.01 dex. In contrast, the largest discrepancies arise in the case of SFHres,  τ, which systematically overestimates the sSFR. This overestimation is mainly driven by the combination of higher inferred SFRs (see Fig. 7) and lower stellar masses returned by the delayed-τ model, both of which bias the sSFR upward. Furthermore, the comparison reveals that the sSFR estimates from SFHglob,  np and SFHglob,  τ are systematically lower than the ground truth (the peak of the residual distribution is below 0 dex), yet likely for different reasons. In the case of SFHglob,  np, we notice that the mass-weighted ages are on average older than the ground truth by 0.02 dex, resulting in lower SFRs and consequently lower sSFR (since M★ is properly recovered). Conversely, in the case of SFHglob,  τ the systematically lower M★ and present SFR (see Fig. 7) result in the underestimation of the sSFR. For quiescent galaxies (log10(sSFRTNG50/yr−1)≤ − 11), SFHres,  np and SFHres,  τ overestimate the sSFR, despite these galaxies having little to no active star formation. This inaccuracy may arise from the limited number of stellar particles in the latest TNG50 time bins, which can produce spurious SFR values at spatial scales, especially for quiescent galaxies. Uncertainties in the stellar ages of individual pixels, linked to the dust-age degeneracy, may further contribute to the overestimation of sSFR in galaxies with minimal star formation (see Fig. 5). On the other hand, the sSFR recovered by both SFHs at global scales may suggest that for passive galaxies, simpler parametric models (e.g., such as τ models) can yield accurate integrated parameters.

Finally, the t★,  mw estimates from SFHres,  np are the closest to the TNG50 ground truth, with a mean offset of μ = 0.01 dex. The SFHglob,  np also performs well, with μ = 0.02 dex. It is worth noting that in both nonparametric models, a potential correlation appears between the residuals and age, where younger ages tend to be overestimated, while older ages are underestimated within 0.1 dex. On the other hand, t★,  mw estimates from SFHres,  τ and SFHglob,  τ show significant deviations, exceeding 0.2 dex. This discrepancy results from the skewed parametric SFHs toward later cosmic times, obscuring the extended SFHs of early-type galaxies and significantly affecting the inferred t★,  mw. These findings align with Suess et al. (2022), Jain et al. (2024), and Euclid Collaboration: Enia et al. (2026), who also reported extremely young inferred stellar ages from SFHs derived on integrated scales with simple parametric SFH models.

5. Discussion

A key takeaway from this work is that nonparametric SFHs provide a more effective means of mitigating the outshining effect from the young stellar populations (< 10 Myr), and therefore improving the accuracy of inferred galaxy properties. However, this improvement should not be overinterpreted, since our results are based on a small sample of galaxies. The primary advantage of spatially resolved SED fitting lies in its ability to recover the internal color variations from the different stellar populations and dust attenuation, capturing the complexity of galaxy evolution that is often missed in global analyses. Certainly, parametric SFH models on spatially resolved scales remain valuable, as they are less computationally intensive, and depending on the functional form can provide a robust reconstruction of SFHs, and yield less biased global galaxy properties (e.g., Giménez-Arteaga et al. 2023, 2024; Jain et al. 2024; Mosleh et al. 2025; Euclid Collaboration: Abdurro’uf et al. 2025).

In this study, we show evidence that nonparametric modeling on resolved and global scales can recover the SFHs and derived properties with higher fidelity, making it a powerful tool for understanding galaxy formation and evolution. In some cases, SFHglob,  np provides a closer match to the TNG50 ground truth than SFHres,  np. This effect likely arises from a combination of factors: residual photometric noise and de-blending uncertainties at the pixel level, stronger sensitivity to outshining in individual regions, and the fact that the integrated light tends to smooth out these localized mismatches. While the resolved fits allow us to recover spatial variations in SFHs, the global fits can sometimes provide a closer match to the average SFH shape, particularly for galaxies with low levels of ongoing star formation.

Our choice to focus on a subset of 25 galaxies from the TNG50-SKIRT Atlas was motivated by the need to balance computational feasibility with sufficient diversity in morphology and star formation properties to test the pixel-by-pixel SED fitting framework. While this limited sample size prevents us from drawing definitive conclusions about the impact of spatial resolution on integrated quantities, it serves as a robust proof of concept for the methodology. The selected galaxies are biased toward relatively massive systems with near-solar metallicities and older stellar populations, reflecting both the stellar mass function of TNG50 at z = 0 and the fact that galaxies with reliable radiative transfer post-processing typically require sufficient stellar particle sampling. Extending this work to include lower-mass, metal-poor galaxies would certainly broaden the dynamic range of stellar population properties, but such an expansion lies beyond the scope of the present study, which is intended as a methodological benchmark in preparation for Euclid.

Furthermore, one persistent limitation of our approach is the dust-age degeneracy discussed in 4.1, which is partly driven by the use of a simplified, parametric attenuation law in the SED fitting, whereas SKIRT models the full radiative transfer and dust geometry. This intentional mismatch reflects a realistic challenge faced in observational studies, where the true dust geometry is unknown and must be approximated by empirical or flexible attenuation models. Thus, the resulting degeneracies are not a flaw of the method but a realistic feature of the problem. A direct validation of the recovered dust attenuation was not feasible here because the current TNG50-SKIRT Atlas only provides diffuse dust attenuation maps. The attenuation in star-forming regions is modeled separately, preventing a consistent reconstruction of the total attenuation. Updated radiative transfer simulations of the TNG50-SKIRT Atlas are currently underway to address this limitation and include the dust emission in the far-infrared, enabling a more accurate assessment of dust attenuation in future work.

The broader implications of these results extend to testing models of galaxy evolution and evaluating the internal consistency of observational datasets. For instance, as highlighted by Fu et al. (2024), in a ΛCDM paradigm of structure formation, galaxies evolve hierarchically within dark matter halos via mergers and star formation. This growth is reflected in a set of interconnected observables, most notably the stellar mass functions, the SFHs, the satellite galaxy abundances, the merger rates, and even other such properties as intra-cluster light, each one capturing different aspects of the galaxy formation process. A physically sound model of galaxy evolution must be able to simultaneously predict and reproduce all these observables in a consistent way. However, this challenge is often complicated by inconsistencies among datasets, driven by differences in observational methods and assumptions. As demonstrated here, Euclid will offer the opportunity to measure nonparametric SFHs, stellar mass functions, as well as intra-cluster light from the same datasets across large areas and at different cosmic epochs, allowing for 1) testing the self-consistency of these data and 2) validating hierarchical models of galaxy formation. The rich and varied datasets from Euclid, including the SFHs put forward here, will help identify systematic tensions among differential and integrated quantities (SFHs versus stellar mass functions, for example), guiding the refinements in both modeling and data calibration. Ultimately, such holistic approaches will provide a critical bridge between observation and theory, enabling the development of galaxy evolution models that are not only descriptive but also predictive, robust, and self-consistent across cosmic time.

6. Summary and conclusions

Euclid will observe approximately 14 000 deg2 of the extragalactic sky, delivering high-resolution imaging that enables spatially resolved studies of nearby galaxies. To assess the feasibility of reconstructing accurate SFHs from Euclid-like data, we presented detailed measurements of SFHs on both global and spatially resolved scales for a sample of 25 TNG50-SKIRT Atlas galaxies using Prospector. We examine the SFHs derived via (1) resolved SED fitting using a nonparametric SFH model (SFHres,  np), (2) resolved SED fitting using a commonly adopted parametric τ-model (SFHres,  τ), (3) global SED fitting with a nonparametric model (SFHglob,  np), and (4) global SED fitting using a parametric SFH (SFHglob,  τ). We applied these models to mock UV-NIR photometry derived from the TNG50 cosmological simulation and processed with the radiative transfer code SKIRT. This approach provides a proof of concept for extracting spatially resolved SFHs of local galaxies with Euclid, highlighting the strengths and limitations of nonparametric SFH modeling in the context of next-generation galaxy surveys. Our main findings and conclusions are the following:

  • The mass surface density is recovered with high accuracy using broadband photometry, thanks to the NIR coverage from Euclid (mean offset μ = 0.01 dex and scatter σ = 0.1 dex from the ground truth values). The SFR surface density and stellar ages are reasonably recovered (μ ≤ 0.04 dex), considering the fact that we only fit broadband photometry without using an age-mass prior. A secondary peak in the stellar age distribution is observed, where Prospector assigns more dust to explain the redder colors in the spatial bins. Finally, we show that photometry-only fits with a linear Z★ prior function do not yield meaningful constraints on Z★ for individual galaxies, although they effectively avoid a downward bias. While the precision remains limited, the accuracy is generally good. We also note that the results for Z★ are highly sensitive to the choice of prior function.

  • In Figs. 7 and 8, we find that the best method to reconstruct the SFHs of nearby galaxies, when using photometry alone, is employing a nonparametric SFH model at spatially resolved scales (SFHres,  np). Indeed, the nonparametric SFH, at resolved scales, is in very good agreement with the ground truth SFH from TNG50. The nonparametric SFH at global scales appears to be a reliable method as well, returning robust stellar mass and age estimates. The flexibility of the nonparametric model allows it to better capture the different variations and amplitudes of star formation across cosmic time, while being less affected by outshining from recent star formation activity.

  • Based on the formation timescales derived from the different SFH models, we find that nonparametric SFHs at global scales (SFHglob,  np) are generally more reliable for reconstructing SFHs when spatial resolution is not available. In contrast, the delayed τ-model, at both resolved (SFHres,  τ) and integrated (SFHglob,  τ) scales, tends to fail in capturing early star formation activity, leading to significant deviations in metrics such as t20 and t50 compared to the ground truth SFHs from TNG50 (see Table 3). This limitation results in systematic underestimation of both the stellar mass and the average age of the stellar populations. Nonetheless, the τ-model remains useful for characterizing recent star formation activity, as UV-bright stars dominate the global SED.

  • We quantified the systematics in three key global properties (i.e., stellar masses, sSFR, and stellar ages) as derived from the different SFH treatments, by comparing them to the ground truth values from TNG50 (see Table 3). Overall, we find that both stellar mass and stellar age are severely underestimated compared to the ground truth values when the simple τ-model is used, while deviations are reduced when using the nonparametric SFH model, especially at resolved scales. In particular, for the global stellar mass, we measure a mean offset of μ = −0.03, μ = −0.12, μ = 0.01, and μ = −0.06 for the SFHres,  np, SFHres,  τ, SFHglob,  np, and SFHglob,  τ, respectively. For the sSFR, we find deviations of μ = 0.01, μ = 0.4, μ = −0.17, and μ = −0.24 for the SFHres,  np, SFHres,  τ, SFHglob,  np, and SFHglob,  τ, respectively. For quiescent galaxies, the difference in sSFR recovered by SFHglob,  τ is reduced, suggesting that simpler parametric models can yield accurate integrated parameters. Finally, for the mass-weighted ages, we estimate a mean offset of μ = 0.01, μ = −0.27, μ = 0.02, and μ = −0.23 for the SFHres,  np, SFHres,  τ, SFHglob,  np, and SFHglob,  τ, respectively.

In summary, spatially resolved nonparametric SFHs provide a more effective approach to mitigating outshining effects, offering improved accuracy in the determination of galaxy stellar properties. This highlights the importance of using spatially resolved analyses to fully capture the complex evolutionary pathways of galaxies and avoid biases inherent in global parametric models. In the context of the Euclid mission, this study demonstrates that photometry-only fits, when combined with nonparametric models and spatial resolution, can enable realistic reconstructions of galaxy SFHs, paving the way for transformative insights into galaxy formation and evolution. Given the complexity and computational demands of spatially resolved SED fitting with a nonparametric SFH, our method is focused on the quality of the inferred galactic properties and can be applied to a small, yet still significant, number of galaxies in the local Universe using deep, high-resolution imaging data from Euclid, complemented by optical observations from other surveys such as UNIONS, DES, and LSST, as well as UV observations from GALEX. The wide sky coverage provided by these surveys will enable the analysis of a diverse range of galaxies, residing in different environments and encompassing a broad range of stellar properties and morphologies.

Acknowledgments

We thank the anonymous referee for the valuable remarks and suggestions that helped us improve the paper. AN gratefully acknowledges the support of the Belgian Federal Science Policy Office (BELSPO) for the provision of financial support in the framework of the PRODEX Programme of the European Space Agency (ESA) under contract number 4000143347. EDC acknowledges funding from the European Union grant WIDERA ExGal-Twin, GA 101158446. This research makes use of ESA Datalabs (datalabs.esa.int/), an initiative by ESA’s Data Science and Archives Division in the Science and Operations Department, Directorate of Science. The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the State Secretariat for Education, Research, and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (www.euclid-ec.org). This work has made use of CosmoHub, developed by PIC (maintained by IFAE and CIEMAT) in collaboration with ICE-CSIC. CosmoHub received funding from the Spanish government (MCIN/AEI/10.13039/501100011033), the EU NextGeneration/PRTR (PRTR-C17.I1), and the Generalitat de Catalunya.

References

  1. Abbott, T. M. C., Adamów, M., Aguena, M., et al. 2021, ApJS, 255, 20 [NASA ADS] [CrossRef] [Google Scholar]
  2. Abdurro’uf, Lin, Y.-T., Wu, P.-F., & Akiyama, M. 2021, ApJS, 254, 15 [CrossRef] [Google Scholar]
  3. Abdurro’uf, Lin, Y.-T., Hirashita, H., et al. 2022a, ApJ, 926, 81 [Google Scholar]
  4. Abdurro’uf, Lin, Y. T., Wu, P. F., & Akiyama, M. 2022b, piXedfit: Analyze spatially resolved SEDs of galaxies, Astrophysics Source Code Library [record ascl:2207.033] [Google Scholar]
  5. Arnouts, S., Le Floc’h, E., Chevallard, J., et al. 2013, A&A, 558, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Baba, J., Saitoh, T. R., & Wada, K. 2013, ApJ, 763, 46 [Google Scholar]
  7. Baes, M., Gebek, A., Trčka, A., et al. 2024, A&A, 683, A181 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Baldry, I. K., Glazebrook, K., Brinkmann, J., et al. 2004, ApJ, 600, 681 [Google Scholar]
  9. Barišić, I., van der Wel, A., Bezanson, R., et al. 2017, ApJ, 847, 72 [Google Scholar]
  10. Barnes, J. E., & Hernquist, L. E. 1991, ApJ, 370, L65 [Google Scholar]
  11. Barnes, J. E., & Hernquist, L. 1996, ApJ, 471, 115 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bell, E. F., & de Jong, R. S. 2001, ApJ, 550, 212 [Google Scholar]
  13. Best, P. N., Kauffmann, G., Heckman, T. M., et al. 2005, MNRAS, 362, 25 [Google Scholar]
  14. Bianchi, L., Conti, A., & Shiao, B. 2014, Adv. Space Res., 53, 900 [Google Scholar]
  15. Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151 [Google Scholar]
  17. Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [NASA ADS] [CrossRef] [Google Scholar]
  18. Bundy, K., Bershady, M. A., Law, D. R., et al. 2015, ApJ, 798, 7 [Google Scholar]
  19. Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682 [NASA ADS] [CrossRef] [Google Scholar]
  20. Camps, P., & Baes, M. 2015, Astron. Comput., 9, 20 [Google Scholar]
  21. Camps, P., & Baes, M. 2020, Astron. Comput., 31, 100381 [NASA ADS] [CrossRef] [Google Scholar]
  22. Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2013, MNRAS, 432, 1862 [NASA ADS] [CrossRef] [Google Scholar]
  23. Carnall, A. C., Leja, J., Johnson, B. D., et al. 2019, ApJ, 873, 44 [Google Scholar]
  24. Casasola, V., Cassarà, L. P., Bianchi, S., et al. 2017, A&A, 605, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Cecchi, R., Bolzonella, M., Cimatti, A., & Girelli, G. 2019, ApJ, 880, L14 [NASA ADS] [CrossRef] [Google Scholar]
  26. Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
  27. Chauke, P., van der Wel, A., Pacifici, C., et al. 2018, ApJ, 861, 13 [NASA ADS] [CrossRef] [Google Scholar]
  28. Cheung, E., Bundy, K., Cappellari, M., et al. 2016, Nature, 533, 504 [NASA ADS] [CrossRef] [Google Scholar]
  29. Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102 [Google Scholar]
  30. Ciesla, L., Elbaz, D., & Fensch, J. 2017, A&A, 608, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Cleveland, W. S., & Devlin, S. J. 1988, J. Am. Stat. Assoc., 83, 596 [CrossRef] [Google Scholar]
  32. Conroy, C. 2013, ARA&A, 51, 393 [NASA ADS] [CrossRef] [Google Scholar]
  33. Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486 [Google Scholar]
  34. Crenshaw, D. M., Kraemer, S. B., & George, I. M. 2003, ARA&A, 41, 117 [NASA ADS] [CrossRef] [Google Scholar]
  35. Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11 [Google Scholar]
  36. Csizi, B., Tortorelli, L., Siudek, M., et al. 2024, A&A, 689, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Cuillandre, J. C., Bolzonella, M., Boselli, A., et al. 2025, A&A, 697, A11 [Google Scholar]
  38. Daddi, E., Cimatti, A., Renzini, A., et al. 2004, ApJ, 617, 746 [NASA ADS] [CrossRef] [Google Scholar]
  39. Daddi, E., Dickinson, M., Morrison, G., et al. 2007, ApJ, 670, 156 [NASA ADS] [CrossRef] [Google Scholar]
  40. Dale, D. A., Anderson, K. R., Bran, L. M., et al. 2020, AJ, 159, 195 [NASA ADS] [CrossRef] [Google Scholar]
  41. Davé, R. 2008, MNRAS, 385, 147 [Google Scholar]
  42. Davé, R., Oppenheimer, B. D., & Finlator, K. 2011, MNRAS, 415, 11 [Google Scholar]
  43. De Vis, P., Jones, A., Viaene, S., et al. 2019, A&A, 623, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Dekel, A., & Burkert, A. 2014, MNRAS, 438, 1870 [Google Scholar]
  45. Devriendt, J. E. G., Guiderdoni, B., & Sadat, R. 1999, A&A, 350, 381 [NASA ADS] [Google Scholar]
  46. Díaz-García, L. A., Cenarro, A. J., López-Sanjuan, C., et al. 2015, A&A, 582, A14 [Google Scholar]
  47. Diemer, B., Sparre, M., Abramson, L. E., & Torrey, P. 2017, ApJ, 839, 26 [Google Scholar]
  48. Dobbs, C., & Baba, J. 2014, PASA, 31, e035 [NASA ADS] [CrossRef] [Google Scholar]
  49. Dome, T., Tacchella, S., Fialkov, A., et al. 2024, MNRAS, 527, 2139 [Google Scholar]
  50. Donnari, M., Pillepich, A., Nelson, D., et al. 2019, MNRAS, 485, 4817 [Google Scholar]
  51. Durán-Camacho, E., & Duarte-Cabral, A. 2025, A&A, submitted [arXiv:2506.13560] [Google Scholar]
  52. Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Elmegreen, B. G. 2011, in EAS Publications Series, eds. C. Charbonnel, & T. Montmerle, EAS Pub. Ser., 51, 19 [Google Scholar]
  54. Euclid Collaboration (Martinet, N., et al.) 2019, A&A, 627, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Euclid Collaboration (Merlin, E., et al.) 2023, A&A, 671, A101 [Google Scholar]
  56. Euclid Collaboration (Abdurro’uf, et al.) 2025, A&A, 702, A72 [Google Scholar]
  57. Euclid Collaboration (Cropper, M., et al.) 2025, A&A, 697, A2 [Google Scholar]
  58. Euclid Collaboration (Jahnke, K., et al.) 2025, A&A, 697, A3 [Google Scholar]
  59. Euclid Collaboration (Kovačić, I., et al.) 2025, A&A, 695, A284 [Google Scholar]
  60. Euclid Collaboration (Mellier, Y., et al.) 2025, A&A, 697, A1 [Google Scholar]
  61. Euclid Collaboration (Corcho-Caballero, P., et al.) 2026, A&A, in press, https://doi.org/10.1051/0004-6361/202554582 [Google Scholar]
  62. Euclid Collaboration (Enia, A., et al.) 2026, A&A, in press, https://doi.org/10.1051/0004-6361/202554576 [Google Scholar]
  63. Euclid Collaboration (Tucci, M., et al.) 2026, A&A, in press, https://doi.org/10.1051/0004-6361/202554588 [Google Scholar]
  64. Fabian, A. C. 2012, ARA&A, 50, 455 [Google Scholar]
  65. Fontana, A., Pozzetti, L., Donnarumma, I., et al. 2004, A&A, 424, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Fontanot, F., De Lucia, G., Monaco, P., Somerville, R. S., & Santini, P. 2009, MNRAS, 397, 1776 [NASA ADS] [CrossRef] [Google Scholar]
  67. Fu, H., Shankar, F., Ayromlou, M., et al. 2024, MNRAS, 532, 177 [NASA ADS] [CrossRef] [Google Scholar]
  68. Gallazzi, A., & Bell, E. F. 2009, ApJS, 185, 253 [NASA ADS] [CrossRef] [Google Scholar]
  69. Gallazzi, A., Charlot, S., Brinchmann, J., White, S. D. M., & Tremonti, C. A. 2005, MNRAS, 362, 41 [Google Scholar]
  70. Galliano, F., Nersesian, A., Bianchi, S., et al. 2021, A&A, 649, A18 [EDP Sciences] [Google Scholar]
  71. García-Benito, R., González Delgado, R. M., Pérez, E., et al. 2017, A&A, 608, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Giménez-Arteaga, C., Oesch, P. A., Brammer, G. B., et al. 2023, ApJ, 948, 126 [CrossRef] [Google Scholar]
  73. Giménez-Arteaga, C., Fujimoto, S., Valentino, F., et al. 2024, A&A, 686, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. González Delgado, R. M., Pérez, E., Cid Fernandes, R., et al. 2014, A&A, 562, A47 [CrossRef] [EDP Sciences] [Google Scholar]
  75. González Delgado, R. M., Cid Fernandes, R., Pérez, E., et al. 2016, A&A, 590, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. González Delgado, R. M., Díaz-García, L. A., de Amorim, A., et al. 2021, A&A, 649, A79 [Google Scholar]
  77. Graves, G. J., & Faber, S. M. 2010, ApJ, 717, 803 [NASA ADS] [CrossRef] [Google Scholar]
  78. Groves, B., Dopita, M. A., Sutherland, R. S., et al. 2008, ApJS, 176, 438 [NASA ADS] [CrossRef] [Google Scholar]
  79. Harrison, C. M. 2017, Nat. Astron., 1, 0165 [NASA ADS] [CrossRef] [Google Scholar]
  80. Harvey, T., Conselice, C. J., Adams, N. J., et al. 2025, MNRAS, 542, 2998 [Google Scholar]
  81. Henden, N. A., Puchwein, E., Shen, S., & Sijacki, D. 2018, MNRAS, 479, 5385 [CrossRef] [Google Scholar]
  82. Hunt, L. K., Annibali, F., Cuillandre, J. C., et al. 2025, A&A, 697, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Ibarra-Medel, H. J., Sánchez, S. F., Avila-Reese, V., et al. 2016, MNRAS, 463, 2799 [NASA ADS] [CrossRef] [Google Scholar]
  84. Iglesias-Navarro, P., Huertas-Company, M., Martín-Navarro, I., Knapen, J. H., & Pernet, E. 2024, A&A, 689, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Iglesias-Navarro, P., Huertas-Company, M., Pérez-González, P., et al. 2025, A&A, 703, A229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111 [Google Scholar]
  87. Iyer, K., & Gawiser, E. 2017, ApJ, 838, 127 [NASA ADS] [CrossRef] [Google Scholar]
  88. Iyer, K., Gawiser, E., Davé, R., et al. 2018, ApJ, 866, 120 [Google Scholar]
  89. Iyer, K. G., Gawiser, E., Faber, S. M., et al. 2019, ApJ, 879, 116 [NASA ADS] [CrossRef] [Google Scholar]
  90. Jain, S., Tacchella, S., & Mosleh, M. 2024, MNRAS, 527, 3291 [Google Scholar]
  91. Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22 [NASA ADS] [CrossRef] [Google Scholar]
  92. Jones, A. P., Köhler, M., Ysard, N., Bocchio, M., & Verstraete, L. 2017, A&A, 602, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Kaushal, Y., Nersesian, A., Bezanson, R., et al. 2024, ApJ, 961, 118 [NASA ADS] [CrossRef] [Google Scholar]
  94. Koposov, S., Speagle, J., Barbary, K., et al. 2022, https://doi.org/10.5281/zenodo.7388523 [Google Scholar]
  95. Kriek, M., & Conroy, C. 2013, ApJ, 775, L16 [NASA ADS] [CrossRef] [Google Scholar]
  96. Labbé, I., Huang, J., Franx, M., et al. 2005, ApJ, 624, L81 [Google Scholar]
  97. Lapiner, S., Dekel, A., Freundlich, J., et al. 2023, MNRAS, 522, 4515 [NASA ADS] [CrossRef] [Google Scholar]
  98. Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints [arXiv:1110.3193] [Google Scholar]
  99. Leitherer, C., Li, I.-H., Calzetti, D., & Heckman, T. M. 2002, ApJS, 140, 303 [Google Scholar]
  100. Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170 [NASA ADS] [CrossRef] [Google Scholar]
  101. Leja, J., Johnson, B. D., Conroy, C., & van Dokkum, P. 2018, ApJ, 854, 62 [NASA ADS] [CrossRef] [Google Scholar]
  102. Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019a, ApJ, 876, 3 [Google Scholar]
  103. Leja, J., Johnson, B. D., Conroy, C., et al. 2019b, ApJ, 877, 140 [NASA ADS] [CrossRef] [Google Scholar]
  104. Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119 [NASA ADS] [CrossRef] [Google Scholar]
  105. López Fernández, R., González Delgado, R. M., Pérez, E., et al. 2018, A&A, 615, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Lower, S., Narayanan, D., Leja, J., et al. 2020, ApJ, 904, 33 [NASA ADS] [CrossRef] [Google Scholar]
  107. Maraston, C., Pforr, J., Renzini, A., et al. 2010, MNRAS, 407, 830 [NASA ADS] [CrossRef] [Google Scholar]
  108. Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113 [NASA ADS] [Google Scholar]
  109. Morishita, T., Abramson, L. E., Treu, T., et al. 2019, ApJ, 877, 141 [NASA ADS] [CrossRef] [Google Scholar]
  110. Mosleh, M., Hosseinnejad, S., Hosseini-ShahiSavandi, S. Z., & Tacchella, S. 2020, ApJ, 905, 170 [NASA ADS] [CrossRef] [Google Scholar]
  111. Mosleh, M., Riahi-Zamin, M., & Tacchella, S. 2025, ApJ, 983, 181 [Google Scholar]
  112. Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206 [Google Scholar]
  113. Nanni, L., Thomas, D., Trayford, J., et al. 2022, MNRAS, 515, 320 [CrossRef] [Google Scholar]
  114. Nanni, L., Thomas, D., Trayford, J., et al. 2023, MNRAS, 522, 5479 [NASA ADS] [CrossRef] [Google Scholar]
  115. Nanni, L., Neumann, J., Thomas, D., et al. 2024, MNRAS, 527, 6419 [Google Scholar]
  116. Narayanan, D., Lower, S., Torrey, P., et al. 2024, ApJ, 961, 73 [NASA ADS] [CrossRef] [Google Scholar]
  117. Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624 [Google Scholar]
  118. Nelson, D., Pillepich, A., Springel, V., et al. 2019, MNRAS, 490, 3234 [Google Scholar]
  119. Nersesian, A., van der Wel, A., Gallazzi, A., et al. 2024, A&A, 681, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  120. Nersesian, A., van der Wel, A., Gallazzi, A. R., et al. 2025, A&A, 695, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  121. Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43 [CrossRef] [Google Scholar]
  122. Noll, S., Burgarella, D., Giovannoli, E., et al. 2009, A&A, 507, 1793 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  123. Ocvirk, P., Pichon, C., Lançon, A., & Thiébaut, E. 2006, MNRAS, 365, 46 [Google Scholar]
  124. Pacifici, C., Iyer, K. G., Mobasher, B., et al. 2023, ApJ, 944, 141 [NASA ADS] [CrossRef] [Google Scholar]
  125. Papovich, C., Dickinson, M., & Ferguson, H. C. 2001, ApJ, 559, 620 [NASA ADS] [CrossRef] [Google Scholar]
  126. Papovich, C., Finkelstein, S. L., Ferguson, H. C., Lotz, J. M., & Giavalisco, M. 2011, MNRAS, 412, 1123 [NASA ADS] [Google Scholar]
  127. Parikh, T., Thomas, D., Maraston, C., et al. 2021, MNRAS, 502, 5508 [NASA ADS] [CrossRef] [Google Scholar]
  128. Paspaliaris, E. D., Xilouris, E. M., Nersesian, A., et al. 2023, A&A, 669, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  129. Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
  130. Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
  131. Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
  132. Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
  133. Pessa, I., Schinnerer, E., Sanchez-Blazquez, P., et al. 2023, A&A, 673, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  134. Pettitt, A. R., Dobbs, C. L., Baba, J., et al. 2020, MNRAS, 498, 1159 [Google Scholar]
  135. Pforr, J., Maraston, C., & Tonini, C. 2012, MNRAS, 422, 3285 [NASA ADS] [CrossRef] [Google Scholar]
  136. Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077 [Google Scholar]
  137. Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196 [Google Scholar]
  138. Pozzetti, L., & Mannucci, F. 2000, MNRAS, 317, L17 [NASA ADS] [CrossRef] [Google Scholar]
  139. Pozzetti, L., Bolzonella, M., Lamareille, F., et al. 2007, A&A, 474, 443 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  140. Saifollahi, T., Voggel, K., Lançon, A., et al. 2025, A&A, 697, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  141. Salim, S., Lee, J. C., Janowiecki, S., et al. 2016, ApJS, 227, 2 [NASA ADS] [CrossRef] [Google Scholar]
  142. Salim, S., Boquien, M., & Lee, J. C. 2018, ApJ, 859, 11 [Google Scholar]
  143. Sánchez-Blázquez, P., Peletier, R. F., Jiménez-Vicente, J., et al. 2006, MNRAS, 371, 703 [Google Scholar]
  144. Sánchez-Blázquez, P., Rosales-Ortega, F. F., Méndez-Abreu, J., et al. 2014, A&A, 570, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  145. Sawicki, M., & Yee, H. K. C. 1998, AJ, 115, 1329 [NASA ADS] [CrossRef] [Google Scholar]
  146. Schawinski, K., Urry, C. M., Simmons, B. D., et al. 2014, MNRAS, 440, 889 [Google Scholar]
  147. Sellwood, J. A., & Carlberg, R. G. 2019, MNRAS, 489, 116 [CrossRef] [Google Scholar]
  148. Semenov, V. A., Kravtsov, A. V., & Gnedin, N. Y. 2021, ApJ, 918, 13 [NASA ADS] [CrossRef] [Google Scholar]
  149. Shapley, A. E., Steidel, C. C., Adelberger, K. L., et al. 2001, ApJ, 562, 95 [NASA ADS] [CrossRef] [Google Scholar]
  150. Shin, E.-J., Tacchella, S., Kim, J.-H., Iyer, K. G., & Semenov, V. A. 2023, ApJ, 947, 61 [Google Scholar]
  151. Silva, L., Granato, G. L., Bressan, A., & Danese, L. 1998, ApJ, 509, 103 [Google Scholar]
  152. Siudek, M., Pucha, R., Mezcua, M., et al. 2024, A&A, 691, A308 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  153. Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24 [Google Scholar]
  154. Skilling, J. 2004, in Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, eds. R. Fischer, R. Preuss, & U. V. Toussaint, AIP Conf. Ser., 735, 395 [NASA ADS] [CrossRef] [Google Scholar]
  155. Skilling, J. 2006, in Bayesian Inference and Maximum Entropy Methods In Science and Engineering, ed. A. Mohammad-Djafari, AIP Conf. Ser., 872, 321 [NASA ADS] [CrossRef] [Google Scholar]
  156. Smith, M. V., van Zee, L., Dale, D. A., et al. 2022, MNRAS, 515, 3270 [NASA ADS] [CrossRef] [Google Scholar]
  157. Snyder, G. F., Cox, T. J., Hayward, C. C., Hernquist, L., & Jonsson, P. 2011, ApJ, 741, 77 [NASA ADS] [CrossRef] [Google Scholar]
  158. Sorba, R., & Sawicki, M. 2015, MNRAS, 452, 235 [NASA ADS] [CrossRef] [Google Scholar]
  159. Sorba, R., & Sawicki, M. 2018, MNRAS, 476, 1532 [NASA ADS] [CrossRef] [Google Scholar]
  160. Springel, V., Di Matteo, T., & Hernquist, L. 2005, ApJ, 620, L79 [Google Scholar]
  161. Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676 [Google Scholar]
  162. Strateva, I., Ivezić, Ž., Knapp, G. R., et al. 2001, AJ, 122, 1861 [CrossRef] [Google Scholar]
  163. Suess, K. A., Kriek, M., Price, S. H., & Barro, G. 2019, ApJ, 885, L22 [NASA ADS] [CrossRef] [Google Scholar]
  164. Suess, K. A., Leja, J., Johnson, B. D., et al. 2022, ApJ, 935, 146 [NASA ADS] [CrossRef] [Google Scholar]
  165. Tacchella, S., Forbes, J. C., & Caplar, N. 2020, MNRAS, 497, 698 [Google Scholar]
  166. Tacchella, S., Conroy, C., Faber, S. M., et al. 2022a, ApJ, 926, 134 [NASA ADS] [CrossRef] [Google Scholar]
  167. Tacchella, S., Finkelstein, S. L., Bagley, M., et al. 2022b, ApJ, 927, 170 [NASA ADS] [CrossRef] [Google Scholar]
  168. Terrazas, B. A., Bell, E. F., Woo, J., & Henriques, B. M. B. 2017, ApJ, 844, 170 [NASA ADS] [CrossRef] [Google Scholar]
  169. Thorne, J. E., Robotham, A. S. G., Davies, L. J. M., et al. 2021, MNRAS, 505, 540 [NASA ADS] [CrossRef] [Google Scholar]
  170. Tojeiro, R., Heavens, A. F., Jimenez, R., & Panter, B. 2007, MNRAS, 381, 1252 [NASA ADS] [CrossRef] [Google Scholar]
  171. Trager, S. C., Faber, S. M., Worthey, G., & González, J. J. 2000, AJ, 120, 165 [NASA ADS] [CrossRef] [Google Scholar]
  172. Trager, S. C., Faber, S. M., & Dressler, A. 2008, MNRAS, 386, 715 [Google Scholar]
  173. Trčka, A., Baes, M., Camps, P., et al. 2022, MNRAS, 516, 3728 [CrossRef] [Google Scholar]
  174. van der Wel, A., Bezanson, R., D’Eugenio, F., et al. 2021, ApJS, 256, 44 [NASA ADS] [CrossRef] [Google Scholar]
  175. van Dokkum, P. G. 2008, ApJ, 674, 29 [Google Scholar]
  176. Walcher, J., Groves, B., Budavári, T., & Dale, D. 2011, Ap&SS, 331, 1 [NASA ADS] [CrossRef] [Google Scholar]
  177. Wellons, S., Torrey, P., Ma, C.-P., et al. 2015, MNRAS, 449, 361 [Google Scholar]
  178. Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, L29 [Google Scholar]
  179. Whitler, L., Stark, D. P., Endsley, R., et al. 2023, MNRAS, 519, 5859 [NASA ADS] [CrossRef] [Google Scholar]
  180. Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879 [NASA ADS] [CrossRef] [Google Scholar]
  181. Worthey, G., Faber, S. M., Gonzalez, J. J., & Burstein, D. 1994, ApJS, 94, 687 [Google Scholar]
  182. Wuyts, S., Labbé, I., Franx, M., et al. 2007, ApJ, 655, 51 [NASA ADS] [CrossRef] [Google Scholar]
  183. Zabel, N., Davis, T. A., Smith, M. W. L., et al. 2021, MNRAS, 502, 4723 [NASA ADS] [CrossRef] [Google Scholar]
  184. Zheng, Z., Wang, H., Ge, J., et al. 2017, MNRAS, 465, 4572 [Google Scholar]
  185. Zhuang, Y., Leaman, R., van de Ven, G., et al. 2019, MNRAS, 483, 1862 [NASA ADS] [CrossRef] [Google Scholar]
  186. Zibetti, S., Charlot, S., & Rix, H.-W. 2009, MNRAS, 400, 1181 [NASA ADS] [CrossRef] [Google Scholar]
  187. Zolotov, A., Dekel, A., Mandelker, N., et al. 2015, MNRAS, 450, 2327 [Google Scholar]
3

The bias in the recovered stellar mass for each method is discussed in Sect. 4.3.

Appendix A: The effect of prior functions on stellar metallicities

In this section, we test the effect of the prior function on the stellar metallicity estimates. In a recent study, Nersesian et al. (2025) showed that when fitting photometry alone, a flat metallicity prior in linear space better recovers the distribution of metallicities measured from the spectrophotometric analysis of the LEGA-C survey (van der Wel et al. 2021). In fact when fitting spectroscopy alone or together with photometry, the results are unaffected by the choice of the Z★ prior (i.e., whether in linear or logarithmic space). In Fig. A.1, we demonstrate that using a flat prior in linear space for Z★ results in measured metallicities that better align with those from our true values from TNG50. For this comparison, all stellar metallicities have been rescaled to the Solar metallicity value of Z⊙ = 0.02. When a flat prior in linear space is used for Z★, the mean offset (μ) from the fiducial Z★ is considerably reduced, from 0.39 dex to 0.15 dex. The scatter (σ) around the mean remains similar yet slightly reduced by 0.02 dex. This improvement in the Z★ measurements does not necessarily imply a better constraint on the stellar metallicities. Nevertheless, Nersesian et al. (2024) showed that stellar metallicities are difficult to be recovered from multiwavelength photometry, even when there is information from both broadbands and narrowbands.

Thumbnail: Fig. A.1. Refer to the following caption and surrounding text. Fig. A.1.

Comparison of the spatially resolved stellar metallicities for 6958 pixel bins in the 25 sample galaxies measured by different prior functions. The true values of Z★ (x-axis) are compared with those from our fiducial run (y-axis), when a flat prior for Z★ is used in logarithmic space (left panel) and in linear space (right panel). The dashed orange line shows the one-to-one relation, while contours enclose 20%, 50%, and 80% of the total data. Top-left: Residual distribution, log10(Z★/Z★,  TNG50). Top-right: Comparison of the true Z★ values (teal distribution) with those recovered from Prospector (orange distribution). The dotted lines correspond to the prior probability distributions, scaled with an arbitrary factor. The mean offset (μ) and scatter (σ) between the two distributions are indicated in each panel. A flat prior in linear space for Z★ returns values closer to the ground truth.

Appendix B: The benefit of Euclid NIR photometry

Similar to Fig. 5, Fig. B.1 presents the distribution of all spatial pixel bins across the 25 galaxies in our sample, but this time derived from SED fitting performed without the Euclid photometric bands. The mean offset (μ) and scatter (standard deviation, σ) of the logarithmic ratio between the best-fit parameters and the ground truth, along with the Spearman rank-order correlation coefficient (ρ) for each stellar property, are summarized in the top-left corner of each panel in Fig. B.1. Points are color-coded according to the average trend of the dust optical depth in the V band ( τ ̂ dust , 2 Mathematical equation: $ \hat{\tau}_{\mathrm{dust,\, 2}} $), calculated with the LOESS method.

Thumbnail: Fig. B.1. Refer to the following caption and surrounding text. Fig. B.1.

Comparisons of the spatially resolved stellar population properties derived from spatially resolved SED fitting with a nonparametric SFH model on seven broadbands. The photometric bands used in this analysis include the GALEX FUV and NUV, and LSST u, g, r, i, and z, in the resolution of GALEX. From left to right, we show the pixel values of the stellar property maps, that is Σ★, ΣSFR, Z★, and t★,  mw. Points are color-coded according to the average trend of the dust optical depth in the V band ( τ ̂ dust , 2 Mathematical equation: $ \hat{\tau}_{\mathrm{dust,\,2}} $), calculated with the LOESS method. Contours enclose 20%, 50%, and 80% of the total data. The dashed orange line shows the one-to-one relation. The Spearman’s correlation coefficient (ρ), measured bias (average offset μ), and scatter (standard deviation σ) relative to the ground truth from TNG50 are indicated in the legend of each panel.

Excluding the Euclid bands introduces a systematic offset of 0.05 dex in Σ★ relative to the ground truth, and an increased scatter around the mean (σ = 0.16 dex), indicating that the Euclid bands play a crucial role in obtaining reliable stellar mass estimates. For ΣSFR, the recovery worsens significantly, with a larger mean offset (μ = 0.25 dex) and increased scatter (σ = 0.41 dex). As seen in Fig. 5, the lower end of the ΣSFR distribution tends to be overestimated, but this effect is even more pronounced in the absence of the Euclid bands. Additionally, we observe a stronger correlation between ΣSFR and τ ̂ dust , 2 Mathematical equation: $ \hat{\tau}_{\mathrm{dust,\, 2}} $, suggesting that the dust-age degeneracy becomes more severe when NIR information is missing.

Furthermore, the exclusion of Euclid photometry results in a completely unconstrained Z★, whereas the recovery of stellar ages remains largely unaffected. We measure an offset (μ = −0.03 dex) and scatter (σ = 0.09 dex), consistent with the results in Fig. 5. However, we again observe the formation of a secondary peak in the stellar age distribution, which is linked to the dust-age degeneracy. In these cases, Prospector infers younger ages, attributing the red colors primarily to dust attenuation rather than the presence of older stellar populations.

All Tables

Table 1.

Free parameters and their associated prior distribution functions in the Prospector physical model.

In the text
Table 2.

Definitions and nomenclature of the different SFH models in our analysis.

In the text
Table 3.

Mean offset and standard deviation of formation times and global properties from the TNG50 reference for each SFH model.

In the text

All Figures

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

SFR as a function of M★ for the TNG50-SKIRT Atlas galaxies at z = 0. Galaxies are separated as star-forming (blue stars) and quiescent (red points) based on their sSFR. The dashed orange line indicates this separation between star formation and quiescence (sSFR = 10−11 yr−1). The black squares indicate the selected sample of 25 galaxies used in our analysis.
In the text
Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Example of a simulated galaxy processed with piXedfit. We show the synthetic images of GALEX, LSST, and Euclid for TNG275545 with the O1 orientation index. Observational effects were introduced to the original synthetic images from the TNG50-SKIRT Atlas database in the form of simulated noise and convolution with the PSF of the corresponding cameras. The FWHM of the PSF as well as the pixel size of the images are 5 . ″ 05 Mathematical equation: $ 5{{\overset{\prime\prime}{.}}}05 $ and 1 . ″ 5 Mathematical equation: $ 1{{\overset{\prime\prime}{.}}}5 $, respectively. The pixels are binned with piXedfit to increase the S/N of the spatially resolved SEDs.
In the text
Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

Examples of the stellar property maps of five TNG50-SKIRT Atlas galaxies, obtained from the spatially resolved SED fitting of GALEX, LSST, and Euclid photometry with Prospector. The maps include the stellar mass surface density (Σ★), the SFR surface density (ΣSFR), stellar metallicity (Z★), and mass-weighted stellar age (t★,  mw). Next to each galaxy’s sub-halo ID, we indicate whether a galaxy is classified as star-forming (SF) or quiescent (Q), based on its global sSFR.
In the text
Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

Examples of residual maps of the stellar properties of two TNG50-SKIRT Atlas galaxies, obtained from the nonparametric SFH configuration. The residual maps include the stellar mass surface density (Σ★), the SFR surface density (ΣSFR), stellar metallicity (Z★), and mass-weighted stellar age (t★,  mw). Each residual map is computed as the logarithmic ratio of the inferred property from Prospector to the ground truth value from TNG50, expressed as log10(Prospector/TNG50).
In the text
Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

Comparisons of the spatially resolved stellar population properties for 6958 pixel bins in the 25 sample galaxies derived from spatially resolved SED fitting with a nonparametric SFH model. The depicted physical properties include the stellar mass surface density (Σ★), the SFR surface density (ΣSFR), the stellar metallicity (Z★), and the mass-weighted stellar age (t★,  mw). The colors indicate the average trend with the dust optical depth in the V band ( τ ̂ dust , 2 Mathematical equation: $ \hat{\tau}_{\mathrm{dust,\, 2}} $), calculated with the LOESS method. Each hexbin includes a minimum of five points. Contours enclose 20%, 50%, and 80% of the total data. The dashed orange line represents the one-to-one relation. The Spearman’s correlation coefficient (ρ), measured bias (average offset μ), and scatter (standard deviation σ) relative to the ground truth from TNG50 are indicated in the legend of each panel.
In the text
Thumbnail: Fig. 6. Refer to the following caption and surrounding text. Fig. 6.

Joint posterior distributions of the main physical properties in our analysis, the best-fit SED, and the recovered SFH for a selected pixel bin of TNG414917 with the O4 orientation index. The corner plot displays the posterior distributions of τ ̂ dust , 2 Mathematical equation: $ \hat{\tau}_{\mathrm{ dust,\, 2}} $, sSFR, t★,  mw, Z★, and Σ★. Contours enclose the 20%, 50%, and 80% of the total data, while the vertical dashed black lines mark the median values for each parameter. The median and the 16th–84th percentile range for each physical property are given at the top of the corresponding posterior distribution. The Spearman’s correlation coefficient, ρ, quantifies the strength of correlation between two posterior distributions. Top right: Comparison of the input photometric data (red points) to the best-fit SED model (blue line). Bottom right: Subpanel illustrating the absolute residuals defined as ΔFν = log10(fit/ground truth). The photometric uncertainties are not visible in the plot because they are smaller than the symbol size. Bottom: Posterior distribution of the SFH.
In the text
Thumbnail: Fig. 7. Refer to the following caption and surrounding text. Fig. 7.

Examples of SFR surface density (ΣSFR) maps and SFH models for ten galaxies in our TNG50-SKIRT Atlas sample. Left panel: ΣSFR map for each galaxy, as estimated with Prospector. Right panel: Comparison of various global SFH models to the ground truth SFH from TNG50. Specifically, we show the global nonparametric SFH inferred from our fiducial spatially resolved run (SFHres,  np; orange), the global parametric SFH inferred from the spatially resolved map (SFHres,  τ; pink), the nonparametric SFH derived from fitting the integrated photometry (SFHglob,  np; yellow), the global parametric SFH based on a simple τ-model (SFHglob,  τ; blue), and the true global SFH from the TNG50 resolved map (SFHres,  TNG50; green).
In the text
Thumbnail: Fig. 8. Refer to the following caption and surrounding text. Fig. 8.

Comparison of the normalized cumulative stellar mass growth curves of all 25 galaxies in our sample. We plot the cumulative stellar mass growth curves from the global nonparametric SFH inferred from our fiducial spatially resolved run (SFHres,  np; orange), the global parametric SFH inferred from the spatially resolved map (SFHres,  τ; pink), the nonparametric SFH derived from fitting the integrated photometry (SFHglob,  np; yellow), the global parametric SFH based on a simple τ-model (SFHglob,  τ; blue), and the true global SFH from the TNG50 resolved map (SFHres,  TNG50; green). The dotted lines indicate the different formation times by which 20%, 50%, and 80% of a galaxy’s stellar mass was formed. Both nonparametric SFHs (SFHres,  np and SFHglob,  np) seem to be in better agreement with SFHres,  TNG50, while most of the parametric SFHs (both SFHres,  τ and SFHglob,  τ) are skewed toward a later cosmic time, thus lacking the contribution from the older stellar populations.
In the text
Thumbnail: Fig. 9. Refer to the following caption and surrounding text. Fig. 9.

Comparison of the formation timescales for the four different SFH models with respect to the ground truth values from the TNG50 simulation. Each marker represents the formation time of a galaxy in our sample, color-coded according to the corresponding SFH model. Star-forming galaxies are indicated with stars, while quiescent galaxies with points. From left to right: Time for 20% (t20), 50% (t50), and 80% (t80) of a galaxy’s stellar mass formation. The four models include the spatially resolved model from pixel-by-pixel SED fitting with a nonparametric (orange) and parametric (pink) SFH, the integrated SED fitting with a nonparametric SFH (yellow), and a simple τ model within Prospector (blue). The measured bias (average offset μ) and scatter (standard deviation σ) relative to the ground truth from TNG50 are indicated in the legends of the panels. The formation times inferred from the spatially resolved and integrated nonparametric SFHs are in best agreement with the ground truth values.
In the text
Thumbnail: Fig. 10. Refer to the following caption and surrounding text. Fig. 10.

Differences in the global stellar properties for the four different SFH models with respect to the ground truth global values from the TNG50 simulation. From left to right: M★, sSFR, and t★,  mw. The properties were derived through SED fitting under four different configurations: spatially resolved fitting with a nonparametric SFHres,  np (orange), spatially resolved SED fitting with a simple τ-model SFHres,  τ (pink), global fitting with a nonparametric SFHglob,  np (yellow), and a global fitting with a simple τ-model SFHglob,  τ (blue). Typical uncertainties for each configuration are shown in the bottom-left corner of each panel. Star-forming galaxies are indicated with stars, while quiescent galaxies with points. The bottom panels show the residual distributions of each stellar property. The measured bias (average offset) and scatter (standard deviation) are indicated in the legends of the panels. The properties inferred from the spatially resolved and integrated nonparametric SFHs are in best agreement with the ground truth values.
In the text
Thumbnail: Fig. A.1. Refer to the following caption and surrounding text. Fig. A.1.

Comparison of the spatially resolved stellar metallicities for 6958 pixel bins in the 25 sample galaxies measured by different prior functions. The true values of Z★ (x-axis) are compared with those from our fiducial run (y-axis), when a flat prior for Z★ is used in logarithmic space (left panel) and in linear space (right panel). The dashed orange line shows the one-to-one relation, while contours enclose 20%, 50%, and 80% of the total data. Top-left: Residual distribution, log10(Z★/Z★,  TNG50). Top-right: Comparison of the true Z★ values (teal distribution) with those recovered from Prospector (orange distribution). The dotted lines correspond to the prior probability distributions, scaled with an arbitrary factor. The mean offset (μ) and scatter (σ) between the two distributions are indicated in each panel. A flat prior in linear space for Z★ returns values closer to the ground truth.
In the text
Thumbnail: Fig. B.1. Refer to the following caption and surrounding text. Fig. B.1.

Comparisons of the spatially resolved stellar population properties derived from spatially resolved SED fitting with a nonparametric SFH model on seven broadbands. The photometric bands used in this analysis include the GALEX FUV and NUV, and LSST u, g, r, i, and z, in the resolution of GALEX. From left to right, we show the pixel values of the stellar property maps, that is Σ★, ΣSFR, Z★, and t★,  mw. Points are color-coded according to the average trend of the dust optical depth in the V band ( τ ̂ dust , 2 Mathematical equation: $ \hat{\tau}_{\mathrm{dust,\,2}} $), calculated with the LOESS method. Contours enclose 20%, 50%, and 80% of the total data. The dashed orange line shows the one-to-one relation. The Spearman’s correlation coefficient (ρ), measured bias (average offset μ), and scatter (standard deviation σ) relative to the ground truth from TNG50 are indicated in the legend of each panel.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.