The universal variability of the stellar initial mass function probed by the TIMER survey

The debate about the universality of the stellar initial mass function (IMF) revolves around two competing lines of evidence. While measurements in the Milky Way, an archetypal spiral galaxy, seem to support an invariant IMF, the observed properties of massive early-type galaxies (ETGs) favor an IMF somehow sensitive to the local star formation conditions. The fundamental methodological and physical di ff erences between both approaches have hampered, however, a comprehensive understanding of IMF variations. We describe here an improved modelling scheme that allows for the first time consistent IMF measurements across stellar populations with di ff erent ages and complex star formation histories. Making use of the exquisite MUSE optical data from the TIMER survey and powered by the MILES stellar population models, we show the age, metallicity, [Mg / Fe], and IMF slope maps of the inner regions of NGC3351, a spiral galaxy with a mass similar to that of the Milky Way. The measured IMF values in NGC3351 follow the expectations from a Milky Way-like IMF, although they simultaneously show systematic and spatially coherent variations, particularly for low-mass stars. In addition, our stellar population analysis reveals the presence of metal-poor and Mg-enhanced star-forming regions that appear to be predominantly enriched by the stellar ejecta of core-collapse supernovae. Our findings showcase therefore the potential of detailed studies of young stellar populations to better understand the early stages of galaxy evolution and, in particular, the origin of the observed IMF variations beyond and within the Milky Way.


Introduction
Much has been debated about the universality of the stellar initial mass function (IMF), which describes the (idealized) mass spectrum of stars at birth.The IMF is a central idea in our understanding and modeling of the baryonic cycle in galaxies, bridging the quantum processes driving stellar evolution to the cosmological scales of galaxy formation (see e.g., Kroupa & Jerabkova 2019, for a recent overview).Therefore, measuring the IMF is a critical and necessary endeavor; however, doing so with observations has proven to be particularly challenging (e.g., Bastian et al. 2010;Smith 2020).
The analysis of nearby systems, where individual stars can be resolved and therefore direct measurements of the IMF are feasible, has revealed a rather invariant IMF behavior.This Milky Way-like standard IMF is characterized by a power-law distribution with a logarithmic slope of Γ 1.3 for stars more massive than ∼0.5 M and a shallower distribution for lower-mass stars (e.g., Miller & Scalo 1979;Kroupa 2001Kroupa , 2002;;Chabrier 2003).
These observational findings cemented the now widely adopted assumption of a universal IMF.
However, the simplicity of a universal IMF faces serious challenges.From a theoretical perspective, it remains unclear as to why, for example, the IMF should or should not depend on the local conditions of the gas from which new stars are formed (Hennebelle & Chabrier 2008;Myers et al. 2011;Hopkins 2012;Krumholz 2014;Chabrier et al. 2014;Fontanot et al. 2018;Guszejnov et al. 2019;Davis & van de Voort 2020;Sharda & Krumholz 2022;Tanvir et al. 2022).Furthermore, even within the Milky Way, systematic deviations from a universal IMF have been evoked in order to explain the observed properties of Galactic globular clusters (Marks et al. 2012).
However, although central to our current understanding of the IMF, massive quiescent ETGs can only provide information about the IMF over a very limited stellar mass range, from the hydrogen-burning limit to ∼1 M , because stars more massive than that have long evolved into dark remnants.Measurements of the high-mass end of the IMF in young stellar populations have been generally limited to emission-line studies, in particular comparing the ionizing flux of massive stars to the overall stellar continuum (e.g., UV flux or optical colors; Hoversten & Glazebrook 2008;Meurer et al. 2009;Lee et al. 2009;Nanayakkara et al. 2017), providing tentative evidence of an IMF deficient in high-mass stars in dwarf galaxies.Complementary, submillimetric measurements (e.g., Romano et al. 2017) of nearby and distant star-forming galaxies suggest that the IMF in these relatively massive objects is biased toward massive stars (Sliwa et al. 2017;Brown & Wilson 2019;Zhang et al. 2018), leaving the possibility for an IMF simultaneously bottomand top-heavy in massive galaxies (e.g., Ferreras et al. 2019).
In this paper, we present an alternative approach to constraining the IMF in young stellar populations based on the analysis of absorption spectra, building upon the success of stellar population synthesis models in revealing the variable nature of the IMF in massive quiescent galaxies.We apply this novel approach to a nearby late-type galaxy (LTG), NGC 3351, drawn from the Time Inference with MUSE in Extragalactic Rings (TIMER) survey (Gadotti et al. 2019), obtaining a spatially resolved two-dimensional map of the IMF in this galaxy.While IMF maps have already been derived for several nearby galaxies (Martín-Navarro et al. 2019, 2021;Barbosa et al. 2021), this is the first time that this has been achieved for a young star-forming system.The layout of this paper is as follows: Sect. 2 presents the TIMER data of NGC 3351.In Sect. 3 we give a detailed explanation of our stellar-population modeling approach.The results of our analysis of NGC 3351 are presented in Sect. 4 and discussed in Sect. 5. Finally, Sect.6 summarizes our findings and lists our concluding remarks.

