© The Authors 2025

1. Introduction

Multi-wavelength photometric and spectroscopic surveys over the past few decades have greatly contributed to our understanding of galaxy formation and evolution. Some major galaxy surveys (including the Sloan Digital Sky Survey, SDSS; York et al. 2000; The 2dF Galaxy Redshift Survey; Colless et al. 2001; Galaxy and Mass Assembly, GAMA; Baldry et al. 2010, Driver et al. 2022) have observed many galaxies and provided important insights into the study of galaxy evolution. The depth and high spatial resolution of the photometric instruments of galaxy surveys, such as the Hubble Space Telescope (HST), the Dark Energy Survey (DES; Dark Energy Survey Collaboration 2016), the Hyper Suprime-Cam SSP Survey (HSC; Aihara et al. 2018), and the James Webb Space Telescope (JWST; Gardner et al. 2023), have enabled detailed observations of galaxies up to high redshifts, probing further back in the history of galaxy evolution.

The current and upcoming galaxy surveys (e.g. the Vera Rubin Observatory Legacy Survey of Space and Time, LSST; Ivezić et al. 2019, the Dark Energy Spectroscopic Instrument, DESI; DESI Collaboration 2016, and the Nancy Grace Roman Space Telescope; Spergel et al. 2015) provide even more data, opening a new era of big data in astronomy. The abundance of astronomical data requires sophisticated techniques for analysing and interpreting them. One of the main ways of measuring the physical properties of galaxies from their light is through so-called spectral energy distribution (SED) fitting (e.g. Sawicki & Yee 1998; Bolzonella et al. 2000; Sawicki 2012; Walcher et al. 2011; Conroy 2013; Pacifici et al. 2023). Major developments have been made in this technique over the last two decades, thanks to advancements in both the modelling of SEDs and the statistical methods of fitting data with the models. In the former aspect, we have seen the inclusion of more complex and physically realistic components in the models. This has enabled the modelling of a wide range of physical processes that occur in galaxies, making it possible to simulate galaxy spectra across a wide range of wavelengths, from X-ray to radio (e.g. Bruzual & Charlot 2003; Conroy et al. 2009; da Cunha et al. 2008; Boquien et al. 2019; Yang et al. 2022). As for statistical techniques, we have seen the emergence of the Bayesian fitting technique with sophisticated posterior sampling methods, such as the Markov chain Monte Carlo (MCMC) and nested sampling (e.g. Acquaviva et al. 2011; Han & Han 2014; Leja et al. 2017; Carnall et al. 2018; Abdurro’uf et al. 2021; Pacifici et al. 2023) and even the applications of machine learning in SED fitting (e.g. Lovell et al. 2019; Dobbels et al. 2020; Gilda et al. 2021; Hahn & Melchior 2022; Euclid Collaboration: Bisigello et al. 2023; Euclid Collaboration: Kovačić et al. 2025).

Despite the abundance of data, most studies on galaxies have focused on the global integrated properties. In the case of spectroscopic observations, a galaxy-integrated spectrum is obtained with spectroscopy (e.g. SDSS, GAMA). In the case of photometric observations, the integrated SED of a galaxy is usually measured within an aperture of a certain radius centred on the galaxy. These analyses have revealed many important properties of galaxy populations across a wide range of redshifts and provided significant insights into our understanding of galaxy structure, formation, and evolution (e.g. Tremonti et al. 2004; Brinchmann et al. 2004; Salucci et al. 2008; Blanton & Moustakas 2009). However, these studies overlook the crucial insights that spatially resolved data can offer into the physical processes within galaxies.

The integral field spectroscopy (IFS) surveys in the past few decades have revolutionized the study of galaxy evolution. Some major IFS surveys are SAURON (de Zeeuw et al. 2002), ATLAS^3D (Cappellari et al. 2011), the Calar Alto Legacy Integral Field Area survey (CALIFA; Sánchez et al. 2012), Mapping nearby Galaxies at Apache Point Observatory (MaNGA; Bundy et al. 2015; Abdurro’uf et al. 2022a), Sydney-AAO Multi-object Integral field spectrograph (SAMI; Croom et al. 2012), and KMOS^3D (Wisnioski et al. 2015). These surveys have enabled the study of galaxies in great detail, thanks to the ability to observe the spectra of spatially resolved regions inside them. One of the important findings from these spatially resolved studies is that key scaling relations among galaxies on a global scale (e.g. the star-forming main sequence, Brinchmann et al. 2004; Noeske et al. 2007; Speagle et al. 2014; Popesso et al. 2023; Euclid Collaboration: Enia et al. 2025; the mass-metallicity relation, Tremonti et al. 2004) likely originate from similar scaling relations on kiloparsec and sub-kiloparsec scales (see e.g. Sánchez 2020).

While the IFS surveys arguably provide the most detailed insights into the properties of each spatially resolved element, to date only a limited number of them have observed a large number of galaxies over a wide area of the sky. Those surveys mostly focus on the local Universe (e.g. CALIFA, MaNGA, and SAMI). The largest IFS survey to date is the MaNGA survey, which observed around 10 000 galaxies at 0.01 < z < 0.15. Unfortunately, IFS surveys tend to have a lower spatial resolution than imaging data, and unlike imaging data, the integral field units usually do not cover the entire optical region of galaxies.

Recent improvements in the depth, spatial resolution, and wavelength coverage of imaging surveys provide opportunities for conducting SED fitting on spatially resolved scales in 2D images of galaxies, namely the spatially resolved SED fitting method. Previous studies have shown the great potential of this method in spatially resolving stellar population properties in galaxies across a wide range of redshifts (e.g. Conti et al. 2003; Zibetti et al. 2009; Welikala et al. 2011; Wuyts et al. 2012, 2013; Viaene et al. 2014; Sorba & Sawicki 2018; Abdurro’uf & Akiyama 2017; Abdurro’uf & Akiyama 2018; Robotham et al. 2022; Abdurro’uf et al. 2022b; Abdurro’uf et al. 2022c; Abdurro’uf et al. 2023; Bellstedt et al. 2024). Wuyts et al. (2012, 2013) applied this method to 0.5 < z < 2.5 star-forming galaxies in the GOODS-South field using imaging data from HST. They studied the radial gradients of rest-frame colours, surface density of stellar mass (Σ_∗), age of the stellar population, and dust attenuation. Abdurro’uf et al. (2022b,c) utilized multi-wavelength imaging data across more than 20 bands, ranging from the far-ultraviolet (FUV) to the far-infrared (FIR), obtained from various telescopes, to create maps of stellar populations and dust properties of galaxies on kiloparsec scales. Recently, some works have employed JWST data to study the stellar populations of galaxies at very high redshifts (up to z = 9, e.g. Abdurro’uf et al. 2023; Giménez-Arteaga et al. 2023). The abundance of high-quality (i.e. deep and high spatial resolution) multi-wavelength imaging data from the current and upcoming surveys opens up avenues for the spatially resolved SED analysis of a large number of galaxies. The wide wavelength range of modern data can provide coverage from the rest-frame UV to the near-infrared (NIR), which is essential for SED fitting across a wide redshift range. Although it does not replace spectroscopy, it offers precious complementary constraints on stellar population properties.

A major step forwards in large galaxy surveys is the European Space Agency’s Euclid mission (Laureijs et al. 2011; Euclid Collaboration: Mellier et al. 2025), launched on July 1st, 2023, which will observe approximately 14 000 deg² of extragalactic sky within the six years of its operation. It will deliver high-quality imaging data for about 1.5 billion galaxies at an optical wavelength with the visible imaging instrument (VIS; Cropper et al. 2016; Euclid Collaboration: Cropper et al. 2025; passband I_E) and in the NIR with the NISP instrument (the Near-Infrared Spectrometer and Photometer, Maciaszek et al. 2022; Euclid Collaboration: Schirmer et al. 2022; Euclid Collaboration: Jahnke et al. 2025; Euclid Collaboration: Hormuth et al. 2025; passbands Y_E, J_E, and H_E). Euclid will provide an unprecedented homogeneous view of the local Universe, re-observing thousands of known galaxies at a high spatial resolution (∼1 kpc), as well as newly discovered low-surface brightness galaxies. The recent early-release observations have highlighted the exceptional potential of Euclid data (e.g. Cuillandre et al. 2025; Hunt et al. 2025; Marleau et al. 2025; Saifollahi et al. 2025).

For a wider wavelength coverage, the Euclid data can be combined with other wide-area and deep imaging surveys in the optical, such as the future Rubin/LSST (Ivezić et al. 2019) that will observe the southern hemisphere, and the UNIONS survey (Ibata et al. 2017), a combined effort by the CFHT, Pan-STARRS, and Subaru telescopes, which observe the northern hemisphere. In combination, these datasets extend from the near-UV to the NIR. For galaxies that are covered by the GALEX Medium Imaging Survey (MIS) or Deep Imaging Survey (DIS), we will incorporate GALEX images to get even wider wavelength coverage extending to FUV. The wide-area coverage will make it possible to analyse many galaxies residing in a wide range of local density environments (from a low-density field to a dense galaxy cluster) that enable a more comprehensive analysis of the effects of internal and external factors in the evolution of galaxies.

In this work, we test a pipeline designed for spatially resolved SED fitting of local galaxies using imaging data primarily from Euclid, complemented by optical data from surveys such as LSST and UNIONS, as well as UV data from GALEX’s medium- and deep-imaging surveys. The pipeline mainly utilizes various functionalities from piXedfit package (Abdurro’uf et al. 2021; Abdurro’uf et al. 2022d). For this testing, we use synthetic images generated from the TNG50 cosmological magneto-hydrodynamical simulations (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2019; Springel et al. 2018). We aim to test to what extent Euclid imaging data (with 4 bands: I_E, Y_E, J_E, H_E) alone can provide constraints on some basic but important stellar population properties (stellar mass, age, and metallicity) and how additional imaging data in the optical and UV can improve the fitting results.

