© The Authors 2023

1. Introduction

Understanding the process by which galaxies assemble their stellar mass through cosmic time is a crucial aspect of galaxy evolution. This process can be modulated either by external (e.g., galaxy merging) or internal (e.g., secular evolution, star formation) mechanisms (Kormendy 2013). These different mechanisms ultimately lead to the wide variety of galaxy properties that we see in the present-day Universe, such as their morphologies and the demographics of their stellar populations. The study of stellar populations through the fossil method (Tinsley 1968), which consists of reproducing the observed spectra with a linear combination of single stellar populations (SSPs) of known ages and metallicities, has been demonstrated to be a powerful tool in unveiling the assembly histories of galaxies (e.g., Conroy 2013; McDermid et al. 2015; Wilkinson et al. 2015; López Fernández et al. 2018; Zhuang et al. 2019; Neumann et al. 2020).

Although the investigation of stellar populations has been used as an approach to address the evolution of galaxies for many years (Sarzi et al. 2005; Rogers et al. 2007), large integral field spectroscopy (IFS or IFU) surveys, such as CALIFA (Sánchez et al. 2012), MaNGA (Bundy et al. 2015), and SAMI (Croom et al. 2012), have further improved these analyses, enabling the simultaneous extraction of spectra from different regions in nearby galaxies, with typical spatial resolutions of ∼1 kpc, which has allowed for the characterization of stellar populations and ionized gas properties across galactic disks. This acknowledges the fact that galaxies are extended and complex objects of which the stellar population and gas properties might change drastically from one specific region to another.

One of the most studied features in the context of spatially resolved stellar populations is the galactocentric radial distribution of their properties such as age and metallicity (see, e.g., Sánchez-Blázquez et al. 2014a; González Delgado et al. 2014a, 2016; Ibarra-Medel et al. 2016; Zheng et al. 2017; García-Benito et al. 2017; López Fernández et al. 2018; Zhuang et al. 2019; Dale et al. 2020; Parikh et al. 2021). These studies have revealed that radial gradients of stellar age and metallicity are primarily connected with the Hubble type of the galaxy (Goddard et al. 2017; Parikh et al. 2021; Smith et al. 2022), with late-type galaxies (LTGs) showing negative age and metallicity gradients (i.e., centers are older and more metal-rich), and early-type galaxies (ETGs) showing a nearly flat age gradient, and a negative metallicity gradient (albeit flatter than that found for LTGs). These works also report a dependence of the gradients on total stellar mass, where more massive galaxies exhibit steeper gradients.

In terms of galaxy evolution, negative age gradients point to an inside-out growth of galaxies, where the inner regions assembled their stellar mass earlier than outer regions. In this regard, Ibarra-Medel et al. (2016) found that LTGs would have a more pronounced inside-out formation mode, compared to ETGs (i.e., a more obvious difference between the formation times of inner and outer regions), and that low-mass galaxies (M_* < 5 × 10⁹) show a large diversity in their radial assembly history.

With respect to the origin of the metallicity gradients, Zhuang et al. (2019) find that they are a reflection of the local stellar mass surface density-metallicity relation (González Delgado et al. 2014b), and conclude that the spatial distribution of stellar populations within a galaxy is primarily the result of the insitu local star formation history (SFH), rather than being shaped by radial migration. However, Neumann et al. (2021) report that the local correlation between stellar metallicity and stellar mass surface density (Σ_*) shows an increased scatter toward the outer parts of galaxies, suggesting that something else, besides Σ_*, drives chemical enrichment at larger galactocentric distances (e.g., gas accretion, outflows). Zhuang et al. (2019) also find that low-mass LTGs show more commonly positive metallicity gradients, consistent with stellar feedback being more efficient at regulating the baryon cycle in the central regions of these galaxies.

Beyond radial gradients of stellar age and metallicity, recent studies have aimed to study the assembly history of the galactic disk in a 2D manner (Guérou et al. 2016; Peterken et al. 2019, 2020; Pinna et al. 2019a,b), by “time-slicing” galaxies across their lifetime (i.e., studying the spatial distribution across the galaxy of stars that were born at a given epoch). However, the limited spatial resolutions of large IFU surveys place limitations on this kind of study. At a spatial resolution of ∼1 kpc, much of the galactic structure is not resolved (e.g., spiral arms typically have widths of a few hundred parsecs; Querejeta et al. 2021). On the other hand, higher resolution measurements (e.g., Guérou et al. 2016) are often limited by small sample sizes.

In this paper, we present what we learn from spatially resolved SFHs of a sample of LTGs from the Physics at High Angular resolution in Nearby GalaxieS (PHANGS¹) survey (Leroy et al. 2021) survey, measured at spatial scales of ∼100 pc, at which different morphological components are clearly resolved and can be studied separately. Thus, we address how the properties of stellar populations vary as a function of the local galactic environment and assess the importance/role of galactic structure in the assembly history of a sample of nearby star forming galaxies. It is worth noticing that while there are other optical IFU surveys with similar spatial resolution that allow the study of the stellar population properties of nearby galaxies, such as TIMER (Gadotti et al. 2019) and Fornax3D (Sarzi et al. 2018; Iodice et al. 2019), they were designed to address specific science questions and their samples are therefore not representative of the population of star-formation main-sequence galaxies in the local Universe. Here we focus primarily on a detailed characterization of the radial trends measured in the stellar population properties of star forming galaxies, and how these trends correlate with global galactic properties. We also show variations in the SFH across the galactic disk of galaxies and study the correlation between kinematic properties and age of their stellar populations.

This paper is structured as follows. In Sect. 2 we present the data and data products used in our analysis. In Sects. 3.1–3.3 we present and discuss the radial distribution of the stellar population properties of the galaxies in our sample, for different local galactic environments. In Sect. 3.4 we show how the kinematic properties of stellar populations change as a function of stellar age, and in Sect. 3.5 we show how the SFHs of galaxies varies across different galactic regions (i.e., not radially integrated). We present a summary and the conclusions of our analysis in Sect. 4, and we validate a revised set of stellar templates to improve fitting to the youngest regions in Appendix E.

2. Data

We use a sample of 19 star-forming galaxies, all of them close to the star-forming main sequence of galaxies (SFMS; e.g., Brinchmann et al. 2004; Daddi et al. 2007; Noeske et al. 2007). These galaxies represent a subsample of the PHANGS survey (Leroy et al. 2021), and are observed as part of the PHANGS-MUSE survey (PI: E. Schinnerer; Emsellem et al. 2022). The PHANGS galaxies have been chosen to constitute a representative set of star-forming disk galaxies where most of the star formation is occurring in the local Universe. They have been selected to be at distances smaller than 20 Mpc to resolve the typical scale of star-forming regions (50 − 100 pc), moderately inclined (i < 60°) to limit the effect of extinction and allow the identification of individual star-forming regions. Their stellar mass limit is defined by log₁₀(M_*/M_⊙)≳9.75, which translates to nearly twice the mass of the LMC, aiming to avoid focusing on low metallicity, low mass galaxies, where the detection of CO becomes increasingly harder (e.g., Schruba et al. 2012, 2017; Bolatto et al. 2013; Cormier et al. 2014). We note that this mass cutoff does not strictly apply to the PHANGS-MUSE subsample because a revision of the fiducial distance estimates (Anand et al. 2021). Table 1 summarizes the properties our sample. We use the global parameters reported by Leroy et al. (2021) and the distance compilation of Anand et al. (2021). The inclination values adopted are those reported by Lang et al. (2020).

2.1. VLT/MUSE

We use data from the PHANGS-MUSE survey. This survey employs the Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2014) optical integral field unit (IFU) mounted on the VLT UT4 to mosaic the star-forming disk of 19 galaxies from the PHANGS sample.

The mosaics combine three to 15 individual MUSE pointings. Each MUSE pointing provides a 1′×1′ field of view sampled at $0\stackrel{″}{.}2$ per pixel, with a typical spectral resolution of ∼2.5 Å FWHM (∼100 km s⁻¹) covering the wavelength range of 4800 − 9300 Å, and with angular resolutions ranging from $\sim 0\stackrel{″}{.}5$ to $\sim 1\stackrel{″}{.}0$ for the targets observed with and without adaptive optics (AO), respectively. The nine galaxies observed with AO are marked with a black dot in the first column of Table 1. The total on-source exposure time per pointing for galaxies in the PHANGS-MUSE Large Program is 43 min. Observations were reduced using recipes from the MUSE data reduction pipeline provided by the MUSE consortium (Weilbacher et al. 2020), executed with ESOREX using the python wrapper developed by the PHANGS team² (Emsellem et al. 2022). The point spread function of the different individual MUSE exposures that are mosaicked together have been homogenized to the largest FWHM within each mosaic (i.e., dataset labeled as “copt” in the public release, see Emsellem et al. 2022, for a detailed description of how the homogenization is done). The total area surveyed by each mosaic ranges from 23 to 441 kpc². Once the data have been reduced, we have used the PHANGS data analysis pipeline (DAP) to derive various physical quantities. Some of these outputs are described in Sect. 2.2. The DAP is described in detail in Emsellem et al. (2022). It consists of a series of modules that perform single stellar population (SSP) fitting and emission line measurements to the full MUSE mosaic.

2.2. Stellar population property maps

The PHANGS-MUSE DAP (Emsellem et al. 2022) includes a stellar population fitting module, which uses a linear combination of SSP templates of known ages, metallicities, and mass-to-light ratios is used to reproduce the observed spectrum. This allows us to infer stellar population properties from an integrated spectrum, such as mass- or luminosity-weighted ages, metallicities, and total stellar masses, together with the underlying SFH.

The full mosaics are initially corrected for Milky Way extinction, using an E(B − V) value obtained from the NASA/IPAC Infrared Science Archive³ (Schlafly & Finkbeiner 2011), and assuming a Cardelli et al. (1989) extinction law. In detail, our spectral fitting pipeline performs the following steps:

First, we apply a Voronoi tessellation (Cappellari & Copin 2003) to bin our MUSE data to a minimum signal-to-noise ratio (S/N) of ∼35, computed at the wavelength range of 5300 − 5500 Å. This value is chosen in order to keep the relative uncertainty in our mass measurements below 15%, even for pixels dominated by a younger stellar population. To do this, we tried different S/N levels to bin a fixed region in our sample, and we bootstrapped our data to have an estimate of the uncertainties at each S/N level.

We use then the Penalized Pixel-Fitting (pPXF) code (Cappellari & Emsellem 2004; Cappellari 2017) to fit the spectrum of each Voronoi bin. We fit the wavelength range 4850 − 7000 Å, in order to avoid spectral regions strongly affected by sky residuals. We further mask the regions around strong emission lines, as well as regions of the spectrum particularly affected by sky line residuals.

We fit our data with a grid of templates consisting of 17 ages, ranging from 6.3 Myr to 15 Gyr, logarithmically-spaced, and four metallicity bins [Z/H] = [ − 0.7, −0.4, 0, 0.22]. The templates are from the E-MILES (Vazdekis et al. 2010, 2012) database, assuming a Chabrier (2003) IMF and Padova isochrones (Girardi et al. 2000). The E-MILES templates were originally computed to start from a minimum age of 63 Myr (for the Padova isochrones). We have complemented this original set of templates with the young extension to the E-MILES SSP models presented in Asa’d et al. (2017), adding five additional age bins, from 6.3 Myr to ∼40 Myr. These young templates have been computed for the same metallicity bins, except for the most metal-rich one, which for the young-extension templates has [Z/H] = 0.41 (instead of 0.22).

This choice of templates is different from that used in Emsellem et al. (2022) for the same data. In Appendix E we present an exhaustive exploration of the parameter space in order to improve the outcome of the spectral fitting of the regions dominated by young (< 400 Myr) stellar populations, and justify this choice. We also refer the reader to Emsellem et al. (2022) for a detailed description of this issue.

We implemented a two-step fitting process. First, we fit our data assuming a Calzetti et al. (2000) extinction law to correct for internal attenuation. We then corrected the observed spectrum using the measured attenuation value before fitting it a second time, including a 12th degree multiplicative polynomial instead of attenuation in this iteration. This two-step fitting process accounts for offsets between individual MUSE pointings. These offsets arise because the different MUSE pointings were not necessarily observed under identical weather conditions, and variations in the sky continuum levels can lead to subtle differences in the flux calibration of individual neighboring MUSE pointing. An analysis of the regions of the mosaic where different pointings overlap revealed that these variations yield differences in the measured stellar extinction levels on the order of ΔE(B − V)∼0.04 mag (see Emsellem et al. 2022, for a detailed description of this issue). Therefore, in the first iteration of the spectral fitting, we measure a reddening value and correct the observed spectra accordingly, and in the second iteration, we use a high-degree multiplicative polynomial to correct for those nonphysical features and homogenize the outcome of the different pointings.

To recover the kinematic properties of the observed spectra, the templates are shifted in velocity space and convolved to match observed spectral features. We fit the first four velocity moments, using a single stellar kinematic component, that is, the same velocity and velocity dispersion, h3 and h4 is applied to all the templates used for the spectral fitting. The spectral resolution of the E-MILES templates is higher than that of the VLT/MUSE data within the wavelength range considered (although at ∼7000 Å they are virtually the same). Thus, templates are convolved to match the spectral resolution of the data, using an appropriate wavelength-dependent kernel, before the fit.

The output of pPXF consists of a vector with the coefficients of the linear combination of templates that best reproduce the observed spectrum. Physically, these weights represent the mass fraction of stars with a given age and metallicity, and they are used to derive stellar mass surface densities, and both light- and mass-weighted ages and metallicities in each spaxel. The stellar mass surface density maps include contributions from live stars and remnants. The total weight for a given age bin, integrating for the different metallicities, also represents the SFH of a given spectrum. For each pixel, the average age and metallicity are computed as follow:

$$\begin{array}{c}\hfill ⟨log\mathrm{age}⟩=\frac{{\mathrm{\Sigma }}_{i}log({\mathrm{age}}_i)w_i}{{\mathrm{\Sigma }}_i{w}_i}\end{array}$$(1)

and

$$\begin{array}{c}\hfill ⟨[\mathrm{Z}/\mathrm{H}]⟩=\frac{{\mathrm{\Sigma }}_i{[\mathrm{Z}/\mathrm{H}]}_i{w}_i}{{\mathrm{\Sigma }}_i{w}_i},\end{array}$$(2)

where age_i and [Z/H]_i correspond to the age and metallicity of each template, and w_i is its corresponding weight in the linear combination. To convert mass-weighted quantities to luminosity-weighted quantities, we use the mass-to-light ratio of each template in the V-band. We compute luminosity-fraction weights (${\mathit{w}}_i^{\mathrm{LumW}}$) of a given template as

$$\begin{array}{c}\hfill w_i^{\mathrm{LumW}}=\frac{w_i}{{(M/L_V)}_i},\end{array}$$(3)

where (M/L_V)_i corresponds to its mass-to-light ratio in the V-band. We can use these luminosity-fraction weights to calculate luminosity-weighted properties, following Eqs. (1) and (2).

Regularization is a standard approach to solve ill-posed problems (see Cappellari 2017) and it is commonly used to reduce the noise in the recovery of the stellar population parameters. However, exhaustive testing showed that using a fixed level of regularization on our dataset leads to strongly biased star-formation histories, with systematic differences across different regions. This is due to the presence of young star forming regions being inhomogeneously distributed across the disk, causing strong differences in the stellar ages from one region to another. Thus, we do not use regularization in our fitting, and instead we use Monte Carlo simulations to estimate the uncertainty in the recovered stellar population parameters. For each spectrum, we perform 20 Monte Carlo iterations, where in each iteration, we perturb the input spectrum assuming a Gaussian noise with a mean of zero and a standard deviation corresponding to the spectral error at each wavelength bin. The uncertainties of stellar population parameters are calculated as the standard deviation of their distributions produced by the Monte Carlo realizations. This is meant as a first-order estimate of the true uncertainties. In Appendix B, we show the distribution of the relative uncertainties in stellar mass surface density, stellar age, and metallicity, for our full sample. We also demonstrate that 20 Monte Carlo realizations is a reasonable compromise between estimating the width of the posterior distribution and computational time to achieve our results. We only run Monte Carlo realizations for the second step of the fitting procedure, once the stellar extinction has been fixed. Hence, no error is computed for the stellar E(B − V). Finally, we have identified foreground stars as velocity outliers in the SSP fitting, and we have masked those pixels for the analysis carried out in this paper.

2.3. Environmental masks

We employ the environmental masks described in Querejeta et al. (2021) to classify the different environments of each galaxy and label them as disks, spiral arms, rings, bars, and centers (see, e.g., Fig. 1). This classification was done using photometric data, mostly from the Spitzer Survey of Stellar structure in Galaxies (S⁴G; Sheth et al. 2010). We refer the reader to that paper for a detailed explanation of how the masks are defined.

Fig. 1. Example of some of the DAP outputs obtained for the mosaic of NGC 1566. Stellar mass surface density (top-left), luminosity-weighted age (top-center), luminosity-weighted [Z/H] (top-right), stellar velocity dispersion (bottom-left), mass-weighted age (bottom-center), and mass-weighted [Z/H] (bottom-right). The spiral arms and bar of the galaxy are marked with blue and green contours, respectively.

In summary, disks and centers are identified via 2D photometric decomposition of 3.6 μm images (see, e.g., Salo et al. 2015). Centers are defined as a central excess of light, independently of its surface brightness profile. If a central excess of light is not identified, the inner galactic region is not labeled as center. Bars and rings properties are defined visually on the NIR images, following Herrera-Endoqui et al. (2015) for S⁴G galaxies. Spiral arms are defined by fitting a log-spiral function to bright regions along the arms on the NIR images, only when the spiral arms are clear features across the galaxy, and their width is empirically determined based on CO emission. We use these environmental masks to probe the properties of stellar populations separately across different galactic environments.

2.4. ALMA CO(2–1) Maps

We compare our measurements to ALMA CO-based measurements of the velocity dispersion. We use the “effective width”-based estimates, which estimate the line width based on the observed peak intensity and line-integrated intensity. These effective width measurements are robust in the presence of noise and are expressed as an equivalent rms velocity dispersion assuming a Gaussian line profile. We expect that these provide a reasonable first-order estimate the local molecular gas velocity dispersion across most of the area in our target galaxies. The primary shortcoming of these measurements will arise in the central parts of galaxies where complex line profiles can mean that the interpretation of the measurement in terms of a single Gaussian line profile represents an oversimplication.