Data
The TIMER project presented in Gadotti et al. (2019) is a MUSE survey of 24 nearby barred galaxies.The exquisite quality and depth of the TIMER data, exemplified in some of its most recent results (e.g., de Lorenzo-Cáceres et al. 2019;Méndez-Abreu et al. 2019;Gadotti et al. 2020;Neumann et al. 2020;Bittner et al. 2021) are ideally suited to measure the subtle effect that a variable IMF has on the integrated spectra of galaxies.These data were obtained through the MUSE integral field unit (Bacon et al. 2010), which cover a wavelength range from ∼4700 Å to ∼9000 Å at a mean resolution of R ∼ 3000.A detailed description of the sample and data properties are given in the survey presentation paper (Gadotti et al. 2019).This paper focuses on the LTG NGC 3351, which has an estimated stellar mass of M = 3.1 × 10 10 M (Muñoz-Mateos et al. 2015) and a distance of 10.1 Mpc 1 .Figure 1 shows a HST color image of NGC 3351 (F438W, F555W, F814W) 2 , where its main morphological features are clearly visible, including the starforming ring and the bar, extending diagonally across the HST image and beyond the MUSE field of view (FOV).In the context of the TIMER survey, Leaman et al. (2019) proposed that the observed (ionized and molecular) outflow in NGC 3351 is powered by the energy and momentum injected by its star forming ring.Moreover, Bittner et al. (2020) reported that the stellar populations in this star forming ring have a relatively high [Mg/Fe] ratio and low metallicity but the origin of such a potentially chemically peculiar gas remains unknown.The relatively high mass and thus luminosity of NGC 3351 and its peculiar stellar population content make it an ideal test case for our pilot study.
Although the spatial resolution of the MUSE data is of ∼10 pc/spaxel, measuring the possible effect of the IMF requires high-signal-to-noise-ratio (S/N) spectra.Therefore, instead of analyzing the stellar population properties on a spaxel-by-spaxel 1 At this distance, the pixel scale of the MUSE spectrograph translates into a spatial scale of ∼10 pc/spaxel. 2Program ID 15654, PI: Janice Lee.A110, page 2 of 15 basis, we used the Voronoi binning code presented in Cappellari & Copin (2003) to spatially combine nearby spaxels in order to achieve a minimum S/N per Å of 100.This binning procedure effectively decreases the spatial scales we can resolve but ensures that the individual spectra to be fit are deep enough for our detailed stellar population analysis.All the results presented in this paper are therefore derived from this Voronoibinned MUSE data cube.

Stellar population synthesis models
Our approach is based on the MILES stellar population synthesis models (Vazdekis et al. 2010).In particular, we use the α-variable models presented in Vazdekis et al. (2015) in order to account for possible variations in the abundance pattern.These models are built by combining the BaSTI set of isochrones (Pietrinferni et al. 2004(Pietrinferni et al. , 2006) ) with the MILES stellar library (Sánchez-Blázquez et al. 2006) and the synthetic stellar library of Coelho et al. (2005).The main advantage of these models is their consistent and semi-empirical treatment of α-element variations, providing intermediate-resolution model spectra from ∼3500 Å to ∼7400 Å (see also Walcher et al. 2009).
The α-variable MILES models cover a range of total metallicity from [M/H] = −2.27 to [M/H] = +0.40,and from 30 Myr to 14 Gyr in age.All α-elements are varied simultaneously and model predictions are given at [α/Fe] = +0.0 and [α/Fe] = +0.4.To model possible IMF variations, we assume in this work a single-power-law IMF parametrization (the so-called unimodal IMF shape in the MILES notation; Vazdekis et al. 1996).Changes in the IMF are parametrized by Γ, the logarithmic slope of the IMF.For reference, a Salpeter IMF (Salpeter 1955) corresponds to Γ = 1.35.We note that other IMF functional forms such as log-normal approximation (e.g., Chabrier 2003) or a broken power law (e.g., Kroupa 2002) have no exact translation into Γ units.
It is worth mentioning that the single power-law IMF description adopted here departs from our previous measurements of the IMF in ETGs.Observations in the Milky Way (e.g., Kroupa 2001;Bastian et al. 2010) and theoretical arguments (e.g., Chabrier et al. 2014;Krumholz 2011Krumholz , 2014) ) support the presence of a characteristic stellar mass in the IMF and therefore strongly disfavor an IMF defined by a single slope across all masses.Moreover, mass-to-light ratio predictions for a single power-law IMF in massive ETGs are inconsistent with dynamical mass measurements (Ferreras et al. 2013;Cappellari et al. 2013).However, given our current lack of knowledge regarding the IMF behavior in young stellar populations beyond the Milky Way, a single power-law IMF parametrization greatly simplifies the interpretation of the obtained results.
Young systems with an extended star formation history (SFH) can, however, contain stellar populations younger than the limit of 30 Myr of the MILES α-variable models.This can be particularly relevant for the spatially resolved IFU data of the TIMER survey, because even if the average age of a galaxy is moderately high, star formation can be restricted to small regions where populations younger than 30 Myr can in fact significantly contribute to the observed spectra.To overcome this potential issue, our analysis complements the α-variable models with the super-young MILES models presented in Asa'd et al. (2017).Fed with the Bertelli et al. (1994) set of isochrones, these superyoung MILES models provide predictions down to 6.3 Myr.
In order to combine both α-variable and super-young MILES models into a regular grid of spectra, we linearly interpolated the latter to match the metallicity sampling of the α-variable models.Moreover, although the effect of a variable abundance pattern becomes weaker as the age of the stellar population decreases (see e.g., Figs. 9 and 10 in Vazdekis et al. 2015), we used the behavior of the α-variable MILES models to expand our treatment of the abundance pattern toward the super-young models.In practice, this is done by calculating the ratio of the [α/Fe] = +0.0 and [α/Fe] = +0.4models at a given metallicity for the youngest α-variable models.This response function is then applied to the super-young models in order to approximate the effect of a variable abundance ratio.We note that, as these super-young models are based on the MILES stars that follow the local [α/Fe]-[Fe/H] relation (see Milone et al. 2011;García Pérez et al. 2021), the correction of the superyoung models depends on the metallicity.Below [M/H] = −0.4,models have effectively [α/Fe] ∼ 0.4 and therefore we use the response functions to calculate model predictions at [α/Fe] ∼ 0.0.Conversely, above [M/H] = −0.4we correct the super-young models to reach [α/Fe] ∼ 0.4.
Briefly, after combining the α-variable and super-young MILES models, our set of single stellar population (SSP) model predictions cover ages from 6.3 Myr to 14 Gyr, with variable total metallicity, [α/Fe], and IMF, with this one parametrized as a single power law with logarithmic slope Γ.