This paper is structured as follows. Section 2 describes the procedures for generating synthetic images of Euclid, LSST, and GALEX. The analysis pipeline and the SED fitting procedures are described in Sect. 3. In Sect. 4 we present our results and further discuss them in Sect. 5, focusing on some factors that drive the SED fitting output. Finally, we provide our conclusions and outlooks in Sect. 6. Throughout this paper, we assume the Chabrier (2003) initial mass function (IMF) with a mass range of 0.1 − 100 M_⊙ and cosmological parameters of Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹.

2. Generating synthetic images of Euclid, LSST, and GALEX

In this section, we first describe the synthetic images from the TNG50-SKIRT Atlas generated using the radiative transfer method (Sect. 2.1). We then describe the procedures for adding some observational factors to the synthetic images in Sect. 2.2.

2.1. Synthetic images from the TNG50-SKIRT Atlas

We use the synthetic images from the TNG50-SKIRT Atlas data products (Baes et al. 2024a,b; Euclid Collaboration: Kovačić et al. 2025) that were generated from the IllustrisTNG cosmological magneto-hydrodynamical simulations (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2019; Springel et al. 2018) using the SKIRT radiative transfer code (Camps & Baes 2015, 2020). The synthetic images were created for galaxies in the TNG50 simulation (Nelson et al. 2019), which has the highest resolution among the three IllustrisTNG runs with a baryonic mass resolution of 8.5 × 10⁴ M_⊙, though covering the smallest cubic volume of 51.7 comoving Mpc on the side. The TNG50 simulation uses the moving-mesh code AREPO (Springel 2010). It includes a wide range of physical processes in the galaxies, including gas heating and cooling, stochastic star formation, stellar evolution, chemical enrichment in the interstellar medium (ISM), seeding and growth of the supermassive black holes (SMBHs), feedback from the supernova and the active galactic nucleus (AGN), and the magnetic fields.

Baes et al. (2024a) generated high-resolution multi-wavelength images of 1160 galaxies from the z = 0 snapshot of the TNG50 simulation and released the dataset publicly as the TNG50-SKIRT Atlas database. Initially, the synthetic images were created in 18 bands (GALEX FUV and NUV, Johnson UBVRI, LSST ugrizy, 2MASS JHK_s, WISE W1 and W2), then Euclid Collaboration: Kovačić et al. (2025) extended the database with additional synthetic images in the four Euclid broad-band filters (I_E, Y_E, J_E, H_E). These images come with stellar mass surface density maps, stellar-mass-weighted metallicity and age maps, and dust mass surface density maps.

The simulated galaxies have a wide range of stellar mass (M_∗) from 10^9.8 to 10¹² M_⊙. To generate synthetic images, the 3D radiative transfer SKIRT code simulates the propagation of stellar light through the dust distribution in galaxies, which directly accounts for the absorption and scattering of light by dust. A detailed description of the radiative transfer post-processing is given in Baes et al. (2024a). In the following, we only briefly overview the synthetic imaging data.

In the radiative transfer analysis, the primary sources of radiation are the stellar particles in the simulated galaxy, which are obtained from the TNG50 database. The stellar particles are assigned the SED of a simple stellar population (SSP) based on their masses, ages, and metallicities. Stellar particles older than 10 Myr are assigned an SED generated from the Bruzual & Charlot (2003, hereafter BC03) stellar population synthesis (SPS) model. The stellar particles younger than 10 Myr are assigned an SED from the MAPPINGS III templates for star-forming regions (Groves et al. 2008). The dust density is determined by assuming a constant dust-to-metal fraction, f_dust = 0.2, adopted from Trčka et al. (2022), which conducted a thorough calibration of the free parameters in radiative transfer simulation using a sample of nearby galaxies from the DustPedia project (Davies et al. 2017). This value is consistent with observational estimates based on nearby galaxy samples (e.g. De Vis et al. 2019; Galliano et al. 2021; Zabel et al. 2021). The diffuse ISM THEMIS model (Jones et al. 2017) is used as the dust grain model at every location in the galaxy. The radiative transfer analysis does not include dust emission modelling because it only focuses on producing synthetic UV to NIR images. Each galaxy is projected at five different viewing angles and 22-band images are generated for each projection with a field of view of 160 × 160 kpc, centred on the galaxy. The synthetic images in all filters have a physical pixel size of 100 pc¹ and pixel value in units of megajanskys per steradian. In this paper, we only use synthetic images in the following filters: FUV, NUV (GALEX), u, g, r, i, z (LSST), I_E, Y_E, J_E, and H_E (Euclid).

2.2. Adding observational effects to the synthetic images

The synthetic images from the TNG50-SKIRT Atlas do not include observational effects, only small random fluctuations generated in the Monte Carlo process by the SKIRT code. To create realistic images that mimic the actual imaging data from the Euclid, LSST, and GALEX surveys, we added some observational effects to the synthetic images. First, we assigned a redshift of z = 0.03 and resampled (i.e. regridded) the images to have the same angular pixel size as that of the cameras: 1$\stackrel{″}{.}$5 for GALEX, 0$\stackrel{″}{.}$1 for the Euclid VIS band I_E, 0$\stackrel{″}{.}$3 for Euclid NISP bands, and 0$\stackrel{″}{.}$2 for LSST. This requires one to calculate the cosmological angular-diameter distance to convert the angular scale on the sky to an actual physical scale and adjust the resampling factor for the images. We used the 2D cubic spline interpolation method to resample the images. We assigned the same redshift to all simulated galaxies in the sample that we use for our testing (to be described in Sect. 3.1).

We adopted a uniform redshift of z = 0.03 in this test for several reasons. First, the TNG50-SKIRT ATLAS images are currently only available for galaxies in snapshot 99 (z = 0) of the TNG50 simulation. As a result, we cannot reliably redshift these images to significantly higher redshifts, where cosmological effects become non-negligible. Moreover, the primary goal of our future work with Euclid data is to spatially resolve stellar populations on sub-kiloparsec scales in local galaxies, limiting our target sample to z ≲ 0.1. Given this low-redshift focus, we consider our current testing set-up appropriate. Redshifting between z = 0 and z = 0.1 does not substantially affect spectral coverage, so we expect the trends in fitting results and the recovery of stellar population properties to remain valid.

The next step is simulating the observational noise and injecting it into the images. The noise includes both the photon shot noise and the sky background noise. For the Euclid and LSST images, we follow the same procedure as Euclid Collaboration: Merlin et al. (2023) and Euclid Collaboration: Martinet et al. (2019). Overall, the noise simulation is designed to achieve a given signal-to-noise ratio (S/N) at a given limiting magnitude by adjusting the instrumental zero points (ZPs). The ZP at a 1 second exposure was estimated using the following equation adapted from Euclid Collaboration: Merlin et al. (2023):

$$\begin{array}{c}\hfill \text{ZP}=2.5{log}_{10}\left[\frac{{(\mathrm{S}/\mathrm{N})}^2\pi r^{2}g}{t_{\exp }}\right]+2m_{\lim }-{\mathrm{SB}}_{\mathrm{bkg}}.\end{array}$$(1)

We adopt the observational parameters listed by Euclid Collaboration: Merlin et al. (2023, see Table 1 therein), which were made to mimic the imaging quality of the Euclid Wide Survey (Euclid Collaboration: Scaramella et al. 2022) and the LSST survey in its early data release. The parameters include the 10σ (i.e. S/N = 10) limiting magnitudes (m_lim) that are expected within an aperture of radius r = 1″, the background surface brightness (SB_bkg), and the exposure time (t_exp). These parameters are summarized in Table 1. The detector gains (in units of electrons/ADU) are g = 3.1 for I_E (Euclid Collaboration: Martinet et al. 2019), g = 2.0 for Euclid NISP bands (Secroun et al. 2018), and g = 1.06 for LSST (O’Connor 2019). The estimated ZP is then used for converting pixel values in the original units of megajanskys per steradian into the observational units; in other words, analog-to-digital units (ADUs). This is done before simulating the photon shot noise, where the pixel values in the noiseless images are perturbed by random noise generated from a Poisson distribution. After perturbing the pixel values, we then convert the images back to megajanskys per steradian.

The sky background noise was simulated as a Gaussian noise with zero mean and a standard deviation that is defined by the depth of the final simulated images. The zero mean reflects that the simulated images are background-subtracted. We adopted the standard deviation values for the background noise from Euclid Collaboration: Merlin et al. (2023) and private communication, and created blank sky maps containing only Gaussian noise. We added these maps to the synthetic images to get the final science images.

In addition to the science images, we needed to generate variance images. For this, we first calculated RMS maps. The noise per pixel is a contribution from the background noise and the photon shot noise from the source. The background noise is reflected in the standard deviation of the Gaussian noise in the empty background maps (RMS_sky), while the photon shot noise is Poissonian in nature. We used the following equation to calculate flux uncertainty of a 1-second exposure of individual pixels,

$$\begin{array}{c}\hfill N=\sqrt{{\mathrm{RMS}}_{\mathrm{sky}}^2+\frac{C}{t_{\exp }}},\end{array}$$(2)

where C and t_exp are the source counts and the exposure time, respectively. The variance images were then created by taking a square of the RMS maps. We have confirmed that our Euclid and LSST images reach similar depths at 10σ as those of the images produced by Euclid Collaboration: Merlin et al. (2023), which are summarized in Table 1. We checked this by looking at the correlations between S/N and magnitudes of 2″-diameter apertures that are randomly drawn across an image stamp, similar to Fig. 5 in Euclid Collaboration: Merlin et al. (2023).

We followed a similar procedure to the one described above to simulate the observational noise of the GALEX images with an exception in the calculation of the background RMS, which we derive from the real observational images. We downloaded four tiles of background-subtracted images for each band from the GALEX website². We specifically choose images from the medium imaging survey (MIS) of GALEX, which has an exposure time of ∼1500 seconds (Bianchi et al. 2014) because we intend to combine Euclid and LSST images with the GALEX images from this type of survey whenever available in our future analyses. To obtain the background RMS of an image tile, first, we remove contamination from the sources within the field by cropping a circular region within the Kron radius around each detected source. The co-ordinates and Kron radii of the sources are taken from the source catalogue associated with the field that we also downloaded from the GALEX website and was generated using SEXTRACTOR (Bertin & Arnouts 1996). We then take the standard deviation of pixel values in each blank sky image. To get the final background RMS, we take an average of the standard deviations from the four fields.