3. Results and discussion

In Sect. 2.2, we described in detail the methodology adopted to fit the spectra of galaxies in our sample with a linear combination of SSPs to derive their SFHs, in a spatially resolved manner. In the following sections, we show how we can use this information to gain insights into galaxy assembly and evolution.

3.1. Radial structure of stellar population properties

In this section, we show the radial distribution of the light- and mass-weighted mean age and metallicity of the stellar populations of PHANGS-MUSE galaxies, at a resolution of ∼100 pc, thus, resolving their different morphological components. While luminosity-weighted quantities are very sensitive to recent star formation events, because young stars dominate the B-band light emission, and thus, encode information about the ongoing galaxy evolution, mass-weighted quantities trace the properties of the long-lived, low-mass stars that dominate the stellar mass budget, and thus, encode information about the mass growth and secular redistribution history of galaxies. Figure 1 shows the stellar mass surface density, velocity dispersion, mean luminosity-weighted (LumW) and mass-weighted (MassW) age, and mean light- and mass-weighted metallicity for a galaxy in our sample (NGC 1566), as an example of our mean stellar population properties maps. Overplotted contours enclose the spiral structure and the bar environments (as defined in Querejeta et al. 2021).

We derive the radial variation of these parameters, averaging the quantities azimuthally in concentric rings. Each radial bin has a width of 0.03 R₂₅, where R₂₅ represents the 25th isophotal magnitude radius in B-band. This radial bin corresponds to a physical size of about 0.4 ± 0.2 kpc for the galaxies in our sample. To define the radial bins consistently, we use the galactic inclinations and position angles reported in Lang et al. (2020), measured from CO(2−1) data. Only those radial bins in which at least one-third of the ring is covered by the MUSE mosaic are considered valid for the radial profiles.

3.1.1. Age radial profiles

Figure 2 shows the LumW age profiles for the 19 galaxies in our sample, for the full field of view (FoV), and each morphological component separately. LumW quantities are strongly biased toward young populations (e.g., Zibetti et al. 2017). Therefore, sudden declines in these profiles are associated with recent local star formation.

Fig. 2. Luminosity-weighted stellar age radial profiles for the galaxies in our sample. Different colors indicate the radial profile measured across different environments, as indicated in the legend of the bottom-right panel. The black line shows the radial profile measured for the entire FoV (i.e., all environments together). The galactocentric distance is measured in units of R25, in order to measure radial distance homogeneously across our sample. The value of R25 (kpc) of each galaxy is indicated in each corresponding panel. The solid gray line shows the best-fit linear gradient for each galaxy.

The first obvious feature here is the overall decreasing trend with radius present for most of the galaxies in our sample. Negative age gradients have been reported in previous studies (e.g., Moustakas et al. 2010; Sánchez-Blázquez et al. 2014a; González Delgado et al. 2014a). The only exceptions are NGC 1087 and NGC 1385, which are both at the low-mass end of our sample (log M_*/M_⊙ of 9.9 and 10.0, respectively). This is consistent with previous findings of a larger diversity in the radial structure of low-mass galaxies (Ibarra-Medel et al. 2016; Smith et al. 2022). This trend can be more clearly appreciated in the top panels of Fig. 3, where low-mass galaxies (blue) show either flat or negative gradients, while higher mass galaxies (green) show mostly negative gradients. The figure also reveals a clear dependence of the normalization of the LumW age radial profile with respect to the total stellar mass of the galaxies, where more massive galaxies are in general older than lower mass galaxies.

Fig. 3. Luminosity-weighted stellar age radial profiles for the galaxies in our sample, considering all environments. Top left panel shows the measured LumW log Age radial profiles, colored by total galaxy stellar mass. Top right panel shows the profiles normalized by their value at 0.3 R25, to better compare their shapes. The bottom panels show the slope (ΔLumW log Age) and intercept of the LumW radial age profile (indicated with a solid gray line in each panel of Fig. 2) as a function of the total stellar mass of galaxies, color coded by total specific star formation rate (sSFR).

The slope of LumW log age radial gradients span a wide range of values. We measure the slope and intercept of the radial profiles using an ordinary least-squares (OLS) fitting routine to fit a linear model to the full FoV profile (i.e., all environments together), but dropping the innermost radial bins (r < 0.06 R₂₅), that are often outliers in the radial trend. The linear models are indicated with a gray solid line in Fig. 2. We acknowledge the caveat that although the linear models generally capture well the observed radial trend, some galaxies show gradients that are intrinsically hard to capture with a simple linear model (e.g., the log (LumW age) radial profiles of NGC 1566 or NGC 3627). In these galaxies, a single slope is not sufficient to describe the observed gradients, and a more appropriate representation of the data would be provided by a compound model with two different slopes. However, a more sophisticated modeling is beyond the scope of the analyses presented here. Figure 3 shows explicitly how the normalization and slope of the LumW age radial profile depend on the galaxy stellar mass and specific SFR.

Now we look at specific galactic structures. Some galaxies show evidence of recent star formation in their inner-most region (R ≲ 0.2 R₂₅), leading to a sudden decrease of stellar age (e.g., NGC 1300, NGC 1365 among others). Although our sample suggests that this rejuvenation occurs preferentially in the central region of barred galaxies, the small number of nonbarred galaxies in our sample makes it hard to conclude exactly the impact of the bar on this phenomenon. Spiral arms, when present, usually show up as younger structures at any given radius, compared to the disk, indicating that star formation is occurring more often in spiral arms. Bars also show a very clear negative profile. Some galaxies show a sharp drop in LumW age at the outer radius of the bar (NGC 1300, NGC 1365, NGC 1672, NGC 7496), consistent with a scenario where the oldest populations dominate in the inner and rounder part of the bar, while young populations are located on more elongated orbits, populating the bar ends, where they spend a larger fraction of their orbiting period (Beuther et al. 2017; Wozniak 2007; Neumann et al. 2020). The impact of bars in the orbital configuration of stars of different ages, together with the suppression of star formation along bars (e.g., Neumann et al. 2019; Fraser-McKelvie et al. 2020), caused by the bar dynamics (i.e., velocity shocks and shear, see, e.g., Zurita et al. 2004; Díaz-García et al. 2020), could be the origin of the consistently higher LumW log Ages in bars, relative to disks, at a fixed galactocentric distance, across most of the radial range covered by both environments (except for the bar ends), in several of our sample galaxies (e.g., NGC 1365, NGC 1433, NGC 1512, NGC 3351). It is also worth noting that many galaxies exhibit a clear flattening of their radial gradients toward larger radii. This break in the radial gradient roughly coincides with the bar radius, when a bar is present.

There have been other works that have reported differences in the stellar age of different morphological components of late-type galaxies. Méndez-Abreu et al. (2021) use data from CALIFA (Sánchez et al. 2012), and also find older bulges that formed earlier in cosmic history than disks. In another work of spectro-photometric decomposition of bulges and disks, Costantin et al. (2022) find that bulge ages show a bimodal distribution, for which they propose a two-wave formation scenario. Similarly, Coccato et al. (2018) separate bulge and disk of a single spiral galaxy (log M_* ≈ 10.7; Leroy et al. 2008), and find that young stars contribute mostly to the disk, with the young contribution increasing with galactocentric radius. They also find a strong negative age gradient in the disk, but no evidence of a gradient in the bulge, consistent with the results that we find in our sample, where some galaxies show a clear age gradient in their center and others do not (although most of them show a clear negative gradient in their disk). Although there is still no consensus about the origin of the stellar age gradients in bulges (See, e.g., Sánchez-Blázquez 2016, and references therein), theoretical models suggest a link between the radial gradient of stellar population properties of bulges and their formation mechanisms (Matteucci et al. 2019). Morelli et al. (2015a) measured the stellar population gradients in the bulges of a sample of spiral galaxies, and found (on average) flat age gradients, albeit with significant dispersion (σ_gradients = 1.3), which could point to a hierarchical merging contribution in the assembly of bulges. On this subject, several authors agree that drops in stellar ages in the central part of galaxies, leading to positive measured age gradients in bulges, are likely caused by the presence of central disks or nuclear rings (Morelli et al. 2008; Bittner et al. 2020). Within this framework, the high angular resolution of our data allows us to clearly determine the presence of a nuclear ring, and to show that its presence indeed often leads to an overall positive LumW age gradient in the innermost galactic region (e.g., NGC 1365 or NGC 3351, among others). Furthermore, the measurement at high-resolution of gradients in stellar population properties in the center of galaxies, together with stellar kinematics could enable a systematic study of the origin of bulges in late-type galaxies (see, e.g., Moorthy & Holtzman 2006; Rosado-Belza et al. 2020; Gadotti et al. 2019, 2020); however, this is beyond the scope of this paper. Lastly, we remind the reader that, unlike other works that compare bulges and disks, defining bulges according to their surface brightness profile or kinematic properties, our definition of “centers” implies only a central excess of light, as explained in Sect. 2.3.

Differences between the LumW log Age of different environments are summarized in Fig. 4. Bars show a clear peak at old ages (∼ + 0.5 dex offset with respect to the median age of each galaxy), being the environment with the oldest LumW ages on average. Centers also peak at old ages, but they show a bimodal distribution, with a second peak toward younger ages (∼ + 0.2 dex), due to the younger stellar populations present in the central region of some galaxies. Rings show a narrower distribution than disks and spiral arms, suggesting a lower dispersion in the age of the stars within this environment (although we acknowledge that this could be partially due to that the Ring distribution is dominated by the outer ring of a few galaxies). Finally, spiral arms show a distribution similar in shape to that of the disk, but displaced toward slightly younger ages. Table 2 shows the median and standard deviation of the LumW log Age distribution within each environment, and the median offset of each environment with respect to the full galaxy measurement for the full sample in the last row (i.e., median and standard deviation of distributions shown in Fig. 4).

Fig. 4. Distribution of the offset of the LumW log Age of the pixels within each environments, with respect to the median galaxy log LumW Age of each galaxy, in bins of 0.04 dex. Different colors indicate the distribution of different environments, as indicated in the legend.

The results from Querejeta et al. (2021) indicate that a spatial resolution of ∼100 pc is sufficient to separate the different morphological components of galaxies. For comparison, in Appendix D we explore how these results would change with lower angular resolution data, by degrading our dataset to a fixed angular resolution of 15 arcsec (1.1 ± 0.3 kpc for the distances of our sample galaxies), at which contamination among different galactic environments is significantly larger.

Figure 5 shows the mass-weighted (MassW) age profiles for the galaxies in our sample. Since MassW quantities are less biased by recent star formation events, as explained above, MassW ages are higher than their LumW counterpart, as the bulk of the mass of a galaxy is dominated by old stellar populations (see, e.g., Cook et al. 2020).

Fig. 5. Mass-weighted stellar age radial profiles for the galaxies in our sample. Different colors indicate the radial profile measured across different environments, as indicated in the legend of the bottom-right panel. The black line shows the radial profile measured for the entire FoV (i.e., all environments together). The galactocentric distance is measured in units of R25, in order to measure radial distance homogeneously across our sample. The value of R25 (kpc) of each galaxy is indicated in each corresponding panel. The solid gray line shows the best-fit gradient for each galaxy.

The most evident feature in the radial MassW age profiles is that environmental differences are basically absent and that the dynamic range in age is drastically reduced, typically confined to values larger than log₁₀ (age yr⁻¹) ≈ 9.5 (∼3 Gyr). This implies that the distribution of older stellar populations is significantly more homogeneous across the galactic structure than for the younger populations. Clear negative slopes persist in the MassW profile meaning that inner regions assembled the bulk of their stellar mass earlier in cosmic history than outer regions (although they are shallower, consistent with findings from Parikh et al. 2021; Neumann et al. 2020). These age gradients have been interpreted as a spatially resolved imprint of the inside-out (Mo et al. 1998) growth of galaxies by several studies (e.g., González Delgado et al. 2014a; Goddard et al. 2017; Parikh et al. 2021). Figure 6 shows the normalization and slope of the MassW age radial profile as a function of galaxy total stellar mass and specific SFR. There is a clear dependence of the normalization of the MassW age profile with the total stellar mass of galaxies, where more massive galaxies show older MassW ages. This can be interpreted as an imprint of the downsizing in galaxy evolution (Cowie et al. 1996; Heavens et al. 2004; Pérez-González et al. 2008, see also Appendix C).

Fig. 6. Mass-weighted stellar age radial profiles for the galaxies in our sample, considering all environments. Top left panel shows the measured MassW log Age radial profiles, colored by total galaxy stellar mass. Top right panel shows the profiles normalized by their value at 0.3 R25, to better compare their shapes. The bottom panels show the slope (ΔMassW log Age) and intercept of the MassW radial age profile (indicated with a solid gray line in each panel of Fig. 5) as a function of the total stellar mass of galaxies, color coded by total specific star formation rate (sSFR).

Finally, it is worth noting that the break of the age gradient at approximately the bar radius discussed above persists for the MassW age gradients, although it is more subtle, due to the smaller dynamic range of the MassW ages. The sharp decrease of age in the central region also persists in the MassW age profile of some galaxies, suggesting that these centers have formed a substantial amount of stellar mass in relatively recent times. In Sect. 3.2, we investigate this further. One possibility is that the MassW age measurement in the center environment is contaminated by the surrounding ring structure in galaxies where both morphological components are present. However, given that we are able to spatially resolve the nuclear rings at the resolution of our dataset (Querejeta et al. 2021) and we observe this pattern in the innermost (thus furthest from the ring) radial bin of some centers, we believe that contamination from the ring is not the primary driver of these trends. Furthermore, a careful inspection reveals that, while low LumW ages in the inner region of these galaxies strongly correlate with the location of the ring environment, the decrease in the MassW age correlates spatially with the center environment, indicating the presence of secularly evolved nuclear disks (see also Bittner et al. 2020). Flatter or even slightly positive radial gradients of stellar population properties in the innermost regions of galaxies have been reported before in the literature (e.g., Sánchez-Blázquez et al. 2014b, for NGC 0628). In Appendix A.1 we explore the correlation between the slope of the stellar age radial gradient and a set of global galaxy properties. Table 3 shows the same information as Table 2, but for the MassW log age distribution rather than the LumW log age distribution.

3.1.2. Metallicity radial profiles

Regarding stellar metallicity, Figs. 7 and 8 show the LumW and MassW stellar metallicity profiles measured for our sample galaxies. The figures show flatter, or generally slightly negative slopes, with the exception of the same two low-mass galaxies mentioned earlier (NGC 1087, NGC 1385). Overall, negative slopes are consistent with results reported previously in the literature (Goddard et al. 2017; Coenda et al. 2020). Moreover, Morelli et al. (2015b) and Coccato et al. (2018) have studied the stellar metallicity gradients in the morphologically decoupled bulges and disks of spiral galaxies, also finding predominantly negative gradients. Interestingly, NGC 1385 shows a strongly positive metallicity gradient. This is consistent with the findings from Zhuang et al. (2019), who reported that low-mass and late morphological types commonly show positive metallicity slopes. In some galaxies, spiral arms show up with a slightly lower metallicity than the disk at a given radius. This is likely an artifact arising from the SSP fitting, and the imperfect masking of the youngest regions (see Appendix E).

Fig. 7. Luminosity-weighted stellar [Z/H] radial profiles for the galaxies in our sample. Different colors indicate the radial profile measured across different environments, as indicated in the legend of the bottom-right panel. The black line shows the radial profile measured for the entire FoV (i.e., all environments together). The galactocentric distance is measured in units of R25, in order to measure radial distance homogeneously across our sample. The value of R25 (kpc) of each galaxy is indicated in each corresponding panel. The solid gray line shows the best-fit gradient for each galaxy.

Fig. 8. Mass-weighted stellar [Z/H] radial profiles for the galaxies in our sample. Different colors indicate the radial profile measured across different environments, as indicated in the legend of the bottom-right panel. The black line shows the radial profile measured for the entire FoV (i.e., all environments together). The galactocentric distance is measured in units of R25, in order to measure radial distance homogeneously across our sample. The value of R25 (kpc) of each galaxy is indicated in each corresponding panel. The solid gray line shows the best-fit gradient for each galaxy.

Zhuang et al. (2019) also show that these metallicity gradients are a consequence of the local [Z/H]–stellar mass surface density (Σ_*) relation, and that they are expected to arise primarily due to the insitu local SFH. Indeed, Neumann et al. (2021) report a tight correlation between [Z/H] and Σ_*, with a higher scatter toward larger galactocentric radii and among lower-mass galaxies, implying that in the outer regions of galaxies (as well as within less massive galaxies), additional mechanisms such as gas accretion or outflows might become relevant in determining the chemical enrichment.

Figures 9 and 10 show explicitly how the slope and intercept of the radial metallicity profiles (LumW and MassW, respectively) depend on total stellar mass and total sSFR of each galaxy. The slope and intercept of the LumW- and MassW–[Z/H] radial profiles are measured in the same fashion as explained above for stellar age profiles. We find a clear trend between the normalization of the stellar metallicity gradients and total stellar mass of each galaxy. This trend is expected from the well studied mass-metallicity relation of galaxies (Foster et al. 2012; Ma et al. 2016), which is thought to arise because more massive galaxies are more likely to retain a larger fraction of their metals due to their stronger gravitational potential, as opposed to low-mass galaxies, in which stellar feedback is able to remove metals from the galactic disk more efficiently (see, e.g., Ma et al. 2016).

Fig. 9. Luminosity-weighted stellar [Z/H] radial profiles for the galaxies in our sample, considering all environments. Top left panel shows the measured LumW [Z/H] radial profiles, colored by total galaxy stellar mass. Top right panel shows the profiles normalized by their value at 0.3 R25, to better compare their shapes. The bottom panels show the slope (ΔLumW [Z/H]) and intercept of the LumW radial [Z/H] profile (indicated with a solid gray line in each panel of Fig. 7) as a function of the total stellar mass of galaxies, color coded by total specific star formation rate (sSFR).