Line-strength behavior
In order to measure the IMF from integrated spectra, we need to probe the contribution to the observed flux of stars with different masses.From a modeling point of view, this is done by mapping (at fixed chemical composition) effective temperature and surface gravity onto stellar mass across a given isochrone.IMF-sensitive features are therefore usually labeled as either temperature-sensitive or gravity-sensitive.However, we note that the net effect of a variable IMF on a given spectral feature is the result of a (nontrivial) combination of its sensitivity to both temperature and surface gravity, modulated by the shape of the isochrone, the intrinsic mass-luminosity relation of the stars, and the adopted IMF.The titanium dioxide molecular band TiO 2 at ∼6200 Å (Worthey et al. 1994) is a paradigmatic example of this complex IMF sensitivity, as for old stellar populations it becomes stronger when increasing the number of dwarf stars (i.e., with bottom-heavier IMFs) even though it is actually more prominent in the atmospheres of giants (see e.g., Fig. 1 in Spiniello et al. 2014).
Therefore, before further presenting our fitting approach, it is worth discussing the behavior of some of the most prominent spectral features within the MUSE wavelength range, particularly for young ages, where their IMF sensitivity remains mostly unexplored.Figure 2 shows how the Mgb, NaD, TiO 2 (Worthey et al. 1994), and H β O (Cervantes & Vazdekis 2009) line-strength indices change as a function of age for three different IMF slopes (at solar metallicity).Ages where the model predictions come from the super-young MILES models are indicated by a dashed line.
Arguably, the most characteristic feature in Fig. 2 is the predicted bump in all indices but H β O for ages of ∼10 Myr.Although this bump roughly coincides with the transition toward the superyoung models, the rapid change in the strengths of the indices is not an artifact of joining the two sets of models but a consequence of the onset of the red supergiant phase (e.g., Mayya 1997).This is shown by the fact that the gray-shaded areas in Fig. 2, which indicate the raw super-young MILES model predictions (i.e., without the metallicity and [α/Fe] corrections detailed above), closely track the behavior of the black lines, both calculated for the same Salpeter IMF.
The importance of including these super-young model predictions in our analysis becomes obvious from Fig. 2, particularly for the case of the TiO 2 feature.In the absence of these models, the expected contribution of red supergiants to the observed spectra would be severely underestimated and we could wrongly interpret a strong TiO 2 absorption as a consequence of a dwarf-rich (i.e., bottom-heavy) IMF, while in reality it results from the presence of very young stellar populations.It is also important to notice that ages in Fig. 2 are shown in logarithmic scale and therefore the actual dominance of red supergiants is effectively restricted to a very short phase peaking at ∼10 Myr.Our fitting procedure will take advantage of this rapid evolution (see Sect. 3.4).
The model predictions represented in Fig. 2 also highlight important differences in the IMF sensitivity of these indices depending on the age of the underlying stellar populations.Both TiO 2 and NaD absorption features have been extensively used to characterize the IMF in old, quiescent galaxies with a rather straightforward interpretation: stronger features translate into steeper IMF slopes.However, this does not hold for younger stellar populations where the opposite trend is expected.This reversal of the IMF dependence is clearly exemplified by the TiO 2 feature, where similarly strong values are predicted for either an old population and a dwarf-dominated IMF or a young population biased towards massive stars.
A final note regarding the effect of a variable IMF on the spectra of young stellar populations.Having young populations implies that more massive stars still contribute to the observed spectra.This has two immediate and profound consequences for the interpretation of the measurements.First, contrary to the old stellar populations found in quiescent galaxies, where only stars less massive than ∼1 M are present, IMF measurements in young populations are age-dependent, because the stellar-mass range probed along the IMF changes as populations become younger.Therefore, changes in the measured slope might be observed -for example, within the same galaxy -not because the IMF itself is varying, but because stellar populations with different ages probe different mass ranges across the IMF, each of them characterized by a different slope.
Second, and related to the point above, in old stellar populations, IMF measurements are limited to a narrow mass range, resulting in a lack of sensitivity to the IMF parametrization.This has greatly simplified the comparison between different observational approaches, and IMF variations are reliably discussed in terms of a single parameter (e.g., IMF slope Γ, dwarf-to-giant ratio F 0.5 , mismatch parameter α, etc.).The spectrum of a young stellar population, however, becomes progressively sensitive to the shape of the IMF, as demonstrated by Fig. 3, where the equivalent width of the IMF-sensitive TiO 2 feature is shown as a function of age for different IMF functional forms.
Figure 3 shows the predicted equivalent width of the TiO 2 as a function of age for three different IMFs and solar metallicity.As in Fig. 2, orange and black lines correspond to a bottomheavy single power law and a Salpeter IMF, respectively.In purple, we show the predictions for a bottom-heavy IMF (i.e., biased toward low-mass stars), but assuming a broken-power-law mass distribution (the so-called bimodal IMF within the MILES notation; Vazdekis et al. 1996).This bimodal IMF is designed to be a generalization of the Kroupa (2002) functional form, with a variable slope for the high-mass end Γ B and a fixed, flat distribution for low-mass stars (below ∼0.4 M ).For Γ B , the bimodal IMF is equivalent to the Milky Way standard.
For increased age, the predicted TiO 2 for the two bottomheavy IMF parametrizations are almost indistinguishable.This is ultimately due to the fact that in old populations, spectra only probe a very narrow stellar mass range, resulting in a lack of sensitivity to the IMF functional form (La Barbera et al. 2013).However, as populations get younger, stars with different masses start to contribute to the observed spectra and line-strength Fig. 3. Optical sensitivity to the shape of the IMF.Black, orange, and purple lines indicate the age sensitivity of the TiO 2 absorption feature for a Salpeter IMF, a bottom-heavy single-power-law IMF, and a bottom-heavy broken-power-law IMF, respectively.For old stellar populations, such as those harbored by massive quiescent galaxies, singlepower-law IMF parametrizations are indistinguishable from brokenpower-law IMF parametrizations because the stellar-mass range probed along the IMF is rather narrow.For young stellar populations, however, different model assumptions on the IMF shape lead to dramatically different model predictions, highlighting the potential use of line-strength indices in constraining the shape of the IMF.Predictions are shown at the native resolution of the MILES models (2.51 Å).
analyses become sensitive to the shape of the IMF.Finally, and in the absence of more flexible and physically motivated functional forms, the trends in Fig. 3 also justify our adoption of a singlepower-law IMF parametrization.Although we acknowledge that such a functional form is inconsistent with dynamical IMF measurements of massive quiescent galaxies (e.g., Ferreras et al. 2013;Lyubenova et al. 2016), it is relatively agnostic about the physics behind IMF variations and has a simple interpretation even for populations with different ages.