After adding the simulated noise to the science images and producing the variance images, we convolve all the images with the point spread functions (PSFs) of the corresponding cameras. For Euclid, we use model PSF images from the Euclid Mission Database models, while for LSST we use custom simulated PSFs from Euclid Collaboration: Merlin et al. (2023) that were created using PhoSim (Peterson et al. 2015). For GALEX, we use empirical PSF images of GALEX generated by Abdurro’uf et al. (2021) that are publicly available³. The PSF full width at half maximum (FWHM) and pixel size of the images are summarized in Table 1.

In this analysis, the convolution with the PSF is performed after adding simulated noise. We have checked that reversing the order of operations–applying PSF convolution before noise injection–yields images with similar characteristics and a minimal overall flux excess (approximately 1%), which does not significantly alter the SED at the pixel scale. We also emphasize that the synthetic images produced here have been cropped around each simulated galaxy, not in the form of wide-field imaging as obtained from observations.

3. Simulating and testing the spatially resolved SED fitting pipeline

Within the Euclid Collaboration, we plan to map the spatially resolved stellar population properties of a large number (100 000 s) of local galaxies at z below 0.1 using upcoming imaging data from Euclid combined with imaging data from other large surveys. For nearby galaxies (within 100 Mpc distance), we plan to include deep imaging data from GALEX surveys to get coverage in the UV. For achieving this goal, we will use piXedfit package (Abdurro’uf et al. 2021; Abdurro’uf et al. 2022d). In other efforts, we will apply SED fitting with stepwise SFH models (Euclid Collaboration: Nersesian et al., in prep.) and a machine learning method (Euclid Collaboration: Kovačić et al. 2025). One of the main objectives of this paper is to test to what extent Euclid-only images can be used to robustly map main stellar population properties using the spatially resolved SED fitting, and what improvements can be achieved by adding imaging data in the optical and UV. Synthetic images introduced in Sect. 2 are used for this analysis. In this section, we will first describe the criteria for selecting our sample of simulated galaxies in Sect. 3.1, and then describe our analysis pipeline inSect. 3.2.

3.1. Sample of simulated galaxies

We select a subset of galaxies from the parent sample of the TNG50-SKIRT Atlas. The sample is selected to cover a wide range of stellar mass (10^9.8–10¹² M_⊙) and specific star-formation rate (sSFR; 10⁻¹³–10⁻⁹ yr⁻¹). To construct the sample, first, we divide the entire mass range into five bins: 10^9.8–10^10.0, 10^10.0–10^10.5, 10^10.5–10^11.0, and 10^11.0–10^11.5, 10^11.5–10^12.0 M_⊙. After that, we further divide the galaxies in each mass bin into five quantile groups with increasing sSFR and then randomly choose one galaxy in each quantile. This results in a sample of 25 galaxies for our analysis in this paper. The distribution of our sample galaxies and the parent sample on the M_∗-sSFR plane are shown in Fig. 1. In addition to this sample selection, for each galaxy, we choose the most face-on configuration among the five viewing angles to minimize the systematic effects of the increasing dust attenuation in an edge-on view. This type of selection has been adopted in multiple spatially resolved studies to minimize projection effects that can otherwise bias the reconstruction of radial profiles from 2D maps of physical properties (e.g. Abdurro’uf & Akiyama 2017; Abdurro’uf & Akiyama 2018; Abdurro’uf et al. 2022b). Given the large number of galaxies that Euclid will observe, limiting the sample to relatively face-on galaxies will not significantly impact the conclusions of structural studies of galaxies. This is because, statistically, there is no expected intrinsic difference in the physical properties of galaxies viewed face-on versus edge-on.

Fig. 1. The Sample of 25 simulated galaxies used for the analyses in this paper. The selected sample and the parent sample from the TNG50-SKIRT Atlas are shown with red squares and blue dots, respectively. The sample is selected to cover the entire mass range in the TNG50-SKIRT catalogue and a wide range of sSFR.

The selected sample galaxies encompass various galaxy types, ranging from star-forming to passive galaxies. However, starburst galaxies are not included, and the imposition of a strict mass cutoff has resulted in the absence of dwarf galaxies.

The synthetic images of two example galaxies shown in Fig. 2 clearly demonstrate the variation in spatial resolution and depth. The GALEX images are the shallowest and have the lowest spatial resolution, while the I_E image is the deepest and sharpest. TNG501725 is a spiral galaxy with prominent spiral arm structures, while TNG414917 is an elliptical galaxy. TNG501725 appears to be clumpy in the GALEX images as these images are dominated by the light from young stars. In the longer-wavelength images, the galaxy appears to be less clumpy and the bulge becomes more prominent. TNG414917 is barely detected in GALEX images but appears to be bright in LSST and Euclid images, as expected for an elliptical galaxy, which has a lack of young stars.

Fig. 2. Examples of simulated images of GALEX, LSST, and Euclid to which observational effects have been added, including spatial resampling (i.e. regridding), simulated noise, and convolution with the PSF of the corresponding cameras. The galaxies here are TNG501725 with O1 orientation index and TNG414917 with O4 orientation index. The original synthetic images are taken from the TNG50-SKIRT Atlas database and are generated from the TNG50 simulation using the SKIRT radiative transfer code. The variation in spatial resolution and depth can be seen in the images with GALEX images being the most blurred and shallowest, and IE images being the sharpest and deepest. The pixel size difference is demonstrated with the red horizontal line (60″ long). The PSF FWHM and pixel size of the images are summarized in Table 1. The colour images are made from the composite of LSST g, r, and i images.

3.2. Analysis pipeline with piXedfit

We used piXedfit (Abdurro’uf et al. 2021; Abdurro’uf et al. 2022d), a state-of-the-art code that provides a comprehensive set of functionalities for analysing spatially resolved SEDs of galaxies. Given a set of multi-wavelength imaging data, piXedfit performs the three following steps:

• image processing, consisting of PSF matching, spatial resampling, reprojection to match the spatial resolution (i.e. PSF size), sampling (i.e. pixel size), and projection of the images based on the world co-ordinate system (WCS) information in their headers

• pixel binning to increase the S/N of the spatially resolved SEDs

• fit those SEDs with models to get the estimates of physical properties on spatially resolved scales.

Detailed descriptions of this code and robustness tests are given in Abdurro’uf et al. (2021). In the following, we test piXedfit using realistic synthetic images of Euclid, LSST, and GALEX, which is the first such test using zoom-in hydrodynamical simulation performed for piXedfit. This work intends to test how good the code is in recovering the maps of stellar population properties of simulated galaxies, given that we know the ground truth.

To investigate the effects of wavelength coverage and variation in the imaging datasets, we performed analyses on three different sets of imaging data cubes: (1) a combination of GALEX, LSST, and Euclid (11 bands) that cover UV to NIR (hereafter called G + L + E data cube); (2) LSST and Euclid (9 bands) that cover optical to NIR (hereafter called L + E); and (3) images that mostly cover NIR (hereafter called E). In the following, we describe briefly the three main steps of our analysis introduced before. We refer the reader to Abdurro’uf et al. (2021) for a more detailed description of each step.

3.2.1. Image processing and generation of photometric data cubes

In the image processing with piXedfit, input multi-wavelength images are matched in spatial resolution (i.e. PSF size), spatial sampling (i.e. pixel size), and projection (i.e. a common WCS). PSF matching is performed to homogenize the spatial resolution of the multi-wavelength imaging data, which is performed by convolving higher-resolution images with pre-calculated kernels that can degrade their spatial resolutions to match those of the lowest-resolution image. We generate convolution kernels based on the PSF images in each filter using the PHOTUTILS package (Bradley et al. 2022). We use the same PSF images as the ones we used for convolving the synthetic images (see Sect. 2.2). After PSF matching, the images are resampled (i.e. regridded) into the same spatial sampling, which is chosen as the largest pixel size among the imaging data being involved. The final PSF size and the pixel size of the post-processed images produced from the image processing depend on the set of imaging data used. The final PSF FWHMs and pixel sizes of the three types of data cubes are summarized in Table 2, showing that Euclid images alone result, obviously, in the highest spatial resolution.

Since all the synthetic images from the TNG50-SKIRT Atlas have the same projection, and we do not assign co-ordinates and WCS information to them, we only resample the images. In an analysis with real imaging data, piXedfit can perform the reprojection and resampling processes simultaneously, which retain the WCS information of the images. The images are then cropped around the target galaxy to produce final stamp images. After that, the region of interest (RoI) of the galaxy, within which the spatially resolved SED fitting will be performed, is defined. First, segmentation maps are obtained using SEP (Barbary 2016), which is called within piXedfit, with the detection threshold (thresh) of 2.0 relative to the background RMS. The segmentation maps in all filters are then merged (i.e. combined) together to get the final RoI of the galaxy. After that, the flux densities of pixels within the RoI are calculated in units of $\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}{\mathrm{\AA }}^{-1}$. This process produces 3D data cubes that are stored in FITS files. We perform the image processing to all the sample galaxies with the three sets of imaging data (i.e. G + L + E, L + E, and Euclid only), resulting in 75 data cubes, 3 data cubes for each sample galaxy. Examples of the maps of flux densities of the galaxy TNG501725 from the 3 types of data cubes are shown in Fig. 3.

Fig. 3. Examples of the maps of flux densities produced from the image processing on the sets of three data cubes: G + L + E (a), L + E (b), and E (c). This example of the galaxy is TNG501725 with the O1 orientation index. The final PSF size and pixel size of the 3 data cubes are summarized in Table 2. The G + L + E data cube has the lowest spatial resolution, trading the high spatial resolution of Euclid and LSST with the wider wavelength coverage from UV to NIR.