Fig. 10. Mass-weighted stellar [Z/H] radial profiles for the galaxies in our sample, considering all environments. Top left panel shows the measured MassW [Z/H] radial profiles, colored by total galaxy stellar mass. Top right panel shows the profiles normalized by their value at 0.3 R25, to better compare their shapes. The bottom panels show the slope (ΔMassW [Z/H]) and intercept of the MassW radial [Z/H] profile (indicated with a solid gray line in each panel of Fig. 8) as a function of the total stellar mass of galaxies, color coded by total specific star formation rate (sSFR).

Finally, there are some peculiarities in the metallicity gradients that are worth mentioning. NGC 3351 exhibits an “S” shape LumW metallicity gradient, which originates from recent star formation occurring in the central region and in the outer ring, leading to enhanced metallicities at these radii. NGC 4535 also has strong star formation at the ends of its bar and in the spiral arms (hence, lower LumW ages), leading to enhanced metallicities in these environments. NGC 7496 shows a sharp radial decrease in metallicity in its bar profile, which is likely an artifact from the spectral fitting (see Appendix E). The LumW metallicity gradient of NGC 1433 shows a strange discontinuous pattern that is likely an artifact due to subtle systematic differences between the MUSE pointings of the mosaic. In Appendix A.2 we explore the correlation between the slope of the stellar metallicity radial gradient and a set of global properties of galaxies. Tables 4 and 5 show the same information as Table 2, for the LumW [Z/H] and MassW [Z/H] distributions, respectively.

3.2. Further insights on the inside-out growth of galaxies

In this section, we explore in more detail what we can learn about the radial structure of the stellar mass assembly of PHANGS galaxies, beyond the mean age radial gradients.

Figure 11 shows the age of the Universe at which each radial bin formed 80% of its current total stellar mass (neglecting the impact of stellar migration, whose contribution in shaping the radial structure of stellar populations is estimated to be of second order; see e.g., Zhuang et al. 2019). A positive trend implies that outer regions assembled their stellar mass later in cosmic time than inner regions. The figure shows that most galaxies follow an overall well-defined positive trend, a direct imprint of the inside-out growth of galaxies. The most remarkable exception is NGC 1385, that shows an almost perfectly flat trend. This, together with its positive age and metallicity radial trends make NGC 1385 a very peculiar system. Another exception is NGC 3627, which also shows a nearly flat MassW age radial profile. NGC 3627 is a member of the interacting group Leo Triplet (Zhang et al. 1993), interactions with nearby galaxies could have triggered episodes of star formation (e.g., Renaud et al. 2019) that might explain its peculiar radial assembly history of stellar mass.

Fig. 11. Radial assembly history of our sample galaxies. The y-axis shows the age of the Universe at which each radial bin assembled 80% of its total current stellar mass. A positive trend in these plots implies that outer regions assembled their stellar mass later in cosmic history than inner regions. The yellow solid lines show the best-fit linear model for each galaxy. The red dashed line show the time at which the full galaxy formed 80% of its total current stellar mass. The gray shaded area represents the uncertainty in our measurement, estimated by propagating the stellar mass uncertainty in the integration of the SFH of each radial bin. The value of R25 (kpc) of each galaxy is indicated in each corresponding panel. The bottom-right panel displays the trends for all galaxies, color-coded by total stellar mass. By construction, the highest possible value for the y-axis is given by the age of our oldest age bin in the SSP fitting (∼15 Gyr, see Sect. 2.2).

Other galaxies, such as NGC 1672 or NGC 1300 show a positive trend across parts of their stellar disk, and then a negative trend in the outer radii. Such a negative trend suggests that in these galaxies, the outer regions have remained relatively quiescent for long periods of time, compared to the inner radii. This does not imply that they lack current star formation, but it means they have not formed a substantial amount of stellar mass in the last ∼5 Gyr. Indeed, Fig. 2 shows that these galaxies do not show higher LumW ages in their outer radii.

Figure 12 shows clearly how the normalization and slope of the of the radial assembly profiles depend on the galaxy stellar mass and specific SFR. We find a trend such that the radial assembly profiles of more massive galaxies (green) are steeper, with values at smaller radii being preferentially located at younger ages (i.e., in the earlier Universe). Less massive galaxies (blue) show flatter gradients, with their inner regions preferentially showing formation times at higher cosmic times.

Fig. 12. Radial assembly history profiles for measured for our sample galaxies, considering all environments. Top left panel shows the measured profiles, colored by total galaxy stellar mass. Top right panel shows the profiles normalized by their value at 0.3 R25, to better compare their shapes. The bottom panels show the slope and intercept of the radial profiles (indicated with a solid yellow line in each panel of Fig. 11) as a function of the total stellar mass of galaxies, color coded by total specific star formation rate (sSFR).

3.3. Searching for evidence of radial mixing in bars

Interaction of stars or gas with nonaxisymmetric structures can lead to a change in angular momentum, ultimately causing their radial migration toward different orbits. Different mechanisms, such as the exchange of angular momentum at the corotation resonance of spiral arms (Sellwood & Binney 2002), or induced by the bar-spiral arms resonance overlap (Minchev et al. 2011) are thought to drive the migration of stars. Radial migration of stars is expected to produce a flattening of radial gradients of the stellar population properties, especially for older stellar populations (e.g., Di Matteo et al. 2013).

Simulations show that bars provide very efficient mechanisms for the radial redistribution of gas and angular momentum (Brunetti et al. 2011; Sormani et al. 2015; Spinoso et al. 2017; Hatchfield et al. 2021); hence, radial gradients of stellar population properties along the bar are expected to be flatter than along the galactic disk due to orbital mixing and stellar radial migration (although the latter is mostly visible from the corotation radius outward, e.g., Friedli et al. 1994; Di Matteo et al. 2013).

A number of studies have tested this prediction using nearby galaxies (Sánchez-Blázquez et al. 2014a; Seidel et al. 2016; Fraser-McKelvie et al. 2019; Neumann et al. 2020), obtaining a variety of results. Sánchez-Blázquez et al. (2014a) see no significant differences between the radial profiles of stellar age and metallicity of barred and unbarred galaxies from the CALIFA sample. Seidel et al. (2016) use data from the BaLROG project (Seidel et al. 2015), and find that stellar metallicity gradients along the bar major axis are considerably flatter than along its minor axis. Fraser-McKelvie et al. (2019) measure flatter stellar age and metallicity gradients along bars, compared to disks in galaxies from the MaNGA sample. Neumann et al. (2020) used data from the TIMER project (Gadotti et al. 2019), and find flatter MassW age and metallicity radial profiles along bars, compared to the disk, but differences are only mild with significant scatter, and no clear trend with bar strength.

Here, we measure MassW age and metallicity radial gradients along pseudo-slits placed across the bar major axis and perpendicular to it. Figure 13 shows the position of the pseudo-slits overplotted on the stellar mass surface density map of NGC 1433. We have measured radial gradients along each of these two pseudo-slits (excluding centers). The radial profiles (and the best-fit gradients, measured using an OLS fitting routine) of MassW age and metallicity measured along each one of these pseudo-slits for NGC 1433 are shown in Fig. 14. We excluded the innermost and outermost radial bins when measuring the radial gradients, in order to avoid potential deviations from the radial trend driven by these extreme values. The gradients along the bar major axis, and their perpendicular counterparts, are calculated as the mean gradient of the two possible directions along each axis from the center (i.e., mean gradient between dark-green (black) and pale-green (gray) pseudo-slits in Fig. 13). The difference between both measurements span a wide range of values, from ∼0, to more than 1 dex, and it is interpreted as the uncertainty in the galaxy gradient (along a given axis).

Fig. 13. Representation of the location of the pseudo-slits along the bar of NGC 1433 (pale and dark green rectangles), and perpendicular to the bar axis (gray and black rectangles). The position of the bar is indicated by the light-green ellipses. The background shows the stellar mass surface density map of the galaxy.

Fig. 14. Mass-weighted age (left) and [Z/H] (right) radial profiles along the bar major axis (dark green), and along its perpendicular direction (black), for NGC 1433. The solid and dashed dark green and black lines represent the gradient measured toward each possible direction from the center. The pale green dashed lines show the best-fitting linear model to the bar major axis gradient (dark green lines), and the silver dashed lines show the best-fitting to the gradients measured in the direction perpendicular to the bar.

Figures 15 and 16 show the slopes of the MassW age and metallicity gradients measured along the bar axis (green), and its perpendicular direction (black), for the 14 barred galaxies in our sample, as a function of bar length (in units of kpc and R₂₅), total stellar mass, and offset from the main sequence of galaxies (ΔMS). NGC 2835 was excluded because the small size of its bar makes the measurement of the gradients unreliable.

Fig. 15. Mass-weighted age gradients measured along bar major axis (green pentagons), and along its perpendicular direction (black squares), as a function of the bar length (in units kpc and R25), total stellar mass, and offset from the global main sequence of galaxies (ΔMS), for 14 barred galaxies in our sample. The red line is calculated following Eq. (4). It is defined to be positive in galaxies that show a flatter gradient along their bar major axes, compared to the bar perpendicular direction, and negative values otherwise.

Fig. 16. Mass-weighted [Z/H] gradients measured along bar major axis (green pentagons), and along its perpendicular direction (black squares), as a function of the bar length (in units kpc and R25), total stellar mass, and offset from the global main sequence of galaxies (ΔMS), for 14 barred galaxies in our sample. The red line is calculated following Eq. (5). It is defined to be positive in galaxies that show a flatter gradient along their bar major axis, compared to the bar perpendicular direction, and negative values otherwise.

The red line is defined as the difference of the absolute value of the gradient measured along the bar and its perpendicular direction, that is,

$$\begin{array}{cc}\hfill \mathrm{\Delta }\mathrm{MassW}\log {\mathrm{Age}}_{\mathrm{disk}-\mathrm{bar}}=& |\mathrm{\Delta }\mathrm{MassW}{\log \mathrm{Age}}_{\mathrm{disk}}|\hfill \\ \hfill & -|\mathrm{\Delta }\mathrm{MassW}{\log \mathrm{Age}}_{\mathrm{bar}}|\hfill \end{array}$$(4)

and

$$\begin{array}{cc}\hfill \mathrm{\Delta }\mathrm{MassW}{[\mathrm{Z}/\mathrm{H}]}_{\mathrm{disk}-\mathrm{bar}}=& |\mathrm{\Delta }\mathrm{MassW}{[\mathrm{Z}/\mathrm{H}]}_{\mathrm{disk}}|\hfill \\ \hfill & -|\mathrm{\Delta }\mathrm{MassW}{[\mathrm{Z}/\mathrm{H}]}_{\mathrm{bar}}|\hfill \end{array}$$(5)

implying that positive values correspond to a flatter profile along the bar, compared to the perpendicular direction, and negative values otherwise.

For ΔMassW age, we do not see evidence of flatter profiles along the bar. On the contrary, we find that, on average, MassW age gradients are steeper along the bar (−0.59 ± 0.33), than along its perpendicular direction (−0.38 ± 0.50), and less than half of the galaxies (6 out of 14) show a flatter MassW age gradient along the bar. We do not find any trend between flatter profiles, and bar length, total stellar mass, or ΔMS.

On the other hand, we observe flatter MassW [Z/H] gradients on average along the bar. We calculate a mean gradient of 0.08 ± 0.31 along bars, and −0.30 ± 0.48 in the perpendicular direction. In terms of individual galaxies, 10 out of the 14 barred galaxies in our sample show a flatter stellar metallicity profile along the bar, without any significant trend with bar length, total stellar mass, or ΔMS.

Finding more negative age profiles along the bar is not surprising. Seidel et al. (2016) also find that age profiles are more negative along the bar semi-major axis, than along its semi-minor axis. This is consistent with chemodynamic simulations (Wozniak 2007) that predict that younger stellar populations will be confined to more elongated orbits, populating the edges of the bar. Neumann et al. (2020) report a similar spatial distribution of stellar populations within the bar of nearby galaxies, with old (> 8 Gyr) stars shaping the inner and rounder part of the bar, intermediate age (2−6 Gyr) stars trapped on more elongated orbits, and an accumulation of young (< 2 Gyr) stars at the bar edges.

In the case of metallicity, we find considerably flatter profiles along bars, compared to their perpendicular direction. This result is consistent with simulations that investigate the bar-driven secular evolution of galaxies, and predict flatter metallicity gradients in barred galaxies (Di Matteo et al. 2013), and that have been observationally confirmed (e.g., Seidel et al. 2016).

In conclusion, we find steeper age gradients, but flatter metallicity gradients along bars, compared to their semi-minor axis. The natural explanation for flatter [Z/H] gradients is orbital mixing of gas along the bar, as gas is more susceptible to nonaxisymmetric structure than stars. Hence, gas from which stars form is homogenized along the bar. Another explanation for flatter [Z/H] gradients along bars could be a flatter surface mass density profile along bars as compared to the disk, such that the local environment in the bar could be more efficient in recycling the gas. However, we acknowledge that due to the limited size of our sample, our results are also consistent with similar age and metallicity gradients along bars, and perpendicular to them.

3.4. Kinematic imprints of different stellar populations across the galactic disk

3.4.1. Previous measurements of the age–velocity dispersion relation

The dynamical evolution of galactic disks is encoded in present-day observations of the kinematic properties of its different stellar populations. It has long been known that older populations of solar neighborhood stars have larger velocity dispersions, compared to younger stars (Roman 1954; Wielen 1977; Freeman 1987). The trend of increasing velocity dispersions with age is known as the age–velocity relation (AVR). Binney et al. (2000) used HIPPARCOS (Perryman et al. 1997) data to quantify the rate at which the velocity dispersion of a coeval group of stars increases with time, and found that it scales with age as τ^0.33, with τ in Gyr. Recent works (Yu & Liu 2018; Mackereth et al. 2019; Tarricq et al. 2021) have performed similar measurements for open clusters and field stars in our Galaxy, using data from Gaia DR2 (Gaia Collaboration 2018).

An alternative explanation to the AVR is that the interstellar medium at high redshift was more turbulent, and thus, older stars retain larger velocity dispersion to present day (Stott et al. 2016). Along this line, Leaman et al. (2017) find that observations of galaxies in the Local Group are consistent with a model in which stars are born with a velocity dispersion close to that of the gas from which they formed, and are then dynamically heated with an efficiency that depends on the galaxy mass, among other factors (see, e.g., Pinna et al. 2018).

Different mechanisms can contribute to the progressive dynamical heating of stars that lead to the AVR, such as interaction with giant molecular clouds (Spitzer & Schwarzschild 1953), interaction with spiral arms (Barbanis & Woltjer 1967; Sellwood & Carlberg 1984; Mackereth et al. 2019), bars (Grand et al. 2016), or even external sources of perturbation (Grand et al. 2016; Pinna et al. 2018). However, measuring the relative importance of these different mechanisms is not straightforward. Studies of the AVR are limited to the Local Group, where individual stars can be resolved. This limitation means that the parameter space of factors that could contribute to the diffusion of stars (e.g., galactic structure) is poorly covered.

3.4.2. Age–velocity dispersion relation in our data

Here, we present an exploration of the AVR measured in nearby galaxies from the PHANGS sample, using the LumW age and velocity dispersion maps derived as explained in Sect. 2. Although the spatial resolution of ∼100 pc is far from being sufficient to resolve individual stars, it is possible to resolve young star-forming regions, and more generally, regions dominated by stellar populations of distinguishable ages. Furthermore, since at this spatial scale the structural components of the galaxy (namely bars, spiral arms, centers, rings, and disks) are clearly resolved, we can also search for changes in dynamical heating of stars as a function of local environment.

Figure 17 shows the AVR (using LumW age) as 2D histograms for three galaxies (NGC 1365, NGC 1566 and NGC 3627) of our sample. These three galaxies are good examples of the general trend that we see in our sample; spiral arms and disk share similar stellar velocity dispersion (σ_*) values, with the former located preferentially at lower ages, and centers and bars dominated by older stellar populations, with higher σ_* values (although some galaxies, such as NGC 1365, can also host young stellar populations in their centers). These positive correlations could, at first glance, be naively interpreted as the AVR. However, in order to properly interpret these trends, it is important to take into account the overall negative σ_* radial gradient set by the underlying stellar mass distribution, following the virial theorem (see, e.g., Bittner et al. 2020). The overall negative σ_* radial profile is clearly visible in the bottom left panel of Fig. 1. The σ_* trend, in combination with the negative LumW age profiles discussed in Sect. 3.1, result in the positive trends seen in Fig. 17.

Fig. 17. 2D histogram that shows the distribution of stellar velocity dispersion (σ*) and log age (luminosity-weighted), for NGC 1365 (left), NGC 1566 (center), and NGC 3627 (right). The color of each [LumW log Age, σ*] bin scales with the number of pixels within that bin. The contours show the one- and two-sigma limits of the distribution of each individual galactic environment, following the color code indicated in the top-right corner of the right of the panel.

Therefore, in order to explore the trend between σ_* and LumW age in a meaningful way, we investigate the radial gradients of σ_* for three different LumW age bins; young (LumW age < 100 Myr), intermediate (100 Myr < LumW age < 600 Myr), and old (LumW age > 3 Gyr). This choice is motivated by the results from Tarricq et al. (2021), who find that the increase in σ_* occurs more rapidly within the first Gyr after stars form, following a much slower increase from this point. To do this measurement, we take the square root of the azimuthally averaged square of the σ_* values in a given radial bin, considering separately spaxels with LumW ages in each one of the three age bins defined.

Figure 18 shows the σ_* radial profiles for the same three galaxies, separately for the three age bins defined above, and across different galactic environments. For reference, we also display the molecular gas σ_mol profile, measured from the CO(2−1) PHANGS-ALMA data.

Fig. 18. Radial profile of stellar velocity dispersion (σ*) for NGC 1365 (left), NGC 1566 (center), and NGC 3627 (right), plotted as a function of galactic structure (colors) and LumW age of the underlying stellar populations. Stellar populations younger than 100 Myr are indicated by diamonds (dotted line). Stellar populations of ages in the range of 100 Myr < log LumW age < 600 Myr are marked with triangles (dashed line). Stars older than 3 Gyr are represented by squares (dashdotted line). Circles (solid line) show the radial profile considering all stellar populations. The gray line shows the radial velocity dispersion gradient measured for the molecular gas, from the PHANGS-ALMA CO(2−1) data.