Limitations of the SSP assumption
The effect of the IMF on the integrated spectra of galaxies is subtle and can be partially mimicked by changes in other stellarpopulation parameters such as the chemical composition or the SFH.To overcome this problem, IMF measurements based on the analysis of absorption spectra have so far been limited to quiescent galaxies whose star-formation histories are reasonably well-represented by SSP models.The SSP assumption simplifies the fitting process by reducing the number of free parameters and isolating the effect of the IMF becomes feasible.
Although it has been central to our current understanding of the stellar populations in galaxies, the SSP assumption faces unavoidable problems when dealing with systems with more complex star-formation histories (e.g., Gallazzi et al. 2005;Serra & Trager 2007;Walcher et al. 2015).The limitations of the SSP approach are illustrated in Fig. 4. The two panels in this figure show a standard index-index diagram with an optimized age-sensitive feature on the vertical axis (H β O , Cervantes & Vazdekis 2009) and an IMF-sensitive feature on the horizontal axis (NaD, Worthey et al. 1994).Each corner of the grid corresponds to a MILES model prediction with a given IMF slope Γ and age, as indicated by the labels.Over-plotted on top of the index-grid, colored symbols indicate the expected H β O and NaD values after combining two SSP models of 0.35 and 14 Gyr with exactly the same Γ = 1.3 Milky Way-like IMF slope.Each of these symbols is color-coded by the fraction of light (left) or mass (right) associated to the young SSP model.
Two main conclusions can be directly drawn from Fig. 4. First, the IMF of composite stellar populations should not be studied under the SSP assumption.In the simplified case illustrated in Fig. 4, when the relative flux emitted by the old and young populations is similar, the best-fitting SSP solution for the IMF slope becomes close to Γ = 3.3, failing to recover the true value of Γ = 1.3.We note that the unreliability of the SSP assumption for composite stellar populations not only applies to IMF variations, but also invalidates measurements of more robust quantities such as metallicities and elemental abundance ratios (see e.g., Appendix B in Seidel et al. 2015).Second, the SSP approach works best for young stellar populations.This, at first, may seem counterintuitive because the concept of SSP has been traditionally used in the context of old, quiescent stellar populations.However, young stellar populations can easily outshine the presence of any old component, as shown in the right panel of Fig. 4. Modeling the flux emitted by young populations is clearly not free from uncertainties and systematic error (see e.g., Leitherer et al. 1999Leitherer et al. , 2014;;Eldridge et al. 2017), but given a set of SSP templates, the importance of accounting for an extended SFH is maximal when old and young stellar populations contribute similarly to the observed flux budget.
These limitations of the SSP approach can be directly tackled by allowing for an extended SFH when fitting the observed absorption spectra of galaxies.This is the foundation of the stellar-population-fitting scheme described in detail in the following section.

Fitting scheme
An ideal stellar population fitting algorithm would be able to retrieve the time evolution of every stellar population property from an integrated spectrum.This would include not only the SFH but also the time evolution of the chemical composition and potentially that of the IMF.However, and despite significant advances in measuring time-varying properties of galaxies from their integrated spectra (see e.g., Cid Fernandes et al. 2005;Ocvirk et al. 2006;Tojeiro et al. 2007;Cappellari 2017;Wilkinson et al. 2017;Johnson et al. 2021), it remains unclear to what extent all that information can be recovered, or whether or not it is even sufficiently encoded in absorption spectra (Ferreras et al. 2023).
As a first step forward, in this work we build upon our approach to measure the IMF described in Martín-Navarro et al. (2019, 2021).Starting with the Voronoi-binned MUSE data described above, our stellar population analysis consists of two steps: First we use the full spectral fitting code pPXF (Cappellari & Emsellem 2004) to measure the kinematics (V and σ) of every spectrum and to model and subtract the ionized gas emission from the observed spectra.Then, fixing the measured kinematics to avoid degeneracies with the stellar population parameters (Sánchez-Blázquez et al. 2011), we rerun pPXF, this time imposing a regularized solution to the best-fitting weights distribution (Cappellari 2017).As detailed in Martín-Navarro et al. (2019), the main goal of this second pPXF run is to estimate the SFH and thus the mean age of each Voronoi-binned spectra.There are three main advantages A110, page 5 of 15 to calculating the age in this way.First, we minimize the age-abundance pattern degeneracy that arises from the Hβ dependence on [C/Fe] (e.g., Conroy & van Dokkum 2012b; La Barbera et al. 2016) and the lack of strong C-sensitive features within the MUSE wavelength range.Second, we robustly measure the age by simultaneously fitting a wide spectral range, from 4700 Å to 6400 Å.Finally, and crucial to our approach, contrary to the standard line-strength analysis where only the zeroth-order moment of the age distribution is measured, by using pPXF we can approximate the full SFH of a given population.
For consistency with the second step explained below, we use pPXF to measure luminosity-weighted quantities.In practice, and related to what it is shown in Fig. 4, we calculate the light fraction contained in stars with different ages.We regularize the pPXF solution in the age-metallicity-IMF space, which minimizes the dependence of the recovered SFH on the assumed IMF.In order to estimate the uncertainty in the measured SFH, we repeat the analysis of each spectrum ten times, adding Gaussian noise to the observed spectra with an amplitude given by the residuals of the best-fitting model.The end result of this first step is therefore that, for each spectra, we measure its kinematics, we remove the gas emission, and we estimate the SFH and its uncertainty.This first part of the fitting scheme is the same as in Martín-Navarro et al. (2019, 2021) and Bittner et al. (2020) and has proven to be robust and efficient at providing detailed stellar population maps based on MUSE data.
During the second step, the goal is to measure the rest of the relevant stellar population parameters, including the slope of the IMF. Figure 4 illustrates the limitations of the standard SSP approach assumed so far in order to measure the IMF from integrated spectra.To overcome this problem, we implemented a more flexible fitting scheme that allows for an extended SFH while assuming a single chemical composition and IMF slope.In practice, the process is as follows.First, we use the SFH derived from pPXF to generate a grid of MILES predictions.We note that this is a fundamental change compared to previous studies, because the set of predictions is not SSP-like, but is based on the measured star formation histories.In this way, by construc-tion, we remove the bias shown in Fig. 4.This combination of an extended SFH and a single value characterizing the stellar population properties is the same as that implemented in most photometric inversion algorithms (e.g., da Cunha et al. 2008;Kriek et al. 2009;Han & Han 2019;Johnson et al. 2021) and is conceptually similar to the scheme followed by Gallazzi et al. (2005), but in our case is optimized to measure the IMF from integrated absorption spectra.
Given these SFH-based MILES model predictions, we then use a Bayesian (Foreman-Mackey et al. 2013) full index fitting (FIF) approach (Martín-Navarro et al. 2019) to measure mean luminosity-weighted stellar population properties.Briefly, this is done by fitting every pixel within the band pass of a set of spectral absorption features after normalizing their continuum.Summarizing the approach followed in this second step, Fig. 5 shows how SSPs with different ages have different equivalent widths and then contribute differently to the total flux of each spectral feature included in our analysis (see details below).The SFH-equivalent prediction is shown as a black line.
We select six prominent spectral features in the MUSE wavelength range, namely Mgb, Fe 5270, Fe 5335, NaD, TiO 1 , and TiO 2 .Modeling the band passes of each of these indices based on the SFH measured during the first step, our Bayesian FIF scheme fits for the following (luminosity-weighted) free parameters: total metallicity [M/H], [Mg/Fe], IMF slope, [Ti/Fe], and [Na/Fe].We note that magnesium is effectively the only α element constrained by our set of indices through the Mgb feature.Therefore, even if MILES models are α-variable, we refer to [Mg/Fe] when describing our results.Moreover, neither [Ti/Fe] nor [Na/Fe] abundances are self-consistently included in the MILES models.In particular, Ti is an α element but does not necessarily track Mg variations (see e.g., Johansson et al. 2012;Conroy et al. 2014;Tang & Worthey 2017).In practice, the inclusion of a variable [Ti/Fe] allows us to approximate the effect of a variable abundance pattern in the behavior of TiO molecular bands as a nuisance parameter rather than constraining the actual [Mg/Fe] value.To account for the impact of [Ti/Fe] and [Na/Fe], we use the response functions of Conroy & van Dokkum (2012b).These response functions are A110, page 6 of 15 calculated for a 13 Gyr stellar population, and as the sensitivity to elemental abundance variations becomes weaker for young populations, the resulting [Na/Fe] and [Ti/Fe] values are not reliable in absolute terms and are not further discussed in this work.Motivated by the behavior of the indices in Fig. 2, we include in our fitting scheme an additional free parameter that freely controls the fraction of flux emitted by a 10 Myr old population.In this way, we can better model the critical contribution of supergiant stars.As noted above, the fact that this phase rapidly peaks at 10 Myr allows us to model it as a pure additional SSP contribution.Finally, we also include a rescaling term in the estimated error in the MUSE data.Assuming uniform prior distributions for all free parameters, we repeat the fitting process ten times, one for each of the measured SFHs (see above).We then combine all the posterior distributions to estimate the best-fitting values and their uncertainties.With all this, our stellar population analysis fits in total for seven free parameters: [M/H], [Mg/Fe], Γ, [Na/Fe], [Ti/Fe], F RSG , and ∆ err .
Figure 6 exemplifies the result of our fitting approach for one of the Voronoi-binned MUSE spectra of NGC 3351.We note that the spectrum shown in Fig. 6 has a relatively young (luminosityweighted) age as a result of an extended SFH and our fitting scheme is able to reproduce the observed strength of all five indices in a regime where the SSP approach no longer holds.The reliability of our approach is summarized in Fig. 7, where the full posterior distribution of metallicity, [Mg/Fe], and IMF slope values is shown for the same spectrum as in Fig. 6.Even in this non-trivial case with a young and extended SFH, our fitting scheme is able to assess the intrinsic degeneracies of the inversion problem.A discussion on potential biases and systematic errors on the recovered stellar population properties, in particular on the IMF values, is included in Appendix A.