Since the galaxy’s RoI is defined by combining segmentation maps in multiple bands, the deepest images control the final region because they usually have a larger segmentation area. Figure 4 shows surface brightness radial profiles of TNG501725 derived from the stamp images of G + L + E data. The deep LSST and Euclid images extend the galaxy’s region up to a radius of around 60 kpc and reach a surface brightness of up to 29 mag arcsec⁻² before sky background noise dominate, while shallower GALEX images only reach up to around 40 kpc with a surface brightness of up to 28 mag arcsec⁻². The galaxy’s RoI extends up to a radius of ∼45 kpc for this particular galaxy, as indicated by the vertical dashed grey line. As we can see from this figure, the galaxy’s RoI is mainly defined by the LSST and Euclid images, and the flux maps of GALEX are extended such that their outskirt regions are noisy, which can also be seen in Fig. 3. In the GALEX map, the missing pixels at the galaxy’s outer regions are caused by negative flux values. Due to the brighter overall surface brightness in the longer wavelengths in this local galaxy, the lowest surface brightness reached in the outskirts of the galaxy’s RoI is deeper in the shorter wavelengths. For TNG501725, the minimum surface brightness of GALEX (FUV and NUV), LSST (u, g, r, i, z), and Euclid (I_E, Y_E, J_E, H_E) are 28.7, 28.9, 29.1, 27.6, 26.9, 26.7, 26.4, 26.9, 26.3, 26.2, and 26.1 mag arcsec⁻², respectively. For all galaxies in the sample, the median values are 28.7, 28.7, 29.1, 27.8, 27.2, 27.0, 26.8, 27.3, 26.7, 26.6, and 26.6 mag arcsec⁻². We note that these values depend on the assumed exposure times in our noise simulations. Moreover, the RoI can also be enlarged by decreasing thresh to include the fainter outskirt region of the galaxies.

Fig. 4. Surface brightness radial profiles of TNG501725 derived from the stamp images of G + L + E data. The vertical dashed grey line indicates the radial extent of the data cube after the galaxy’s RoI is defined, roughly corresponding to a detection threshold of 2 above the background RMS of Euclid images.

Given the nature of spatially resolved analyses that involve multi-band imaging data from various telescopes with differing spatial resolutions, there is a concern that PSF matching may alter a galaxy’s surface brightness. To assess this, we examined the effect of PSF convolution on surface brightness radial profiles and found that the impact is minimal, even when convolving to the broadest resolution of the NUV band.

We also investigated the effect on peak surface brightness (i.e. the central value) and found no significant overall change due to PSF convolution. However, a slight systematic effect emerges when convolving to the NUV resolution, where the peak surface brightness in the convolved images is decreased by ∼0.2 dex compared to the original images. This effect is not present in the FUV and NUV images themselves but is observed in the other bands.

3.2.2. Pixel binning

The fluxes of pixels usually have low S/N and cannot provide sufficient constraints on the model SEDs. piXedfit provides a pixel binning feature that can maximize the S/N of SEDs on spatially resolved scales by binning neighbouring pixels that have similar SED shapes. This scheme of piXedfit can retain important information on the pixel level while binning pixels together, which degrades the spatial resolution. This is because it takes into account the similarity of SED shape among the pixels such that only pixels that have similar SED shape are binned together. The pixel binning scheme takes four input parameters: (1) the S/N thresholds that can be set for multiple filters; (2) reduced χ² value for evaluating the similarity of SED shape; (3) the minimum diameter of spatial bins; and (4) reference filter for sorting pixels based on their brightness. We refer the reader to Abdurro’uf et al. (2021) for a more detailed description.

We experiment to get optimum binning parameters that can give a sufficient S/N without losing spatial resolution and detailed spatial structures of our sample galaxies. We converge to the following binning parameters:

• The minimum S/N of 5 in all filters other than FUV, NUV, and u;

• Maximum reduced χ² = 5 for the evaluation of the similarity in SED shape among pixels;

• The smallest bin’s diameter of 5 pixels;

• I_E is chosen as the reference filter for sorting the pixel brightness.

We do not set the S/N threshold for FUV, NUV, and u filters because pixels in those images have low S/N such that setting constraints on them can result in bigger overall bin sizes as compared to the binning map obtained with the current setting. This can cause the loss of detailed spatial information in the galaxies. We apply these binning parameters to all 75 photometric data cubes from our three analyses (see Sect. 3.2.1). The total number of spatial bins from these data cubes is 219 158.

Binning maps and S/N plotted as a function of pixels and spatial bins for TNG501725, obtained with the three data cubes, are shown in Fig. 5. As we can see from the plots, the pixel binning by piXedfit can achieve the S/N thresholds. The G + L + E data cube has the least number of spatial bins (343) mainly due to the larger pixel size of this as compared to the other two data cubes. The L + E and E data cubes have a total number of spatial bins of 8391, and 9044, respectively. The smaller number of spatial bins in the former, although they both have the same pixel size, is caused by the slightly lower average S/N per pixel of the LSST images.

Fig. 5. Examples of the pixel binning maps of TNG501725 obtained with the three data cubes: G + L + E (a), L + E (b), and Euclid only (c). In each row, the three panels from left to right show the binning map, S/N radial profile of pixels, and S/N radial profile of spatial bins. The pixel binning scheme of piXedfit can achieve the S/N threshold of 5 (indicated by the dashed black lines) for the g, r, i, z, IE, YE, JE, and HE filters.

3.2.3. SED fitting set-up

After performing pixel binning to the data cubes, we proceed with SED fitting. For each galaxy, piXedfit fits the SED of spatial bins in the galaxy and obtains the estimates of the physical properties of the underlying stellar populations. Since our synthetic imaging data only covers UV to NIR, we only include the stellar emission modelling in our SED fitting. We turn off the modelling of dust emission because we do not have data in MIR and FIR that can constrain the models. We also turn off the nebular emission modelling in this test analysis because the broad-band photometry is expected to be weakly affected by the nebular emission, especially for our sample of local galaxies that are forming stars only moderately. We will implement this modelling in future analyses with real data. A detailed description of the SED modelling and fitting method is given in Abdurro’uf et al. (2021). In the following, we briefly outline the priors that we assume in our SED modelling andfitting.

We assume the Chabrier (2003) IMF, the Calzetti et al. (2000) dust attenuation law, and a parametric star-formation history (SFH) model in the form of a double power law. It has been shown in Abdurro’uf et al. (2021) that this SFH model has sufficient flexibility to be able to reconstruct the shape of the true SFH when tested with mock SEDs of simulated galaxies from the TNG100 simulation. The SED modelling has six free parameters and three dependent parameters. We summarize these parameters along with the assumed priors in Table 3. The SFR is defined as the recent instantaneous rate determined from the SFH at the time of observation. We have verified that this SFR is consistent with the value obtained by averaging over a specific timescale (e.g. 100 Myr), which is a common approach in the literature.

For the fitting method, we use the Bayesian inference method with the random dense sampling of parameter space (RDSPS) as the posterior sampling method. This method samples the parameter space by generating random values for each parameter that are uniformly distributed within the assumed range. We choose this method because it is much faster than the MCMC method⁴. We refer the reader to Abdurro’uf et al. (2021, Sect. 4.2.2 therein) for more detailed descriptions of this method. In this fitting method, M_∗ is not a free parameter. For each model, M_∗ is determined by the best model normalization derived by minimizing the χ² (see e.g. Abdurro’uf & Akiyama 2017, their Eq. (6)). The mass-weighted age t_M, which is a representative age of a stellar population, is directly calculated from the model’s SFH using the following equation:

$$\begin{array}{c}\hfill t_{\mathrm{M}}={\displaystyle \frac{{\int }_{t_{\mathrm{form}}}^{t_z}t\mathrm{SFR}(t)\mathrm{d}t}{{\int }_{t_{\mathrm{form}}}^{t_z}\mathrm{SFR}(t)\mathrm{d}t}},\end{array}$$(3)

where t_z is the age of the Universe at the galaxy’s redshift and t_form is the cosmic time when star formation is started in the galaxy. We refer to this quantity interchangeably as ‘age’ or ‘mass-weighted age’ in the paper. It is worth mentioning that the assumed flat priors induce a ‘hidden’ non-flat prior to the mass-weighted age, peaking around 4–5 Gyr. This age distribution closely resembles the observed distribution of spatially resolved stellar population age distribution in local galaxies by MaNGA (Abdurro’uf et al. 2022a; Sánchez et al. 2022; Lacerda et al. 2022). It is also a distribution found in the priors of a variety of previous studies, as shown in the comparative study of Leja et al. (2019).

For each data cube, we performed multiple SED fitting runs by varying two factors: the adopted SPS models and the priors. We consider two SPS models: the flexible stellar population synthesis (FSPS; Conroy et al. 2009; Conroy & Gunn 2010) and the 2016 version of Bruzual & Charlot (2003), hereafter BC16. The former is the native SPS model used in piXedfit. With this, we can learn possible systematic effects caused by the SPS model on the fitting results. In the FSPS model, we use the Padova isochrones (Girardi et al. 2000; Marigo & Girardi 2007; Marigo et al. 2008) and stellar spectral library from the Medium-resolution Isaac Newton Telescope library of empirical spectra (MILES; Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011; Röck et al. 2016)⁵.

For generating model SEDs based on the BC16 model that assumes the double power-law SFH, we used BAGPIPES (Carnall et al. 2018). This code applies the same isochrones and stellar spectral library as those used in FSPS, but a different IMF from Kroupa & Boily (2002). However, the expected difference in the estimated M_∗ between Kroupa & Boily (2002) and Chabrier (2003) IMFs is very small (within ∼0.03 dex, see e.g. Speagle et al. 2014). In both the SPS models, each SED is characterized by a constant metallicity.

For the variation in priors, we defined two settings: the first corresponds to the classical flat priors described in Table 3, and the second adds a special prior in the form of mass-Z and mass-age relations (hereafter called mass-Z-age prior). We have constructed the latter prior from the combination of empirical relations taken from the MaNGA survey (Bundy et al. 2015) and theoretical predictions from the TNG50 simulation. We describe this in more detail in Appendix A. The extra prior is particularly relevant when the wavelength coverage is limited and physical information available in the data is scarce, as will be further discussed later.

4. Results

4.1. Examples of the maps of stellar population properties