The figures show a clear age gradient in the σ_* value at a given galactocentric distance, with younger populations having lower velocity dispersions than older ones. This effect is very clear across disk and spiral arms in general, but the age gradient is generally nonexistent in bars. Out of the 14 barred galaxies in our sample, the only noticeable exceptions are NGC 1365 and NGC 5068, with a very clear age gradient present also in the bar. It is interesting to note that these two galaxies represent the high- and low-stellar mass extremes of our sample (log M_*/M_⊙ = 10.99 and 9.40, respectively). Therefore, this similarity is likely not connected to their total stellar mass. On the other hand, we do not observe significant differences between spiral arms and disks, at any given radius. Finally, the molecular gas from which young stars form shows a nearly flat radial σ_mol profile, with very low (≲20 km s⁻¹) velocity dispersion values, and a rapid increase toward the inner radii. Overall, these findings suggest that (i) the dynamical heating of stellar populations occurs at any given radius, increasing the velocity dispersion from ≲20 km s⁻¹, up to ≳50 km s⁻¹ on timescales of hundreds of Myr; (ii) assuming that molecular gas traces the velocity dispersion of stars at birth (Leaman et al. 2017), the increase of velocity dispersion takes place in shorter timescales within the bar, as we do not see clear difference in the σ_* of the different populations in bars for most barred galaxies. However, we acknowledge that this result could be driven by the lack of young stellar populations in the bar hindering the measurement of a clear age gradient. (iii) Spiral arms do not play a primary role in the heating of stellar populations, as we observe this phenomenon across the entire disks, including galaxies that do not contain distinct spiral arms. Alternatively, spiral arms could heat stellar populations in the radial and azimuthal directions, contributing only mildly to the line-of-sight velocity dispersion on galaxies with relatively low inclinations (as those in our sample).

To illustrate these trends across our full sample, Fig. 19 shows the σ_* values measured for stellar populations within the three different age bins defined at a fixed radius (0.25−0.30 R₂₅) for the disk environment. We choose this specific range of radii because it is covered by the disk of most of our sample galaxies. For those few galaxies that do not show these three stellar ages present within the previously defined radius range (NGC 1087, NGC 1433, NGC 1512, NGC 3351, NGC 4303), we opted by choose a different radial bin of the same length (0.05 R₂₅), located at the minimum possible radius in which the three populations are present in the disk environment. The figure reveals a positive trend between σ_* and LumW age of the underlying stellar populations that exists for most of the galaxies in our sample. However, the differences in σ_* between different age bins are often comparable or even smaller than the standard deviation of the distribution of σ_* values within a given age bin (shown as the error bars in Fig. 19). Although this could, in principle, imply that the measured σ_* differences are not statistically significant, the fact that we find the positive slope of σ_* with LumW age to be nearly ubiquitous in our sample and that error bars within a given age bin are intrinsically overestimated due to the underlying gradient (introducing an additional spread of ∼5 − 10 km s⁻¹ in a 0.05 R₂₅ bin), points to the fact that these differences in σ_* across stellar populations of different ages are indeed a real feature in our data. Furthermore, while it is true that the light (and thus, measurement of σ_*) in the pixels with young LumW ages is dominated by the spectra of young stars, these pixels also likely have a subdominant (in light) old component with higher velocity dispersion, that would bias the σ_* measurements toward higher values, reducing the measured differences. This effect would become stronger in pixels with relatively higher light contributions from old stars.

Fig. 19. Mean stellar velocity dispersion (σ*) measured at a galactocentric radius of 0.25−0.30 R25 in the disk environment of all our sample galaxies for stellar populations with luminosity-weighted ages within the three different age bins defined in the main text. We note that for visualization purposes, the x-axis position of each data point is located at the edge of the age bin it represents. In each panel, the leftmost triangle represents the mean (area-weighted) velocity dispersion of all populations younger than 100 Myr. The central square marks the mean velocity dispersion of populations with 100 Myr < LumW age < 300 Myr, and the rightmost triangle shows the same quantity for populations older than 3 Gyr. The errorbar indicates the standard deviation of the σ* values of the pixels within a given age bin. The bottom right panel shows the trend for all galaxies, colored by total stellar mass.

It is also worth noting that the finite spectral resolution of VLT/MUSE limits our ability to measure low (< 40 km s⁻¹) σ_* values depending on the S/N in each Voronoi bin. As described in Sect. 2.2, we bin our data to a S/N of at least 35, though single spaxels can have higher S/N, especially in the high surface brightness innermost (< 0.3 R₂₅) regions of galaxies. This suggest that the results in Fig. 19 are robust, but the low σ_* values in Fig. 17 should be treated with cauthion, and likely explains some of the scatter at low σ_* values in Fig. 18, especially at large radii in the young disk.

Finally, the bottom right panel of Fig. 19 shows that there is not any clear trend between the change of velocity dispersion with stellar age and total stellar mass of the galaxies.

The scenario in which young stellar populations get heated progressively toward higher vertical heights (i.e., increasing σ_*, z), producing a continuous gradient between age and σ_*, z (at a given radius), is also consistent with the results from Bovy et al. (2012), where the authors report that the vertical mass distribution of the disk of our Galaxy is consistent with a single structural component, rather than a combination of a ‘t‘hin” and a “thick” disk. This diffusion would also naturally lead to higher galactic latitudes being populated exclusively by old stars, while galactic latitudes closer to the galactic plane could be populated by a wider variety of stellar populations, a trend that has been indeed observed in our Galaxy (Casagrande et al. 2016).

It is worth noting that because of the generally low inclination of the galaxies in our sample (< 60°), our setting allows us to measure the velocity dispersion in the vertical direction primarily (i.e., σ_*, z). Although different works have shown that the diffusion of stars occurs similarly in all directions (Yu & Liu 2018; Mackereth et al. 2019; Tarricq et al. 2021), there is debate in the literature about the relative strength of the diffusion toward different directions. Yu & Liu (2018), Mackereth et al. (2019) and Sharma et al. (2021) find that the diffusion in the vertical direction is stronger than in the azimuthal and radial directions. More recently, Tarricq et al. (2021) found that the diffusion in the radial and angular directions are stronger than in the vertical direction. In these works, the rate at which the velocity dispersion of young stars increases with time is quantified by the index of the power law (β) between age and σ_*, such that σ_* ∝ τ^β, following Binney et al. (2000). A similar quantification scheme is not applicable to our data. This is because we define age bins significantly broader than the time-resolution of the measured SFH in order to maximize the radial coverage at a given age, and decrease the statistical uncertainty in our measurements. The exact offset between the young, intermediate, and old populations depends on the definition of the bins, the galactocentric distance, local galactic environment, as well as systematic uncertainties associated with our spectral fitting. Nevertheless, our analysis provides useful insights on the impact of galactic structure on the progressive dynamical heating of young stellar populations, and the speed at which it occurs in star-forming galaxies beyond our Local Group.

Finally, these findings suggest it could be worthwhile to consider a refinement of our spectral fitting methodology currently implemented in the PHANGS-MUSE DAP, to include two independent kinematic components rather than just one. Zhang et al. (in prep.) will present an exploration of this dimension in more detail, characterizing the bias in the derived stellar population parameters yielded by the different kinematic conditions of young stellar populations.

3.5. Star formation histories across the galactic disks

In this section, we present how different galactic regions have formed stars through the lifetime of a galaxy. Interpreting the present-day spatially resolved SFHs as the real assembly history of a given galactic region neglects the impact of mechanisms that could cause stars to travel far away from their birthplace. However, the negative age and metallicity profiles, and the general inside-out formation found for these galaxies point to an evolution dominated primarily by insitu star formation, with radial migration being a second-order effect (Ibarra-Medel et al. 2016; Zhuang et al. 2019; Neumann et al. 2021).

Figures 20 and 21 show the time-averaged SFR in five time bins, calculated as the total stellar mass formed within that time window (according to the derived SFH), divided by the length of the time-bin, for NGC 1566, and NGC 3627. There is a number of interesting features worth noting. The most obvious one is the lack of SFR in the second-youngest age bin (25 Myr–100 Myr). The low SFR in this time bin is a consequence of the artifact in the SFHs described in Appendix C, i.e., the dip in the SFH at log age ∼7.7 Myr. Hence, we avoid physically interpreting this feature.

Fig. 20. Time-averaged SFR in five different time bins, calculated from the spatially resolved SFH, for NGC 1566. The panels show how different regions of the galactic disk have assembled their stellar mass in different moments of cosmic history. The blue and green contours mark the position of the spiral arms and the bar, respectively.

Fig. 21. Time-averaged SFR in five different time bins, calculated from the spatially resolved SFH, for NGC 3627. The panels show how different regions of the galactic disk have assembled their stellar mass in different moments of cosmic history. The blue and green contours mark the position of the spiral arms and the bar, respectively.

The central region of NGC 1566 seems to have been quenched (i.e., stopped forming stars) after the oldest age bin, and restarting star formation activity only very recently. A possible scenario could be that bar-driven secular evolution has enabled star formation in the center, by supplying gas to the innermost region. However, as discussed in Sect. 3.2, the lack of nonbarred galaxies in our sample makes it difficult to quantify the role of the bar in the reactivation of star formation in centers. On the other hand, NGC 3627 seems to have been forming stars in its center continuously (this is also reflected in its inner radial bin, in Fig. 11). In our sample, we see a variety of behaviors, such as a drastic decrease of star formation in the inner regions, from early times to present day (e.g., NGC 4254), a continuous star formation activity, or a reactivation of star formation in recent times. These findings are consistent with Querejeta et al. (2021), who also report a wide range of depletion times for centers, possibly due to bars fuelling gas episodically to the center, cycling between star forming and quiescent phases.

Another clear feature is how galactic structure (sp. arms) is evident in the youngest age bins, but it is progressively less prominent toward older ages, to become essentially nonexistent for old stars. This could be interpreted as a consequence of the dynamical heating discussed in Sect. 3.4 disrupting the galactic substructure, leading to a smooth distribution. Along this line, Neumann et al. (2020) use magneto-hydrodynamical cosmological simulations from the Auriga project (Grand et al. 2017) and find that while older stars (> 8 Gyr) follow a smooth distribution with a central concentration, younger stars (< 4 Gyr) are strongly coupled to the galactic structure. However, we can not rule out the scenario in which the AVR is driven by a more turbulent interstellar medium in the past, rather than a progressive increase of velocity dispersion with time.

The diffusion of galactic structure can be clearly appreciated in the time-averaged SFR maps of NGC 1566, particularly for its spiral arms. While in the youngest bin, stars are located preferentially along the spiral arms (marked in blue), in the third age bin (150−700 Myr), the spiral pattern is still clearly visible, although due to the diffusion of stars, the spiral arms are considerably broader. It is worth mentioning that while in Sect. 3.4 we presented evidence for dynamical heating in the vertical direction (i.e., increase of σ_*, z with time), it is the dynamical heating in the perpendicular directions (i.e., radial and angular directions) that is responsible of the broader spiral arms we see in the third age bin, compared to the youngest bin. We also note that the diffusion of stars occurs everywhere across the galactic disk, not only across the spiral arms. However, the existence of a clear spiral pattern provides the opportunity of observing the process of stellar diffusion, and the subsequent “washing out” of galactic structure across discrete age bins. A quantification of this diffusion from these SFR maps is challenging, due to the generally low spiral arm-disk contrast in the third age bin map, and due to artifacts in the maps (pointing-to-pointing jumps, or “holes” in the maps). A sophisticated approach to quantify the time-scales of this broadening is thus, beyond the scope of this paper.

In order to visually identify this diffusion phenomenon in the time-averaged SFR maps, we require the galaxy to have clear spiral arms, and that the spiral arms had been persisting in time, not formed recently by, for instance, a local instability sheared by differential rotation (Goldreich & Lynden-Bell 1965). Thus, it is not surprising that this effect is more prominent in some galaxies than in other. In the particular case of NGC 3627, the spiral arms-disk contrast is very low in the third age bin, which could be indicating a more recent origin of its spiral arms, possibly connected to a disruption driven by interactions with other galaxies from the Leo triplet (NGC 3623, and NGC 3628). Thus, this approach can provide a powerful tool to study the nature and prevalence of spiral arms across time.

4. Summary

In this paper, we have presented a detailed analysis of the stellar populations’ properties of PHANGS-MUSE galaxies, including investigating the radial structure of age and metallicity, radial assembly history, dynamical heating of young stellar populations, and time-averaged SFRs across the galactic disks of 19 massive star-forming galaxies. Our main conclusions are the following:

• We find generally negative mass-weighted (MassW) and luminosity-weighted (LumW) age profiles. Low-mass galaxies NGC 1385 and NGC 1087 are exceptions, showing nearly flat profiles. We find that the slope of the radial gradients correlates reasonably well with global quantities, such that galaxies with higher offset from the main sequence or higher molecular gas content show flatter LumW age profiles. Spiral arms host younger (lower LumW ages) populations than other galactic environments at any given radius.

• Galaxies have mostly negative metallicity profiles and show some degree of correlation of the slope of the metallicity gradient with total stellar mass and SFR, such that more massive (and star-forming) galaxies exhibit steeper metallicity gradients. The intercept of the metallicity radial profiles also correlates with total mass, such that higher mass galaxies show overall higher metallicites.

• We measure a generally well-defined radial assembly history, in which inner regions assembled their stellar mass earlier than outer regions, pointing to an “inside-out” growth. This is consistent with previous findings in the literature, where late-type galaxies have been found to show a clear inside-out growth mode.

• We find flatter metallicity gradients along galactic bars, compared to their minor axis. This is consistent with the expectation from bar-driven secular evolution inducing radial migration within bars. However, we cannot robustly establish that metallicity gradients are flatter along bars due to our small sample size.

• We investigated the progressive increase of velocity dispersion of young stellar populations, and we find a clear gradient in velocity dispersion with age, such that at any given radius, younger stellar populations (< 100 Myr) have in general lower velocity dispersion than intermediate age [100 Myr, 600 Myr] and old (> 3 Gyr) populations. This gradient is can be interpreted as an imprint of the mechanisms dynamically heating young stellar populations, such as interactions with molecular clouds or nonaxisymmetric galactic substructure. We find this process to occur regardless of the presence of spiral arms. The age gradient in velocity dispersion is not present in bars, and we interpret this as dynamical heating being more efficient in the bar region.

• We have used the spatially resolved SFHs to reconstruct time-averaged SFR maps and study how different regions of galaxies assembled their mass. We find a wide variety of features, such as centers continuously forming stars, or quenched early in the galaxy lifetime. Younger populations are coupled to galactic structural components, but older stars follow a much more homogeneous and centrally concentrated distribution. Finally, in some galaxies we see a clear increase in the width of spiral arms with increasing age of the stellar populations probed, and interpret this as a sign of the diffusion of young stellar populations in the radial and angular directions.

Overall, the derived SFHs tell us a story in which galaxies assembled their stellar mass in an inside-out mode, leading to negative age and metallicity profiles, and where more massive galaxies assembled their stellar mass earlier in cosmic history than less massive galaxies, consistent with the expectation from previous studies. Furthermore, thanks to the high spatial resolution achieved by PHANGS-MUSE, for the first time, we have been able to investigate the properties of stellar populations across multiple galactic environments separately, covering a significant part of the stellar disk of nearby star-forming galaxies, being also able the measure distinct kinematic properties for the different stellar populations present in our sample galaxies.

Acknowledgments

We thank the anonymous referee for their insightful comments which helped to improve the paper. This work was carried out as part of the PHANGS collaboration. Based on observations collected at the European Southern Observatory under ESO programmes 094.C-0623 (PI: Kreckel), 095.C-0473, 098.C-0484 (PI: Blanc), 1100.B-0651 (PHANGS-MUSE; PI: Schinnerer), as well as 094.B-0321 (MAGNUM; PI: Marconi), 099.B-0242, 0100.B-0116, 098.B-0551 (MAD; PI: Carollo) and 097.B-0640 (TIMER; PI: Gadotti). This paper also makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.01161.S, ADS/JAO.ALMA#2015.1.00925.S, ADS/JAO.ALMA#2015.1.00956.S and ADS/JAO.ALMA#2017.1.00886.L. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). I.P. acknowledges funding by the European Research Council through ERC-AdG SPECMAP-CGM, GA 101020943. T.G.W. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 694343). J.M.D.K. gratefully acknowledges funding from the ERC under the European Union’s Horizon 2020 research and innovation programme via the ERC Starting Grant MUSTANG (grant agreement number 714907). COOL Research DAO is a Decentralized Autonomous Organization supporting research in astrophysics aimed at uncovering our cosmic origins. R.S.K. and S.C.O.G. acknowledge support from DFG via the collaborative research center “The Milky Way System” (SFB 881; project ID 138713538; sub-projects B1, B2 and B8), from the Heidelberg cluster of excellence EXC 2181“STRUCTURES” (project ID 390900948), funded by the German Excellence Strategy, from ERC via the Synergy Grant “ECOGAL” (grant 855130), and from the German Ministry for Economic Affairs and Climate Action for funding in project “MAINN” (funding ID 50OO2206). K.K. gratefully acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the form of an Emmy Noether Research Group (grant number KR4598/2-1, PI: Kreckel). P.S.B. acknowledges support through the project grants PID2019-107427GB-C31 from the Spanish Ministry of Science, Innovation and Universities (MCIU). F.B. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No.726384/Empire). M.Q. acknowledges support from the Spanish grant PID2019-106027GA-C44, funded by MCIN/AEI/10.13039/501100011033. J.N. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 694343) and the Science and Technology Facilities for support through the Consolidated Grant Cosmology and Astrophysics at Portsmouth, ST/S000550/1.

References

Appendix A: Slope of radial gradients of stellar population properties versus global galactic parameters

A.1. Stellar age

Here we investigate if the gradients measured for the LumW log Age radial profiles correlate with global galaxy properties. Figure A.1 shows the slope of the LumW-age gradients, as a function of several global galaxy parameters: total stellar mass (log M_*), total SFR (log SFR), offset from the global main sequence of galaxies as defined in Leroy et al. (2019), (ΔMS), total molecular gas mass-to-stellar mass ratio (M_mol/M_*), total atomic gas mass-to-stellar mass ratio (M_HI), total gas (molecular + atomic)-to-stellar mass ratio M_gas/M_*, and molecular-to-atomic gas mass ratio (M_mol/M_HI), with their corresponding Pearson correlation coefficient (ρ). We also include inclination and distance, to ensure that the different gradients are not related to these nonintrinsic properties. The uncertainty in the Pearson coefficient has been calculated by bootstrapping; for each correlation, we repeatedly chose 100 subsamples of datapoints (allowing for repetitions), each subsample with the same size than the original sample, and calculated their ρ coefficient. We adopt the standard deviation of the set of 100 coefficients as the uncertainty for ρ.