Results
Figure 8 shows the measured age, metallicity, [Mg/Fe], and IMF slope maps for NGC 3351.The age map corresponds to the mean luminosity-weighted values measured during the first step of the fitting process as described above and clearly reveals the star-forming nuclear ring of NGC 3351 (see Leaman et al. 2019;Bittner et al. 2020).Towards the outskirts, the stellar populations of NGC 3351 become predominantly old.The prominent dust lines shown in Fig. 1 do not have a clear impact on the age map, which is in agreement with the lack of (dust-corrected) Hα emission (Bittner et al. 2020).Such dust lanes are known to lack star formation, which is possibly due to the effects of shear on the molecular clouds (e.g., Athanassoula 1992;Kim & Stone 2012;Emsellem et al. 2015;Neumann et al. 2019).
The metallicity map is characterized by the presence of the metal-enriched nuclear disk.Almost perpendicularly, the bar of NGC 3351 is clearly visible as an elongated structure that extends to the edge of the MUSE field of view.In addition, the metallicity map reveals the presence of what appear to be scattered metal-poor regions along the star-forming nuclear ring (Bittner et al. 2020;Pessa et al. 2023).These are starforming regions with varying mass, dust attenuation level, and A110, page 7 of 15 The bottom right panel of Fig. 8 shows the IMF slope map of NGC 3351.Reassuringly, the measured IMF slope values are close to the Milky Way standard (Γ = 1.3), with a slight trend following changes in local metallicity.In this regard, the metalrich stellar population of the bar (Neumann et al. 2020) is characterized by relatively steep IMF slopes.This relation between IMF slope and metallicity is similar to what has already been reported in massive ETGs (e.g., Martín-Navarro et al. 2015a, 2021;Parikh et al. 2018;Sarzi et al. 2018;Zhou et al. 2019).On top of that, the edges of the nuclear star-forming ring seem to exhibit relatively steep IMF slopes, in particular in the area around the youngest star forming knots.We note that these bins with a steep IMF slope do not perfectly correspond to the CP regions, suggesting that they are not the result of a trivial systematic effect.Finally, it is worth emphasizing that the steps described above in order to account for an extended SFH and the presence of young stars, in particular supergiants, are mandatory in order to obtain reliable IMF measurements.To exemplify the importance of these steps, in Appendix B we show the (incorrect) IMF slope map resulting from an analysis of the spectra of NGC 3351 following the same approach as for quiescent ETGs.
Figure 9 shows how the measured IMF slopes change as a function of metallicity (top panel), [Mg/Fe] (middle panel), and stellar surface mass density (bottom panel).The trend with metallicity already suggested by the maps in Fig. 8 becomes evident in the top panel, and is consistent with IMF measurements in ETGs (Martín-Navarro et al. 2015b;Parikh et al. 2018;Sarzi et al. 2018;Zhou et al. 2019).Although the scatter increases as age decreases, IMF variations seem to be dependent on the local metallicity over the whole age range probed in our maps.Moreover, these metal-rich and therefore chemically evolved regions also tend to exhibit lower [Mg/Fe] values and thus a subtle relation between Γ and [Mg/Fe] emerges in the middle panel of Fig. 9. Finally, the measured Γ values correlate with the local stellar mass density, which is not surprising given the fact that the evolution of both quantities is expected to be tightly connected (e.g., Sánchez et al. 2017;Zhuang et al. 2019).A correlation between stellar mass density and IMF slope has also been reported in ETGs (La Barbera et al. 2019;van Dokkum et al. 2017;Martín-Navarro et al. 2021).