After SED fitting is performed on all the spatial bins in a galaxy, 2D maps of stellar population properties can be constructed. The properties include stellar mass surface density (Σ_∗), SFR surface density (Σ_SFR), mass-weighted age, and metallicity (Z). Examples of the maps of stellar population properties of TNG501725 obtained from the analyses with the three data cubes are shown in Fig. 6. The analyses with L + E and E data cubes yield maps that can resolve the spiral arms and bulge structures in great detail, in contrast to the G + L + E due to its worse resolution. These structures appear in all the maps of properties, not only in the Σ_∗ and Σ_SFR. From the maps, especially the ones obtained with L + E data cube, we can see that overall, bulge and spiral arms have higher Σ_∗, Σ_SFR, metallicity, and younger stellar populations than the interarm regions. The maps of the rest of our sample galaxies are presented in Appendix C.

Fig. 6. Examples of the maps of the stellar population properties of TNG501725 obtained from the spatially resolved SED fitting on the three data cubes: G + L + E (a), L + E (b), and E (c). The maps derived from the analysis with the L + E and E data cubes alone have a high spatial resolution and reveal detailed spatial variation in the physical properties compared to those obtained from the analysis with the G + L + E data cubes.

4.2. Recovering spatially resolved stellar population properties

In this section, we analyse the robustness of our analysis pipeline in recovering the true stellar population properties on spatially resolved scales. We investigate the effects of three important factors on the fitting results: wavelength coverage of the data cube, SPS model, and assumed priors. For the prior variations, we compare between flat (i.e. uniform) and the mass-Z-age priors. First, in this section, we directly compare the inferred parameters of all spatial bins in the sample galaxies derived from SED fitting and the true values from the TNG50 simulations. We focus on four main properties in this paper: Σ_∗, Σ_SFR, age, and Z.

In Figs. 7, 8, and 9, we show the comparisons between our fitting results and the ground truth from the analyses with the three imaging data cubes. The contours show the distributions of all spatial bins in the sample galaxies. We show the results obtained with the FSPS and BC16 models, with and without the mass-Z-age prior. To make quantitative comparisons, we show the histograms of the logarithmic ratio between the best-fit parameters and the ground truth in the figures and calculate the mean offset (μ), scatter (i.e. standard deviation; σ), and the Spearman rank-order correlation coefficient (ρ). ρ is a nonparametric measure of the monotonicity of a relationship between two datasets and can indicate the significance of a correlation. The derived parameters are summarized in Table 4 for the fitting that applies the mass-Z-age prior and Table 5 for the one that does not apply the special prior.

Fig. 7. Comparisons of the spatially resolved stellar population properties derived from SED fitting on the G + L + E data cubes (ordinate) and the true properties (abscissa). These include all spatial bins from the whole sample of galaxies. We focus on four main properties in this test analysis: Σ∗, Z, age, and ΣSFR. The top row (a) shows results obtained with the FSPS model, while the bottom row (b) shows results using the BC16 model. The blue and red colours represent results obtained with and without applying the mass-Z-age prior.

Fig. 8. Similar to Fig. 7 but for results of SED fitting on the L + E data cubes.

Fig. 9. Similar to Fig. 7 but for results of SED fitting on the E data cube.

Not too surprisingly, from Figs. 7, 8, and 9, we can see that Σ_∗ is recovered very well in all the analyses with the three types of data cubes, corroborated by small offset (μ ≲ 0.06 dex), small scatter (σ ≲ 0.12 dex), and high Spearman rank-order correlation coefficient (ρ ≳ 0.97). This is the case for both SPS models and whether or not the mass-Z-age prior is applied. This result suggests that stellar mass can be well recovered on spatially resolved scales using photometry data that cover rest-frame NIR. This is encouraging and shows the great potential of Euclid for mapping Σ_∗ of local galaxies.

Star formation rate density can be recovered reasonably well with the G + L + E data cube but poorly with the L + E and E data cubes. Notably, most data points cluster along the one-to-one line for L + E data cube. A critical threshold emerges at an Σ_SFR of approximately 10⁻³ M_⊙ yr⁻¹ kpc⁻², where SFR recovery becomes markedly differentiated–performing well above this density and poorly below it, particularly for L + E and E data cubes. This differentiation is evidenced by the mean (μ) and standard deviation (σ) of the logarithmic ratio between true and inferred Σ_SFR as shown in Tables 4 and 5, with values for spatial bins above 10⁻³ M_⊙ yr⁻¹ kpc⁻² noted after slash symbol.

These results show that our pipeline can recover spatially resolved SFR well with G + L + E data irrespective of how active star formation in the spatial region is, while with L + E and E data cube, it can recover SFR well only in star-forming regions and fails in passive regions. The critical SFR density of 10⁻³ M_⊙ yr⁻¹ kpc⁻² is relatively low, even by local galaxy standards, a threshold typically exceeded in spiral arm regions. The less accurate Σ_SFR estimates derived from L + E and E data cubes are expected, primarily due to their lack of UV coverage. UV data are crucial as indicators of SFR and dust attenuation. Although L + E data offers a marginal improvement by capturing the 4000 Å break, a valuable age indicator, Euclid-only data provides the most limited constraints on SFR, dust properties, and stellar age.

Without the mass-Z-age prior, the inferred metallicities and ages are systematically underestimated even in the G + L + E data cube. The biases may have their origin in the traditional flat priors (combined with parameter-estimates based on marginalized posterior distributions rather than best-fits), since a flat distribution in log₁₀(Z) contains more low-Z than high-Z models, and a flat distribution of ${\widehat{\tau }}_2$ contains more low values of $\exp (-{\widehat{\tau }}_2)$ than high values. It is also worth recalling the impact of the non-flat distribution of the mass-weighted ages.

Once the prior is applied, age is recovered fairly well with the G + L + E and L + E data cubes, but barely recovered with the E data cube. The trends obtained with the two SPS models are similar regardless of the data cubes being used, which suggests a weak effect of the SPS model variation. In the analysis with the G + L + E data, the results obtained by applying the special prior have a smaller systematic offset (μ ∼ −0.07) as compared to those obtained without applying the priors (μ ∼ −0.13). Interestingly, the fitting results without applying the priors do not deviate significantly from the one-to-one line, as is the case for Z. The systematic offset becomes slightly larger in the analysis with the L + E data. The deviation from the true value is more significant for the older stellar populations than the younger ones. We also note the similarity between the results obtained for age with and without priors in the analysis with the L + E data. This might indicate that the special priors cannot improve the estimate of age further for this dataset that only covers the rest-frame optical-to-NIR, and missing the UV data, which is an important constraint for age. As we can see from Fig. 9, the E data cube struggles to recover age, as shown by the broadly flat distribution of the fitting results across a wide range of the true values, flatter than the results with L + E data (Fig. 8). This is understandable because of the absence of UV and optical data. Although it has I_E filter, its wavelength range does not cover the Balmer break, which is an important age indicator.

For metallicity, the best ground truth recovery is achieved using the G + L + E data cube and applying the mass-Z-age prior to the fitting process. The variation in the SPS model does not seem to affect the results, unlike the priors that play a big role. With the help of the mass-Z-age prior, the result is encouraging, with a small absolute offset of ≲0.05 dex, a scatter of around 0.08 dex, and ρ greater than 0.6. Removing GALEX and LSST data, but still applying the prior, results in a slightly larger offset than that achieved with the G + L + E data cube. The L + E gives μ ∼ −0.12, σ ∼ 0.11, and ρ ∼ 0.6 with the two SPS models, while the E data cube gives slightly larger offset and scatter but similar ρ. These results show the significant influence of priors in the SED fitting with the Bayesianmethod.

In summary, it is important to emphasize the joint role of priors and data. For example, the comparison of offsets in Tables 4 and 5 demonstrates that the priors are essential for accurately recovering age and metallicity. However, incorporating LSST data and, subsequently, GALEX data also plays a significant role. Specifically, for the BC16 case, applying the special priors improves the age offsets from −0.13 to −0.11, −0.13 to −0.11, and −0.13 to −0.06 for the E, L + E, and G + L + E data cubes, respectively. The most significant improvement is naturally achieved when all datasets are used. When only E data are used, the inferred age and metallicity are mainly driven by the priors.

To get a better look at the influences of priors on the fitting results, we plot the distributions of derived age and metallicities on the mass-Z-age prior in Fig. A.2. Overall, we can conclude that the influence of priors becomes more prominent when fewer data points (covering a narrower wavelength range) are used, as expected.

When using only the E data with flat priors, the estimated ages exhibit a very narrow distribution centred around ${log}_{10}({\mathrm{Age}}_{\mathrm{true}}[\mathrm{Gyr}])\sim 0.6-0.7$. This corresponds to the peak of the mass-weighted age distribution of the SPS model libraries used in the SED fitting procedure. This result confirms that with limited data and without imposing the mass-Z-age prior, the ‘hidden’ mass-weighted age prior can bias the fitting results.

4.3. Offset maps: Locating the good and poor fitting results

Next, we analyse the maps of spatially resolved offsets between the fitting results and the true values to investigate the spatial regions where good and poor fitting results are located and try to get more insights into the contributing factors to the discrepancies. We show offset maps of a spiral galaxy TNG501725, obtained from the analyses with the G + L + E, L + E, and E data cubes. We only show results from SED fitting with the FSPS model in this example because the BC16 model produces offset maps that are similar overall. The mass-Z-age prior is applied.

The offset maps show that our SED fitting slightly overestimates Σ_∗ in the spiral arm regions of TNG501725, though the median offset is small. This is also the case for other spiral galaxies in our sample. In contrast to Σ_∗, the offset map of age is roughly flat across the galaxy’s region in all data types. Both good estimates of age (i.e. nearly zero offset) and overestimation can be seen in the spiral arms. For metallicity, the G + L + E data give good estimates in almost the whole region, except the outer shell where Z is overestimated. However, the offset maps are roughly flat for the results with L + E and E data cubes, showing an overall underestimation of Z. In contrast to Σ_∗, SFR is underestimated in the spiral arm regions and overestimated in the interarm regions. The colour map shows that the absolute offset is larger for SFR than other properties, exceeding 0.4 dex (the range set for consistency with other maps) for results using L + E and E data cubes.

4.4. Recovering radial profiles of stellar population properties