Fig. A.1. Slope of the LumW-age gradient, as a function of (from left to right, and from top to bottom) total stellar mass (log M*), total SFR (log SFR), offset from the global main sequence of galaxies (ΔMS), total molecular gas mass-to-stellar mass ratio (Mmol/M*), total atomic gas mass-to-stellar mass ratio (MHI), total gas (molecular + atomic)-to-stellar mass ratio Mgas/M*, and molecular-to-atomic gas mass ratio (Mmol/MHI), galaxy inclination and distance. The Pearson correlation coefficient for each parameter is indicated in the bottom-right corner of each panel.

The highest ρ correlations are with M_mol/M_* (ρ = 0.53) and ΔMS (ρ = 0.43) implying that galaxies with higher molecular gas content and a global enhancement of SFR also have a less negative (i.e., flatter) LumW-age radial gradient. This is consistent with the findings from Ellison et al. (2018) and Pessa et al. (2021), where the authors find that a global enhancement of SFR is reflected primarily on the SFR of the inner galactic regions. In this case, since the LumW-age is very sensitive to recent star-formation, a more negative gradient implies an older and less star-forming central region, with respect to the outer radii.

Figure A.2 provides plots of the MassW-age gradient against the same quantities as Fig. A.1. The slopes of the MassW-age gradients do not correlate strongly neither with ΔMS, M_mol/M_*, nor with any of the parameters explored. However, we see a trend in which galaxies with a high total gas content (M_gas/M_*) span a wide range of gradients, while galaxies with a low gas content exhibit exclusively flatter gradients. This trend could indicate that galaxies with a low gas content have consumed their gas supply for fueling star formation, leading to younger populations in their centers, hence, flatter MassW-age profiles. Even though the level of correlation of M_mol/M_* with LumW-age is not particularly high (ρ = 0.37), the trend seems to be opposite to the one with MassW-age. This is not in contradiction with the scenario suggested for the origin of the trend between the MassW-age gradient and M_mol/M_*, since a low total gas content could certainly imply a lower total present-day SFR (and hence, less SFR in the central region), but it does not imply a lack of SFR on longer time-scales, to which MassW-age is sensitive (as opposed to LumW-age, that is strongly biased toward recent SFR).

Fig. A.2. Slope of the MassW-age gradient, as a function of (from left to right, and from top to bottom) total stellar mass (log M*), total SFR (log SFR), offset from the global main sequence of galaxies (ΔMS), total molecular gas mass-to-stellar mass ratio (Mmol/M*), total atomic gas mass-to-stellar mass ratio (MHI), total gas (molecular + atomic)-to-stellar mass ratio Mgas/M*, and molecular-to-atomic gas mass ratio (Mmol/MHI), galaxy inclination and distance. The Pearson correlation coefficient for each parameter is indicated in the bottom-right corner of each panel.

In summary, we find generally negative age radial gradients, with LumW profiles being steeper than the MassW ones. Galaxies with higher molecular gas content and that are experiencing a global enhancement of SFR exhibit flatter LumW age profiles.

A.2. Stellar metallicity

Figures A.3 and A.4 show the LumW- and MassW- stellar metallicity gradients, as a function of the same global properties that we explored for the age gradients. We do not see any obvious trends between metallicity gradients and these global properties. Goddard et al. (2017) report steeper negative metallicity profiles for more massive galaxies. In our sample, we find that the highest ρ values are indeed obtained with total stellar mass and total SFR, however, possibly due to our low-number statistics, we cannot robustly confirm such a trend. In general, it seems that low-mass galaxies (log (M_*/M_⊙)≲10) span a more diverse range of gradients, while higher mass galaxies exhibit stellar metallicity gradients that correlate better with total mass. In summary, we find negative stellar metallicity radial profiles, with more massive galaxies exhibiting overall higher metallicities and steeper gradients.

Fig. A.3. Slope of the LumW-[Z/H] gradient, as a function of (from left to right, and from top to bottom) total stellar mass (log M*), total SFR (log SFR), offset from the global main sequence of galaxies (ΔMS), total molecular gas mass-to-stellar mass ratio (Mmol/M*), total atomic gas mass-to-stellar mass ratio (MHI), total gas (molecular + atomic)-to-stellar mass ratio Mgas/M*, and molecular-to-atomic gas mass ratio (Mmol/MHI), galaxy inclination and distance. The Pearson correlation coefficient for each parameter is indicated in the bottom-right corner of each panel.

Fig. A.4. Slope of the MassW-[Z/H] gradient, as a function of (from left to right, and from top to bottom) total stellar mass (log M*), total SFR (log SFR), offset from the global main sequence of galaxies (ΔMS), total molecular gas mass-to-stellar mass ratio (Mmol/M*), total atomic gas mass-to-stellar mass ratio (MHI), total gas (molecular + atomic)-to-stellar mass ratio Mgas/M*, and molecular-to-atomic gas mass ratio (Mmol/MHI), galaxy inclination and distance. The Pearson correlation coefficient for each parameter is indicated in the bottom-right corner of each panel.

Appendix B: Uncertainties in age, metallicity, and stellar mass surface density

In this section we show the distributions of the random uncertainties for the physical quantities derived through our full spectrum fitting approach, as described in Sec. 2.2, for the full sample. Additionally, we show that performing 20 Monte Carlo iterations to estimate the uncertainties in these quantities for each Voronoi bin is a reasonable compromise in terms of characterizing the dispersion of the distribution of values produced by the Monte Carlo realizations.

Figure B.1 shows the distribution of the relative errors of Σ_*, LumW Age, MassW Age, LumW [Z/H] and MassW [Z/H]. The relative uncertainty in stellar mass surface density is typically below 10%, while the uncertainties in ages and metallicities are typically below 20%. The mass-weighted quantities exhibit a long wing toward larger values, although this is not surprising, as mass-weighted quantities are dominated by fainter stars that are often outshone by younger and brighter stars.

Fig. B.1. Probability density function of the relative uncertainties estimated for stellar mass surface density (Σ*), LumW Age, MassW Age, LumW [Z/H], and MassW [Z/H], for our full sample, that is, all available pixels.

To test how many Monte Carlo realizations are required to appropriately estimate the width of the produced distribution, we have chosen five representative Voronoi bins (one for each galactic environment), and we have computed the dispersion of the results of the Monte Carlo iteration, for different numbers of iterations, from 5 to 100.

Figure B.2 shows the dispersion of the distribution as a function of the number of Monte Carlo iterations, for stellar mass surface density, age, and metallicity, for these five Voronoi bins. In most cases, the width of the dispersion does not change significantly beyond 20 iterations, although in some cases there is some level of variability until ∼50 iterations. Nevertheless, 20 realizations is a reasonable number to estimate the dispersion of the results of the Monte Carlo iterations, and hence, the magnitude of the random uncertainties of our measurements. The figure also shows clearly that the uncertainty in the MassW quantities is significantly larger than that for LumW quantities when young stars are present (e.g., in spiral arms).

Fig. B.2. Dispersion of the resulting posterior distribution of stellar mass velocity dispersions (Σ*), ages and metallicity from the Monte Carlo realizations, as a function of the number of Monte Carlo iterations, for 5 representative Voronoi bins. The color code indicates the galactic environment of each Voronoi bin. Solid lines correspond to mass-weighted quantities, and dashed lines to their light-weighted counterparts. The vertical black line indicates the adopted number (20) of Monte Carlo iterations to estimate the uncertainties in the analyses presented in this paper.

Appendix C: Downsizing in galaxy evolution

Another interesting topic that we can probe within our sample is the downsizing paradigm of galaxy evolution (see, e.g., Cowie et al. 1996; Heavens et al. 2004; Thomas et al. 2005; Pérez-González et al. 2008). This phenomenon describes how more massive galaxies assembled their stellar mass earlier in cosmic history than less massive galaxies.

C.1. Signal of downsizing in the mosaic-integrated spectra of galaxies

Figure C.1 shows, in the top panels, the total stellar mass formed in a given age bin, and in the bottom panels, the average SFR in a given age bin, defined as the total stellar mass formed, divided by the time-separation between each bin, and the previous one, for two galaxies in our sample (NGC 1365 and NGC 1566). In the case of the first bin, the time-separation is calculated with respect to the present time (lookback time = 0).

Fig. C.1. Star formation history (top) and star formation rate history (bottom) derived for NGC 1365 (left) and NGC 1566 (right). The best-fitting linear model for the SFH curve of each galaxy is shown as a red dashed line, and its slope and intercept are shown in the top-left corner of the top panels. The shaded areas mark sudden decreases in the SFR at ages log (age/yr) ∼7.7 and log (age/yr) ∼9.5, found to exist in all galaxies in some degree, deemed as an artifact of the computed SFHs.

The SFHs (and time-averaged SFRs) of most of the galaxies in our sample show a sudden decrease in the formation of stars at ages log (age/yr) ∼7.7 and log (age/yr) ∼9.5 (indicated as the shaded area in Fig. C.1). We acknowledge that this is very likely an artifact of the derivation of the SFH. This decrease occurs everywhere across the galaxy disk, and we have been not able to identify a clear cause for it. Similar systematic problems at these ages have been identified in the literature. Peterken et al. (2020) report a correlation between the weights assigned for two similar age bins, finding evidence that it is related to the “UV upturn” in old stellar populations (Yi 2008). Although we do not report here an identical problem (we do not see a sudden increment in the SFH, as they do), SFH derivations are subject to systematics that depend on the spectral library used (e.g., Walcher et al. 2011; Martins 2021), as well as the code used to derive it (e.g., San Roman et al. 2019). Finding the exact cause of this artifact is beyond the scope of this paper, but we exclude it from the interpretation of our results.

The SFHs derived are well represented by almost perfectly linear power-laws (see top panels of Fig. C.1). This is partially by construction, as we designed our age grid to be log-spaced, and the time-averaged SFR is nearly flat across cosmic time (neglecting the artifacts described earlier). However, deviations from this linearity reflect real features of the mass assembly process of these galaxies. Galaxies that assembled most of their stellar mass earlier in cosmic history show steeper power-laws (i.e., higher ΔM). On the other hand, galaxies that have assembled a substantial amount of their current stellar mass later in cosmic history are characterized by flatter SFHs (i.e., lower ΔM).

Figure C.2 shows the slope of the SFH (ΔM), measured using an OLS fitting routine, as a function of the same global parameters explored in the previous sections. We find clear correlations between ΔM and several of the global properties, including total stellar mass (ρ = 0.61), atomic gas fraction (ρ = −0.56), and total gas fraction (ρ = −0.64).

Fig. C.2. Slope of the SFH of PHANGS galaxies, as a function of (from left to right, and from top to bottom) total stellar mass (log M*), total SFR (log SFR), offset from the global main sequence of galaxies (ΔMS), total molecular gas mass-to-stellar mass ratio (Mmol/M*), total atomic gas mass-to-stellar mass ratio (MHI), total gas (molecular + atomic)-to-stellar mass ratio Mgas/M*, and molecular-to-atomic gas mass ratio (Mmol/MHI), galaxy inclination and distance. The Pearson correlation coefficient for each parameter is indicated in the bottom-right corner of each panel.

Despite the artifacts present in our SFHs, we believe that trends in the assembly history of galaxies are robust, where more massive galaxies, with lower gas content assembled their stellar mass earlier in cosmic history, than galaxies with lower stellar masses, and higher gas content.

This reflects the downsizing effect in galaxy evolution, and is consistent with findings from previous works (e.g., Ibarra-Medel et al. 2016; Peterken et al. 2020). An alternative explanation for this trend is that the ratio of present-day to past SFR is higher for low mass galaxies than for high mass galaxies Thus, we conclude that the derived SFHs tell a consistent story of the stellar mass assembly of the galaxies in our sample, with strong imprints of an inside-out formation, and with higher mass galaxies assembling their stellar mass on average earlier in cosmic history than lower mass galaxies.

C.2. Downsizing in individual galactic environments

Recent works have studied the assembly history of galaxies across their different galactic environments, by performing a photometric decomposition to isolate bulges and disks, and investigating their stellar populations’ properties separately (e.g., Costantin et al. 2021, 2022; Johnston et al. 2022). They find that bulges are generally older than their surrounding disks, pointing to an inside-out growth of galaxies. Costantin et al. (2022) also finds that the downsizing paradigm is present in individual galactic environments, such that older disks and older bulges are also more massive than their younger counterparts, while (Johnston et al. 2022) find that bulges in more massive galaxies (≳10¹⁰ M_⊙) form earlier than those in lower mass galaxies.

Motivated by these results, we explore whether we can recover the downsizing paradigm in individual environments. To do so, we have simplified our environmental label, grouping spiral arms, outer rings, and stellar bars together with the “Disk” component, and inner rings together with the “Center” component. Hereafter, we refer to these combination of morphological structures as combined disk and combined center, respectively. The left and middle panel of Fig. C.3 show the median MassW age of these combined morphological components as a function of total stellar mass, and stellar mass of each combined morphological component, respectively, for all galaxies in our sample. The figure shows that while the downsizing paradigm is clearly preserved in disks, this is not the case for centers. Furthermore, we do not see a correlation between the age of bulges with the total stellar mass either.

Fig. C.3. Downsizing paradigm in separated galactic environments. Left: correlation between the median MassW log Age and the total galaxy stellar mass for the (combined) disk and center environments, for each galaxy in our sample. The total galaxy stellar mass corresponds to the sum of the masses of the (combined) disk and center environments. Middle: median MassW log Age of the (combined) disk and center environments as a function of total stellar mass in each galactic environment. Right: Median MassW log Age of the (combined) disk environment as a function of the median MassW log Age age of the (combined) center environment for our sample galaxies. Galaxies that present a nuclear ring are indicated with an additional orange marker.

Finding such different behavior of centers, compared to the properties of bulges reported in Costantin et al. (2022) and Johnston et al. (2022) is not surprising. There are a number of reasons that explain this apparent disagreement between this work and these seemingly similar studies. On one hand, there are inherent differences in the sample selection. For instance, Johnston et al. (2022) uses a sample of mostly unbarred S0 galaxies, while Costantin et al. (2022) uses galaxies at redshift 0.14 < z ≤ 1, and it is known that the abundance of galaxies with stellar bars steadily decreases toward redshift 1 (Simmons et al. 2014; Melvin et al. 2014). This implies that their sample is considerably less affected by the presence of bars than our sample, composed of mostly barred galaxies (14 out of the 19 PHANGS-MUSE galaxies exhibit bars). This is a key difference, as bars can efficiently transport gas from the outer galactic regions to the innermost regions of galaxies (Lu et al. 2022). Hence, galactic centers in the (Johnston et al. 2022) and Costantin et al. (2022) samples are less likely to present bar-driven secularly built nuclear disks/rings. Moreover, if bars are present, a pure bulge-disk decomposition is most likely too simplified to capture and isolate the properties of the bulge (e.g., Gadotti 2009).

Beyond differences in the sample and in the methodology to separate morphological components, the spatial resolution of the analyses also plays a relevant role. For example, at the spatial resolution of MaNGA (in the case of Johnston et al. 2022), even if a small secularly built central substructure exists, it would probably be unresolved and mixed with the stellar disk in the decomposition. Finally, the definition of bulges is by itself ambiguous. In this work, we define centers as a central excess of light, independently of their exact surface brightness profile or dynamical properties (see Sec. 2.3). This definition differs from that used in other works, which implies a disparity in the comparison.

Finally, the right panel of Fig. C.3 shows the median MassW age of the (combined) disk as a function of the median MassW age of the (combined) center component. We do not find any clear systematic difference between them, contrasting with previous works that find systematically older bulges. However, this is not unexpected considering the fundamental differences between this and the recent studies discussed, described above. In fact, 5 out of the 7 galaxies where the disk is older than the center show clear bars and inner rings, which further highlights how relevant the sample differences are. This result is consistent with the radial profiles shown in Sec. 3.1, where some galaxies exhibit clear evidence of young stellar populations in their innermost regions.

To summarize, a direct comparison of our results with recent studies that perform photometric decompositions to investigate the assembly history of bulges and disks is difficult due to the underlying sample and resolution differences. We recover the downsizing paradigm in disks, but not in centers, and we do not find any evidence of centers being systematically older than disks. These apparent discrepancies with previous works can be explained by the fundamental dissimilitude between these different studies.

Appendix D: Decreasing the angular resolution of our dataset

In this section, we show how our results depend on the angular resolution of our data. To do this, we have convolved the MUSE mosaics to a fixed angular resolution of 15 arcsecs, and derived the stellar populations’ property maps following the same methodology described in Sec. 2.2. Then, we recomputed radial profiles and overall distributions of the stellar population properties for the convolved dataset, across different galactic environments, using the same environmental masks as for our fiducial analysis (see Sec. 2.3). An angular resolution of 15 arcsecs corresponds to physical scales from 0.4 kpc to 1.4 kpc at the distances of our sample galaxies, with a mean value and dispersion of 1.1 ± 0.3 kpc.

Figures D.1 and D.2 are identical to Figs. 2 and 4, but using now the convolved dataset. Figure D.1 shows smoother radial profiles. Differences among environments are in general, more subtle than at higher resolution, with rings becoming essentially impossible to distinguish from the center anymore. On the other hand, spiral arms are overall still distinguishable from the stellar disk. This is not unexpected, as Querejeta et al. (2021) found that spiral arms of galaxies from the PHANGS sample exhibit typical width of ∼1.5 kpc. Thus, we can marginally resolve them at a resolution of ∼1 kpc. Other works at ∼1 kpc resolution have reported differences between arms and their surrounding disk. For instance, Sánchez-Menguiano et al. (2017) find subtle differences in the oxygen abundances of spiral arms with respect to the interarm region using a sample of galaxies from the CALIFA dataset.