Is the IMF universal in LTGs?
For quiescent galaxies, with typical ages of around 10 Gyr, only stars with masses of ∼0.5 M contribute to the observed spectra.Therefore, the observational evidence supporting the nonuniversality of the IMF in these systems only applies to its low-mass end (e.g., Conroy et al. 2017).On the other hand, IMF measurements of resolved stellar populations are able to constrain the shape of the IMF across a wide range of stellar masses and this is commonly done by characterizing the slope of the IMF Γ as a function of stellar mass.This is illustrated in Fig. 10 where red filled squares are measurements compiled in Scalo (1998), yellow squares are globular clusters from Piotto & Zoccali (1999), and purple squares correspond to the bulge of the Milky Way by Holtzman et al. (1998) and Zoccali et al. (2000).Horizontal green dashed lines in Fig. 10 indicate the average IMF slope of the solar neighborhood as measured by Kroupa (2002), while the shaded green area is the associated uncertainty.These measurements in green are commonly adopted as the Milky Way standard.
Bridging the gap between resolved and unresolved IMF measurements, colored circles in Fig. 10 represent the best-fitting IMF slope of each Voroni bin in NGC 3351 as a function of the mean stellar mass of the underlying stellar population.We note that the horizontal axis in Fig. 10 does not correspond to the surface stellar mass density of the different bins but to the mean mass of the stars still alive and therefore still contributing to the observed spectra.In practice, m was calculated as the age-dependent arithmetic mean between the maximum and minimum stellar mass in the BaSTI set of isochrones.For old stellar populations, this mean m is roughly equal to ∼0.5 M , which is within the mass range explored by IMF studies in quiescent galaxies.The young populations within NGC 3351 allow us to expand IMF measurements from integrated spectra toward more massive stars.This approximation to the mean stellar mass is independent of the adopted IMF parametrization but it does depend on the exact combination of spectral indices used in the analysis.Attempting a more accurate characterization of the A110, page 8 of 15 exact mass range probed by each absorption spectrum is beyond the scope of this paper, but will be addressed in the future.
The agreement between our IMF measurements and those based on resolved stars in the Milky Way is remarkable.Not only is the average value consistent with observations of the Milky Way but the observed scatter also matches the dispersion around the solar neighborhood.In addition, and mimicking the behavior observed in the Milky Way, we recover flatter IMF slope values when probing lower mass (i.e., older) stars.We note that in our case there is a limit in the minimum stellar mass we can measure with our approach, as the average stellar mass of any population older than ∼10 Gyr remains approximately constant at ∼0.5 M (e.g., La Barbera et al. 2013).Nevertheless, and despite these important methodological differences, our IMF slope measurements closely follow expectations for the Milky Way.
Going beyond the comparison with the Milky Way, the colorcoding of our measurements shown in Fig. 10 quantifies the systematic behavior of the IMF map of NGC 3351, as the measured slope correlates with the surface stellar mass density of each bin.As noted above, IMF variations in NGC 3351 are not purely stochastic around a mean value of Γ = 1.3 but appear spatially coherent, tracing some of the morphological structures of the galaxy and changes in the stellar population properties (Fig. 9), in agreement with IMF measurements in massive ETGs as noted above (e.g., Martín-Navarro et al. 2015b, 2021;van Dokkum et al. 2017;La Barbera et al. 2019).
Finally, it is worth reflecting on some of the implicit assumptions and limitations of our approach.In particular, characterizing the slope of the IMF with a single power law leads to well-known issues with the predicted stellar mass-to-light ratio because of the unrealistically high contribution of very low-mass stars and stellar remnants to the mass budget M/L (see e.g., Ferreras et al. 2013).Therefore, IMF slope measurements presented here should not be naïvely translated into M/L values.Moreover, the interpretation of the measured Γ as a proxy for the local shape of the IMF (i.e., for the slope of the IMF over a range of stellar masses) becomes less robust if there is a change in the underlying IMF slope over the explored mass range.This could be particularly problematic for young stellar populations where stars with very different masses contribute to the observed spectrum.To assess the dependence of our results on the adopted IMF slope, we repeated our stellar population analysis using a broken-power-law IMF parametrization.The results of this test are included in Appendix C and are very similar the trends shown in Fig. 10.Nevertheless, neither of the two IMF functional forms is realistic enough and more flexible forms should be implemented in the future.measured slope of a single-power-law IMF, simply because older populations probe lower stellar masses of the IMF.Despite these caveats and until newer generations of stellar population models and fitting approaches are available, the adoption of a singlepower-law IMF parametrization greatly simplifies the analysis and interpretation of the observed spectra.In conclusion, our inference relies on the assumption that the IMF is fully sampled.
For young stellar populations in spatially resolved nearby galax-  1998) and Zoccali et al. (2000).Green dashed lines and shaded areas are the solar neighborhood standard measured by Kroupa (2002).Our IMF slope measurements, color-coded by the local surface stellar mass density, follow a behavior similar to the resolved IMF studies, although the scatter in our measurements is not purely stochastic and a systematic dependence on the stellar density is found.
ies, an incomplete sampling of the IMF could result in important biases.In our case, this would only apply to a handful of Voroni bins and therefore our results are robust, but the importance of this effect cannot be overlooked in younger nearby galaxies.