In galaxies, dynamical timescales, gas properties relevant to star formation, the history of stellar feedback, and the influence of neighbouring galaxies, all depend on radius. The radial profiles of stellar population properties are often used for describing the local distribution of galaxy properties because their variations with distance from the galactic centre can inform us about the growth history of the galaxies. For instance, the negative gradient of the age radial profile and the positive gradient of specific SFR is thought to be the consequence of the growth history of the galaxy structure that happens in an inside-out manner (e.g. Abdurro’uf et al. 2022b; Abdurro’uf et al. 2023). Here we check the robustness of the radial profiles of Σ_*, age, and Z. For simplicity, we only analyse the radial profiles from the fitting results with the FSPS model and apply the mass-Z-age prior because this fitting set-up gives the best recovery of the ground truth as shown in Sect. 4.2.

The radial profiles of the systematic offset of all galaxies in our sample are shown in Fig. 11. As we can see from the figure, the absolute offset varies with radius in the majority of the galaxies in our sample. The offset is generally smaller in the central region than in the outskirts, which is reasonable given the overall higher S/N in the central region. The Σ_* offset is centred around zero and has a scatter that increases with radius. However, the overall absolute offset is smaller than 0.3 dex across the radial extent. No galaxy exhibits Σ_* values that systematically deviate from the true values over the entire galaxy’s radial extent (see, for example, Fig. 10). In contrast, some galaxies show consistent deviations from the ground truth in age and Z throughout their entire regions. In those galaxies, the offset radial profile tends to be flat (i.e. constant across the galaxy’s region). Despite the full deviation, the overall offset is mostly within ±0.3 dex.

Fig. 10. Examples of spatially resolved distributions of the systematic offsets (in logarithmic scale) between the parameters derived from SED fitting and the true values. Offset maps of a galaxy TNG501725, are shown here for example. The first, second, and third rows are for analyses with G + L + E, L + E, and E data cubes, respectively. Only results obtained with the FSPS model and applying the mass-Z-age prior are shown.

Fig. 11. The Systematic offsets between estimated and input properties (μ) as a function of the radial distance from the galactic centre. The radial profiles of the offsets, as obtained from fits with the FSPS models and with the mass-Z-age prior, are shown for all sample galaxies with three types of data cubes. Columns from left to right show results from the analysis with the G + L + E, L + E, and E data cubes. Each galaxy’s radial extent is divided into ten radial bins equally spaced and radial profiles are then derived by taking the median of μ in circular annuli. The radius is normalized to the half-mass radius that is taken from the catalogue of the TNG50-SKIRT Atlas. The two horizontal dashed grey lines mark ±0.3 dex offset.

4.5. Recovering the integrated stellar mass

Previous studies have found some indications that global M_∗ derived from fitting to the integrated SED of galaxies could be underestimated as compared to the integrated M_∗ derived from spatially resolved SED fitting (e.g. Zibetti et al. 2009; Sorba & Sawicki 2015, 2018). The underestimation of M_∗ is thought to be due to the domination of light from young stellar populations that make up the integrated SED of a galaxy. This outshining effect can bias the fitting of integrated SED towards inferring younger and less massive stellar population properties, especially if the assumed parametric SFH model is too simple and not flexible enough to capture the star-formation activity that occurred farther back in time. This effect becomes more significant at high redshifts, as young stars mostly dominate the stellar populations in galaxies. Spatially resolved SED fitting can alleviate this effect because it can directly disentangle young and old stellar populations in the process. Here, we compare integrated M_∗ derived from the two methods (i.e. spatially resolved and integrated SED fitting) and the true values from TNG50. The integrated SEDs are obtained by summing up fluxes of all pixels within the galaxy’s RoI, then SED fitting is performed on them. The set-up for SED modelling and all the assumed priors in the fitting process are the same as those used in the spatially resolved SED fitting (see Sect. 3.2.3 and Table 3). For the true M_∗, we adopt the values from the TNG50-SKIRT catalogue, which are the total mass of stellar particles in the subhalo of thegalaxies.

We show the comparisons between the two estimates of the global stellar masses and the ground truth in Fig. 12. For M_∗ derived from SED fitting, we show the six variations in the fitting, consisting of two SPS models and three sets of imaging data cubes. These variations are performed to both the integrated and spatially resolved SED fitting. From this figure, we can see that overall, integrated M_∗ is well recovered by all the methods, which are shown by the very small systematic offset and scatter in the one-to-one relations. We also compare the integrated M_∗ derived from the integrated and spatially resolved SED fitting. Considering the μ values, stellar masses from integrated SED fitting are only imperceptibly lower than those from resolved SED fitting. Within the uncertainties, there are no noticeable differences between the two estimates, which suggests that our integrated SED fitting is not affected by the outshining effects of young stellar populations in the galaxies, in contrast to the results of previous studies (e.g. Zibetti et al. 2009; Sorba & Sawicki 2015, 2018). Other previous studies have also found consistent stellar mass derived from spatially resolved SED fitting and integrated one (e.g. Wuyts et al. 2012; Bellstedt et al. 2024), consistent with our finding. This can be caused in part by two factors. First, our sample galaxies are forming stars only moderately, and their stellar populations are not dominated by young stars. Therefore, the outshining effect is expected to be minor. The outshining effect is expected to be more significant in galaxies that are actively forming stars, which are more prevalent in high-redshift samples (e.g. Sorba & Sawicki 2018). Second, our SFH model in the form of the double power-law function and the assumed priors for the associated parameters are flexible enough to capture the SFH of the galaxies, especially those that occurred farther back in time, which formed the old stars. Abdurro’uf et al. (2021) has shown that this SFH model is flexible enough to reconstruct the true SFH of TNG100 galaxies in a fitting test with mock SEDs.

Fig. 12. The Robustness of integrated M∗ derived from the spatially resolved SED fitting and the integrated SED fitting compared to the true M∗ from the TNG50 simulation. Left panel: Comparison between M∗ derived from spatially resolved SED fitting (i.e. by summing up M∗ of pixels within the ROI of the galaxies) and the true values. Middle panel: Comparison between M∗ derived from integrated SED fitting and the true values. Right panel: Comparison between integrated M∗ from spatially resolved and integrated SED fitting.

5. Discussion

In our tests, we have made variations in some of the main and basic factors involved in SED fitting, including the SPS model, the wavelength coverage of the photometric data, and priors. In the following, we further discuss the influences of these three factors on the fitting results. We also put our results in the context of findings from previous studies in the literature.

5.1. The influence of SPS model

The variation in the SPS model is made to accommodate the SPS model used in generating the synthetic images by the TNG50-SKIRT Atlas, which is BC16 (see Sect. 2.1), while piXedfit uses FSPS by default. We use the two models in our experiment for two reasons: (a) by applying the same BC16 model, the systematic effects caused by the difference in the SPS model can be eliminated; and (b) we get the opportunity to compare the results obtained using the two SPS models. As we can see from Figs. 7, 8, and 9 (Sect. 4.2), there is no significant difference between the results obtained with the FSPS and BC16 models, whether or not the mass-Z-age prior is applied. The influence of the SPS model is weak on the derived M_∗, Z, and age in our SED fitting tests. Therefore, using two of the most commonly adopted SPS models in the literature, we demonstrate that the uncertainty in our description of stellar populations in nature has a negligible effect on our SED fitting process.

5.2. The influence of wavelength coverage of the imaging data

Euclid covers the NIR and red optical wavelengths, but it is well known that data at shorter wavelengths, preferably covering the Balmer break and the UV, hold sensitive complementary information on stellar populations. This is particularly important when dealing with spatially resolved data where S/N is typically relatively low. The combination of wavelengths is helpful to break degeneracies among age, Z, and dust attenuation (e.g. Worthey 1994; Kodama & Arimoto 1997; Conroy 2013), although more passbands are not necessarily sufficient for this purpose (e.g. Pforr et al. 2012). Therefore, it is essential to test to what extent Euclid imaging data alone can provide constraints on the main stellar population properties and what improvements can be gained by adding more imaging data in the optical and UV.

From Figs. 7, 8, and 9, we can see that M_∗ measurement is very stable against the variation in the wavelength coverage of photometric data. M_∗ is recovered very well using any of the three data cubes as Euclid covers the NIR wavelength. On the other hand, age and Z are affected more significantly by the presence or absence of short wavelength data. Age is poorly constrained when using the E data cube with the classical ‘flat’ priors: the estimated ages cluster tightly around 4 − 5 Gyr, corresponding to the peak of the mass-weighted age distribution of the SPS models for these priors. However, even with the implementation of informative priors (i.e. mass-Z-age), ages remain more poorly constrained compared to metallicities. In the fitting that implements the informative priors, the systematic offset (i.e. underestimation) in Z measurement is enhanced with the exclusion of GALEX (μ ∼ 0.11 dex) and GALEX+LSST (μ ∼ 0.15 dex), while the strong proportionality between the inferred and true values is stable against the data variation.

5.3. The influence of the mass-Z-age prior

The prior is one of the important components in the Bayesian SED fitting, providing additional information and constraints in the fitting process. In the absence of highly constraining data (e.g. spectroscopy or multi-wavelength photometry over a wide wavelength range), the role of priors can be significant in getting reliable fitting results. There are multiple forms of priors in SED fitting, and they can be assumed explicitly (e.g. ranges in the grid of free parameters) or implicitly (e.g. dust attenuation law, and the form of analytic SFH). There have been only a few studies that adopted empirical scaling relations for a prior in the Bayesian SED fitting. Leja et al. (2019) implemented a prior in the form of mass-Z relation in their integrated SED fitting of high-z galaxies. They adopted the empirical mass-Z relation observed by Gallazzi et al. (2005). Similar scaling relations also exist among stellar population properties on kpc scales in the local galaxies (see review by Sánchez 2020). This provides an opportunity to use them for informative priors in spatially resolved SED fitting.

Applying the classical flat priors has resulted in very poor constraints on stellar population parameters, underestimating metallicity and age, and overestimating dust attenuation.

The implementation of the mass-Z-age in our work has showcased, instead, its potential. As we can see from Figs. 7, 8, and 9, the role of the informative priors is significant in recovering age and Z. The most significant effect can be seen in the measurement of Z with the G + L + E data cube, where fitting with priors can significantly improve the fitting results.