Fig. D.1. Luminosity-weighted stellar age radial profiles for the galaxies in our sample, measured after convolving our MUSE mosaics to a fixed angular resolution of 15 arcsecs, and computing the stellar population properties for these convolved cubes. Different colors indicate the radial profile measured across different environments, as indicated in the legend of the bottom-right panel. The black line shows the radial profile measured for the entire FoV (i.e., all environments together). The galactocentric distance is measured in units of R25, in order to measure radial distance homogeneously across our sample. The value of R25 (kpc) of each galaxy is indicated in each corresponding panel. The solid gray line shows the best-fit gradient for each galaxy. The bottom-right panel displays the entire FoV trends for all galaxies, color-coded by total stellar mass.

Fig. D.2. Distribution of the offset of the LumW log Age of the pixels within each environments, with respect to the median galaxy log LumW Age of each galaxy, in bins of 0.04 dex, obtained for the MUSE mosaics convolved to a fixed angular resolution of 15 arcsecs. Different colors indicate the distribution of different environments, as indicated in the legend.

However, beyond environmental differences, the ability to distinguish particularly old or young regions with respect to their surroundings is drastically reduced. This is also illustrated by Fig. D.2. In contrast with the fiducial version of the figure, which shows extended tails that reach very young (<  − 1.0 dex with respect to the median galaxy age) and old (∼1.0 dex with respect to the median galaxy age) ages, the age distributions in the convolved dataset are confined to values much closer to the median galaxy age, mostly within −0.5 dex and +0.5 dex from this value, implying that although environments might be still distinguishable (as in the case of spiral arms), contamination due to the lower resolution reduces drastically the dynamical range of ages we can measure, in addition of making environments harder to differentiate (or totally washing out their stellar population features, as is the case for inner rings).

Appendix E: Improving the fitting of young stellar populations

In Emsellem et al. (2022) is reported the presence of low metallicity values (LumW [Z/H] < −1.3) in a few regions encompassing very young stellar clusters (LumW age < 400 Myr). Such low metallicity values would be inconsistent with an internal and progressive chemical enrichment of the interstellar medium (e.g., Ho et al. 2017). This suggests that the fitting process converges toward a misleading local minimum, the bluest available stellar template, constrained by the youngest age bin (30 Myr) of the implemented template library. This issue has already been reported in several studies (e.g., Carrillo et al. 2020; Bittner et al. 2020), who similarly reported unexpected low metallicities in young regions.

In this section, we perform a broader exploration of the parameter space, aiming to define a better approach to fit the spectra of these problematic regions and reduce the artifacts in their derived stellar population properties. This exploration includes the usage of different libraries of templates (Sec. E.1), the subtraction of the stellar continuum before fitting the observed spectra (Sec. E.2), accounting for nebular emission (Sec. E.3), and varying the metallicity grid of the templates (Sec. E.4).

E.1. Impact of different sets of templates

The stellar library, together with the initial mass function, and stellar evolution models (isochrones) are the basic ingredients to generate the spectra of SSPs. Variations in these three ingredients will lead to differences in the resulting SSPs, that ultimately can produce systematic differences in the measured properties of a galaxy, or galactic region, from its integrated spectrum (see, e.g., Ge et al. 2019). In this section, we explore how our results vary as a function of the set of SSP templates used to fit our data. In addition to the fiducial E-MILES (Vazdekis et al. 2016) templates used in Emsellem et al. (2022), which include 13 log-spaced age bins of [0.03, 0.05, 0.08, 0.15, 0.25, 0.40, 0.60, 1.0, 1.75, 3.0, 5.0, 8.5, 13.5] Gyr, and 6 metallicity bins of [Z/H] = [-1.49, -0.96, -0.35, +0.06, +0.26, +0.4], for these tests we consider the following set of templates:

• Charlot & Bruzual (2007, priv. comm.) – updated version 2016, to which we will refer as “CB07”. These templates are constructed using the empirical stellar library STELIB (Le Borgne et al. 2003), and use the stellar evolution prescription of Marigo & Girardi (2007) to model the TP-AGB evolution of low- and intermediate mass stars. We employ five [Z/H] bins: [-1.7, -0.7, -0.4, 0, 0.4], and 25 age bins ranging from 1 Myr to 20 Gyr, of which 9 are younger than our fiducial age limit of 30 Myr.

• E-MILES + young extension presented in González Delgado et al. (2005, 2014a), to which we will refer as “E-MILES + G-D”. This combined set of templates consists of our fiducial grid, plus six extra age bins, ranging from 1 Myr to 20 Myr, with four [Z/H] bins: [-0.71, -0.4, 0, 0.22]. These younger templates are only available at a spectral resolution of ∼6 Å, lower than our fiducial E-MILES templates (FWHM ∼2.51 Å)

• E-MILES + newly computed young MILES SSP models, presented in Asa’d et al. (2017), to which we will refer as “E-MILES young”. These templates are an extension of the MILES library toward younger ages, available for a Padova isochrone. The combined library consist of 17 age bins, ranging from 6.3 Myr to 15 Gyr, and five [Z/H] bins, ranging from -1.33 to 0.22 for templates older than 60 Myr, and from -1.33 to 0.41 for templates younger than 60 Myr.

Finally, the tests are performed in the central region of NGC 3351, which shows a prominent star-forming ring, and provides a perfect testing sample to evaluate the impact of the different templates/setting on the resultant fitting of young stellar populations. The quantities that we consider to evaluate the improvement of the quality of the solution are LumW age, LumW metallicity, E(B-V), and stellar mass surface density (Σ_*).

Figure E.1 shows the maps of these four quantities, for each one of the set of templates tested here. It is clear that the addition of younger templates to the age grid is not sufficient to overcome the issue of extremely metal-poor young regions, as the low metallicity values in the inner ring of NGC 3351 persist for all SSP templates. To quantify variations, in terms of these four quantities, we consider their distributions within the ring (between the black and red ellipses in the top-right panel of Fig. E.1), and at radii larger than the ring (i.e., disk, outside the red ellipse). Specifically, we compare the 5^th percentile of the ring distribution, with the median of the disk distribution, for each one of these four quantities. The reason for this difference is that young clusters within the ring are young, often “erroneously” characterized by low metallicity, low extinction, and low stellar masses. Therefore, the low-value tail of these ring distributions encodes information about these problematic regions. On the other hand, the distribution beyond the ring is relatively smooth, and its median is a robust metric for characterizing it.

Fig. E.1. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (second row), stellar E(B-V) (third row), and stellar mass surface density (bottom row) maps derived for the central star-forming ring of NGC 3351. Each column shows the maps obtained using the set of templates indicated at the top of each column. The maps obtained using our fiducial set of templates are shown in the left column. The area outside the red ellipse in the top-left panel is used to probe disk stellar population properties, and the area between the red and the black ellipses is used to probe stellar population properties of the young star-forming ring (see text).

Figure E.2 shows the ring and disk distributions for each of these four quantities, for each set of templates used to produce Fig. E.1. We define:

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{X}={\mathrm{X}}_{p5}^{\mathrm{Ring}}-{\mathrm{X}}_{p50}^{\mathrm{Disk}}\end{array}$$(E.1)

Fig. E.2. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (second row), stellar E(B-V) (third row), and stellar mass surface density (bottom row) normalized distributions for the star-forming ring (black) and disk (red), delimitated in the top-left panel of Fig. E.1. Each column shows the distributions obtained using the set of templates indicated in the legends of the top-row panels. The distributions obtained using our fiducial set of templates are shown in the left column.

for X in [log age, [Z/H], E(B-V), log Σ_*], where ${\mathrm{X}}_{p5}^{\mathrm{Ring}}$ stands for the 5^th percentile of the ring distribution for the quantity X, and ${\mathrm{X}}_{p50}^{\mathrm{Disk}}$ stands for the 50^th percentile (median) of the disk distribution for the same quantity.

Hence, in terms of these four quantities, a good solution should yield to: (i) lower (negative) Δlog age, that is, star-forming clusters are expected to be younger than the surrounding disk, (ii) higher (positive or ∼0) Δ[Z/H], or, in other words, young clusters should not be associated with very metal-poor metallicities, compared to the disk, (iii) not strongly negative ΔE(B-V) values, since holes in the E(B-V) map associated with these young regions are unphysical, and finally, (iv) positive (or ∼0) ΔlogΣ_*, due to the expected smooth underlying stellar mass surface density gradient.

In the following, we show the output maps generated with these four set of templates, under different considerations, and at the end of the section, we present the final comparison in terms of Δlog age, Δ[Z/H], ΔE(B-V), and ΔlogΣ_*.

E.2. Is the stellar continuum shape driving the low-metallicity feature?

Low-metallicity SSPs have a continuum bluer than their high-metallicity counterparts, at a fixed stellar age. Thus, a natural feature to explore the origin of this problem is the stellar continuum shape. One possible explanation for the young and metal-poor feature is that the solution is being driven by the continuum shape, as the software might be struggling to find blue enough templates to reproduce the observed spectra of these young regions (i.e., age-Z degeneracy, see, e.g., Worthey 1994). To test if the shape of the stellar continuum is responsible for the metal-poor feature, we have repeated our measurements after subtracting the continuum of both, the templates and the observed spectra of each Voronoi bin. We have used the specutils⁴ python package to calculate the continuum shape for a given spectrum, iteratively fitting a 5^th order polynomial to account for the stellar continuum, removing outlier absorption and emission features in each iteration. This step is repeated until no more pixels are identified as outliers and rejected from the fit.

Figure E.3 shows the maps obtained with the four template sets under study when the stellar continuum is subtracted before performing the fit. The extinction is calculated before subtracting the continuum, therefore this presents no change with respect to the fiducial case. The stellar mass surface density maps can not be computed once the continuum is subtracted, and hence, the last row has been omitted from the figure. We highlight two aspects from Fig. E.3. First, the solution in age and metallicity is much noisier, compared to the fiducial case. This is because the unsupervised continuum subtraction is not perfect, and creates artifacts in the spectra, often associated with absorption features or residual emission lines perturbing the definition of the continuum level. Nevertheless, the ring is still recovered as a young and metal-poor feature for 3 out of the 4 libraries investigated, and we conclude that the shape of the continuum is not the primary driver of the low metallicity solutions for young regions.

Fig. E.3. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (middle row), and stellar E(B-V) (bottom row) maps derived for the central star-forming ring of NGC 3351, once the stellar continuum has been removed from the observed spectra and the templates. Each column shows the maps obtained using the set of templates indicated at the top of each column. The maps obtained using our fiducial set of templates (after removing the stellar continuum) are shown in the left column. We note that the E(B-V) maps are identical to those in Fig. E.2, as the stellar extinction is measured before subtracting the stellar continuum.

E.3. Could nebular emission be the responsible?

High-energy ultraviolet photons emitted from nearby hot (and young) stars can ionize the surrounding interstellar medium, leading to nebular emission. Nebular emission is composed of nebular line emission and nebular continuum. While the former is primarily produced by radiative recombination processes and emission from specific line transitions, the latter is a continuous emission spectrum that consists of free-free (Bremsstrahlung), free-bound (recombination continuum), and two-photon emission (Byler et al. 2017).

The strength of emission from these two components depends on both the ionizing radiation field and the metallicity of the gas. The amount of nebular emission thus varies from galaxy to galaxy, and can evolve with cosmic time. As stated in Sec. 2.2, for our analysis we have masked emission lines, as they make measurements of the underlying absorption features unreliable. However, nebular continuum emission might also have a pronounced impact on the SSP fitting. In this regard, Cardoso et al. (2022) very recently showed that ignoring the presence of nebular continuum emission in the spectra can lead to important biases in the stellar population properties derived through full spectral fitting.

Figure E.4 shows the spectrum of a solar metallicity SSP at an age of 1 Myr, with and without the inclusion of nebular emission. The nebular component is computed with FSPS (Flexible Stellar Population Synthesis; Conroy et al. 2009; Conroy & Gunn 2010), a python package that generates spectra and photometric predictions for arbitrary stellar populations, using Cloudy (Ferland et al. 2013; Byler et al. 2017; Byler 2018) to calculate the emission produced by ionized gas.

Fig. E.4. Example of the impact of nebular continuum emission on the SED of a 1 Myr old solar metallicity SSP. The sudden decrease of nebular emission (green line) at ∼3646 Å and ∼8207 Å correspond to the Balmer and Paschen breaks, respectively. The vertical red dashed lines enclose the wavelength range considered for the full spectral fitting [4850,7000] Å.

In order to test the influence of nebular continuum emission in our fitting, we first choose two representative regions within the young ring in the center of NGC 3351, one of them is found to be extremely metal-poor, while the other is relatively metal-rich in our fiducial fitting. For each region, we extract a high-SN spectrum, and we use these spectra for further tests. We will refer to these two young regions as MR (metal-rich) and MP (metal-poor). Figure E.5 shows these regions, and the stellar metallicity map determined with E-MILES templates as background. The MR and MP regions are indicated with an orange circle and a blue ellipse, respectively.

Fig. E.5. Two very young regions within the central star-forming ring of NGC 3351 with different inferred metallicity values. The orange circle encloses the metal-rich (MR) region, and the blue ellipse encloses the metal-poor (MP) region. The background map shows the luminosity-weighted stellar metallicity in the central region.

We use FADO (Fitting Analysis using Differential evolution Optimization; Gomes & Papaderos 2017), a SED-fitting software (similar to pPXF) that aims for recovering the best-fitting SFH, while consistently accounting for the observed nebular emission characteristics of a star-forming (SF) galaxy, to fit the spectra of these two regions, and quantify the contribution expected from nebular emission.

Figure E.6 shows that both spectra contain some degree of nebular continuum emission contribution. This is expected, as both of them correspond to young SF regions. However, region MP hosts a much larger nebular continuum contribution than MR, while also being much younger. Moreover, the metal-poor peaks in Fig. E.5 spatially correlate well with peaks in the Hα flux surface density map. This suggests that not accounting properly for nebular continuum contribution in the spectra of the youngest regions could be responsible for the low metallicity obtained for the MP spectrum. Indeed, in the wavelength range of interest, the nebular continuum can be represented, to zero^th order, as a flat spectrum (see Fig. E.4). This boost of the continuum yields weaker absorption features (i.e., shallower relative to the continuum level). This effect can be further quantified by measuring different line-strength indices for different contributions of nebular continuum emission.

Fig. E.6. Best SED-fitting obtained with FADO for the MR and MP regions, defined in Fig. E.5. The orange and blue lines represent the MR and MP observed spectra, respectively. The black line show the best-fitting model, and the gray and green curves show separately the stellar and nebular continuum contribution, respectively. It is clear that the MP region (bottom) requires a more significant contribution of nebular continuum than the MR (top) region.

Figure E.7 shows a set of line-strength indices (Burstein et al. 1984; Faber et al. 1985) as a function of stellar metallicity measured for the CB07 templates, at a fixed age of 1 Myr. It also shows the measurement at solar metallicity obtained under different assumptions for the nebular continuum contribution, from an escape fraction of 1 (i.e., no nebular continuum), to an escape fraction of 0 (i.e., the full contribution of nebular emission as computed using FSPS).

Fig. E.7. Set of line-strength indices commonly used as stellar metallicity tracers, measured for the SED of a 1 Myr SSP model from the CB07 library, for the full range of metallicity values available. For solar metallicity, we show the measurement for different levels of nebular continuum contribution. It is clear that a stronger nebular continuum contribution (lower escape fraction) yields a line-strength index measurement consistent with lower metallicities (with TiOCaH being the only exception, bottom right panel). Also, we note that the most metal-poor metallicity bin often appears as a strong outlier in the line-strength-[Z/H] relation.

For most line-strength indices that correlate with stellar metallicity (neglecting the most metal-poor bin), the contribution of nebular emission causes a shift in the index, consistent with a more metal-poor metallicity. Therefore, in the next tests, we explore the impact that including nebular continuum emission in the SSP fitting has in the output maps.

It is worth noting that although FADO does account for nebular emission, it is currently not implemented in python, and thus, we opt for different approaches to account for this contribution using pPXF. To include nebular emission in the fitting, we investigate two different options:

• Include nebular contribution in the templates, as computed by FSPS: this approach has the advantage of including an accurate description of the nebular emission in the fitting. However, it has the drawback that it is only useful for libraries that already include very young templates (i.e., ≲5 Myr). We will refer to this approach as “Neb. A” in the following text.

• Tie the nebular continuum contribution to the emission line measurement: as they are originating from the same source, it is expected that the amount of nebular continuum scales with the flux of emission lines (e.g., Byler et al. 2017). Specifically, we model the nebular continuum emission as a flat spectrum with its flux equal to 0.5% of the Hβ flux, and we subtract this contribution before proceeding with the SSP fit. As this is a very simplistic representation of the nebular contribution, we do not perform additional corrections (i.e., we do not correct by extinction), which could, in principle, lead to an underestimation of the nebular continuum emission. Despite the simplicity of this approach, it offers the advantage of being applicable to any of the libraries of templates used. We will refer to this approach as “Neb. B” in the following text.

E.3.1. Impact of neglecting nebular continuum contribution in synthetic spectra

Here we use a set of mock spectra built from the CB07 library to assess how problematic it can be to not account for existing nebular emission during the full spectrum fitting. To evaluate this, we define our mock-spectra library to harbor a fixed contribution of 100 solar masses of an old (∼10 Gyr) and metal-poor ([Z/H] = -0.7) component, a fixed contribution of 20 solar masses of a intermediate age (∼1 Gyr) and medium metallicity ([Z/H] = -0.4) component, and a variable contribution (30 values covering the range from 0 to 1.5 solar masses) of a young (2 Myr) and solar metallicity component. We have chosen such a configuration aiming for a realistic representation of the combination of multiple stellar populations, in which most of the mass is contributed by old and metal-poor stars (formed from unprocessed gas), and subsequent stellar generations formed from a recycled (and thus, metal-enriched) interstellar medium. Finally, young stars are subdominant in terms of mass, but their light strongly influences the spectrum of the combined population. We build this mock library using young templates that include nebular continuum contribution (computed using FSPS), and then we fit them to recover their original SFH, removing the nebular contribution from the templates used in the fitting. Since we want to explore exclusively the impact of neglecting the contribution of existing nebular emission in the mock spectra, we assume an infinite S/N regime, in order to remove noise-driven differences.