Chemical enrichment of the star-forming ring
As described above, the star-forming ring of NGC 3351 seems to contain a series of metal-poor regions.Emsellem et al. (2022) argued that the chemical enrichment of such low-metallicity regions would be inconsistent with the mixing timescales of the interstellar medium (e.g., Ho et al. 2017), andPessa et al. (2023) performed a thorough analysis of the possible model systematic errors that could be biasing the fitting of these young stellar populations.Our stellar population approach based on the analysis of individual absorption features is different from the full spectral fitting scheme followed by Emsellem et al. (2022) and Pessa et al. (2023) and yet leads to similar low-metallicity bestfitting values.
Our extension of the α-enhanced MILES models toward young ages allows us not only to measure the total metallicity content of these regions, but also to be sensitive to variations in the abundance ratio.In Sect. 4 we note that the measured [Mg/Fe] is not the same across all the star-forming and metalpoor regions, as there are four high-[Mg/Fe] knots clearly visible in the map of NGC 3351.To understand the origin of these four CP regions, Fig. 11  the youngest star-forming regions can dilute the strength of the observed absorption features, which would in turn lead to unrealistically low metallicities.Significant improvements have indeed been achieved over recent years by simultaneously and consistently fitting the stellar absorption and gas emission spectra of galaxies (e.g., Gomes & Papaderos 2017;Cardoso et al. 2019Cardoso et al. , 2022)).We tested this hypothesis by repeating our stellar population analysis but including the presence of a nonstellar continuum as an additional free parameter.As expected, this test resulted in more metal-rich populations once a possible additive continuum is included, but still the recovered metallicity values remain low.Moreover, the presence of such nebular continuum cannot easily explain the relative difference between the strength of Fe and Mg lines (i.e., the apparent [Mg/Fe] enhancement).
The fact that there is a clear decoupling between the metallicity and the [Mg/Fe] behavior across the star-forming ring of NGC 3351 may hint towards an alternative, physically motivated explanation.While magnesium is rapidly produced in core-collapse supernovae, the release of iron-peak elements is delayed as Type Ia supernovae are the main polluters (e.g., Kobayashi et al. 2020).The differences between [Mg/Fe] and metallicity maps could therefore be understood as different stages of the chemical evolution process, where the CP regions are just too young for Type Ia supernovae to have enriched the interstellar medium with iron.In this scenario, the metalpoor but [Mg/Fe] ∼ 0 knots would be more chemically evolved regions, although still powered by relatively metal-poor gas (e.g., Maiolino & Mannucci 2019).
The comparison with the photometric analysis of Calzetti et al. (2021) sheds further light on the chemical enrichment of the star-forming ring of NGC 3351.The four CP regions visible in our [Mg/Fe] plane correspond perfectly to those regions where the multiband photometric analysis reveals the presence of very young ( 10 Myr) star-formation bursts, which are presently too young for Type Ia supernovae to have exploded.In agreement with these findings, our metal-poor regions with solar [Mg/Fe] values are located in regions with a more extended SFH over the last ∼300 Myr.These differences between the photometric star-formation timescales of the different regions, although calculated completely independently, coincide with the chemical timescales derived from our stellar population analysis.
Finally, it is worth noting that, although the luminosityweighted age of the CP regions is slightly older than the age of the metal-poor and [Mg/Fe] ∼ 0 knots (see Sect. 4), this does not invalidate the proposed scenario.On the contrary, according to our measured SFHs, up to 55% of the flux of the CP regions is contributed by stars younger than 10 Myr, while for the metalpoor and [Mg/Fe] ∼ 0, this percentage drops to 35%.It is precisely the extended SFH of these long-lasting bursts that leads to younger mean ages on average.Our results are in notable agreement with the findings of Calzetti et al. (2021) and highlight the limitations of using mean values to characterize the age of stellar populations with complex and extended SFHs.
As a cautionary note, however, our results regarding the potential chemical peculiarities of the star-forming knots of NGC 3351, despite showing encouraging consistency with the independent analysis of Calzetti et al. (2021), are sensitive to significant model systematic uncertainties and limitations.In particular, the sensitivity of integrated spectra to variations in the chemical composition, and in particular to the abundance pattern of stars becomes progressively weaker in younger stellar populations.Moreover, from a modeling perspective, the theoretical corrections of Coelho et al. (2005) applied to the Vazdekis et al. (2015) SSP models are limited to stellar atmospheres cooler than 7000 K, which is insufficient to properly cover the whole temperature range explored in this work.A more reliable modeling of the chemical composition of young stellar populations (e.g., Percival et al. 2009) would therefore be a desirable next step, combined with the analysis of the gas phase chemistry (e.g., Díaz et al. 2007;Moustakas et al. 2010).

Summary and conclusions
Making use of the tools and ideas that have revealed the nonuniversality of the IMF in massive ETGs, we developed a stellar population fitting approach capable of analyzing the absorption spectra of galaxies with young populations and complex SFHs.We measured the spatially resolved age, metallicity, [Mg/Fe], and IMF slope Γ maps of the LTG NGC 3351, which is similar to the Milky Way in terms of mass and is characterized by the presence of a nuclear disk and star-forming ring.Our main results can be summarized as follows: - chemical signature of core-collapse supernova enrichment, which is in agreement with (independent) spectroscopic and photometric age determinations.In summary, our results highlight the feasibility of measuring the properties of young stellar populations in detail based on their integrated absorption spectra.When analyzed in the context of understanding the origin of the observed variations in the IMF, these young populations are a natural bridge between the universality of the Milky Way and the systematic changes measured in the centers of massive ETGs.In this regard, the IMF map of NGC 3351 allows a fine and complete characterization of the IMF across the inner regions of the galaxy, contrary to the relatively scarce resolved IMF measurements.It is only because of our large number of independent IMF measurements that we are able to detect the subtle IMF variations that appear to be mostly related to its low-mass end.Although observations in the Milky Way seem to favor some sort of universal IMF shape, it is possible that our ability to detect systematic deviations is hampered by the limited number of measurements.
Looking ahead, a detailed modeling of the absorption spectra of young stellar populations will play a capital role in understanding the high-redshift Universe in the JWST and ELT era.In order to maximize the amount of information extracted from these observations, models and analysis tools have to become flexible enough to account for the complex and potentially unfamiliar behavior of the young Universe.The methodology presented here is a step in this direction, but several aspects need to be improved, from a better treatment of noncanonical ingredients to a self-consistent measurement of the time evolution of the stellar population properties that takes into account the potentially time-varying nature of the IMF.