To further test the performance of our SED fitting pipeline, we perform additional fitting tests using simple mock SEDs generated using FSPS without preprocessing, unlike the synthetic images that are generated using the radiative transfer process. We present the results in Appendix B. Overall, we obtain trends similar to those seen in the tests using the synthetic images. M_∗ is recovered very well with the three types of data cubes and seems to be independent of priors, while Z and age are recovered relatively well with the three data cubes only if the mass-Z-age prior is applied.

6. Conclusions and outlooks

We tested a pipeline for spatially resolved SED fitting of local galaxies as part of Euclid science missions on the local Universe. We aimed to test to what extent Euclid-data only can be used to measure spatially resolved stellar populations of local galaxies, and what improvements can be obtained by including imaging data in the optical and UV from other telescopes. We used synthetic images representative of Euclid, GALEX, and LSST data based on the TNG50-SKIRT Atlas (Baes et al. 2024a,b), which were generated using a 3D Monte Carlo radiative transfer process.

Our pipeline uses functionalities provided by piXedfit (Abdurro’uf et al. 2021; Abdurro’uf et al. 2022d) for processing the imaging data (i.e. matching the spatial resolution and sampling of the various passbands), binning neighbouring pixels to achieve chosen S/N thresholds in multiple bands, and fitting the spatially resolved SEDs to infer the local properties of the stellar populations. In this work, we ran our pipeline on TNG50-SKIRT galaxies with modest levels of current star formation to test its robustness in recovering three main properties on spatially resolved scales: stellar mass surface density (Σ_∗), age, and metallicity (Z). We constructed three types of data cubes: GALEX + LSST + Euclid, LSST + Euclid, and Euclid-only (G + L + E, L + E, and E, respectively), and ran SED fitting on each of them with varying SPS models (FSPS and BC16) and priors. We implemented a new informative prior in the form of the mass-age and mass-Z relations (i.e. mass-Z-age prior), adopted from a combination of empirical and simulated data (see Appendix A).

The mass surface density can be recovered very well using the three data cubes, regardless of the variations in the SPS model and priors (i.e. whether or not the mass-Z-age prior is implemented). The systematic offset (μ) and scatter (σ) of the logarithmic ratio between the inferred and true values are minimal, ≲0.05 (absolute value) and ∼0.1 dex, respectively. This suggests that the Euclid-only data are sufficient to robustly map Σ_∗ of mildly star-forming local galaxies. A slight overestimation of Σ_∗ seems to occur in the regions that are more actively forming stars or contain younger stellar populations than the rest of the galaxy’s region, such as in the spiral arms.

Not surprisingly, the added prior can improve the Z and age retrieved from the simulated data. In the analysis with the G + L + E data cube, Z can be recovered reasonably well (μ ∼ 0.03 and σ ∼ 0.1 dex) when the fitting implements these priors, while the inferred Z is significantly offset by around 0.6 dex when these priors are not implemented. In the analysis with the L + E and E data cubes, fitting that implements the added priors results in an offset of around −0.12 dex, and −0.15 dex, respectively. Despite the offset, there is a good proportionality between the values inferred from SED fitting and the true ones. Without the priors, Z is offset by around 0.5 dex. For the age, the new priors can significantly improve the fitting result for the simulated G + L + E data cube but less so for the other two data cubes. In the analysis with the E data cube, age is poorly recovered. This is mainly due to the lack of coverage in the optical and UV, which is essential for constraining the stellar population age.

In future analyses, we plan to further test our ability to recover stellar population parameters by exploring a wider wavelength range, including medium and far-infrared data. With these data, we will further assess our capacity to constrain dust attenuation and the extinction law. Additionally, we will investigate the importance of using single-metallicity models in the analysis of more complex populations, and if needed, allow for more flexibility in star-formation histories with non-parametric descriptions. Spectroscopic measurements will be incorporated where available and help assess the retrieved properties in real, rather than simulated, galaxies.

Using our spatially resolved SED fitting pipeline, we aim to map the stellar population properties of a large number (100 000 s) of galaxies in the local Universe using imaging data from Euclid, combined with optical images from other telescopes, especially LSST and UNIONS surveys. The depth and high spatial resolution of the imaging data, combined with the wide-area coverage of the surveys, make it possible to perform our analysis for many galaxies that have a wide range of global properties and reside in various local density environments, from a low-density field to a dense galaxy cluster. We will produce comprehensive and detailed sets of stellar mass profiles and stellar population gradients across various galactic scales, from centres to peripheries. This will provide constraints on in-situ and ex-situ stellar populations, and the role of feedback processes and mergers as a function of galaxy mass (e.g. Tortora et al. 2010, 2011). Euclid is poised to uncover thousands of previously undiscovered dwarf low-surface brightness galaxies and the ultra-faint peripheral regions of large galaxies. Consequently, we will enhance and validate our pipeline designed for determining the stellar populations of these faint objects. This will allow us to comprehensively study the internal properties of galaxies and the influencing factors, both internally and externally, for a deeper understanding of galaxy evolution.