Figure E.8 shows the real luminosity-weighted mean age, luminosity-weighted mean metallicity, total mass, and extinction values as a function of their recovered values, accounting for nebular emission (left column) and neglecting nebular emission (right column) in the fitting. It is remarkable that when existing nebular emission is neglected in the fitting, we recover an anticorrelation between the real metallicity and its measured values. This effect is similar to the trend revealed by the line-strength indices; the nebular continuum contribution reduces the equivalent width of absorption features, mimicking lower metallicities. Another interesting effect is that the stronger continuum is interpreted as a higher total stellar mass, when the templates lack this contribution. This is a direct consequence of the higher level of the continuum due to the nebular continuum emission. However, in our data, we see that the stellar mass surface density is underestimated in the young regions (see, e.g., Fig. E.1), that is, the observed bias in stellar mass in regions expected to have a strong nebular continuum emission contribution is in the opposite direction than that found in this test. A higher extinction is not surprising, as the recovered low-metallicity of the young templates will lead to a bluer spectrum. Thus, a higher extinction value is needed to match the spectral shape of the observed spectrum.

Fig. E.8. Mock-spectra test designed to investigate the impact on the outcome of the SSP fitting when nebular continuum emission is present in the data, but neglected in the analysis. We created a synthetic sample of 30 linear combinations of SSP, with a constant contribution of old (metal-poor), and intermediate age (intermediate metallicity) SSPs, and a varying contribution of a young (metal-rich) SSP that include nebular continuum contribution (computed using FSPS, i.e., Neb. A approach). Both columns show the real age (top row), metallicity (second row), and stellar mass (third row) of each mock spectra in the y-axis, and the recovered value in the x-axis. The bottom panels show an histogram of the recoved E(B-V) values, and the vertical dashed red line shows the real E(B-V) applied to the mock spectra (0.2). The left column shows the results obtained when nebular continuum emission is accounted for in the templates used to fit the data, and the right column shows the same results when nebular continuum emission is neglected from the fitting. It is remarkable that when existing nebular emission is neglected in the fitting, we obtain an anticorrelation between the real metallicity and its recovered values.

E.3.2. Effect of accounting for nebular continuum contribution in the output maps

Figures E.9 and E.10 show the output maps obtained when nebular emission is included in the SSP fitting using the two approaches described earlier (Neb. A and Neb. B, respectively). Despite the impact that nebular continuum emission is expected to have on the spectra of SSPs, particularly in young populations, its addition to the fitting does not lead to a visually significant improvement of the low-metallicity feature.

Fig. E.9. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (second row), stellar E(B-V) (third row), and stellar mass surface density (bottom row) maps derived for the central star-forming ring of NGC 3351. Each column shows the maps obtained using the set of templates indicated at the top of each column, accounting for nebular continuum emission using the approach Neb. A (see main text). The maps obtained using our fiducial set of templates are shown in the left column.

Fig. E.10. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (second row), stellar E(B-V) (third row), and stellar mass surface density (bottom row) maps derived for the central star-forming ring of NGC 3351. Each column shows the maps obtained using the set of templates indicated at the top of each column, accounting for nebular continuum emission using the approach Neb. B (see main text). The maps obtained using our fiducial set of templates (with the modifications described) are shown in the left column.

E.4. Do we trust our full metallicity grid?

In the previous paragraphs, it has become clear that the issue of the young metal-poor feature issue persists, even when we add younger templates to the age grid, remove the stellar continuum of the spectra, or account for nebular continuum emission.

In the next set of tests, we explore if this metal-poor feature could be due to the lack of many hot stars in the stellar libraries used to build the models, that makes the predictions uncertain (Villaume et al. 2017; Conroy et al. 2018). A hint for potential a mismatch between the metal-poor templates and the true galaxy spectra is provided by Fig. E.7, which shows that the most metal-poor bin is often an outlier in the line-strength index versus [Z/H] correlation. Furthermore, metallicity values such as -1.49 (MILES), -1.7 (CB07) and -1.33 (MILES young extensions) are extremely low, and unlikely to be significantly present in the stellar disk of massive star-forming galaxies, where the interstellar medium is being constantly recycled and enriched by star formation and its subsequent feedback (see, e.g., Ho et al. 2017). Thus, we repeat our measurements, removing these extremely low metallicity bins from our grid. We note that the MILES+G-D library already lacks metallicities lower than −1.0 for ages younger than 30 Myr. We will refer to the libraries lacking their most metal-poor metallicity bin as NoMP.

Figure E.11 shows the output maps after removing the most metal-poor metallicity bin from each stellar library. The young metal-poor feature is drastically improved after implementing this change in the metallicity prior, meaning that the most metal-poor values get a significant weight in the fitting of the spectra of the young regions, when they are present.

Fig. E.11. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (second row), stellar E(B-V) (third row), and stellar mass surface density (bottom row) maps derived for the central star-forming ring of NGC 3351. Each column shows the maps obtained using the set of templates indicated at the top of each column, removing the lowest metallicity bins from the metallicity grid (see main text). The maps obtained using our fiducial set of templates (with the modifications described) are shown in the left column.

With this in mind, and for completeness, we explore once again if the inclusion of nebular contribution has a positive impact on the fitting, after removing the most metal-poor bin of each library, by examining the output maps obtained under this combination of settings.

Figures E.12 and E.13 show the maps obtained when nebular emission is included in the fit, together with the removal of low-metallicity templates, using the Neb. A and Neb. B descriptions, respectively. The metallicity maps show weaker low-metallicity features under the Neb. B + NoMP approach, compared to the Neb. A + NoMP one. Nevertheless, none of these tests exhibit a visually significant improvement compared to the NoMP scenario. Thus, it seems that the inclusion of nebular continuum emission does not significantly improve the quality of the output maps, even when the most metal-poor templates have been removed.

Fig. E.12. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (second row), stellar E(B-V) (third row), and stellar mass surface density (bottom row) maps derived for the central star-forming ring of NGC 3351. Each column shows the maps obtained using the set of templates indicated at the top of each column, accounting for nebular continuum emission using the approach Neb. A, and removing the lowest metallicity bins from the metallicity grid (see main text). The maps obtained using our fiducial set of templates (with the modifications described) are shown in the left column.

Fig. E.13. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (second row), stellar E(B-V) (third row), and stellar mass surface density (bottom row) maps derived for the central star-forming ring of NGC 3351. Each column shows the maps obtained using the set of templates indicated at the top of each column, accounting for nebular continuum emission using the approach Neb. B, and removing the lowest metallicity bins from the metallicity grid (see text). The maps obtained using our fiducial set of templates (with the modifications described) are shown in the left column.

E.5. How the different physical conditions in which young stars are embedded bias the results of our SSP fitting?

Young stellar populations are embedded in deeper layers of gas and dust, remnants of the cloud from which they formed, than old stellar populations. In this work, we have assumed a single extinction value for all the stellar populations present in a given region, parametrized by a Calzetti et al. (2000) extinction law. Although young stars disperse molecular clouds on relatively short time scales of a few Myr (Chevance et al. 2020), young stellar populations in star-forming regions are likely affected by greater extinction values than their older counterparts. The impact of this different extinction across coexisting stellar populations on the measured SFHs has not been determined yet. To evaluate the effect that this feature would have on the SSP fitting, we have repeated the test with mock data described in Appendix E.3.1, but this time, applying the extinction curve (E(B-V) = 0.2) only to the SED of the youngest stellar population.

Figure E.14 is analogous to Fig. E.8, except that in this case, each input spectrum is composed by constant a linear combination of unextincted old (∼10 Gyr) and intermediate age (∼1 Gyr) stellar population, and a variable contribution of a dust-extincted young (2 Myr) stellar population (see Appendix E.3.1 for details). The figure shows that even when the templates used to build the mock spectra are identical to the templates used to measure their stellar properties (left column), the measured age, metallicity, mass, and extinction are severely biased respect to the real values. Moreover, the bias in the metallicity is particularly complex, being always toward more metal-poor values, but showing two clearly different regimes. Lower contributions of young stars lead to an anticorrelation between the measured and the true metallicity, whereas for higher contribution of young stars, we see a nearly linear correlation between the measured metallicity and the true metallicity. Nevertheless, the dynamical range of the measured metallicity is small, implying that overall, an extinction curve applied exclusively to young (and metal-rich) stellar populations lead to apparently lower metallicities.

Fig. E.14. Analogous to Fig. E.8, but in this case, dust extinction (E(B-V) = 0.2) has been applied only to the SED of the youngest stellar population. The left column shows the results obtained when nebular continuum emission is accounted for in the templates used to fit the data, and the right column shows the same results when nebular continuum emission is neglected from the fitting. The figure shows that even when the templates used to build the mock spectra are identical to the templates used to measure their stellar properties, the measured age, metallicity, mass, and extinction values are severely biased respect to the real values.

We stress here that although the specific magnitude of the biases found for the stellar population properties driven by this effect might depend on the assumptions made in this test (e.g., SSP models used to build our mock data), the general trend is preserved under different assumptions. Thus, this effect could indeed contribute to the unphysically low stellar metallicities measured for very young regions.

Finally, implementing a fitting routine that accounts for this feature (e.g., using a different extinction value for each kinematic component, in the case of pPXF) is beyond the scope of this paper. However, this test shows that addressing this issue would lead to more reliable SFHs, and thus improve the utility of these spectral fitting tools.

E.6. Deciding on best fitting results

In the previous sections, we have shown the output log LumW age, LumW [Z/H], E(B-V), and log Σ_* maps obtained using different spectral libraries and approaches to fit our data in order to obtain a visualization of the resulting maps that these different approaches yield. We parametrize quantitatively the quality of these results, in terms of Δlog age, Δ[Z/H], ΔE(B-V), and ΔlogΣ_*, introduced earlier.

Figure E.15 shows the values calculated for these quantities, for each one of the sets of templates used, under each different approach. The fiducial case for the four libraries is indicated in gray. The two implementations of nebular emission are marked in red and blue. The tests removing the most metal-poor templates are marked in yellow, and finally, the tests that include a combination of including nebular emission, and removing the most metal-poor templates are indicated with green and pink backgrounds.

Fig. E.15. Values measured for Δlog age, Δ[Z/H], ΔE(B-V), and ΔlogΣ*, following Eq. E.1, from the distribution of each one of these four quantities within the star-forming ring of NGC 3351, and in its surrounding disk, using each one of the sets of templates described in Sec. E. We find that the MILES young templates, with the most metal-poor metallicity bin removed from the grid offers a good compromise between improving the young and metal-poor feature, reaching younger ages than our fiducial set of templates, and not leading to noticeable additional artifacts neither in the stellar mass surface density, nor in the stellar extinction map. See discussion in Sec. E.6.

The figure clearly shows that the highest (closest to zero) Δ[Z/H] values (2nd row) are achieved when the most metal-poor templates are removed from the metallicity grid (yellow, green, and pink background). In particular, the CB07 templates yield the highest Δ[Z/H] value among the set of templates tested. This is visually consistent with what can be seen in Fig. E.11. However, the last row shows that the CB07 templates also yield low stellar mass surface density values in the young regions (i.e., low ΔlogΣ_*). This low-mass feature is also clearly visible in Fig. E.11. Furthermore, in some cases, the CB07 templates also yield particularly low E(B-V) values in the young regions. Out of the remaining set of templates (MILES + G-D and MILES young, with metal-poor metallicity bins removed), both show similarly low Δlog age, and high Δ[Z/H], without any noticeable artifacts neither in E(B-V) nor in stellar mass surface density. Since the MILES + G-D templates suffer from a low spectral resolution in the youngest ages, and this different spectral resolution is not properly treated by the DAP in its current state, we finally choose the MILES young templates (NoMP) for future analysis of the stellar populations.

Despite choosing the best set of templates according to our criteria, among those libraries tested, we acknowledge that the fitting of the youngest regions is still suboptimal, as the metal-poor metallicities persist for these regions. We consider these metallicities unreliably, and therefore mask them in further analysis. We have tested and adopted an extinction-corrected Hα⁵ surface density emission threshold for our masking scheme, as it traces ongoing star-formation and thus, the presence of young stellar populations. Specifically, we masked out from the metallicity maps those pixels in which we measure logHα_corr > 35 erg s⁻¹ pc⁻²⁶. We find that this value successfully identifies the youngest and most problematic regions. Figure E.16 shows the extinction-corrected Hα and LumW [Z/H] maps of the central region of NGC 3351 after applying our masking scheme. The low metallicity pixels with high Hα emission values are thus removed from our sample. Although this masking scheme can bias the measured mean metallicity in environments with significant presence of young regions (e.g., sp. arms and rings), it is necessary in order to keep our analysis free from this unphysical values.

Fig. E.16. Stellar metallicity (left) map of the center of NGC 3351, after masking pixels with an extinction-corrected Hα surface density (right) emission higher than 10−35 erg s−1 pc−2. This value successfully characterizes the youngest and most problematic regions.

There are ongoing efforts to further improve the SSP fitting in the regions dominated by young stellar populations. It is known that stellar populations of different ages are also characterized by different kinematic conditions, with older stars being dynamically hotter than younger stars (see, e.g., Tarricq et al. 2021). In this regard, Zhang et al. (in prep.) will present a detailed exploration on how these differences in the kinematic conditions across different stellar populations can bias the recovered SFHs. This will potentially lead to improvements in the methodology currently employed in the PHANGS-MUSE pipeline, by including multiple stellar kinematic components to reproduce the observed spectra. An additional possibility to significantly improve our SSP fitting is using available photometric data from the PHANGS-HST survey to constrain bluer wavelength ranges, valuable to shed light on the properties of young stellar populations (see, e.g., González Delgado et al. 2005). However, the inclusion of additional photometric data to the spectral fitting technique is not trivial. It requires further development of the currently available tools, and thus, it is beyond the scope of the analysis presented here.

All Tables

All Figures

Fig. 1. Example of some of the DAP outputs obtained for the mosaic of NGC 1566. Stellar mass surface density (top-left), luminosity-weighted age (top-center), luminosity-weighted [Z/H] (top-right), stellar velocity dispersion (bottom-left), mass-weighted age (bottom-center), and mass-weighted [Z/H] (bottom-right). The spiral arms and bar of the galaxy are marked with blue and green contours, respectively.

In the text

Fig. 2. Luminosity-weighted stellar age radial profiles for the galaxies in our sample. Different colors indicate the radial profile measured across different environments, as indicated in the legend of the bottom-right panel. The black line shows the radial profile measured for the entire FoV (i.e., all environments together). The galactocentric distance is measured in units of R25, in order to measure radial distance homogeneously across our sample. The value of R25 (kpc) of each galaxy is indicated in each corresponding panel. The solid gray line shows the best-fit linear gradient for each galaxy.

In the text

Fig. 3. Luminosity-weighted stellar age radial profiles for the galaxies in our sample, considering all environments. Top left panel shows the measured LumW log Age radial profiles, colored by total galaxy stellar mass. Top right panel shows the profiles normalized by their value at 0.3 R25, to better compare their shapes. The bottom panels show the slope (ΔLumW log Age) and intercept of the LumW radial age profile (indicated with a solid gray line in each panel of Fig. 2) as a function of the total stellar mass of galaxies, color coded by total specific star formation rate (sSFR).

In the text

	Fig. 4. Distribution of the offset of the LumW log Age of the pixels within each environments, with respect to the median galaxy log LumW Age of each galaxy, in bins of 0.04 dex. Different colors indicate the distribution of different environments, as indicated in the legend.
In the text

Fig. 5. Mass-weighted stellar age radial profiles for the galaxies in our sample. Different colors indicate the radial profile measured across different environments, as indicated in the legend of the bottom-right panel. The black line shows the radial profile measured for the entire FoV (i.e., all environments together). The galactocentric distance is measured in units of R25, in order to measure radial distance homogeneously across our sample. The value of R25 (kpc) of each galaxy is indicated in each corresponding panel. The solid gray line shows the best-fit gradient for each galaxy.

In the text

Fig. 6. Mass-weighted stellar age radial profiles for the galaxies in our sample, considering all environments. Top left panel shows the measured MassW log Age radial profiles, colored by total galaxy stellar mass. Top right panel shows the profiles normalized by their value at 0.3 R25, to better compare their shapes. The bottom panels show the slope (ΔMassW log Age) and intercept of the MassW radial age profile (indicated with a solid gray line in each panel of Fig. 5) as a function of the total stellar mass of galaxies, color coded by total specific star formation rate (sSFR).

In the text

Fig. 7. Luminosity-weighted stellar [Z/H] radial profiles for the galaxies in our sample. Different colors indicate the radial profile measured across different environments, as indicated in the legend of the bottom-right panel. The black line shows the radial profile measured for the entire FoV (i.e., all environments together). The galactocentric distance is measured in units of R25, in order to measure radial distance homogeneously across our sample. The value of R25 (kpc) of each galaxy is indicated in each corresponding panel. The solid gray line shows the best-fit gradient for each galaxy.

In the text

Fig. 8. Mass-weighted stellar [Z/H] radial profiles for the galaxies in our sample. Different colors indicate the radial profile measured across different environments, as indicated in the legend of the bottom-right panel. The black line shows the radial profile measured for the entire FoV (i.e., all environments together). The galactocentric distance is measured in units of R25, in order to measure radial distance homogeneously across our sample. The value of R25 (kpc) of each galaxy is indicated in each corresponding panel. The solid gray line shows the best-fit gradient for each galaxy.

In the text

Fig. 9. Luminosity-weighted stellar [Z/H] radial profiles for the galaxies in our sample, considering all environments. Top left panel shows the measured LumW [Z/H] radial profiles, colored by total galaxy stellar mass. Top right panel shows the profiles normalized by their value at 0.3 R25, to better compare their shapes. The bottom panels show the slope (ΔLumW [Z/H]) and intercept of the LumW radial [Z/H] profile (indicated with a solid gray line in each panel of Fig. 7) as a function of the total stellar mass of galaxies, color coded by total specific star formation rate (sSFR).

In the text

Fig. 10. Mass-weighted stellar [Z/H] radial profiles for the galaxies in our sample, considering all environments. Top left panel shows the measured MassW [Z/H] radial profiles, colored by total galaxy stellar mass. Top right panel shows the profiles normalized by their value at 0.3 R25, to better compare their shapes. The bottom panels show the slope (ΔMassW [Z/H]) and intercept of the MassW radial [Z/H] profile (indicated with a solid gray line in each panel of Fig. 8) as a function of the total stellar mass of galaxies, color coded by total specific star formation rate (sSFR).

In the text