Fig. 1 .
Fig. 1.HST color image of NGC 3351.The star-forming ring is clearly visible in this color-composite image in the inner regions of NGC 3351, while the presence of a long bar also becomes clear as an elongated structure that runs diagonally across the HST image, extending beyond the MUSE FOV indicated with a white square (1 × 1 ).

Fig. 4 .
Fig. 4. H β O vs. NaD index-index diagram.Gray grids show the MILES predictions at the native resolution of the models (2.51 Å) and solar metallicity for different ages and IMF slopes, as indicated by the labels.Blue solid lines highlight the age variation for a Salpeter IMF.Dashed blue lines and filled circles track the range of expected values for a composite stellar population of two SSP models (14 Gyr and 0.35 Gyr), both of them with a Milky Way-like IMF.Symbols are color-coded according to the light (left panel) and mass (right panel) fraction associated to the young SSP model.This simplified yet realistic example shows how the best-fitting SSP solution fails to recover the true value of Γ = 1.3 when applied to composite stellar populations, which could result in unrealistically high IMF slope values (Γ ∼ 3.3).

Fig. 5 .Fig. 6 .
Fig. 5. Spectral dependence on the SFH.Colored lines illustrate the age dependence of the six spectral features included in our analysis.The solid black line shows the SFH-equivalent prediction resulting from weighting each SSP model according to the measured SFH.We note that the continuum of each model has been normalized to better represent the flux of the different SSPs.However, during the fitting process the normalization is only applied to the SFH-equivalent model prediction (i.e., the black line).

Fig. 7 .
Fig.7.Full posterior distribution for the three main stellar population parameters measured during the second step, namely total metallicity, [Mg/Fe], and IMF slope Γ, for the observed spectrum shown in Fig.6.The green solid lines indicate the median of the distribution, while dashed vertical lines mark the 16th and 84th percentiles.Our fitting approach is able to capture and successfully break the degeneracies of the inversion problem even when fitting stellar populations characterized by an extended SFH.

Fig. 8 .
Fig. 8. Best-fitting stellar population maps for NGC 3351.The top left panel shows the age map, with a clearly visible nuclear star-forming ring characterized by the presence of young stellar populations.In the top right panel, the total metallicity map of NGC 3351 shows the metal enrichment of both the nuclear disk and the bar.Regions of apparently low metallicity are visible across the star-forming ring.The [Mg/Fe] map (bottom left) exhibits a behavior similar (but inverted) to the metallicity map, with the nuclear disk and the bar appearing as low [Mg/Fe] components.Finally, the bottom right panel shows the IMF slope map where the typical values are close to what is observed in the Milky Way (Γ = 1.3).IMF variations, although mild, partially track changes in both metallicity and stellar density.For reference, a horizontal white line of 10 or 0.5 kpc at the distance of NGC 3351 is shown in each panel.

Fig. 9 .
Fig. 9. Measured slope of the IMF as a function of metallicity (top panel), [Mg/Fe] (middle panel), and stellar surface mass density (bottom panel).In each panel, symbols are color coded according to age (top), stellar surface mass density (middle), and total metallicity (bottom).Measured IMF values appear to correlate with the three quantities, with the IMF becoming increasingly steep for more metal-rich, less [Mg/Fe]-enhanced, and denser regions in NGC 3351.

Fig. 10 .
Fig. 10.Measured slope of the IMF as a function of the (logarithmic) stellar mass for a compilation of Milky Way measurements (filled squares) and for our measurements of NGC 3351 (colored circles).Red squares are measurements from Scalo (1998), purple ones are globular cluster measurements from Piotto & Zoccali (1999), and finally yellow symbols are IMF measurements of the Milky Way bulge presented in Holtzman et al. (1998) andZoccali et al. (2000).Green dashed lines and shaded areas are the solar neighborhood standard measured byKroupa (2002).Our IMF slope measurements, color-coded by the local surface stellar mass density, follow a behavior similar to the resolved IMF studies, although the scatter in our measurements is not purely stochastic and a systematic dependence on the stellar density is found.

Fig. 11 .
Fig. 11.Star-forming ring spectra.The average spectrum of the CP regions of NGC 3351 is shown in blue, along with the average spectrum of the star-forming ring in black.The depth of the metallicitysensitive absorption features (indicated with gray shaded regions) is clearly weaker in the CP spectrum.These relatively weak absorption features lead to the observed low metallicity and high [Mg/Fe] values.
The nuclear and morphological structures of NGC 3351 are clearly visible on its stellar population maps.The nuclear disk stands as a relatively young and metal-rich component, with a [Mg/Fe] ∼ 0. Similarly, the bar of NGC 3351 hosts metal-enriched and low [Mg/Fe] stars, a clear indication of chemically evolved populations.-The measured values of the IMF slope in NGC 3351 are consistent with an underlying Milky Way-like IMF.However, embedded in this universal behavior, the slope of the IMF changes in a systematic and spatially coherent manner.In particular, the bar seems to be characterized by a relatively steep IMF slope compared to the surrounding inner parts of the disk.Moreover, the outskirts of the star-forming nuclear ring show evidence of steep IMF values as well.-Regions with peculiar chemical properties are measured within the nuclear star-forming ring.These regions, associated with star-forming bursts, are characterized by the presence of metal-poor populations.Interestingly, the measured [Mg/Fe] values over these regions show a bimodal distribution, with either [Mg/Fe] ∼ 0 or [Mg/Fe] ∼ 0.3.We interpret this decoupling between total metallicity and [Mg/Fe] as the A110, page 11 of 15 Fig. 2. Line strength, age, and IMF sensitivity.From left to right and top to bottom, we show the age dependence of the Mgb, H β O , Nad, and TiO 2 spectral features, respectively.Solid lines correspond to the MILES α-variable models and dashed ones indicate their (adapted) super-young extension.The gray shaded area in the background shows the unmodified super-young model predictions for a Milky Way-like IMF (see main text for more details).Orange lines are predictions for a bottom-heavy IMF (i.e., steeper than the Milky Way standard), blue lines for a top-heavy (i.e., flatter) IMF, and black lines for the reference Salpeter IMF.The rapid increase in the equivalent width of the indices (but H β O ) at ages of ∼10 Myr is caused by the sudden appearance of red supergiant stars.Predictions are shown at the native resolution of the MILES models (2.51 Å).