Acknowledgments

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 (http://www.euclid-ec.org). Co-funded by the European Union. Views and opinions expressed are however, those of the author(s) only and do not necessarily reflect those of the European Union. Neither the European Union nor the granting authority can be held responsible for them. CT acknowledges the INAF grant 2022 LEMON. JHK acknowledges grant PID2022-136505NB-I00 funded by MCIN/AEI/10.13039/501100011033 and EU, ERDF. This project makes use of the MaNGA-Pipe3D data products. We thank the IA-UNAM MaNGA team for creating this catalogue, and the Conacyt Project CB-285080 for supporting them.

References

Appendix A: Constructing informative priors in the forms of scaling relations on spatially resolved scales

A prior can be constructed based on empirical data from previous observations that are deemed reliable. In the study of spatially resolved properties of local galaxies, it has been known that Σ_∗ correlates strongly with other key physical properties, such as SFR surface density, Z, stellar population age, and gas-phase metallicity (see review by Sánchez 2020). These empirical scaling relations can be used as a prior in SED fitting. Previous studies have applied this on a global scale. For example, Leja et al. (2019) adopted the mass-Z relation observed by Gallazzi et al. (2005) as a prior in their SED fitting. We extend this attempt to spatially resolved scales.

To construct a prior in the form of mass-age and mass-Z relations (i.e. mass-Z-age prior), first, we combine observational data taken from the MaNGA survey and simulated data from the TNG50-SKIRT database. The MaNGA data is in the form of 2D maps of stellar population properties of 10 220 local (0.01 < z < 0.15) galaxies that are produced by the pipe3D pipeline (Sánchez et al. 2022; Lacerda et al. 2022).⁶ We retrieve the maps of all the galaxies and combine their spaxels, which totals 16 250 999. For the simulated data, we combine pixels from all maps of 1 160 galaxies produced by the TNG50-SKIRT Atlas. For simplicity and to limit the number of pixels, we only use the maps from one viewing angle: “O1”. The total of collected pixels is 488 235 803. These are original simulated maps, without post-processing as done for the synthetic imaging data (see Sect. 2.1).

Using these data, we examine the scaling relations between Σ_∗, age, and Z in both observations and simulations. These mass-Z-age relations are shown in Fig. A.1. Tight relationships are evident from the density contours and the median profiles of both simulation and empirical data, especially for the mass-Z and mass-age relations. The age-Z relation exhibits greater scatter than the other two relations. Its median profiles are derived from combined median profiles of mass-Z and mass-age relations rather than from the data directly. This is to illustrate how the mass-Z and mass-age priors behave on the Z-age plane. Further analysis of these scaling relations is beyond the scope of this paper and we defer this for future work.

The deviation in the median profiles of the mass–age and mass–Z relations from MaNGA and TNG50 primarily occurs in low Σ_∗ regions, below ∼ 10⁷ M_⊙ kpc⁻². In higher mass density regions, the relations from MaNGA and TNG50 appear consistent, albeit with slightly different normalizations. The contours of MaNGA data reveal that most spaxels lie above this density threshold, and the majority of our sample galaxies exhibit spatially resolved mass densities above this threshold (see Figures C.1–C.3).

From these relations, we construct priors by further taking the median between the MaNGA and TNG50 median profiles for mass bin above Σ_* = 10⁶ M_⊙ kpc⁻² and then fit the median profile with a second-degree polynomial function. The prior is defined to be a Gaussian function along the y axis, centred at the median profile (for each point of Σ_∗) and standard deviation of 0.3 dex.

In Fig. A.2, we have also explored the impact of these priors on the results of the paper shown in Figs. 7, 8, and 9. In the upper row, the location of the retrieved properties relative to the prior, shown in blue when the mass-Z-age prior is implemented, is comparable to the location seen for the TNG50 models in Fig. A.1, and suggests that the prior is efficient without being the sole driver of the results when G+L+E data are used (the retrieved parameters are not exactly located along the ridge-line of the prior). In the lower two panels, the results are closer to the prior, indicating the role of the prior relative to other constraints is becoming predominant.

Fig. A.1. Constructing informative priors in the forms of mass-Z-age relations based on observational data from the MaNGA survey and simulations from TNG50. The median profiles are shown by solid black lines, which we use for priors. The dashed black lines indicate ±0.2 around the median profiles. The rightmost panel shows age-Z relation. Its median profiles are derived from combined median profiles of mass-Z and mass-age relations to illustrate how the mass-Z and mass-age priors behave on the Z-age plane.

Fig. A.2. Distributions of derived surface mass density, Z, and age values, compared with the mass–Z–age priors. From top to bottom, results for G + L + E, L + E and E-only data cubes are shown. The solid lines show the median prior, while the dashed lines show the ±0.2 dex range. The colour scheme for SED fitting results is similar to Figs. 7, 8, and 9.

Appendix B: SED fitting tests using simple mock SEDs

The complexity of the radiative transfer process in generating synthetic images might distort the SEDs of individual pixels and change them from being pure spectra of the model of a composite stellar population (CSP). This complexity arises from the intricate spatial configuration of stars and dust within galaxies, which shapes the composite SED of each pixel in the synthetic images produced through radiative transfer modelling. In a dust-free environment, the pixel’s SED is simply the sum of the spectra from the stellar populations within that pixel. However, when dust is present, the situation becomes more complex: the amount of attenuation affecting each spectrum varies depending on the distribution of dust particles (and the resulting cumulative optical depth) along the line of sight to each source.

This can make our fitting tests on the synthetic images non-trivial. To further test our SED fitting pipeline, especially our new priors, we perform additional SED fitting tests using simple mock SEDs generated from CSP models. We generate mock SEDs using FSPS. Instead of generating random parameters, the properties of mock data are adopted from the inferred properties of spatial bins that are randomly chosen from our simulated galaxy sample. For this, we use results from the analysis on the G + L + E data, in which the fitting uses the FSPS model and applies the mass-Z-age prior. We randomly chose 40 spatial bins from each galaxy to take their properties and produce 1 000 mock SEDs. To check the effect of wavelength coverage (by the photometric data), we make three data types (G + L + E, L + E, and E) by simply removing the photometry associated with the filters. We assign S/N of 5 in all filters. The flux is perturbed by adding Gaussian noise with a standard deviation equal to Flux/5. The flux uncertainty is then computed based on the perturbed flux value.

We run fitting to the mock SEDs using the same set-ups as those for the tests on synthetic images, described in Sect. 3.2.3. To check the effect of the mass-Z-age prior, we fit each of the mock SEDs twice, one applies the priors and the other does not. The results, obtained with G + L + E, L + E, and E data cubes, are shown in Fig. B.1. As we can see from this figure, M_∗ is recovered very well, and the M_∗ estimate does not seem to be affected by whether or not the priors are applied. This is in part due to the coverage of NIR from Euclid.

Fig. B.1. Results of SED fitting tests using simple mock SEDs generated with FSPS. Three sets of mock SEDs are made: G + L + E, L + E, and E, and the results of fitting on them are shown in the top, middle, and bottom rows, respectively. We also fit the mock data with and without applying the mass-age and mass-Z priors and the results are shown in blue and red colours, respectively. The histograms in the insets show ratios (in logarithmic scale) between the best-fit and true values.

Without the aid of the mass-Z-age priors, age and Z are barely recovered. Similar to the case with synthetic TNG images, both parameters are underestimated. The worst-case scenario occurs when only Euclid bands are used. In this case, ages cluster around 4-5 Gyr (i.e. centred around the peak of the mass-weighted age distribution in the SPS models), and metallicities are concentrated around 0.5 − 0.6 Z_⊙. For data cubes with a wider wavelength coverage, age estimation improves, while metallicities show no further improvement. In all data cubes, to counterbalance the underestimated ages and metallicities, larger dust attenuations (not shown) are generated in the SED fitting procedure.

Instead, age and Z are recovered relatively well with the three datasets and by applying the priors. Z is recovered equally well with the three datasets, in contrast to the results of fitting tests to the synthetic images, where L + E and E data cubes give ∼0.1 dex offset (see Figs. 8 and 9). On the other hand, using the priors, age is recovered well above ∼5 Gyr with the 3 data cubes, although a slight offset is present in the results with E data. This trend is also different from the results obtained with synthetic images, where age is barely recovered with L + E and E data cube.

Appendix C: Maps of the spatially resolved stellar population properties

Maps of the spatially resolved stellar population properties of 24 galaxies in our sample, derived from the analyses with the three data cubes (G + L + E, L + E, and E) are shown in Figs. C.1, C.2, and C.3, respectively. The maps of the other galaxy TNG501725 are shown in Fig. 6.

Fig. C.1. Maps of the spatially resolved stellar population properties of 24 galaxies in our sample obtained from analysis on G + L + E data cubes. The maps include stellar mass surface density (Σ∗; a), age (b), and metallicity (c).

Fig. C.2. Same as Fig. C.1 but for the results obtained from the analysis on L + E data cubes.

Fig. C.3. Same as Fig. C.1 but for the results obtained from the analysis on E data cubes.

All Tables

All Figures

	Fig. 1. The Sample of 25 simulated galaxies used for the analyses in this paper. The selected sample and the parent sample from the TNG50-SKIRT Atlas are shown with red squares and blue dots, respectively. The sample is selected to cover the entire mass range in the TNG50-SKIRT catalogue and a wide range of sSFR.
In the text

Fig. 2. Examples of simulated images of GALEX, LSST, and Euclid to which observational effects have been added, including spatial resampling (i.e. regridding), simulated noise, and convolution with the PSF of the corresponding cameras. The galaxies here are TNG501725 with O1 orientation index and TNG414917 with O4 orientation index. The original synthetic images are taken from the TNG50-SKIRT Atlas database and are generated from the TNG50 simulation using the SKIRT radiative transfer code. The variation in spatial resolution and depth can be seen in the images with GALEX images being the most blurred and shallowest, and IE images being the sharpest and deepest. The pixel size difference is demonstrated with the red horizontal line (60″ long). The PSF FWHM and pixel size of the images are summarized in Table 1. The colour images are made from the composite of LSST g, r, and i images.

In the text

Fig. 3. Examples of the maps of flux densities produced from the image processing on the sets of three data cubes: G + L + E (a), L + E (b), and E (c). This example of the galaxy is TNG501725 with the O1 orientation index. The final PSF size and pixel size of the 3 data cubes are summarized in Table 2. The G + L + E data cube has the lowest spatial resolution, trading the high spatial resolution of Euclid and LSST with the wider wavelength coverage from UV to NIR.

In the text

	Fig. 4. Surface brightness radial profiles of TNG501725 derived from the stamp images of G + L + E data. The vertical dashed grey line indicates the radial extent of the data cube after the galaxy’s RoI is defined, roughly corresponding to a detection threshold of 2 above the background RMS of Euclid images.
In the text

Fig. 5. Examples of the pixel binning maps of TNG501725 obtained with the three data cubes: G + L + E (a), L + E (b), and Euclid only (c). In each row, the three panels from left to right show the binning map, S/N radial profile of pixels, and S/N radial profile of spatial bins. The pixel binning scheme of piXedfit can achieve the S/N threshold of 5 (indicated by the dashed black lines) for the g, r, i, z, IE, YE, JE, and HE filters.

In the text

Fig. 6. Examples of the maps of the stellar population properties of TNG501725 obtained from the spatially resolved SED fitting on the three data cubes: G + L + E (a), L + E (b), and E (c). The maps derived from the analysis with the L + E and E data cubes alone have a high spatial resolution and reveal detailed spatial variation in the physical properties compared to those obtained from the analysis with the G + L + E data cubes.

In the text

Fig. 7. Comparisons of the spatially resolved stellar population properties derived from SED fitting on the G + L + E data cubes (ordinate) and the true properties (abscissa). These include all spatial bins from the whole sample of galaxies. We focus on four main properties in this test analysis: Σ∗, Z, age, and ΣSFR. The top row (a) shows results obtained with the FSPS model, while the bottom row (b) shows results using the BC16 model. The blue and red colours represent results obtained with and without applying the mass-Z-age prior.

In the text

	Fig. 8. Similar to Fig. 7 but for results of SED fitting on the L + E data cubes.
In the text

	Fig. 9. Similar to Fig. 7 but for results of SED fitting on the E data cube.
In the text

Fig. 10. Examples of spatially resolved distributions of the systematic offsets (in logarithmic scale) between the parameters derived from SED fitting and the true values. Offset maps of a galaxy TNG501725, are shown here for example. The first, second, and third rows are for analyses with G + L + E, L + E, and E data cubes, respectively. Only results obtained with the FSPS model and applying the mass-Z-age prior are shown.

In the text

Fig. 11. The Systematic offsets between estimated and input properties (μ) as a function of the radial distance from the galactic centre. The radial profiles of the offsets, as obtained from fits with the FSPS models and with the mass-Z-age prior, are shown for all sample galaxies with three types of data cubes. Columns from left to right show results from the analysis with the G + L + E, L + E, and E data cubes. Each galaxy’s radial extent is divided into ten radial bins equally spaced and radial profiles are then derived by taking the median of μ in circular annuli. The radius is normalized to the half-mass radius that is taken from the catalogue of the TNG50-SKIRT Atlas. The two horizontal dashed grey lines mark ±0.3 dex offset.

In the text

Fig. 12. The Robustness of integrated M∗ derived from the spatially resolved SED fitting and the integrated SED fitting compared to the true M∗ from the TNG50 simulation. Left panel: Comparison between M∗ derived from spatially resolved SED fitting (i.e. by summing up M∗ of pixels within the ROI of the galaxies) and the true values. Middle panel: Comparison between M∗ derived from integrated SED fitting and the true values. Right panel: Comparison between integrated M∗ from spatially resolved and integrated SED fitting.

In the text

Fig. A.1. Constructing informative priors in the forms of mass-Z-age relations based on observational data from the MaNGA survey and simulations from TNG50. The median profiles are shown by solid black lines, which we use for priors. The dashed black lines indicate ±0.2 around the median profiles. The rightmost panel shows age-Z relation. Its median profiles are derived from combined median profiles of mass-Z and mass-age relations to illustrate how the mass-Z and mass-age priors behave on the Z-age plane.

In the text

	Fig. A.2. Distributions of derived surface mass density, Z, and age values, compared with the mass–Z–age priors. From top to bottom, results for G + L + E, L + E and E-only data cubes are shown. The solid lines show the median prior, while the dashed lines show the ±0.2 dex range. The colour scheme for SED fitting results is similar to Figs. 7, 8, and 9.
In the text

Fig. B.1. Results of SED fitting tests using simple mock SEDs generated with FSPS. Three sets of mock SEDs are made: G + L + E, L + E, and E, and the results of fitting on them are shown in the top, middle, and bottom rows, respectively. We also fit the mock data with and without applying the mass-age and mass-Z priors and the results are shown in blue and red colours, respectively. The histograms in the insets show ratios (in logarithmic scale) between the best-fit and true values.

In the text

	Fig. C.1. Maps of the spatially resolved stellar population properties of 24 galaxies in our sample obtained from analysis on G + L + E data cubes. The maps include stellar mass surface density (Σ∗; a), age (b), and metallicity (c).
In the text

	Fig. C.2. Same as Fig. C.1 but for the results obtained from the analysis on L + E data cubes.
In the text

	Fig. C.3. Same as Fig. C.1 but for the results obtained from the analysis on E data cubes.
In the text