Fig. 11. Radial assembly history of our sample galaxies. The y-axis shows the age of the Universe at which each radial bin assembled 80% of its total current stellar mass. A positive trend in these plots implies that outer regions assembled their stellar mass later in cosmic history than inner regions. The yellow solid lines show the best-fit linear model for each galaxy. The red dashed line show the time at which the full galaxy formed 80% of its total current stellar mass. The gray shaded area represents the uncertainty in our measurement, estimated by propagating the stellar mass uncertainty in the integration of the SFH of each radial bin. The value of R25 (kpc) of each galaxy is indicated in each corresponding panel. The bottom-right panel displays the trends for all galaxies, color-coded by total stellar mass. By construction, the highest possible value for the y-axis is given by the age of our oldest age bin in the SSP fitting (∼15 Gyr, see Sect. 2.2).

In the text

Fig. 12. Radial assembly history profiles for measured for our sample galaxies, considering all environments. Top left panel shows the measured profiles, colored by total galaxy stellar mass. Top right panel shows the profiles normalized by their value at 0.3 R25, to better compare their shapes. The bottom panels show the slope and intercept of the radial profiles (indicated with a solid yellow line in each panel of Fig. 11) as a function of the total stellar mass of galaxies, color coded by total specific star formation rate (sSFR).

In the text

	Fig. 13. Representation of the location of the pseudo-slits along the bar of NGC 1433 (pale and dark green rectangles), and perpendicular to the bar axis (gray and black rectangles). The position of the bar is indicated by the light-green ellipses. The background shows the stellar mass surface density map of the galaxy.
In the text

Fig. 14. Mass-weighted age (left) and [Z/H] (right) radial profiles along the bar major axis (dark green), and along its perpendicular direction (black), for NGC 1433. The solid and dashed dark green and black lines represent the gradient measured toward each possible direction from the center. The pale green dashed lines show the best-fitting linear model to the bar major axis gradient (dark green lines), and the silver dashed lines show the best-fitting to the gradients measured in the direction perpendicular to the bar.

In the text

Fig. 15. Mass-weighted age gradients measured along bar major axis (green pentagons), and along its perpendicular direction (black squares), as a function of the bar length (in units kpc and R25), total stellar mass, and offset from the global main sequence of galaxies (ΔMS), for 14 barred galaxies in our sample. The red line is calculated following Eq. (4). It is defined to be positive in galaxies that show a flatter gradient along their bar major axes, compared to the bar perpendicular direction, and negative values otherwise.

In the text

Fig. 16. Mass-weighted [Z/H] gradients measured along bar major axis (green pentagons), and along its perpendicular direction (black squares), as a function of the bar length (in units kpc and R25), total stellar mass, and offset from the global main sequence of galaxies (ΔMS), for 14 barred galaxies in our sample. The red line is calculated following Eq. (5). It is defined to be positive in galaxies that show a flatter gradient along their bar major axis, compared to the bar perpendicular direction, and negative values otherwise.

In the text

Fig. 17. 2D histogram that shows the distribution of stellar velocity dispersion (σ*) and log age (luminosity-weighted), for NGC 1365 (left), NGC 1566 (center), and NGC 3627 (right). The color of each [LumW log Age, σ*] bin scales with the number of pixels within that bin. The contours show the one- and two-sigma limits of the distribution of each individual galactic environment, following the color code indicated in the top-right corner of the right of the panel.

In the text

Fig. 18. Radial profile of stellar velocity dispersion (σ*) for NGC 1365 (left), NGC 1566 (center), and NGC 3627 (right), plotted as a function of galactic structure (colors) and LumW age of the underlying stellar populations. Stellar populations younger than 100 Myr are indicated by diamonds (dotted line). Stellar populations of ages in the range of 100 Myr < log LumW age < 600 Myr are marked with triangles (dashed line). Stars older than 3 Gyr are represented by squares (dashdotted line). Circles (solid line) show the radial profile considering all stellar populations. The gray line shows the radial velocity dispersion gradient measured for the molecular gas, from the PHANGS-ALMA CO(2−1) data.

In the text

Fig. 19. Mean stellar velocity dispersion (σ*) measured at a galactocentric radius of 0.25−0.30 R25 in the disk environment of all our sample galaxies for stellar populations with luminosity-weighted ages within the three different age bins defined in the main text. We note that for visualization purposes, the x-axis position of each data point is located at the edge of the age bin it represents. In each panel, the leftmost triangle represents the mean (area-weighted) velocity dispersion of all populations younger than 100 Myr. The central square marks the mean velocity dispersion of populations with 100 Myr < LumW age < 300 Myr, and the rightmost triangle shows the same quantity for populations older than 3 Gyr. The errorbar indicates the standard deviation of the σ* values of the pixels within a given age bin. The bottom right panel shows the trend for all galaxies, colored by total stellar mass.

In the text

	Fig. 20. Time-averaged SFR in five different time bins, calculated from the spatially resolved SFH, for NGC 1566. The panels show how different regions of the galactic disk have assembled their stellar mass in different moments of cosmic history. The blue and green contours mark the position of the spiral arms and the bar, respectively.
In the text

	Fig. 21. Time-averaged SFR in five different time bins, calculated from the spatially resolved SFH, for NGC 3627. The panels show how different regions of the galactic disk have assembled their stellar mass in different moments of cosmic history. The blue and green contours mark the position of the spiral arms and the bar, respectively.
In the text

Fig. A.1. Slope of the LumW-age gradient, as a function of (from left to right, and from top to bottom) total stellar mass (log M*), total SFR (log SFR), offset from the global main sequence of galaxies (ΔMS), total molecular gas mass-to-stellar mass ratio (Mmol/M*), total atomic gas mass-to-stellar mass ratio (MHI), total gas (molecular + atomic)-to-stellar mass ratio Mgas/M*, and molecular-to-atomic gas mass ratio (Mmol/MHI), galaxy inclination and distance. The Pearson correlation coefficient for each parameter is indicated in the bottom-right corner of each panel.

In the text

Fig. A.2. Slope of the MassW-age gradient, as a function of (from left to right, and from top to bottom) total stellar mass (log M*), total SFR (log SFR), offset from the global main sequence of galaxies (ΔMS), total molecular gas mass-to-stellar mass ratio (Mmol/M*), total atomic gas mass-to-stellar mass ratio (MHI), total gas (molecular + atomic)-to-stellar mass ratio Mgas/M*, and molecular-to-atomic gas mass ratio (Mmol/MHI), galaxy inclination and distance. The Pearson correlation coefficient for each parameter is indicated in the bottom-right corner of each panel.

In the text

Fig. A.3. Slope of the LumW-[Z/H] gradient, as a function of (from left to right, and from top to bottom) total stellar mass (log M*), total SFR (log SFR), offset from the global main sequence of galaxies (ΔMS), total molecular gas mass-to-stellar mass ratio (Mmol/M*), total atomic gas mass-to-stellar mass ratio (MHI), total gas (molecular + atomic)-to-stellar mass ratio Mgas/M*, and molecular-to-atomic gas mass ratio (Mmol/MHI), galaxy inclination and distance. The Pearson correlation coefficient for each parameter is indicated in the bottom-right corner of each panel.

In the text

Fig. A.4. Slope of the MassW-[Z/H] gradient, as a function of (from left to right, and from top to bottom) total stellar mass (log M*), total SFR (log SFR), offset from the global main sequence of galaxies (ΔMS), total molecular gas mass-to-stellar mass ratio (Mmol/M*), total atomic gas mass-to-stellar mass ratio (MHI), total gas (molecular + atomic)-to-stellar mass ratio Mgas/M*, and molecular-to-atomic gas mass ratio (Mmol/MHI), galaxy inclination and distance. The Pearson correlation coefficient for each parameter is indicated in the bottom-right corner of each panel.

In the text

	Fig. B.1. Probability density function of the relative uncertainties estimated for stellar mass surface density (Σ*), LumW Age, MassW Age, LumW [Z/H], and MassW [Z/H], for our full sample, that is, all available pixels.
In the text

Fig. B.2. Dispersion of the resulting posterior distribution of stellar mass velocity dispersions (Σ*), ages and metallicity from the Monte Carlo realizations, as a function of the number of Monte Carlo iterations, for 5 representative Voronoi bins. The color code indicates the galactic environment of each Voronoi bin. Solid lines correspond to mass-weighted quantities, and dashed lines to their light-weighted counterparts. The vertical black line indicates the adopted number (20) of Monte Carlo iterations to estimate the uncertainties in the analyses presented in this paper.

In the text

Fig. C.1. Star formation history (top) and star formation rate history (bottom) derived for NGC 1365 (left) and NGC 1566 (right). The best-fitting linear model for the SFH curve of each galaxy is shown as a red dashed line, and its slope and intercept are shown in the top-left corner of the top panels. The shaded areas mark sudden decreases in the SFR at ages log (age/yr) ∼7.7 and log (age/yr) ∼9.5, found to exist in all galaxies in some degree, deemed as an artifact of the computed SFHs.

In the text

Fig. C.2. Slope of the SFH of PHANGS galaxies, as a function of (from left to right, and from top to bottom) total stellar mass (log M*), total SFR (log SFR), offset from the global main sequence of galaxies (ΔMS), total molecular gas mass-to-stellar mass ratio (Mmol/M*), total atomic gas mass-to-stellar mass ratio (MHI), total gas (molecular + atomic)-to-stellar mass ratio Mgas/M*, and molecular-to-atomic gas mass ratio (Mmol/MHI), galaxy inclination and distance. The Pearson correlation coefficient for each parameter is indicated in the bottom-right corner of each panel.

In the text

Fig. C.3. Downsizing paradigm in separated galactic environments. Left: correlation between the median MassW log Age and the total galaxy stellar mass for the (combined) disk and center environments, for each galaxy in our sample. The total galaxy stellar mass corresponds to the sum of the masses of the (combined) disk and center environments. Middle: median MassW log Age of the (combined) disk and center environments as a function of total stellar mass in each galactic environment. Right: Median MassW log Age of the (combined) disk environment as a function of the median MassW log Age age of the (combined) center environment for our sample galaxies. Galaxies that present a nuclear ring are indicated with an additional orange marker.

In the text

Fig. D.1. Luminosity-weighted stellar age radial profiles for the galaxies in our sample, measured after convolving our MUSE mosaics to a fixed angular resolution of 15 arcsecs, and computing the stellar population properties for these convolved cubes. Different colors indicate the radial profile measured across different environments, as indicated in the legend of the bottom-right panel. The black line shows the radial profile measured for the entire FoV (i.e., all environments together). The galactocentric distance is measured in units of R25, in order to measure radial distance homogeneously across our sample. The value of R25 (kpc) of each galaxy is indicated in each corresponding panel. The solid gray line shows the best-fit gradient for each galaxy. The bottom-right panel displays the entire FoV trends for all galaxies, color-coded by total stellar mass.

In the text

	Fig. D.2. Distribution of the offset of the LumW log Age of the pixels within each environments, with respect to the median galaxy log LumW Age of each galaxy, in bins of 0.04 dex, obtained for the MUSE mosaics convolved to a fixed angular resolution of 15 arcsecs. Different colors indicate the distribution of different environments, as indicated in the legend.
In the text

Fig. E.1. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (second row), stellar E(B-V) (third row), and stellar mass surface density (bottom row) maps derived for the central star-forming ring of NGC 3351. Each column shows the maps obtained using the set of templates indicated at the top of each column. The maps obtained using our fiducial set of templates are shown in the left column. The area outside the red ellipse in the top-left panel is used to probe disk stellar population properties, and the area between the red and the black ellipses is used to probe stellar population properties of the young star-forming ring (see text).

In the text

Fig. E.2. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (second row), stellar E(B-V) (third row), and stellar mass surface density (bottom row) normalized distributions for the star-forming ring (black) and disk (red), delimitated in the top-left panel of Fig. E.1. Each column shows the distributions obtained using the set of templates indicated in the legends of the top-row panels. The distributions obtained using our fiducial set of templates are shown in the left column.

In the text

Fig. E.3. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (middle row), and stellar E(B-V) (bottom row) maps derived for the central star-forming ring of NGC 3351, once the stellar continuum has been removed from the observed spectra and the templates. Each column shows the maps obtained using the set of templates indicated at the top of each column. The maps obtained using our fiducial set of templates (after removing the stellar continuum) are shown in the left column. We note that the E(B-V) maps are identical to those in Fig. E.2, as the stellar extinction is measured before subtracting the stellar continuum.

In the text

	Fig. E.4. Example of the impact of nebular continuum emission on the SED of a 1 Myr old solar metallicity SSP. The sudden decrease of nebular emission (green line) at ∼3646 Å and ∼8207 Å correspond to the Balmer and Paschen breaks, respectively. The vertical red dashed lines enclose the wavelength range considered for the full spectral fitting [4850,7000] Å.
In the text

	Fig. E.5. Two very young regions within the central star-forming ring of NGC 3351 with different inferred metallicity values. The orange circle encloses the metal-rich (MR) region, and the blue ellipse encloses the metal-poor (MP) region. The background map shows the luminosity-weighted stellar metallicity in the central region.
In the text

Fig. E.6. Best SED-fitting obtained with FADO for the MR and MP regions, defined in Fig. E.5. The orange and blue lines represent the MR and MP observed spectra, respectively. The black line show the best-fitting model, and the gray and green curves show separately the stellar and nebular continuum contribution, respectively. It is clear that the MP region (bottom) requires a more significant contribution of nebular continuum than the MR (top) region.

In the text

Fig. E.7. Set of line-strength indices commonly used as stellar metallicity tracers, measured for the SED of a 1 Myr SSP model from the CB07 library, for the full range of metallicity values available. For solar metallicity, we show the measurement for different levels of nebular continuum contribution. It is clear that a stronger nebular continuum contribution (lower escape fraction) yields a line-strength index measurement consistent with lower metallicities (with TiOCaH being the only exception, bottom right panel). Also, we note that the most metal-poor metallicity bin often appears as a strong outlier in the line-strength-[Z/H] relation.

In the text

Fig. E.8. Mock-spectra test designed to investigate the impact on the outcome of the SSP fitting when nebular continuum emission is present in the data, but neglected in the analysis. We created a synthetic sample of 30 linear combinations of SSP, with a constant contribution of old (metal-poor), and intermediate age (intermediate metallicity) SSPs, and a varying contribution of a young (metal-rich) SSP that include nebular continuum contribution (computed using FSPS, i.e., Neb. A approach). Both columns show the real age (top row), metallicity (second row), and stellar mass (third row) of each mock spectra in the y-axis, and the recovered value in the x-axis. The bottom panels show an histogram of the recoved E(B-V) values, and the vertical dashed red line shows the real E(B-V) applied to the mock spectra (0.2). The left column shows the results obtained when nebular continuum emission is accounted for in the templates used to fit the data, and the right column shows the same results when nebular continuum emission is neglected from the fitting. It is remarkable that when existing nebular emission is neglected in the fitting, we obtain an anticorrelation between the real metallicity and its recovered values.

In the text

Fig. E.9. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (second row), stellar E(B-V) (third row), and stellar mass surface density (bottom row) maps derived for the central star-forming ring of NGC 3351. Each column shows the maps obtained using the set of templates indicated at the top of each column, accounting for nebular continuum emission using the approach Neb. A (see main text). The maps obtained using our fiducial set of templates are shown in the left column.

In the text

Fig. E.10. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (second row), stellar E(B-V) (third row), and stellar mass surface density (bottom row) maps derived for the central star-forming ring of NGC 3351. Each column shows the maps obtained using the set of templates indicated at the top of each column, accounting for nebular continuum emission using the approach Neb. B (see main text). The maps obtained using our fiducial set of templates (with the modifications described) are shown in the left column.

In the text

Fig. E.11. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (second row), stellar E(B-V) (third row), and stellar mass surface density (bottom row) maps derived for the central star-forming ring of NGC 3351. Each column shows the maps obtained using the set of templates indicated at the top of each column, removing the lowest metallicity bins from the metallicity grid (see main text). The maps obtained using our fiducial set of templates (with the modifications described) are shown in the left column.

In the text

Fig. E.12. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (second row), stellar E(B-V) (third row), and stellar mass surface density (bottom row) maps derived for the central star-forming ring of NGC 3351. Each column shows the maps obtained using the set of templates indicated at the top of each column, accounting for nebular continuum emission using the approach Neb. A, and removing the lowest metallicity bins from the metallicity grid (see main text). The maps obtained using our fiducial set of templates (with the modifications described) are shown in the left column.

In the text

Fig. E.13. Luminosity-weighted age (top row), luminosity-weighted [Z/H] (second row), stellar E(B-V) (third row), and stellar mass surface density (bottom row) maps derived for the central star-forming ring of NGC 3351. Each column shows the maps obtained using the set of templates indicated at the top of each column, accounting for nebular continuum emission using the approach Neb. B, and removing the lowest metallicity bins from the metallicity grid (see text). The maps obtained using our fiducial set of templates (with the modifications described) are shown in the left column.

In the text

Fig. E.14. Analogous to Fig. E.8, but in this case, dust extinction (E(B-V) = 0.2) has been applied only to the SED of the youngest stellar population. The left column shows the results obtained when nebular continuum emission is accounted for in the templates used to fit the data, and the right column shows the same results when nebular continuum emission is neglected from the fitting. The figure shows that even when the templates used to build the mock spectra are identical to the templates used to measure their stellar properties, the measured age, metallicity, mass, and extinction values are severely biased respect to the real values.

In the text

Fig. E.15. Values measured for Δlog age, Δ[Z/H], ΔE(B-V), and ΔlogΣ*, following Eq. E.1, from the distribution of each one of these four quantities within the star-forming ring of NGC 3351, and in its surrounding disk, using each one of the sets of templates described in Sec. E. We find that the MILES young templates, with the most metal-poor metallicity bin removed from the grid offers a good compromise between improving the young and metal-poor feature, reaching younger ages than our fiducial set of templates, and not leading to noticeable additional artifacts neither in the stellar mass surface density, nor in the stellar extinction map. See discussion in Sec. E.6.

In the text

	Fig. E.16. Stellar metallicity (left) map of the center of NGC 3351, after masking pixels with an extinction-corrected Hα surface density (right) emission higher than 10−35 erg s−1 pc−2. This value successfully characterizes the youngest and most problematic regions.
In the text