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

A&A, 684, A110 (2024) https://doi.org/10.1051/0004-6361/202348060 c The Authors 2024 Astronomy & Astrophysics The universal variability of the stellar initial mass function probed by the TIMER survey Ignacio Martín-Navarro1,2 , Adriana de Lorenzo-Cáceres1,2, Dimitri A. Gadotti3,4, Jairo Méndez-Abreu1,2 , Jesús Falcón-Barroso1,2 , Patricia Sánchez-Blázquez5, Paula Coelho6, Justus Neumann7, Glenn van de Ven8, and Isabel Pérez9,10 1Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, 38205 La Laguna, Tenerife, Spain e-mail: [email protected] 2Departamento de Astrofísica, Universidad de La Laguna, 38200 La Laguna, Tenerife, Spain 3European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching bei München, Germany 4Centre for Extragalactic Astronomy, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK 5Departamento de Física de la Tierra y Astrofísica & IPARCOS, UCM, 28040 Madrid, Spain 6Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Universidade de São Paulo, R. do Matão 1226, São Paulo, Brazil 7Max-Planck Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany 8Department of Astrophysics, University of Vienna, Turkenschanzstrasse 17, 1180 Wien, Austria 9Departamento de Física Teórica y del Cosmos, Campus de Fuente Nueva, Edificio Mecenas, Universidad de Granada, 18071 Granada, Spain 10 Instituto Carlos I de Física Teórica y Computacional, Facultad de Ciencias, 18071 Granada, Spain Received 25 September 2023 /Accepted 11 December 2023 ABSTRACT 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. However, the fundamental methodological and physical differences between the two approaches have hampered a comprehensive understanding of IMF variations. Here, we describe an improved modeling scheme that, for the first time, allows consistent IMF measurements across stellar populations with different ages and complex star-formation histories (SFHs). 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 starforming regions that appear to be predominantly enriched by the stellar ejecta of core-collapse supernovae. Our findings therefore showcase the potential of detailed studies of young stellar populations to provide the means 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. Key words. Galaxy: formation – galaxies: evolution – galaxies: fundamental parameters – galaxies: star formation – galaxies: stellar content 1. 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.5Mand a shallower distribution for lower-mass stars (e.g., Miller & Scalo 1979;Kroupa 2001,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). Beyond the limited range of star-formation conditions probed by the Milky Way environment, early-type galaxies (ETGs) are ideal systems to assess the universality of the IMF. Their old stellar populations are relatively simple (e.g., Faber 1973;Tinsley & Gunn 1976;Vazdekis et al. 1997; Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. This article is published in open access under the Subscribe to Open model.Subscribe to A&A to support open access publication. A110, page 1 of 15 Martín-Navarro, I., et al.: A&A, 684, A110 (2024) Trager et al. 2000;Thomas et al. 2005) and their high luminosities and well-behaved radial light distributions allow precise analyses of their internal dynamics (e.g., Jeans 1922; Schwarzschild 1979;van den Bosch et al. 2008;Cappellari 2008;Zhu et al. 2018). Both stellar population (e.g., van Dokkum & Conroy 2010;Spiniello et al. 2012,2014,2015; Conroy & van Dokkum 2012a;Smith et al. 2012;Ferreras et al. 2013;La Barbera et al. 2013;Tang & Worthey 2017; Lagattuta et al. 2017;Rosani et al. 2018;Lonoce et al. 2021, 2023) and dynamics-based (e.g., Treu et al. 2010;Thomas et al. 2011;Auger et al. 2010;Cappellari et al. 2012;Dutton et al. 2012;Wegner et al. 2012;Tortora et al. 2013,2014;Läsker et al. 2013;Posacki et al. 2015;Li et al. 2017;Corsini et al. 2017; Sonnenfeld et al. 2019) IMF measurements in massive ETGs systematically point toward an excess of low-mass stars in these objects compared to expectations when a Milky Way-like IMF is assumed. While there are systematic differences between the two approaches (Smith & Lucey 2013;Smith 2014;Lyubenova et al. 2016;Davis & McDermid 2017;Martín-Navarro et al. 2019; Lu et al. 2023), these results suggest that the extreme condition under which stars in massive ETGs formed altered the fragmentation of molecular clouds, leading to the observed IMF variations. 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 ∼1M, 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), NGC3351, 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. 2presents the TIMER data of NGC3351. In Sect. 3we give a detailed explanation of our stellar-population modeling approach. The results of our analysis of NGC3351 are presented in Sect. 4and discussed in Sect. 5. Finally, Sect. 6summarizes our findings and lists our concluding remarks. 2. Data The TIMER project presented in Gadotti et al. (2019) is a MUSE survey of 24 nearby barred galaxies. The exquisite 10h44m04s00s43m56s 11°44' 43' 42' R.A. Dec MUSE FOV Fig. 1. HST color image of NGC3351. The star-forming ring is clearly visible in this color-composite image in the inner regions of NGC3351, 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 (10×10). 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×1010 M(Muñoz-Mateos et al. 2015) and a distance of 10.1Mpc1. Figure 1shows a HST color image of NGC3351 (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 NGC3351 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 ∼10pc/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 1At this distance, the pixel scale of the MUSE spectrograph translates into a spatial scale of ∼10pc/spaxel. 2Program ID 15654, PI: Janice Lee. A110, page 2 of 15 Martín-Navarro, I., et al.: A&A, 684, A110 (2024) 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. 3. Stellar population modeling 3.1. 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,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 30Myr to 14Gyr 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 2011,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 30Myr 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 30Myr 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.4 models 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.4 we 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.3Myr to 14 Gyr, with variable total metallicity, [α/Fe], and IMF, with this one parametrized as a single power law with logarithmic slope Γ. 3.2. 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 TiO2at ∼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 2shows how the Mgb, NaD, TiO2 (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. 2is the predicted bump in all indices but HβOfor ages of ∼10Myr. 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 A110, page 3 of 15 Martín-Navarro, I., et al.: A&A, 684, A110 (2024) 0.00 1.00 2.00 3.00 4.00 Mgb [Å] Bottom-heavy IMF ( =2.8) Top-heavy IMF ( =0.3) Salpeter IMF ( =1.3) MILES super young ( =1.3) 3.00 4.00 5.00 6.00 7.00 8.00 HO [Å] 0.01 0.10 1.00 10.00 Age [Gyr] 1.00 2.00 3.00 4.00 5.00 NaD [Å] 0.01 0.10 1.00 10.00 Age [Gyr] 0.03 0.05 0.08 0.10 0.12 0.15 TiO2 [mag] 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 TiO2spectral 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 ∼10Myr is caused by the sudden appearance of red supergiant stars. Predictions are shown at the native resolution of the MILES models (2.51Å). 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 TiO2feature. 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 TiO2absorption 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. 2are shown in logarithmic scale and therefore the actual dominance of red supergiants is effectively restricted to a very short phase peaking at ∼10Myr. Our fitting procedure will take advantage of this rapid evolution (see Sect. 3.4). The model predictions represented in Fig. 2also highlight important differences in the IMF sensitivity of these indices depending on the age of the underlying stellar populations. Both TiO2and 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 TiO2 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 ∼1M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 F0.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 TiO2feature is shown as a function of age for different IMF functional forms. Figure 3shows the predicted equivalent width of the TiO2as 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 ΓBand a fixed, flat distribution for low-mass stars (below ∼0.4M). For ΓB, the bimodal IMF is equivalent to the Milky Way standard. For increased age, the predicted TiO2for 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 A110, page 4 of 15 Martín-Navarro, I., et al.: A&A, 684, A110 (2024) 0.01 0.10 1.00 10.00 Age [Gyr] 0.02 0.04 0.06 0.08 0.10 TiO2[mag] Broken power-law IMF ( B=3.3) Single power-law IMF ( =2.3) Salpeter IMF ( =1.3) Fig. 3. Optical sensitivity to the shape of the IMF. Black, orange, and purple lines indicate the age sensitivity of the TiO2absorption 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. 3also 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. 3.3. 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βOand NaD values after combining two SSP models of 0.35 and 14Gyr 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. 1999,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. 3.4. 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 Martín-Navarro, I., et al.: A&A, 684, A110 (2024) 123456 NaD [Å] 3 4 5 6 7 8 HO [Å] =1.3 =3.3 0.35 Gyr 1 Gyr 10 Gyr Light-weighted 0.0 0.2 0.4 0.6 0.8 1.0 Young light fraction 123456 NaD [Å] 3 4 5 6 7 8 HO [Å] =1.3 =3.3 0.35 Gyr 1 Gyr 10 Gyr Mass-weighted 0.0 0.2 0.4 0.6 0.8 1.0 Young mass fraction Fig. 4. HβOvs. 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 (14Gyr 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). 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 4illustrates 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 construction, 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, Fe5270, Fe5335, NaD, TiO1, and TiO2. 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 Martín-Navarro, I., et al.: A&A, 684, A110 (2024) 516 518 [nm] 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 Normalized flux Mgb 525 528 [nm] Fe 5270 532 535 [nm] Fe 5335 588 590 [nm] Na D 595 598 [nm] TiO1 620 625 [nm] TiO2 5 10 Age [Gyr] SFH equivalent 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). 516 518 [nm] 0.80 0.85 0.90 0.95 1.00 1.05 Normalized flux Mgb Data Model 525 528 [nm] Fe5270 532 535 [nm] Fe5335 588 590 [nm] NaD 595 598 [nm] TiO1 620 625 [nm] Age = 4.5 Gyr TiO2 Fig. 6. Data (in black) and best-fitting models (green) for one of the Voronoibinned data spectra of NGC3351. Each panel shows the central band pass of the five spectral features used in our analysis: Mgb, Fe5270, Fe 5335, NaD, TiO1, and TiO2. Each green line shows a best-fitting model resulting from sampling the full posterior distribution, taking into account the uncertainty in the recovered SFH. The mean luminosityweighted age is indicated in the rightmost panel (TiO2), although the prediction for the index strength is calculated based on the measured SFH (see text for more details). calculated for a 13Gyr 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 10Myr 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 10Myr 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], FRSG, and ∆err. Figure 6exemplifies 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. 6has 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. 4. Results Figure 8shows the measured age, metallicity, [Mg/Fe], and IMF slope maps for NGC3351. 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 NGC3351 (see Leaman et al. 2019; Bittner et al. 2020). Towards the outskirts, the stellar populations of NGC3351 become predominantly old. The prominent dust lines shown in Fig. 1do 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 NGC3351 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 Martín-Navarro, I., et al.: A&A, 684, A110 (2024) [M/H] = 0.16+0.03 0.03 0.03 0.06 0.09 0.12 [Mg/Fe] [Mg/Fe] = 0.08+0.02 0.02 0.10 0.15 0.20 0.25 [M/H] 0.6 1.2 1.8 2.4 0.03 0.06 0.09 0.12 [Mg/Fe] 0.6 1.2 1.8 2.4 = 1.62+0.47 0.54 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. formation history as suggested by their photometric properties (Calzetti et al. 2021). Interestingly, four [Mg/Fe]-enhanced ([Mg/Fe]>0.22) regions are clearly visible in the map of NGC3351 shown in Fig. 8. These chemically peculiar (CP) regions, seemingly characterized by metal-poor and [Mg/Fe]- enhanced stellar populations, correspond to the youngest starforming knots identified by Calzetti et al. (2021), with ages of between ∼5 and ∼10Myr. Beyond these CP regions, the [Mg/Fe] map shows rather smooth variations, where [Mg/Fe] values tend to decrease towards the center of NGC3351 and along the direction of the bar. This coupling between the behavior of the [Mg/Fe] and the total metallicity maps indicates a long-lasting chemical enrichment and thus star formation activity in both the nuclear stellar disk and the bar of NGC3351. The bottom right panel of Fig. 8shows the IMF slope map of NGC3351. 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 Bwe show the (incorrect) IMF slope map resulting from an analysis of the spectra of NGC3351 following the same approach as for quiescent ETGs. Figure 9shows 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. 8becomes 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). 5. Discussion 5.1. Is the IMF universal in LTGs? For quiescent galaxies, with typical ages of around 10Gyr, only stars with masses of ∼0.5M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.5M, which is within the mass range explored by IMF studies in quiescent galaxies. The young populations within NGC3351 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 Martín-Navarro, I., et al.: A&A, 684, A110 (2024) 10" 0.5 kpc 3 6 9 12 Age [Gyr] 10" 0.5 kpc 0.2 0.1 0.0 0.1 0.2 0.3 [M/H] 10" 0.5 kpc 0.00 0.06 0.12 0.18 0.24 0.30 [Mg/Fe] 10" 0.5 kpc 0.6 1.2 1.8 2.4 3.0 (IMF slope) Fig. 8. Best-fitting stellar population maps for NGC3351. 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 NGC3351 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 1000 or 0.5kpc at the distance of NGC 3351 is shown in each panel. 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 ∼10Gyr remains approximately constant at ∼0.5M(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 NGC3351, 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 Cand 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. Furthermore, as suggested by Fig. 10, a change in the age of the observed spectrum can naturally lead to a change in the A110, page 9 of 15