# Mild radial variations of the stellar IMF in the bulge of M31.

F. La Barbera1, A., Vazdekis2,3, I. Ferreras2,3,4, A. Pasquali5
1INAF-Osservatorio Astronomico di Capodimonte, sal. Moiariello 16, Napoli,
80131, Italy
2Instituto de Astrofísica de Canarias, Calle Vía Láctea s/n, E-38205 La
Laguna, Tenerife, Spain
3Departamento de Astrofísica, Universidad de La Laguna (ULL), E-38206 La
Laguna, Tenerife, Spain
4Department of Physics and Astronomy, University College London, Gower Street,
London WC1E 6BT, UK
5Astronomisches Rechen-Institut, Zentrum für Astronomie, Universität
Heidelberg, Mönchhofstr. 12-14, D-69120 Heidelberg, Germany
E-mail<EMAIL_ADDRESS>(FLB)

###### Abstract

Using new, homogeneous, long-slit spectroscopy in the wavelength range from
$\sim 0.35$ to $\sim 1\,\mu$m, we study radial gradients of optical and near-
infrared (NIR) IMF-sensitive features along the major axis of the bulge of
M31, out to a galactocentric distance of $\sim 200$ arcsec ($\sim 800$ pc).
Based on state-of-the-art stellar population synthesis models with varying Na
abundance ratio, we fit a number of spectral indices, from different chemical
species (including TiO’s, Ca, and Na indices), to constrain the low-mass
($\lesssim 0.5$ $M_{\odot}$) end slope (i.e. the fraction of low-mass stars)
of the stellar IMF, as a function of galactocentric distance. Outside a radial
distance of $\sim 10$”, we infer an IMF similar to a Milky-Way-like
distribution, while at small galactocentric distances, an IMF radial gradient
is detected, with a mildly bottom-heavy IMF in the few inner arcsec. We are
able to fit Na features (both NaD and $\rm NaI8190$), without requiring
extremely high Na abundance ratios. $\rm[Na/Fe]$ is $\sim 0.4$ dex for most of
the bulge, rising up to $\sim 0.6$ dex in the innermost radial bins. Our
results imply an overall, luminosity-weighted, IMF and mass-to-light ratio for
the M31 bulge, consistent with those for a Milky-Way-like distribution, in
contrast to results obtained, in general, for most massive early-type
galaxies.

###### keywords:

galaxies: stellar content – galaxies: fundamental parameters – galaxies:
formation – galaxies: elliptical and lenticular, cD

††pubyear: 2021††pagerange: Mild radial variations of the stellar IMF in the
bulge of M31.–E

##  1 Introduction

The Initial Mass Function (IMF), i.e. the distribution of the masses of stars,
at birth, in a stellar system, is a key ingredient of astrophysics, as it sets
the overall mass-scale of galactic systems, controls the intensity of the
stellar feedback processes and drives the chemical enrichment abundance
patterns. While the IMF has long been assumed to be “universal”, and the same
as in the solar neighborhood, a number of observational studies have found
evidence of IMF variations for both star-forming and quiescent systems (see
Hopkins, 2018, for a recent review). Understanding the origin of these
variations is one of the main challenges that nowadays the astronomical
community is facing on.

In early-type galaxies (hereafter ETGs), Cappellari et al. (2012, 2013) have
found a systematic increase in the stellar $\rm M/L$ with galaxy mass, based
on detailed dynamical models of the galaxies’ kinematics, with a Kroupa IMF
normalization at low velocity dispersion ($\sigma\sim 80$ km s-1)
transitioning to a Salpeter, i.e. bottom-heavier than the Milky-Way, one at
$\sigma\sim 260$ km s-1. Independent dynamical studies have achieved similar
conclusions (see, e.g., Thomas, Maraston & Johansson, 2011; Wegner et al.,
2012; Cappellari et al., 2012; Dutton, Mendel & Simard, 2012; Tortora et al.,
2013). Strong gravitational lensing over galaxy scales can also be exploited
to constrain the stellar $\rm M/L$. While strong lensing studies of low-mass
spheroids excluded a Salpeter IMF in favour of a Milky-Way-like distribution
(Ferreras et al., 2005, 2008, 2010), other works have detected a systematic
variations towards higher stellar $\rm M/L$ with increasing mass (Auger et
al., 2010; Treu et al., 2010; Barnabé et al., 2011), with some exceptions
(Smith & Lucey, 2013; Smith, Lucey & Conroy, 2015; Leier et al., 2016). In
general, most lensing and dynamical studies have consistently pointed to a
scenario where the IMF is either bottom- or top-heavier (both implying a
higher $\rm M/L$) than Kroupa/Chabrier, in massive ETGs. Indeed, the analysis
of spectral features sensitive to the presence of low-mass stars, such as the
Na I doublet feature at $\lambda\lambda 8183,8195$ Å (Faber & French, 1980;
Schiavon et al., 1997b, hereafter NaI8200), have shown that the above M/L
variations are likely due to a bottom-heavy IMF, i.e. an enhanced relative
contribution of dwarf versus giant stars, in most massive ETGs (see, e.g.,
Cenarro et al. 2003; van Dokkum & Conroy 2010; Ferreras et al. 2013; La
Barbera et al. 2013; Spiniello et al. 2014; but see Alton, Smith, Lucey 2017,
2018). Unfortunately, all these studies lacked spatial resolution within
galaxies, constraining only the “integrated” (light-weighted) IMF, on
apertures encompassing most of the galaxy total light. Since, overall, both
velocity-dispersion, elemental abundance ratios, age, and stellar metallicity
increase as a function of galaxy mass in ETGs (see, e.g., Gallazzi et al.
2021, and references therein), studies of integrated galaxy properties are not
able to fully disentangle the effect of different possible drivers of IMF
variations (La Barbera, Ferreras, & Vazdekis, 2015).

Therefore, subsequent studies have focused on radial gradients of the IMF in
ETGs (Martín-Navarro et al., 2015a, b; La Barbera et al., 2016, 2017;
Zieleniewski et al., 2017; van Dokkum et al., 2017; Parikh et al., 2018; Sarzi
et al., 2018; Domínguez et al., 2019; La Barbera et al., 2019), finding that a
bottom heavy IMF is only present in the central regions of the most massive
galaxies. This result seems to be also consistent with a two-phase scenario of
galaxy formation (Oser et al., 2010). Although radially resolved studies have
allowed us to enlarge significantly the parameter space where IMF variations
can be constrained, the problem is far from being sorted out. Based on
spatially-resolved measurements of the TiO2 IMF-sensitive features for ETGs in
the CALIFA survey, Martín-Navarro et al. (2015b) found that metallicity,
rather than local velocity dispersion, might be the primary local driver of
IMF variations. However, La Barbera et al. (2019, hereafter LB19) found that
in most massive ETGs, mostly brightest cluster galaxies, regions with high
metallicity do not necessarily exhibit a bottom-heavy IMF, while Martín-
Navarro et al. (2019) showed that in a disk-dominated galaxy, IMF variations
do not mimic metallicity changes.

The bulges of spiral galaxies might allow us to further expand the parameters’
space where IMF variations can be measured through spectral features. At least
some bulges have high metallicity and old ages in their central regions,
similar to those of massive ETGs (see, e.g., fig. 3 of Thomas & Davies 2006),
but lower velocity dispersion (a proxy for galaxy mass). Because of their
lower mass, bulges set tighter upper limits to possible IMF variations, and
thus offer a good benchmark to pinpoint possible systematics in the analysis
of IMF-sensitive spectral features. The bulge of M31, due to its proximity,
gives a unique opportunity to obtain high S/N spectroscopy, at low
observational cost, and apply the same methodology used so far to constrain
the stellar IMF of ETGs. Indeed, the stellar populations in the bulge are old,
enhanced in alpha elements, and have super-solar metallicity in the center,
similar to massive ETGs (Saglia et al., 2010, 2018). On the contrary, the
bulge’s velocity dispersion is low ($\sim 150$ km s-1).

The stellar IMF of the M31 bulge has been the subject of an intense debate
over the years, since the 1970s. Spinrad & Taylor (1971) measured a very
strong $\rm NaI8190$ feature, explaining it through a dwarf-dominated IMF,
with a high mass-to-light ratio ($\rm M/L=44$). A similar result was obtained
by O’Connell (1976), but with lower M/L (between 2 and 15). On the contrary,
based on NIR spectroscopy targeting the FeH Wing-Ford band (Schiavon et al.,
1997b) and the K-band CO absorption, Whitford (1977) and Cohen (1978)
concluded that the bulge IMF is similar to that of the solar-neighborhood. The
debate has continued over the years (see, e.g., Carter et al. 1986; Delisle &
Hardy 1992), but it is only with the advent of modern instrumentation and
improved stellar population models that we are now in the position to make a
significant leap forward.

Indeed, recent studies have shown that the bulge of M31 has a much lighter IMF
than what suggested in the 1970s. Conroy & van Dokkum (2012b, hereafter
CvD12b) analyzed spectra for the innermost few arcsec region of M31.
Performing spectral fitting in the ranges $0.4<\lambda<0.55$ $\mu$m and
$0.8<\lambda<1.02$ $\mu$m with stellar population models including the effect
of non-solar abundance ratios, CvD12b found evidence for an IMF normalization
between Kroupa and Salpeter ($\rm M/L_{r}\sim 6.2$). Zieleniewski et al.
(2015, hereafter Z15) analyzed NIR spectral features (NaI8200, CaT, and FeH)
at six radial positions in the bulge of M31 (out to a radial distance of
700”), performing a qualitative comparison of observed line-strengths to SSP
model predictions from Conroy & van Dokkum (2012a, hereafter CvD12a). Z15
concluded that M31 is well described by a Chabrier IMF at all radial
positions, hence favouring a lighter IMF in the center, with respect to
CvD12b. Both CvD12b and Z15 concluded that the strong NaI8200 absorption in
the center of M31 is likely due to Na over-abundance, rather than IMF (as
originally proposed by Spinrad & Taylor 1971), with $\rm[Na/Fe]$ as high as
$\sim 1$ dex. The NaI8200 radial gradient would then be explained as a strong
$\rm[Na/Fe]$ radial gradient (see Z15). However, an extremely high value of
$\rm[Na/Fe]$ is not without problems. In the center of the most massive ETGs,
LB19 found $\rm[Na/Fe]$ abundances not higher than $\sim 0.6$–$0.7$ dex, even
in the most metal-rich regions. Bensby et al. (2017), analyzing stars in the
Galactic bulge, have found, for individual stars, $\rm[Na/Fe]$ abundances
below $\sim 0.3$ dex, even at the highest metallicity probed ($\rm[Fe/H]\sim
0.5$ dex), far below the above estimate for the center of M31. We point out
that $\rm[Na/Fe]$ has been measured indirectly by CvD12b (i.e. based on model
predictions for the effect of $\rm[Na/Fe]$ on other spectral features), as
they lacked the region around NaD, the most sensitive feature to $\rm[Na/Fe]$
in optical galaxies’ spectra. Instead, Z15 used $\rm NaD$ and $\rm NaI8190$
measurements from different sources, with $\rm NaD$ from Davidge (1991, 1997).
Therefore, the question remains open, about the origin of Na absorption in the
center of M31.

The present work presents the most detailed, homogeneous, study of IMF-
sensitive spectral features along the major axis of the M31 bulge, from the
optical throughout NIR spectral range. To this effect, we have acquired new,
dedicated, high-S/N, spectroscopy along the major axis of the bulge, with the
OSIRIS (Optical System for Imaging and low-Intermediate-Resolution Integrated
Spectroscopy) spectrograph at the GRANTECAN (GTC) telescope, ensuring a
continuous, homogeneous, spectral coverage from $\lambda\sim 0.35\,\mu$m, to
$\lambda\sim 1\,\mu$m. Contrary to previous works, we target all IMF-sensitive
features up to $\lambda\sim 1\mu$m, including Mg, Fe, TiO, Na (both NaD and
$\rm NaI8190$), and Ca absorptions, and analyze them based on state-of-the-art
stellar population models that take the coupled effect of IMF and Na abundance
explicitly into account, applying the same approach that we have used so far
to infer the stellar IMF of massive ETGs (see, e.g., LB19 and references
therein).

The layout of the paper is as follows. In Sec. 2, we describe the new OSIRIS
data for the M31 bulge, the kinematics extracted from the new spectroscopy, as
well as the radially binned spectra used in our analysis. Sec. 3 describes the
stellar population models used to analyze the spectra. The analysis is
presented in Sec. 4, including the definition of spectral indices used in this
work (Sec. 4.1), the methods used to constrain IMF and other stellar
population properties (Sec. 4.2 and 4.3), as well as the comparison of model
and best-fitting spectral indices (Sec. 4.4). Sec. 5 shows the main results of
the present work, and in particular, the IMF and $\rm[Na/Fe]$ radial profiles
for the bulge of M31 (Sec. 5.2 and 5.3, respectively). We discuss the results
in Sec. 6. Summary and conclusions follow in Sec. 7. Throughout the present
work, we adopt a distance of 785 kpc from the MW to M31 (McConnachie et al.,
2005), implying a conversion scale of $\rm\sim 3.8$ pc/arcsec.

##  2 New UVRI spectroscopy for the M31 bulge

###  2.1 Observations

We obtained new long-slit spectroscopy for the bulge of M31 on August and
September 2017, using the OSIRIS instrument at the Nasmyth-B focus of the Gran
Telescopio CANARIAS (GTC), at Roque de los Muchachos Observatory. OSIRIS is
equipped with a mosaic of two 2k$\times$4k red-optimized CCDs, separated by a
narrow gap of 9.4”. Observations were carried out with the U-, V-, R-, and I-
R2500 grisms, using a 7.4’-long slit of width 0.4”, resulting in a spectral
coverage from $3460$ to $10150$ Å with a uniform spectral resolution of $\sim
38$ km s-1 (FWHM; as measured from arc lamps and sky emission lines). The M31
pointings were centered at RA=00:42:44.57 and DEC=41:16:05.7, with the slit
aligned along the bulge major axis, at PA=48 deg (see Saglia et al. 2010,
hereafter S10). The target was centered on CCD#2, 55” apart from the gap,
allowing us to observe the bulge up to a distance of $\sim 180^{\prime\prime}$
($\sim 260^{\prime\prime}$) along the East (West) direction. We adopted a 2x
(1x) binning along the spatial (dispersion) direction, resulting in a spatial
scale of $0.254^{\prime\prime}$/pixel, and a dispersion of $\sim 0.3$, $\sim
0.4$, $\sim 0.5$, and $\sim 0.7$ Å per pixel for the UVRI grisms,
respectively. For each grism, we performed two observations of the bulge, each
followed by a pointing on a blank sky region, centered at RA=00:44:27.20 and
DEC=40:49:54.8. For the VRI grisms, each observation consisted of three
dithered exposures of 150sec each, with a dithering pattern of $\sim$15”,
resulting in a total on-target exposure time of 900sec. For the U grism, we
performed five exposures, giving a total integration time of 1500sec. The
average seeing was 1, 1.2, 1.1, and $1^{\prime\prime}$ (FWHM) for the UVRI
grisms, respectively.

###  2.2 Data reduction

The data were reduced using dedicated FORTRAN and IRAF scripts written by the
authors. For each CCD, all frames were bias subtracted, trimmed, and flat-
fielded using twilight sky frames. For the I grism, wavelength calibration was
performed using sky lines only, while for the V and R grisms, we combined sky
lines with HgAr+Ne+Xe arc lamp lines. For the U grism, wavelength calibration
was performed with HgAr+Xe calibration lamp lines. For the R and I grisms, we
estimated the accuracy of the wavelength calibration by measuring the position
of sky lines in the calibrated frames, finding an rms (across the spatial
direction and the spectral range of each grism) of $\sim 4$ km s-1. Flux
calibration was performed using data for a spectrophotometric standard star
observed with the same instrumental setup as for the science exposures. For
each grism, the flux standard was observed right before the observations of
M31. For the R and I grisms, a standard star was also observed at different
positions on both OSIRIS CCDs, in order to map spatial variations of the
response function across the entire field of view. We found variations of up
to 10$\%$ at the edges of each spectral range, among the response functions of
CCD#2 and CCD#1. These response functions were interpolated across the spatial
direction and used to further improve the flux calibration.

Sky subtraction was performed by correcting each M31 observation with the
corresponding sky frame. Sky lines were rescaled to match those in the M31
frame by using the software SKYCORR (Noll et al., 2014). To this aim, each M31
frame was binned along the spatial direction, and a one-dimensional spectrum
was extracted from the M31 and the corresponding sky frames. For each bin, we
ran SKYCORR to derive the corresponding sky-lines rescaling factors.
Interpolating the scaling factors across the spatial direction provided us
with a rescaled 2D sky spectrum, that was subtracted off from the given M31
frame. Particular care was taken in the subtraction of the NaD sky lines at
$\lambda\lambda\sim 5890,5896$ Å. The (prominent) NaD stellar absorption in
M31 is potentially contaminated by dust absorption due to the ISM in M31 and
our Galaxy (see App. B for details). Since the ISM absorption from our Galaxy
overlaps, in wavelength, to the NaD sky lines, rescaling the NaD in the sky
spectra with those in the M31 frames may lead to underestimate the true
intensity of the NaD sky lines. Since we found that the intensity of the NaD
sky lines was stable among different observations, we decided to not apply any
rescaling. Instead, we modeled the NaD lines in each sky frame with a
combination of Gaussian functions, applied a shift to match the position of
the lines in the M31 frame (as a function of the spatial direction along the
slit), and subtracted them off. Before sky subtraction, all frames were
corrected for bad pixels and cosmic rays, by performing a linear interpolation
of affected pixels along the slit direction.

For the R- and I- grism, each sky-subtracted frame was corrected for telluric
lines, using the software MOLECFIT (Smette et al., 2015; Kausch et al., 2015).
MOLECFIT produces a theoretical transmission model by fitting selected
regions, with prominent telluric lines, in the object spectrum itself. We ran
MOLECFIT on a spectrum extracted in the innermost region of M31, verifying
that results were unchanged when considering spectra extracted at different
positions along the bulge.

For each observation, the sky-corrected M31 frames of both CCDs were joined
together, rectified (i.e. corrected for spatial variations of the photometric
center of the galaxy along the dispersion direction) and combined using the
IRAF task IMCOMBINE. At each reduction step, and for each science frame, a
variance map was produced and updated during the reduction process to account
for different sources of uncertainty, including that on sky subtraction111 To
estimate the uncertainty on the correction of NaD sky lines, we repeated the
subtraction using different sky frames, and included the corresponding error
budget in the variance maps of sky-subtracted frames.. Variance maps were
finally combined in the same way as for the science data.

###  2.3 Radially binned spectra

For each grism, we extracted radially binned spectra along both sides of the
slit. The central bin was set to be 1.5′′-wide around the photometric center
of the galaxy, while for the other bins we adopted a minimum width of
$0.75^{\prime\prime}$ ($\sim 3$ pixels). When required, the bin width was
adaptively increased outwards, to ensure a minimum signal-to-noise ratio, $\rm
S/N_{min}$. The S/N was computed in the U-grism spectral range, from 3900 to
4200 Å (corresponding to the CaH+K absorption lines), as the S/N was found to
be higher at redder wavelengths. Since we did not perform absolute flux
calibration, the binned spectra from different grisms were just joined by
applying suitable multiplicative factors, derived from overlapping spectral
regions.

We extracted two sets of binned spectra. First, we set $\rm S/N_{min}=25$
222We also repeated the analysis with $\rm S/N_{min}=15$, finding consistent
results to those for $\rm S/N_{min}=25$. to derive the kinematics of the
bulge, i.e. radial profiles of rotation velocity $\rm V_{rot}$, velocity
dispersion, $\rm sigma$, and higher moments of the line-of-sight (LOS)
velocity distribution, $\rm H_{3}$ and $\rm H_{4}$, respectively. As detailed
in App. E, the kinematics was extracted with the software pPXF (Cappellari &
Emsellem, 2004; Cappellari, 2017), performing spectral fitting on different
spectral regions, and combining results into final radial profiles. In App. E,
we show that the new kinematics is in good agreement with that obtained by
previous studies.

At each position along the spatial direction, the rotation velocity profile
was interpolated and used to correct the two-dimensional spectra of M31 to the
restframe. To perform the stellar population analysis (see below), a second
set of radially binned spectra was then extracted with $\rm S/N_{min}=70$
333We repeated the analysis with $\rm S/N_{min}=90$, finding consistent
results with those for $\rm S/N_{min}=$70\. Note that adopting $\rm
S/N_{min}=70$ ensured a S/N ratio larger than $\sim 100$ Å-1 for all the IMF
sensitive features analyzed in this work, allowing both IMF and abundance
ratios to be properly constrained (see, e.g., La Barbera et al. 2013,
hereafter LB13)., in order to perform the stellar population analysis, as
detailed below.

Fig. 1 plots, as an example, some of the $\rm S/N_{min}=70$ radially binned
spectra of M31, showing the excellent quality of the new OSIRIS data.

Figure 1: The figure plots, as an example, some of the M31 radially binned
spectra used for the stellar population analysis (with $\rm S/N_{min}=70$; see
the text). The distance to the center along the slit, $\rm R$, is encoded with
different colours, from blue (North-East side of the slit) through red (South-
West), as shown by the inset coloured bar. From bottom to top, the plot shows
the M31 spectra at distances of $R\sim$ -205, -27, 0 (i.e. the central
spectrum), +28, and +152 arcsec ($\sim-787$, -103, 0, 106 and 578 pc,
respectively) from the galaxy center. The bluest and reddest spectra
correspond to the largest galactocentric distances probed with our stellar
population analysis. The spectra have been normalized to one in the spectral
range from 4800 to 5400 Å, and arbitrarily shifted for displaying purposes.
The IMF-sensitivy features analyzed in the present work are marked with black
arrows.

##  3 Stellar population models

To analyze the spectra of M31, we rely on Na–EMILES and $\alpha$–enhanced
MILES stellar population models.

The Na–EMILES models are a dedicated version of the EMILES models, covering a
range of $\rm[Na/Fe]$ abundance ratios (see LB19 and references therein).
EMILES simple stellar population (SSP) models cover the spectral range from
$0.35$ to $5\,\mu$m, based on different empirical stellar libraries, namely
MILES in the optical range (Sánchez-Blázquez et al., 2006), up to
$\lambda\\!\sim\\!7410$ Å, Indo-US (Valdes et al., 2004) and CaT (Cenarro et
al., 2001a) out to $\lambda\\!\sim\\!8950$ Å (Vazdekis et al. 2012), and the
IRTF stellar library (Cushing et al., 2005; Rayner et al., 2009) at redder
wavelengths (see Röck et al. 2016 for details). The spectral resolution is
kept constant with wavelength (at FWHM=2.5 Å) for all libraries, except for
IRTF, having a constant $\sigma$=60 km s-1 (see figure 8 of Vazdekis et al.
2016). Note that the above spectral libraries follow the abundance pattern of
the Milky-Way, (i.e. they are approximately scaled-solar at solar metallicity,
and significantly $\alpha$–enhanced at metallicity below $\rm[Z/H]$ about
$-0.3$ dex). Hence, EMILES SSPs should be considered as “base” (not scaled-
solar) models. The Na–EMILES models are computed for two sets of scaled-solar
theoretical isochrones, namely the ones of Girardi et al. (2000) (Padova00)
and those of Pietrinferni et al. (2004) (BaSTI), the latter having cooler
temperatures at the low-mass end (see Vazdekis et al. 2015, and references
therein, for details). To cover a range in $\rm[Na/Fe]$, we apply theoretical
differential corrections for $\rm[Na/Fe]$ overabundance to each individual
stellar spectrum in the empirical libraries, and construct the Na–EMILES SSPs
based on scaled-solar isochrones (see La Barbera et al. 2017, hereafter LB17,
and LB19). Note that this approach differs from that used to construct
$\alpha$–enhanced MILES models (hereafter $\alpha$–MILES; see V15), where
corrections for non-solar $[\rm\alpha/{\rm Fe}]$ abundance ratios were applied
directly to the model SSPs, rather than to individual stars. Also, the
$\alpha$–enhanced models only cover the optical (MILES) spectral range and are
computed for BaSTI isochrones only, while Na–enhanced models are computed over
the optical plus NIR spectral range for both Padova00 and BaSTI isochrones.

Models are computed for different values of age, from $1.0$ to $\sim 14$ Gyr,
and total metallicity, $\rm[Z/H]$. For BaSTI SSPs (either Na–EMILES or
$\alpha$–enhanced MILES), five values of total metallicity are considered,
i.e. $\rm[Z/H]=\\{$ $-0.66$, $-0.35$, $-0.25$, $0.06$, $+0.26\\}$ 444Note that
BaSTI models are also computed for $\rm[Z/H]$$=0.4$. However, given the lower
quality of these models they are not used in the present analysis., while for
Padova00, we rely on models having $\rm[Z/H]=\\{$ $-0.71$, $-0.4$, $0$,
$+0.22\\}$, respectively. The Na–EMILES models are computed for
$\rm[Na/Fe]$$=${$0,0.3,0.6$, $0.9$} dex, while $\alpha$–enhanced models are
given for $[\rm\alpha/{\rm Fe}]$$=0$ (solar scale), and $[\rm\alpha/{\rm
Fe}]$$=+0.4$ dex.

Na–EMILES and $\alpha$–enhanced MILES models are computed for different IMF
shapes, and in particular for two power-law distributions, i.e a single power-
law and a low-mass tapered single power-law IMF, also referred to as
“unimodal” and “bimodal” distributions, respectively (see Vazdekis et al.
2003, V15, and Vazdekis et al. 1996). The single power-law and low-mass
tapered IMFs are defined by their logarithmic slopes, $\Gamma$ and
$\Gamma_{b}$, respectively. The lower and upper mass-cutoffs of the IMF are
set to $0.1$ and $100$ M⊙, respectively. We note that the difference between a
single power-law and a low-mass tapered IMF is that the latter is smoothly
tapered towards masses below $\sim 0.5\,M_{\odot}$; hence, varying the slope
$\Gamma_{b}$ changes the ratio of dwarf-to-giant stars in the IMF through its
overall normalization. While this approach is different with respect to a
change of the low-mass slope (e.g. Conroy & van Dokkum 2012b), we emphasize
that a low-mass tapered parametrization is suitable for our purposes, as most
IMF-sensitive features are only sensitive to the dwarf-to-giant ratio in the
IMF (e.g. LB13, La Barbera et al. 2016, hereafter LB16). Moreover, a low-mass
tapered IMF has been shown to provide mass-to-light ratios more consistent
with dynamical constraints (Lyubenova et al., 2016), and to be able to
describe both optical and NIR IMF-sensitive features in galaxy spectra (see
LB16). Therefore, in the present work, we consider only models with a low-mass
tapered distribution, that do also allow us to perform a direct comparison to
previous works for massive galaxies (see Sec. 6). The Na–EMILES and
$\alpha$–MILES models are computed for the following low-mass tapered IMF
slopes, $\rm\Gamma_{b}$=
{$0.3,0.5,0.8,1.0,1.3,1.5,1.8,2.0,2.3,2.5,2.8,3.0,3.3,3.5$}. For
$\Gamma_{b}=1.3$, the low-mass tapered IMF closely approximates the Kroupa
(2001) Universal IMF.

##  4 Stellar population analysis

###  4.1 Definition of spectral indices

The wide spectral range provided by the OSIRIS UVRI grisms allows us to probe
a wide set of optical and NIR spectral features. Following the same approach
as in our previous works (e.g. LB13, LB17, and LB19), we select a set of
optical and NIR spectral indices, and constrain the stellar population
properties of M31 by comparing their observed line-strengths to model
predictions. We consider the age-sensitive Balmer indices, $\rm H\beta_{o}$
and $\rm H\gamma_{F}$, the total metallicity indicator $\rm[MgFe]^{\prime}$,
as well as the IMF-sensitive features $\rm TiO1$, $\rm TiO2_{SDSS}$, $\rm
aTiO$, and $\rm Mg4780$, the two Na indices, $\rm NaD$ and $\rm NaI8190$, and
the Calcium triplet lines, $\rm Ca1$, $\rm Ca2$, and $\rm Ca3$, with the
combined $\rm CaT$=$\rm Ca1\\!+\\!Ca2\\!+\\!Ca3$ (Cenarro et al., 2001a). We
also include other indices in the analysis, such as $\rm Ca4227r$, $\rm
C4668$, $\rm CN2$, $\rm Fe4531$, Fe5015, $\rm Mgb5177$, $\rm Mg1$, $\rm Mg2$,
which are mostly sensitive to abundance ratios (see table 1 of LB15), and for
which we adopt the same central passband and pseudo-continua definitions as in
the Lick system (Trager et al., 1998). The index definitions for $\rm TiO1$,
$\rm H\gamma_{F}$, and $\rm NaD$ are also the same as in Trager et al. (1998),
while $\rm H\beta_{o}$ is the optimized $\rm H\beta$ index defined by
Cervantes & Vazdekis (2009). The $\rm TiO2_{SDSS}$ is defined as in LB13,
being a modified version of $\rm TiO2$ from Trager et al. (1998). The total-
metallicity indicator
$\rm[MgFe]^{\prime}$$\rm=[Mgb5177\cdot(0.72Fe5270+0.28Fe5335)]^{1/2}$ is a
combined Mgb and Fe index, defined by Thomas et al. (2003a) to be insensitive
to $\rm[Mg/Fe]$ abundance ratio (see also V15). Finally, the $\rm Mg4780$ is
from Serven (2005), $\rm aTiO$ from Spiniello et al. (2014), while $\rm
NaI8190$ is defined as in CvD12a, with some modifications as described in
LB17.

We note that (i) the spectral indices $\rm H\beta_{o}$, $\rm H\gamma_{F}$,
$\rm Mgb5177$, and Fe5015, are corrected for contamination from emission
lines, as detailed in App. A; (ii) for each radial bin, when fitting observed
and model line-strengths (see Sec. 4.2 and 4.3 below), models are first
smoothed to match the sigma 555We refer to sigma here to indicate the root
square of the sum in quadrature of the velocity dispersion of a given spectrum
and the instrumental resolution. of the given bin 666This approach maximizes
the information we can extract from the data, as it does not require any
smoothing of the spectra to bring all of them to the same sigma. However, in
all Figures of the present paper, to compare observed and model line-strengths
among different radial bins, we correct all line-strengths to the same sigma
of $\sigma_{0}=$150 km s-1. The correction is given as the index variation
between a reference SSP model smoothed to $\sigma_{0}$ and to the sigma of a
given bin. As a reference model we use, for each radial bin, the best-fitting
model from method A (see Sec. 4.3).; (iii) NaD was corrected for dust
contamination, as detailed in App. B.

###  4.2 Constraining age, metallicity, and $\rm[Mg/Fe]$

As a first step in the analysis, we estimate the age, metallicity, and
$\rm[Mg/Fe]$ abundance ratio. We adopt different methods:

i.

We use 1SSP model predictions from EMILES base models with Padova00 isochrones
(hereafter, EMILES iP), for a Kroupa-like IMF (i.e. a low-mass tapered IMF
with $\rm\Gamma_{b}$$=1.3$). Age and metallicity are estimated by fitting $\rm
H\beta_{o}$ and $\rm[MgFe]^{\prime}$ simultaneously (see e.g. LB13). Note that
the effect of $\rm[Mg/Fe]$ on $\rm H\beta_{o}$ is rather small (see V15),
while $\rm[MgFe]^{\prime}$ is independent of $\rm[Mg/Fe]$ (Thomas et al. 2004;
V15). A proxy for $\rm[Mg/Fe]$ is obtained from $\rm Mgb5177$ and $\rm\langle
Fe\rangle$ and then converted to $\rm[Mg/Fe]$ with a similar approach as
described in LB13 and V15.

ii.

Same as at point i, but using EMILES base models with BaSTI isochrones
(hereafter BASTI iT).

iii.

We fit simultaneously $\rm H\beta_{o}$, $\rm[MgFe]^{\prime}$, $\rm Mgb5177$,
and $\rm\langle Fe\rangle$ using $\alpha$–MILES SSP models.

iv.

Age and metallicity are obtained by averaging out results from different sets
of spectral indices, fitting simultaneously also the IMF slope (see Sec. 4.3
below).

Note that we use different methods/models in order to account for possible
systematic effects, and/or the effect of uncertainties on stellar population
models.

###  4.3 Constraining the stellar IMF

We constrain the stellar IMF following the same approach as in LB19 (and
references therein). For each spectrum, we minimize the expression,

$\displaystyle\rm\chi^{2}(Age,{\rm[Z/H]},\Gamma_{b},[X/Fe])=\rm\sum_{i}\left[\frac{E_{obs,i}-E_{mod,i}}{\sigma_{E_{obs,i}}}\right]^{2}$
(1)

where the index $i$ runs over a selected set of spectral features (see below);
$\rm E_{mod,i}$ are line-strength predictions for Na–EMILES stellar population
models; $\rm E_{obs,i}$ and $\sigma_{E_{obs,i}}$ are observed line-strengths
and their uncertainties; $\rm[X/Fe]$ are the elemental abundance ratios (e.g.
$\rm[Na/Fe]$) included in the fitting procedure. For each fitting case (see
below), uncertainties on best-fitting parameters, {Age, ${\rm[Z/H]}$,
$\rm\Gamma_{b}$, $\rm[X/Fe]$}, are obtained from $N=1000$ bootstrap
iterations, where the fitting is repeated after shifting observed line
strengths according to their uncertainties.

In order to obtain robust results, we consider different fitting methods, by
(i) changing the set of fitted indices, (ii) changing the method to account
for non-solar abundance ratios (see below); (iii) exploring either 1SSP or a
combination of two SSP models, and/or imposing additional constraints to the
(luminosity-weighted) age of the best-fitting model; (iv) using models based
on different isochrones, i.e. BaSTI and Padova00 (see Sec. 3). The options
adopted for different methods are summarized in Tab. 1. The rationale behind
options (ii) and (iii) is the following:

Abundance ratios -

In our “reference” approach (method A), we correct line-strengths to solar
scale as a function of metallicity (see col. 6 of Tab. 1). The corrections are
based on observed trends of line-strengths, at fixed galaxy velocity
dispersion, for SDSS ETGs’ spectra. As shown in LB16 and LB19, this approach
was able to match well the observed line indices of massive ETGs at different
galacto-centric distances. Since we found that the same approach does not
allow us to match all observed indices for M31 (see below), we also consider
different options (see col. 5 of Tab. 1), where additional elemental abundance
ratios, $\rm[X/Fe]$, are included in the fitting procedure. To this effect,
the sensitivity of a given index to $\rm[X/Fe]$ is estimated from CvD12a
stellar population model predictions, and an additional term is included in
Eq. 1 (see LB13 and LB15 for details). Note that this method is not applied to
$\rm[Na/Fe]$, whose effect is already taken into account by Na–EMILES models.
In general, when including a specific elemental abundance in the fitting, we
also include indices with a prominent sensitivity to it 777For instance, in
method E, we also include $\rm[Ca/Fe]$ abundance ratio in the fitting
procedure, since CaT lines are sensitive to both IMF and Ca abundance. The
$\rm[Ca/Fe]$ is constrained by adding $\rm Ca2$ as well as $\rm Ca4227$ to the
list of indices. Since $\rm Ca4227$ is anticorrelated to $\rm[C/Fe]$, we
include also $\rm[C/Fe]$ as a further fitting parameter (see Tab. 1), and
constrain the $\rm[C/Fe]$ by fitting also the $\rm C4668$ band.. Note that
when applying empirical corrections and including $\rm[X/Fe]$’s in the fitting
procedure, one should interpret the $\rm[X/Fe]$’s as “residual” abundance
ratios (i.e. not accounted for by the empirical correction, as discussed in
LB13).

Constraints to age -

As discussed in LB13, Balmer lines have some sensitivity to IMF. The $\rm
H\beta_{o}$ index decreases with IMF, so that for a bottom-heavier
distribution one tends to infer younger ages, and vice versa. However, this
IMF sensitivity turns out to be model dependent, as CvD12a models do not show
a dependence of Balmer lines on IMF (see V15). Moreover, the $\rm H\beta_{o}$
decreases with $\rm[Mg/Fe]$, and is sensitive to other abundance ratios (e.g.
$\rm[C/Fe]$). In order to overcome these issues, for each spectrum of M31, we
also derive a luminosity-weighted age estimate with spectral fitting, using
the software pPXF. For some fitting methods (see col. 7 of Tab. 1), the pPXF
luminosity-weighted age is also added as extra constraint (i.e. an extra term)
in Eq. 1 (see eq. 1 of LB19). To explore the effect of different star-
formation histories, we also use both 1SSP and 2SSP models, where the effect
of a second young component is taken into account (see methods J and K in Tab.
1).

Table 1: Summary of the different methods used to infer the IMF slope. Column 1 lists the labels used for different methods. Column 2 reports the number of SSPs in the fitting procedure, while column 3 refers to the models’ isochrones, with label “iP” (“iT”) for Padova00 (BaSTI) isochrones (see Sec.3). Columns 4 and 5 give the list of indices and the elemental abundance ratios included in the fitting. Columns 6 and 7 flag the case where empirical corrections to line-strengths and age constraints from spectral fitting were applied (see the text). Method | #SSPs | isochrones | Spectral indices |  $\rm[X/Fe]$’s | Empirical | Age
---|---|---|---|---|---|---
|  |  |  |  | corrections | constraints
(1) | (2) | (3) | (4) | (5) | (6) | (7)
A | 1 | iP |  $\rm H\beta_{o}$, $\rm[MgFe]^{\prime}$, $\rm Mg4780$, $\rm TiO1$ | none | yes | yes
|  |  |  $\rm TiO2$, $\rm aTiO$, $\rm NaD$, $\rm NaI8190$, $\rm CaT$ |  |  |
B | 1 | iP | same as A w/o TiO’s | none | yes | yes
C | 1 | iP | same as A w/o CaT | none | yes | yes
D | 1 | iP | same B | none | yes | no
E | 1 | iP |  same as A, $\rm Ca2$, $\rm Ca4227$, $\rm C4668$ |  $\rm[C/Fe]$, $\rm[Ca/Fe]$ | yes | yes
F | 1 | iP |  same as E w/o $\rm TiO2$, $\rm aTiO$ |  $\rm[C/Fe]$, $\rm[Ca/Fe]$ | yes | yes
G | 1 | iT | same as E |  $\rm[C/Fe]$, $\rm[Ca/Fe]$ | yes | yes
H | 1 | iP |  same as E w/o $\rm H\beta_{o}$ |  $\rm[C/Fe]$, $\rm[Ca/Fe]$ | yes | yes
I | 1 | iP |  $\rm H\beta_{o}$, $\rm[MgFe]^{\prime}$, $\rm TiO1$, $\rm TiO2$, $\rm aTiO$, |  $\rm[Ca/Fe]$, $\rm[Ti/Fe]$, $\rm[O/Fe]$, $\rm[C/Fe]$, | yes | no
|  |  |  $\rm Mg4780$, $\rm NaD$, $\rm NaI8190$, $\rm Ca1$, |  $\rm[N/Fe]$, $\rm[Mg/Fe]$, $\rm[Si/Fe]$ |  |
|  |  |  $\rm Ca2$, $\rm CaT$, $\rm Ca4227$, $\rm Mg1$, $\rm Mg2$ |  |  |
|  |  |  $\rm C4668$, $\rm CN2$, $\rm Mgb5177$, $\rm Fe4531$ |  |  |
J | 2 | iP | same as B | none | yes | no
K | 2 | iP |  same as I + $\rm H\gamma_{F}$ | none | yes | yes

###  4.4 Best-fit spectral indices

Fig. 2 compares observed and best-fit line-strengths for the bulge of M31, as
a function of galactocentric distance $\rm R$, in units of arcsec. Fig. 3
shows the same comparison but zooming into an innermost region of about $\pm
20$” from the center. To illustrate the quality of the fits, the Figures show
a selected set of spectral indices, including all the IMF-sensitive features
analyzed in this work, i.e. the TiO’s (panels d and e), Na indices (panels g
and h), $\rm aTiO$ (panel f), $\rm Mg4780$ (panel c), and the $\rm CaT$ (panel
i), as well as the age and metallicity indicators $\rm H\beta_{o}$ and
$\rm[MgFe]^{\prime}$ in panels a and b, respectively. In order to compare
different radial bins, all the line-strengths plotted in Figs. 2 and 3 have
been corrected 888 Note that this correction is only performed to display
observed and model indices for all radial bins in the same figure. For each
radial bin, we fit observed line-strengths with predictions of models smoothed
at the sigma of the given bin. to the same sigma 999The correction is
estimated using results from method A. In a given radial bin, the SSP model
corresponding to the best-fitting parameters (age, metallicity, IMF, and
$\rm[X/Fe]$) is smoothed at $\sigma=150$ km s-1 and at the actual sigma of the
bin. The correction is performed by computing the difference of the model
line-strengths measured for both sigma values, and adding this difference to
the observed line-strength.  of 150 km s-1. Observed line-strengths, with
error bars, are plotted in black, while best-fitting indices are plotted with
different colours, for a representative set of fitting methods from those
listed in Tab. 1.

\-

Method A (red line) is our basic approach, where we fit all IMF sensitive
features, together with $\rm H\beta_{o}$ and $\rm[MgFe]^{\prime}$, applying an
empirical correction to line-strengths for the effect of abundance ratios (see
Sec. 4.3). Fig. 2 shows that most of the best-fitting indices (red line) are
consistent with the observed ones within the error bars, although a
significant offset is seen for $\rm TiO1$, $\rm aTiO$, and $\rm CaT$. For the
TiO’s and CaT, the best-fit does not match the rapid increase of observed
indices in the bulge central region (Fig. 3).

-

Method E (blue line) shows the case where some residual abundance ratios, i.e.
$\rm[C/Fe]$ and $\rm[Ca/Fe]$, are included in the fitting. While including
$\rm[Ca/Fe]$ removes most of the offset for $\rm CaT$ (see above), as in the
case of method A, the best-fit solution does not provide a good match of the
radial trends in the innermost radial bins.

-

Method I (green line) fits all abundance ratios by combining Na–EMILES models
with responses from CvD12a 101010The effect of $\rm[Na/Fe]$ is already
accounted for by Na–MILES. Hence, no response from CvD12a is used in this
case. (see Tab. 1). We see that best-fitting line-strengths match extremely
well all the observed line-strengths, but for the innermost radial bin, where
models (green curve) do not match the peak of $\rm TiO1$ and $\rm TiO2$ in the
center (see Fig. 3). Method K (magenta curve) shows that a younger component
(i.e. a 2SSP fit) does not improve the matching.

In summary, Figs. 2 and 3 show that we can fit reasonably well all observed
indices of M31, at all radial position, but for the central bin, where we do
not match the high values of $\rm TiO1$ and $\rm TiO2$. As discussed in App.
D, the mismatch is likely due to the models’ extrapolation in the very high
metallicity regime, that predicts both TiO indices to be independent of
metallicity up to $\rm[Z/H]$$\sim 0.4$ and above. Indeed, if we modify our
best-fitting method I to allow for both $\rm TiO1$ and $\rm TiO2$ to increase
with $\rm[Z/H]$ (see App. D for details), we are able to fit all observed
line-strengths of M31 in the innermost bin, as shown by the the green
asterisks in Figs. 2 and 3. Some deviation remains only for $\rm Mg4780$ in
the central bin (see green asterisk in panel c of Fig. 3), though the
discrepancy is not significant within the error bars ($<1.5$ sigma level).
Note also that $\rm Mg4780$ shows some asymmetric behaviour with radius, with
marginally lower values at negative $\rm R$, though differences are not
significant within the error bars 111111 These differences might be explained
by the fact that negative and positive values of $\rm R$ have been observed
with different OSIRIS CCDs (see also Sec. 2.2 and Sec. 5.2). . For the present
purposes, the key point is that, regardless of the fitting quality of
individual methods, the fitting results, and in particular the IMF
determination for M31, are very robust, as shown in the following sections.

Figure 2: Observed indices for the bulge of M31 as a function of
galactocentric distance. Different panels correspond to different spectral
indices. All indices (black dots with error bars) have been corrected to a
sigma of 150 km s-1 (see the text). Black error bars account for different
sources of uncertainties on observed line-strengths, as detailed in Sec. 4.
Grey error bars account for the uncertainty on the empirical correction
(methods A and E; see the text). Note that grey error bars are not shown for
$\rm H\beta_{o}$, $\rm[MgFe]^{\prime}$, and Na indices, as no correction is
applied to these indices. For all the other indices but $\rm CaT$, grey error
bars are almost indistinguishable from the black ones. Solid lines with
different colours show best-fitting results for a representative set of
different fitting methods (i.e. methods A, E, I, and K, see labels in panel b,
and Sec. 5.2). The green asterisks show best-fitting results for the innermost
radial bin when treating the dependence of $\rm TiO1$ and $\rm TiO2$ on
metallicity as free fitting parameters (see the text).  Figure 3: Same as Fig.
2 but zooming into the inner region within $\sim 20$” from the bulge center.

##  5 Results

###  5.1 Age, metallicity, and $\rm[Mg/Fe]$ radial profiles

Fig. 4 plots the radial profiles of age, metallicity, and $\rm[Mg/Fe]$ for the
bulge of M31, comparing the different methods we adopted to derive these
quantities, as detailed in Sec. 4.2. The right panels in the Figure zoom into
the inner bulge region (within a radial distance of $\sim 20$” from the
center).

Our OSIRIS spectroscopy shows that the bulge of M31 is characterized by old
stellar populations, with ages older than $\sim 11$ Gyr, at all radial
positions, but for the innermost region, within $\sim 6$” ($\sim 23$ pc),
where the age is as young as $\sim 5$–$7$ Gyr. For what concerns metallicity,
some offset exists among results from different models. For radial distances
$\gtrsim 20$” from the center, metallicity is slightly subsolar (around
solar), with $\rm[Z/H]$$\sim-0.2$ ($\rm[Z/H]$$\sim 0$) for EMILES
($\alpha$–MILES) models. In the innermost region (within $\sim 10$ arcsec),
where age is younger, metallicity is as high as $\sim 0.3$–$0.4$ dex
(depending on the method). The $\rm[Mg/Fe]$ is almost constant ($\sim
0.2$–$0.24$ dex), at all radii probed in this work (out to $\sim 200$”),
decreasing to $\sim~{}0.1$–$0.15$ dex (depending on the method used to
estimate $\rm[Mg/Fe]$) only in the innermost radial bin. The presence of a
(mildly) positive $\rm[Mg/Fe]$ gradient is also confirmed by the analysis of
the $\rm Mgb5177$–Fe index-index diagram, as shown in App. C, where we also
discuss differences between results from EMILES and $\alpha$–MILES models,
respectively. The fact that the innermost radial bins are characterized by a
young age and lower $\rm[Mg/Fe]$ points to a more extended star-formation
history in the central region (within $\sim 20$ pc) of the bulge. The magenta
curves in the bottom panels of Fig. 4 show that the $\rm[Mg/Fe]$ is
essentially independent of the assumed IMF, with a bottom-heavy distribution
($\rm\Gamma_{b}$$=2.8$) giving consistent results 121212The fact that
$\rm[Mg/Fe]$ is independent of IMF is relevant as we use the $\rm[Mg/Fe]$
estimates obtained from $\alpha$–MILES models with a Kroupa-like IMF, as input
to our fitting procedure to constrain the stellar IMF (Sec. 4.3). to those for
a Kroupa-like IMF, in agreement with what already found by LB17.

Figure 4: Age (top), metallicity (middle), and $\rm[Mg/Fe]$ (bottom) as a
function of galactocentric distance for the bulge of M31. The panels on the
right zoom into the innermost region of the bulge, within $\sim 20$” from the
center. The bottom (top) axes show the radial distance $\rm R$, in units of
arcsec ($\rm pc$). Black, red, blue, and green curves show results obtained
with different stellar population models, as well as by averaging results from
different methods to constrain the stellar IMF (green curve), as labeled in
the top–left panel (see the text). In the bottom panels, magenta curves and
symbols show $\rm[Mg/Fe]$ estimates obtained from $\alpha$–MILES models with a
bottom-heavy ($\rm\Gamma_{b}$$=2.8$), rather than a Kroupa-like, IMF. Error
bars denote 1-sigma uncertainties.

###  5.2 The IMF radial gradient

Fig. 5 shows the IMF radial profile for the bulge of M31. Results from
different methods to infer the IMF slope, $\rm\Gamma_{b}$, are shown with
different colours (as labeled in the top–left of the Figure; see Sec. 4.3 and
Tab. 1). The final IMF profile is obtained by averaging results from different
methods, and is shown as a thick black curve with error bars. We find a mild,
but significant, IMF radial gradient for the bulge of M31, with
$\rm\Gamma_{b}$ being consistent with a Kroupa-like distribution
($\rm\Gamma_{b}$$\sim 1.3$) in the outermost regions probed by OSIRIS ($\rm
R\gtrsim 100$”), and a “mildly” bottom-heavy distribution
($\rm\Gamma_{b}$$\sim 2.5$) in the central bins. For comparison, massive ETGs
can have a value of $\rm\Gamma_{b}$ as high as $3$, or even higher, in their
central regions (see Sec. 6). Note that although some offset exists among
different methods to infer the IMF, likely because of different
systematics/model uncertainties when using different sets of spectral indices,
the IMF radial trend is qualitatively the same for all methods. This points to
the importance of using a wide set of absorption features, from different
chemical species, when constraining the IMF (see LB13 and LB17). We note that
also modifying method I to allow for a possible dependence of both TiO indices
on $\rm[Z/H]$ (see App. D), does not change significantly our conclusions. In
this case, we get an IMF slope of $\rm\Gamma_{b}$$\sim 2.3$ at $\rm R=0$,
compared to $\rm\Gamma_{b}$$\sim 2.4$ for method I (see Fig. 5).

The IMF profile of M31 suggests the presence of three different radial
regions, marked with grey thick segments in Fig. 5, i.e. (i) the radial range
from $\sim 100$ to $\sim 200$”, with an average value of $\rm\Gamma_{b}$$\sim
1.4$ (see the grey dotted thick segments), where the IMF slope is fully
consistent with a Kroupa-like distribution (i.e. $\rm\Gamma_{b}$$\sim 1.3$,
plotted as a grey horizontal dot-dashed line); (ii) the range from $\sim 10$
to $\sim 100$”, where the IMF is marginally above a Kroupa-like distribution,
with $\rm\Gamma_{b}$$\sim 1.8$ (see grey dashed thick segments); and (iii) the
innermost region, with $\rm|R|<10$”, where the IMF is mildly bottom-heavy,
with $\rm\Gamma_{b}$$\sim 2.5$ (see the grey solid thick segment). Note that
the IMF profile appears to be slightly asymmetric, with values of
$\rm\Gamma_{b}$ higher for negative (relative to positive) values of $\rm R$.
The effect is more pronounced at larger galactocentric distances. For $\rm
R\gtrsim 100$”, the median value of $\rm\Gamma_{b}$ is $\sim 1.3$, while for
$\rm R\lesssim-100$”, we obtain a median value of $\rm\Gamma_{b}$$\sim 1.75$.
We note that the difference is small, given the quoted error bars, and mimics
the small asymmetry seen in some of the observed line-strength profiles. In
particular, as shown in Fig. 2, while the radial profile of $\rm NaI8190$ is
symmetric, the $\rm NaD$ tends to be lower for $\rm R<-100$”, hence favouring
lower $\rm[Na/Fe]$, and thus higher $\rm\Gamma_{b}$, than for $\rm R>100$”.
This behaviour might be due to (i) intrinsic small differences in the stellar
IMF among the two sides of the slit; (ii) a result of some residual
contamination of NaD from interstellar absorption (see App. B); (iii) the fact
that negative and positive values of $\rm R$ are observed with different
OSIRIS CCDs, and thus may be affected by different (small) systematic effects
(e.g. flux calibration). Note that some (small) asymmetries are also seen in
the radial profiles of $\rm TiO1$, $\rm aTiO$, and $\rm Mg4780$. While it is
hard to establish if such asymmetries are real, or result from data reduction,
this issue does not affect at all the main conclusions of our analysis. Also,
we emphasize that it is far more difficult to constrain the IMF slope when it
is close to $\sim 1.3$ compared to the range with, e.g.,
$\rm\Gamma_{b}$$\gtrsim 2.3$, as indices have a much weaker dependence on IMF
slope for low $\rm\Gamma_{b}$ values (see e.g. LB13).

Figure 5: The IMF slope, $\rm\Gamma_{b}$, for the bulge of M31 is plotted as a
function of the radial distance $\rm R$, in units of arcsec ($\rm pc$), along
the bottom (top) horizontal axis. Negative and positive values of $R$ refer to
opposite directions along the OSIRIS slit (see Sec. 2.3). Different colours
are the IMF profiles obtained with different fitting methods to spectral
indices (see labels in the top–left of the Figure and Tab. 1). The median-
combined profile of $\rm\Gamma_{b}$ is shown with a black thick line and
circles with error bars, corresponding to 1–sigma standard errors. The
horizontal grey dot-dashed line marks the value of $\rm\Gamma_{b}$$=1.3$, for
which the low-mass tapered IMF approximates the Kroupa IMF. Dotted, dashed,
and solid thick segments, plotted in grey, define different radial regions,
where the IMF is consistent with a Kroupa-like distribution ($100\leq|R|\leq
200$”), marginally above it ($10\leq|R|\leq 100$”), mildly bottom-heavy
($|R|<10$”), respectively (see the text). The y-axis on the right side of the
plot shows the mass fraction of low-mass (<$0.5\,M_{\odot}$) stars in the IMF,
as a function of $\rm\Gamma_{b}$, as defined in LB13.

Fig. 6 (top) plots the stellar mass-to-light ratio in r band, $\rm
M_{\star}/L_{r}$, as a function of radial distance to the center of the M31
bulge (see black dots with error bars). For each radial bin, the $\rm
M_{\star}/L_{r}$ is obtained by averaging results from all different fitting
methods (see above). The $\rm M_{\star}/L_{r}$ is approximately flat outside
$60$”, with $\rm M_{\star}/L_{r}$$\sim 3.9$, increasing up to $\sim 6$ in the
innermost few arcsec region. Such an increase results from the radial
variation of age, metallicity, and IMF. To single out the IMF effect, the
bottom panel of Fig. 6 plots the “mass-excess” factor, $\rm\alpha_{r}$, i.e.
the actual $\rm M_{\star}/L_{r}$ normalized to that predicted for a Kroupa-
like IMF, keeping the other stellar population parameters (i.e. age and
metallicity) unchanged. A value of $\rm\alpha_{r}\sim 1$ ($\sim 1.54$)
corresponds to a Kroupa (Salpeter) IMF, as shown by the horizontal grey dashed
lines in the bottom panel. Overall, the radial behaviour of $\rm\alpha_{r}$
reflects that of $\rm\Gamma_{b}$, with a mildly bottom-heavy IMF (Salpeter, or
“slightly” above Salpeter), in the innermost few arcsec, and a distribution
consistent with a Kroupa-like IMF in the outer radial bins. Note that from
$\rm\Gamma_{b}$$=1.3$ to $1.8$, the $\rm M_{\star}/L_{r}$ increases by only
$\sim 12$% (assuming an SSP model with age of 12 Gyr and solar metallicity).
Hence, it is difficult to define different radial ranges for $\rm
M_{\star}/L_{r}$ and $\rm\alpha_{r}$ as for $\rm\Gamma_{b}$ (see the grey
segments in Fig. 5).

Figure 6: Stellar mass-to-light ratio, $\rm M_{\star}/L_{r}$ (top), and mass-
excess factor, $\rm\alpha_{r}$ (bottom), for the bulge of M31, as a function
of galactocentric distance (i.e. the same as in Fig. 5 for IMF slope). The
$\rm\alpha_{r}$ is obtained by normalizing the $\rm M_{\star}/L_{r}$ by the
expected value for a Kroupa-like IMF (see the text). The black circles with
error bars (1-sigma confidence levels) are obtained by averaging results from
different fitting methods (Sec. 4.3). The red dots correspond to estimates
obtained by CvD12b for the nuclear region (within a radius of 4”), while the
light-green region in the top panel marks the range of values for $\rm M/L$
found by S10, under the assumption of a fixed Kroupa IMF (see Sec. 6.1). In
the bottom panel, the dashed grey horizontal lines mark the values of
$\rm\alpha_{r}$ for a Kroupa ($=1$) and Salpeter ($=1.54$) IMF, as labeled.

###  5.3 [Na/Fe] radial gradients

Since we rely on Na–enhanced EMILES stellar population models, our fitting
procedure provides also an estimate of $\rm[Na/Fe]$ at different radial
positions. The $\rm[Na/Fe]$ radial profile for the bulge of M31 is shown in
Fig. 7, with different methods plotted with different colours, and the
combined $\rm[Na/Fe]$ shown in black, as for $\rm\Gamma_{b}$ in Fig. 5.
Overall, the $\rm[Na/Fe]$ profile of M31 tends to be rather flat at all
positions, but for the innermost region ($\lesssim 5$”). For $R>5$”, the
median value of $\rm[Na/Fe]$ amounts to $0.41\pm 0.06$ dex, where the error
bar reflects the scatter among different radial bins. The maximum value of
$\rm[Na/Fe]$ is reached in the central bin, with $\rm[Na/Fe]$$=0.59\pm 0.05$
dex. For comparison, massive ETGs can have $\rm[Na/Fe]$ in the range
$0.6$–$0.9$ dex, in their central regions (see, e.g., LB19; and Sec. 6). We
point out that, as discussed in LB17, the effect of $\rm[Na/Fe]$ and IMF are
coupled in Na–MILES models, so that the effect of $\rm[Na/Fe]$ on Na indices
is stronger for a bottom-heavy, relative to a Kroupa-like, distribution. This
allows us to match the strong Na line-strengths of M31, in the innermost
regions, without requiring a very bottom-heavy IMF, nor an extremely high
$\rm[Na/Fe]$.

Figure 7: The same as Fig. 5 but plotting $\rm[Na/Fe]$ abundance ratio, rather
than IMF slope.

##  6 Discussion

###  6.1 Age, metallicity, and $[\rm\alpha/{\rm Fe}]$

Our age and metallicity radial profiles (Fig. 4) are qualitatively consistent
with those of S10 (see also Saglia et al. 2018), who found that, overall, the
bulge of M31 is characterized by old ages ($>10$ Gyr) and solar metallicity,
with younger ages ($\sim 8$–10 Gyr) and supersolar metallicities (up to
$\rm[Z/H]$$\sim$0.4 dex) at galactocentric distances smaller than a few arcsec
(see fig. 12 of S10). We also find a small, positive, radial gradient of
$\rm[Mg/Fe]$ (see App. C), not detected by S10, likely because of the
uncertainties on $[\rm\alpha/{\rm Fe}]$ (see, e.g., their fig. 12). Since the
$\rm[Mg/Fe]$ can be interpreted as a proxy of the star-formation timescale in
a stellar system (Thomas et al. 2005; de La Rosa et al. 2011; but see Fontanot
et al. 2017, 2018), the presence of young ages and lower $\rm[Mg/Fe]$ in the
innermost radial bins of M31 consistently point to a more extended/recent
star-formation in the bulge center. For the first time, we also find a (mild)
radial variation of the stellar IMF in the bulge of M31, with a Salpeter-like
($\rm\alpha_{r}\sim 1.7\pm 0.2$) IMF normalization (Fig. 6) within the
innermost few arcsec, decreasing to a Kroupa-like normalization at larger
radii (see Sec. 5.2).

Our results for the central radial bins of M31 can be compared to those of
CvD12b, who constrained the stellar IMF within a region of $4$” from the
center of M31, using a different set of models, and a different approach (i.e.
full spectral fitting, rather than index fitting). CvD12b found an IMF
normalization between Kroupa and Salpeter for the nuclear region of M31. Fig.
6 plots, as red dots, the $\rm M_{\star}/L_{r}$ and $\rm\alpha_{r}$ derived by
CvD12b. Our $\rm\alpha_{r}$ tends to be higher, and only marginally consistent
(at about 2 sigma level), with that of CvD12b. Nonetheless, the $\rm
M_{\star}/L_{r}$ is consistent with that of CvD12b, most likely because of the
IMF parametrization adopted in the present work (i.e. a low-mass tapered
distribution, rather than a multi-component power-law, as in CvD12b), as well
as the young age (implying lower $\rm M_{\star}/L_{r}$) that we estimate in
the innermost radial bins. CvD12b found also an extremely high $\rm[Na/Fe]$
abundance ratio ($\sim 1$ dex), in disagreement with our results, of
$\rm[Na/Fe]$ $\sim 0.4$ dex, rising up to $\sim 0.6$ dex (at most) in the
inner few arcsec (see Sec. 5.3). Note that the data used by CvD12b for M31 did
not include the $\rm NaD$ feature, which is crucial, indeed, to estimate
$\rm[Na/Fe]$ and thus disentangle the effect of $\rm[Na/Fe]$ and IMF on Na
features (see LB17).

###  6.2 Mass-to-light ratios

S10 analysed spectral indices for the bulge of M31, using stellar population
models by Maraston (2005). Under the assumption of a constant Kroupa IMF, they
predicted R-band M/L ratios in the range of 4 to 4.5. This range,
corresponding to $\rm 4.3\lesssim M/L_{r}\lesssim 4.9$ 131313Based on EMILES
SSP models, we adopt a conversion factor of $\rm(M/L_{R})/(M/L_{r})\sim
0.925$. , is shown as a green shaded region in Fig. 6 (top panel). To perform
a more direct comparison, we also computed from our M/L profile, under the
assumption of circular symmetry, the luminosity-weighted value of $\rm
M/L_{r}$ for the entire bulge, out to the maximum galactocentric distance
($\sim 200$”) covered by our analysis. We find $\rm\langle
M/L_{r}\rangle=4.1^{+0.6}_{-0.4}$. The predictions of S10 are fully consistent
with our $\rm\langle M/L_{r}\rangle$ estimate. However, since we find a
radially varying IMF, our M/L’s change significantly with radius, from $\sim
3.8$ (at R$\sim 200$”) up to $\sim 6$ (in the innermost radial bins). Note
that Saglia et al. (2018) also found, under the assumption of a Kroupa IMF,
V-band M/L ratios in the range of 4.5 to 5 (see, e.g., their figs. 21–22),
corresponding to $3.9\lesssim M/L_{r}\lesssim 4.4$, still consistent with our
$\rm\langle M/L_{r}\rangle$, and consistent with our M/L estimates for
galactocentric distances $\gtrsim 10^{\prime\prime}$. S10 concluded that there
is a good agreement between their estimated M/L and those from the dynamical
model of Widrow, Perrett, Suyu (2003, hereafter W03), with $\rm M/L_{R}\sim
4.5$ (i.e. $\rm M/L_{r}\sim 4.9$). Note again that this value is consistent
with our average M/L. However, since the dynamical models of W03 assume a
radially constant M/L, the comparison should be taken with some caution. The
same applies when comparing our findings to those of Díaz et al. (2018). The
authors constructed a detailed model of the M31 bulge, including the
contribution from different dynamical structures, such as a “classical” bulge,
a boxy/peanut bulge, and a thin bar, concluding that the M/L of the two
(classical and boxy/peanut) bulge components is consistent with that expected
for a Chabrier IMF ($\rm M/L_{r}\sim 4$), under the assumption of a radially
constant M/L. Note that, indeed, our $\rm\langle M/L_{r}\rangle$ is very
similar to that one would infer for a Kroupa-like IMF (i.e. $M/L_{r}\sim 4$),
as the innermost bins, where we infer a bottom-heavier IMF, give only a minor
contribution to the integrated light (mass) of the bulge.

Dutton et al. (2013) used gravitational lensing to constrain the stellar IMF
normalization (i.e. the stellar M/L) in the bulges of massive spiral galaxies,
finding evidence for a Salpeter-like normalization. We point out that this
result is not in disagreement with our findings for the bulge of M31, as the
bulges analyzed by Dutton et al. (2013) are in the velocity dispersion range
of $200$–$250$ km s-1, while the sigma of the M31 bulge is $150$–$160$ km s-1.
Considering the relation between IMF slope and velocity dispersion of ETGs,
for a velocity dispersion of 160 km s-1, the expected $M/L_{r}$ is $3.5\pm 1$
(for a low-mass tapered IMF; see fig. 12 of LB13), fully consistent with our
luminosity-weighted estimate for the bulge of M31.

The presence of a radially varying IMF, and thus the increase of $\rm M/L_{r}$
in the center of the bulge, might also affect dynamical models of the M31
nucleus (at $R<1$”), consisting of an eccentric stellar disk orbiting a
central black hole (BH). Assuming a negligible contribution ($<1\%$) of the
bulge component to the nucleus, Peiris & Tremaine (2003) performed a detailed
dynamical model, obtaining a BH mass estimate of $\sim 1.0\times 10^{8}$ $\rm
M_{\odot}$, by a factor of $\sim 2$ larger (but still within the
uncertainties) compared to the estimate of $6\times 10^{7}\,\rm M_{\odot}$
from the correlation of BH mass and bulge dispersion of nearby galaxies
(Tremaine, 2002). Indeed, an increase of $\rm M_{\star}/L_{r}$ by a factor of
$\sim 2$ in the nucleus, as predicted by our results (Fig. 6), would still
imply a negligible contribution of the bulge mass to the nuclear region, and
thus not decrease significantly the BH mass estimate.

Indeed, dynamical models allowing for radial M/L variations would help to
further scrutinize our results, and test if radial IMF variations are fully
consistent with the bulge kinematics, and how/if they affect the inferred
properties (such as the BH mass) of the M31 nucleus. We postpone this analysis
to a forthcoming contribution.

###  6.3 Bottom-heavy or Na-enhanced?

The radial behaviour of IMF-sensitive features for the bulge of M31 has been
also analyzed by Zieleniewski et al. (2015, hereafter Z15), who found a
negative radial gradient for $\rm NaI8190$ and $\rm CaT$, decreasing from a
radial distance of a few arcsec out to a few hundred arcsec, in agreement with
our results (see panels h and i in Fig. 2). Z15 compared the radial behaviour
of $\rm NaI8190$ and $\rm NaD$ for M31, with SSP model predictions from
Vazdekis et al. (2012, hereafter V12) (i.e. an earlier version of EMILES
models) and CvD12a model predictions. They concluded that CvD12a models
predict a strong dependence of $\rm NaI8190$ on IMF, compared to $\rm NaD$,
while for V12 models, also $\rm NaD$ depends significantly on IMF (see also
LB13). For this reason, CvD12a models would favour a MW-like IMF at all radial
positions in M31, with an extreme $\rm[Na/Fe]$ abundance ratio in the center
(up to 1 dex); while V12 models would suggest a bottom-heavy IMF in this
region, and a Salpeter-like IMF for the rest of the bulge. We note that the
comparison performed by Z15 suffered from the lack of predictions with varying
$\rm[Na/Fe]$ for V12 models, and relied on $\rm NaI8190$ and $\rm NaD$
measurements drawn from different sources (with $\rm NaD$ from Davidge 1991,
Davidge 1997). Indeed, using Na-enhanced stellar population models, and an
homogeneous spectroscopic data-set, our analysis shows that for most of the
bulge radial extent, all optical and NIR spectral features of M31 are well
matched with a Kroupa-like IMF, and an $\rm[Na/Fe]$ abundance ratio of $\sim
0.4$ dex.

###  6.4 M31 bulge vs. massive ETGs

Based on the $\rm TiO2$ spectral index, Martín-Navarro et al. (2015b,
hereafter MN15b) found a tight correlation between $\rm[Z/H]$ and IMF slope
for ETGs in the CALIFA spectroscopic survey (Sánchez et al., 2012), pointing
to stellar metallicity as a possible (local) driver of IMF variations in
massive ETGs, in agreement with what also found by van Dokkum et al. (2017).
However, in LB19, we showed that very massive ETGs, mostly brightest cluster
galaxies (BCGs), do not follow the IMF-metallicity relation. Martín-Navarro et
al. (2019) also showed that in a disk-like structure, IMF variations do not
follow metallicity changes. Hence, metallicity does not seem to be the only
culprit of IMF variations in galaxies. Our results for M31 further support
this conclusion. Indeed, in the bulge innermost radial bins, where stellar
metallicity is very high (up to $\rm[Z/H]$$\sim 0.4$), we find only a mildly
bottom-heavy IMF ($\rm\Gamma_{b}$$\sim 2.5$; see Fig. 5), while based on the
IMF-metallicity relation (see, e.g., fig. 2 of MN15b), one would expect
$\rm\Gamma_{b}$$\sim 3$ or above. Note that based on $\rm TiO1$ and $\rm
TiO2$, one would actually infer a very bottom-heavy IMF also for the innermost
bins of M31, as these indices increase very steeply in the innermost few
arcsec (see Fig. 3). However, as discussed in App. D, the increase of $\rm
TiO2$ (and $\rm TiO1$) with $\rm[Z/H]$ in the center of M31, is most likely
due to a genuine dependence of the TiO’s on metallicity in the high-
metallicity regime, rather than an IMF variation. This highlights the
importance of using a wide set of spectral features, throughout a large
wavelength baseline, when constraining the IMF.

In order to gain further insights into the origin of IMF variations, in Fig. 8
we compare the estimates of $\rm\Gamma_{b}$ for the bulge of M31 with those
for very massive ETGs from LB19, as a function of galactocentric distance (top
panel), and logarithmic stellar mass density ($\rm log(\Sigma)$; see bottom
panel). The data for M31 (red crosses) cover a region out to $\sim 0.8$ kpc,
and thus can be compared only with the innermost data-points for galaxies in
LB19. To perform a direct comparison, we have computed, under the assumption
of circular symmetry, the luminosity-weighted value of $\rm\Gamma_{b}$ for the
entire bulge, out to the maximum galactocentric distance ($r\sim 200$”)
covered by our analysis. The red dot in the top panel of the Figure, with
error bars, shows the luminosity-weighted IMF slope, i.e.
$\rm\Gamma_{b}$$=1.7\pm 0.4$. Note that the integrated IMF slope is only
marginally above the Kroupa-like value of $\rm\Gamma_{b}$$=1.3$, as the
innermost bins, where we infer a bottom-heavier IMF, give only a minor
contribution to the integrated light (mass) of the bulge. Fig. 8 shows that
the integrated IMF slope for the bulge is significantly lower compared to that
for massive ETGs, within the same galactocentric distance. Hence, observing a
spiral galaxy like M31 with the same spatial resolution as for ETGs, one would
infer an IMF very similar to a Kroupa-like distribution. Also, although
metallicity is very high in the center of the bulge, the IMF slope is $\sim
2.5$, while all massive ETGs tend to have higher $\rm\Gamma_{b}$($>2.6$). Note
also that different radial bins for the M31 bulge cover a wide range in
stellar mass density (see bottom panel of the Figure), allowing a direct
comparison with different radial bins for massive ETGs. At fixed local
density, the M31 spectra reveal a significantly lower $\rm\Gamma_{b}$ compared
to ETGs. The reason why the IMF is not very bottom-heavy in the center of M31,
despite the high metallicity, might be its “low” velocity dispersion ($\sim
160$ km s-1), compared to that of massive galaxies ($\sim 300$ km s-1), or
likely related to the formation mechanism itself of the bulge. From the
relation between IMF slope and galaxy velocity dispersion of LB13 (see their
fig. 12), ETGs with $\sim 160$ km s-1 are expected to have an “integrated”
$\rm\Gamma_{b}$ consistent with a Kroupa-like distribution, as we observe,
indeed, for the bulge of M31. Note that, as shown in Tamm et al. (2012, see
their fig. 3), the nucleus gives a prominent contribution to the M31 surface
brightness profile only at radial distances <$1--2$”, while we find IMF radial
variations over a larger radial range ($<10$”; see, e.g., Fig. 5). Hence, the
IMF variations are not driven by the nucleus of M31, but do actually reflect
the properties of the bulge itself. Indeed, as shown in Blaña et al. (2017),
the bulge of M31 consists of a “classical” and a “box/peanut” (BP) component,
the latter contributing by $\sim$ 2/3 of the bulge total stellar mass. The BP
component is believed to form later from the disk material (possibly with a
“normal” IMF), through the buckling instability of the bar. Therefore, one may
speculate that any genuine IMF variation in the classical bulge might have
been actually diluted by the stellar material formed in the BP component.
Detailed dynamical models, including the effect of a varying IMF, might help
to address this point. The comparison shown in Fig. 8 reinforces the idea that
there is no single parameter, such as stellar metallicity or local density,
able to explain IMF variations in stellar systems (see LB19).

Figure 8: The IMF slope, $\rm\Gamma_{b}$, of M31 is compared to the IMF radial
profiles of massive ETGs, as a function of radial distance to galactocentric
distance (top), and logarithmic stellar mass density (bottom). In both panels,
red crosses plot the IMF slope values for different radial bins of M31, while
other symbols, with different colours, plot the IMF radial profiles for the
massive “XSG” ETGs from LB19 (see labels on the upper–right). In the top
panel, the red dot with error bars is the expected M31 IMF slope obtained by
mimicking a circular aperture with average radius of 0.3 kpc, including the
contribution from all radial bins (see the text for details). Note that the
circularized value for M31 can be compared to the the innermost radial points
of massive ETGs. The y-axis on the right of each panel shows the fraction of
low-mass (<$0.5\,M_{\odot}$) stars corresponding to a given value of
$\rm\Gamma_{b}$ (see LB13). In the top panel, the slope value for a Kroupa-
like IMF is marked by the green horizontal dashed line.

##  7 Summary and Conclusions

We have acquired new, long-slit spectroscopy with OSIRIS at GTC, along the
major axis of the bulge of M31, obtaining an homogeneous, continuous spectral
coverage in the wavelength range from $\sim 0.35$ to $\sim 1\mu$m, with a
constant spectral resolution of $\sim 38$ km s-1(FWHM). Our dataset is the
most homogeneous obtained so far (in terms of spectral coverage and
resolution) along the bulge major axis. We fit a wide set of spectral indices,
including Ca, Na, Mg, and TiO spectral features, with state-of-the-art stellar
population models (Na–MILES), that take into account the coupled effect of
varying IMF and Na abundance ratios.

Our analysis shows that the bulge of M31 has an IMF consistent with a Kroupa-
like, or “slightly” above Kroupa, IMF (with IMF slope $\rm\Gamma_{b}$$<2$) at
all radial positions but in the innermost 10” (corresponding to $40$ pc),
where a mildly bottom-heavy IMF (with a normalization close to the Salpeter
distribution) is required. We are able to match both the $\rm NaD$ and $\rm
NaI8190$ spectral indices simultaneously, without requiring extremely high
$\rm[Na/Fe]$, in contrast with previous studies. The best-fitting $\rm[Na/Fe]$
is $\sim 0.4$ dex, rising up to $\sim 0.6$ dex, in the inner few arcsec. The
Na–enhanced models allow us to fit extremely well all IMF-sensitive spectral
indices of M31, as well as a large set of Lick-like indices, once the effect
of other individual abundance ratios is taken into account (using CvD12a SSP
models). Compared to very massive ETGs (with a velocity dispersion $\gtrsim
300$ km s-1), the bulge of M31 has a less bottom-heavy IMF at fixed stellar
mass density, and a less bottom-heavy IMF in the center, despite having
similarly high stellar metallicity. Significant IMF variations in the bulge
are concentrated within a much smaller region ($<100$ pc) compared to most
massive ETGs, where the IMF changes occur at galactocentric distances larger
than $\sim 1$ kpc. The integrated, luminosity-weighted, IMF slope is
consistent with that expected for a Milky-Way-like distribution, and
consistent with lower-mass ETGs, having the same velocity dispersion of the
bulge ($\sigma\sim 150$–$160$ km s-1). The IMF radial gradient of the M31
bulge might originate from the presence of different dynamical structures
(such as a classical and a box/peanut bulge), and IMF radial gradients might
actually help to constrain the formation and evolution scenarios of different
sub-systems. In the future, this point should be addressed with detailed
dynamical models, including the effect of a radially varying IMF. Our results
reinforce the idea that there is no single driver of IMF variations in
galaxies, and we cannot entirely ascribe IMF variations to a single parameter,
such as stellar metallicity, or observed mass density. The IMF gradients we
observe at $z\sim 0$ likely result from a complex multi-component physical
process, related to the details of the physical properties of the gas at the
time of formation, as well as the assembly history itself of a galactic
system.

## Acknowledgments

We thank the anonymous referee for his/her helpful comments. Observations of
M31 were carried out with the Gran Telescopio CANARIAS (GTC), proposal ID:
GTC39-17A (P.I. A. Vazdekis). The authors thank dr. S. Geier for the kind help
he provided for the preparation of the phase2 material, and to perform the
observations of M31. FLB acknowledges the Instituto de Astrofísica de Canarias
for the kind hospitality when this project started. FLB and AV acknowledge
support from grant PID2019-107427GB-C32 from Ministerio de Ciencia e
Innovación. FLB and AP acknowledge financial support from the INAF PRIN
1.05.01.85.11.

## Data Availability

The E-MILES SSP models are publicly available at the MILES website
(http://miles.iac.es). The Na-MILES models are available at the same website
(under the section "Other predictions/data"). The CvD12a SSP models are
available upon request to the authors (see
https://scholar.harvard.edu/cconroy/projects). Data for the LB19 sample
(programmes 092.B-0378, 094.B-0747, 097.B-0229; PI: FLB) can be downloaded
from the ESO archive (http://archive.eso.org). Raw spectroscopic data for the
bulge of M31 (programme GTC39-17A) are available from the GTC Public Archive
(http://gtc.sdc.cab.inta-csic.es/gtc/jsp/searchform.jsp). The reduced
radially-binned spectra of M31 are available upon request to FLB.

## References

  * Alton, Smith, Lucey (2017) Alton, P. D., Smith, R. J., Lucey, J. R., 2017, MNRAS, 468, 1594
  * Alton, Smith, Lucey (2018) Alton, P. D., Smith, R. J., Lucey, J. R., 2018, MNRAS, 478, 4464
  * Auger et al. (2010) Auger, M. W., Treu, T., Bolton, A. S., Gavazzi, R., Koopmans, L. V. E., Marshall, P. J., Moustakas, L. A., Burles, S., 2010, ApJ, 724, 511
  * Blaña et al. (2017) Blaña, M., Wegg, C., Gerhard, O., Erwin, P., Portail, M., Opitsch, M., Saglia, R., Bender R., 2017, MNRAS, 466, 4279
  * Barnabé et al. (2011) Barnabé, M., Czoske, O., Koopmans, L.V.E., Treu, T., Bolton, A.S., 2011, MNRAS, 415, 2215
  * Bensby et al. (2017) Bensby, T., et al., 2017, A&A, 605, 89
  * Delisle & Hardy (1992) Delisle, S., Hardy, E., 1992, AJ, 103, 711
  * Díaz et al. (2018) Díaz, M. B.; Gerhard, O., Wegg, C., Portail, M., Opitsch, M., Saglia, R., Fabricius, M., Erwin, P., Bender, R., 2018, MNRAS, 481, 3210
  * Domínguez et al. (2019) Domínguez Sánchez, H., Bernardi, M., Brownstein, J. R., Drory, N., Sheth, R. K., 2019, MNRAS, 489, 5612
  * Cappellari & Emsellem (2004) Cappellari, M., Emsellem, E., 2004, PASP, 116, 138
  * Cappellari et al. (2012) Cappellari M., et al., 2012, Natur, 484, 485
  * Cappellari et al. (2013) Cappellari M., et al., 2013, MNRAS, 432, 1862
  * Cappellari (2017) Cappellari, M., 2017, MNRAS, 466, 798
  * Carter et al. (1986) Carter, D., Visvanathan, N., Pickles, A.J., 1986, ApJ, 311, 637
  * Cassisi et al. (2000) Cassisi, S., Castellani, V., Ciarcelluti, P., Piotto, G., Zoccali, M., 2000, MNRAS, 315, 679
  * Cenarro et al. (2001a) Cenarro, A. J., Cardiel, N., Gorgas, J., Peletier, R. F., Vazdekis, A., Prada, F., 2001a, MNRAS, 326, 959
  * Cenarro et al. (2003) Cenarro, A. J., Gorgas, J., Vazdekis, A., Cardiel, N., Peletier, R. F.,
  * Cervantes & Vazdekis (2009) Cervantes, J. L., Vazdekis, A., 2009, MNRAS, 392, 691
  * Cohen (1978) Cohen, J. G., 1978, ApJ, 221, 788
  * Conroy & van Dokkum (2012a) Conroy, C., van Dokkum, P., 2012a, ApJ, 747, 69 (CvD12a)
  * Conroy & van Dokkum (2012b) Conroy, C., van Dokkum, P., 2012b, ApJ, 760, 71
  * Cushing et al. (2005) Cushing, M. C., Rayner, J. T., Vacca, W. D., 2005, ApJ, 623, 1115
  * Davidge (1991) Davidge, T. J., 1991, AJ, 101, 884
  * Davidge (1997) Davidge, T. J., 1997, AJ, 113, 985
  * Dutton, Mendel & Simard (2012) Dutton, A. A., Mendel, J. T., Simard, L., 2012, MNRAS, 422, 33
  * Dutton et al. (2013) Dutton, A. A., et al., 2013, MNRAS, 428, 3183
  * de La Rosa et al. (2011) de La Rosa I. G., La Barbera F., Ferreras I., de Carvalho R. R., 2011, MNRAS, 418, L74
  * Faber & French (1980) Faber S. M., French H. B., 1980, ApJ, 235, 405
  * Ferreras et al. (2005) Ferreras, I., Saha, P., Williams, L. L. R., 2005, ApJ, 623, 5
  * Ferreras et al. (2008) Ferreras, I., Saha, P., Burles, S., 2008, MNRAS, 383, 857
  * Ferreras et al. (2010) Ferreras, I., Saha, P., Leier, D., Courbin, F., Falco, E. E. 2010, MNRAS, 409, L30
  * Ferreras et al. (2013) Ferreras, I., La Barbera, F., de la Rosa, I. G., Vazdekis, A., de Carvalho, R. R., Falcón-Barroso, J., Ricciardelli, E., 2013, MNRAS, 429, L15 (F13)
  * Fontanot et al. (2017) Fontanot F., De Lucia G., Hirschmann M., Bruzual G., Charlot S., ZibettiS., 2017, MNRAS, 464, 3812
  * Fontanot et al. (2018) Fontanot F., La Barbera F., De Lucia G., Pasquali A., Vazdekis A., 2018, MNRAS, 479, 5678
  * Gallazzi et al. (2021) Gallazzi, A. R., Pasquali, A., Zibetti, S., La Barbera, F., MNRAS, 2021, MNRAS, in press, arXiv201004733G
  * Girardi et al. (2000) Girardi, L., Bressan, A., Bertelli, G., Chiosi, C., 2000, A&AS, 141, 371
  * Hopkins (2018) Hopkins A. M., 2018, PASA, 35, 39
  * Jeong et al. (2013) Jeong, H., Yi, S. K., Kyeong, J., Sarzi, M., Sung, E.-C., Oh, K., 2013, ApJS, 208, 7J
  * Kausch et al. (2015) Kausch, W., et al., 2015, A&A, 576, 78
  * Kroupa (2001) Kroupa, P., 2001, MNRAS, 322, 231
  * La Barbera et al. (2013) La Barbera, F., Ferreras, I., Vazdekis, A., de la Rosa, I. G., de Carvalho, R. R., Trevisan, M., Falcón-Barroso, J., Ricciardelli, E., 2013, MNRAS, 433, 3017 (LB13)
  * La Barbera, Ferreras, & Vazdekis (2015) La Barbera F., Ferreras I., Vazdekis A., 2015, MNRAS, 449, L137
  * La Barbera et al. (2016) La Barbera F., Vazdekis A., Ferreras I., Pasquali A., Cappellari M., Martín-Navarro I., Schönebeck F., Falcón-Barroso J., 2016, MNRAS, 457, 1468 (LB16)
  * La Barbera et al. (2017) La Barbera F., Vazdekis A., Ferreras I., Pasquali A., Allende Prieto, C., Rc̈k, B., Aguado, D. S., Peletier, R. F., 2017, MNRAS, 464, 3597 (LB17)
  * La Barbera et al. (2019) La Barbera F., et al., 2019, MNRAS, 489, 4090 (LB19)
  * Leier et al. (2016) Leier D., Ferreras I., Saha P., Charlot S., Bruzual G., La Barbera F., 2016, MNRAS, 459, 3677
  * Lyubenova et al. (2016) Lyubenova, M., Martín-Navarro, I., van de Ven, G., et al., 2016, MNRAS, 463, 3220
  * Maraston (2005) Maraston, C., 2005, MNRAS, 362, 799
  * Martín-Navarro et al. (2015a) Martín-Navarro I., La Barbera F., Vazdekis A., Falcón-Barroso J., Ferreras I., 2015a, MNRAS, 447, 1033 (MN15a)
  * Martín-Navarro et al. (2015b) Martín-Navarro I., et al., 2015b, ApJ, 806, L31 (MN15b)
  * Martín-Navarro et al. (2019) Martín-Navarro I., et al., 2019, MNRAS, in press (arXiv:1903.10514)
  * McConnachie et al. (2005) McConnachie, A. W., Irwin, M. J., Ferguson, A. M. N., Ibata, R. A., Lewis, G. F., Tanvir, N., 2005, MNRAS, 356, 979
  * Noll et al. (2014) Noll, S., Kausch, W., Kimeswenger, S., Barden, M., Jones, A. M., Modigliani, A., Szyszka, C., Taylor, J., 2014, A&A, 567, 25
  * O’Connell (1976) O’Connell, R. W., 1976, ApJ, 206, 370O
  * Opitsch et al. (2018) Opitsch, M., Fabricius, M. H.; Saglia, R. P., Bender, R., Blaña, M., Gerhard, O., 2018, A&A, 611, 38O
  * Oser et al. (2010) Oser L., Ostriker J. P., Naab T., Johansson P. H., Burkert A., 2010, ApJ, 725, 2312
  * Parikh et al. (2018) Parikh, T., et al., 2018, MNRAS, 477, 3954P
  * Pietrinferni et al. (2004) Pietrinferni, A., Cassisi, S., Salaris, M., Castelli, F., 2004, ApJ, 612, 168
  * Peiris & Tremaine (2003) Peiris, H. V., Tremaine, S., 2003, ApJ, 599, 237
  * Tremaine (2002) Tremaine, S., 2002, ApJ, 574, 740
  * Rayner et al. (2009) Rayner, J. T., Toomey, D. W., Onaka, P. M., Denault, A. J., Stahlberger, W. E., Vacca, W. D., Cushing, M. C., Wang, S., 2003, PASP, 115, 362
  * Röck et al. (2016) Röck, B., Vazdekis, A., Ricciardelli, E., Peletier, R. F., Knapen, J. H., Falcón-Barroso, J., 2016, A&A, 589, 73
  * Saglia et al. (2010) Saglia, R. P., et al., 2010, A&A, 509, A61 (S10)
  * Saglia et al. (2018) Saglia, R. P., Opitsch, M., Fabricius, M. H., et al. 2018, A&A, 618, 156
  * Sánchez-Blázquez et al. (2006) Sánchez-Blázquez, P., et al., 2006, MNRAS, 371, 703
  * Sánchez et al. (2012) Sánchez, S. F., Kennicutt, R. C., Gil de Paz, A., et al., 2012, A&A, 538, A8
  * Sarzi et al. (2006) Sarzi, M., Falcón-Barroso, J., Davies, R. L., et al., 2006, MNRAS, 366, 1151
  * Sarzi et al. (2018) Sarzi, M., Spiniello, C., La Barbera, F., Krajnović, D., van den Bosch, R.,
  * Schiavon et al. (1997b) Schiavon, R., Barbuy, B., Rossi, S. C. F., Milone, A., 1997, ApJ, 479, 902
  * Serven (2005) Serven, J., Worthey, G., Briley, M. M., 2005, ApJ, 627, 754
  * Smith & Lucey (2013) Smith, R. J., Lucey, J. R., 2013, MNRASL, 434, 1964
  * Smith, Lucey & Conroy (2015) Smith, R. J., Lucey, J. R., Conroy, C., 2015, MNRASL, 449, 3441
  * Smette et al. (2015) Smette, A., et al., 2015, A&A, 576, 77
  * Spiniello et al. (2014) Spiniello, C., Trager, S. C., Koopmans, L. V. E., Conroy, C., 2014, MNRAS, 438, 1438
  * Spinrad & Taylor (1971) Spinrad H., Taylor B. J., 1971, ApJS, 22, 445
  * Tamm et al. (2012) Tamm, A., Tempel, E., Tenjes, P., Tihhonova, O., Tuvikene, T., 2012, A&A, 645, 4
  * Thomas et al. (2003a) Thomas, D., Maraston, C., Bender, R., 2003, MNRAS, 339, 897
  * Thomas et al. (2004) Thomas D., Maraston C., Korn A., 2004, MNRAS, 351, L19
  * Thomas et al. (2005) Thomas, D., Maraston, C., Bender, R., Mendes de Oliveira, C., 2005, ApJ, 621, 673
  * Thomas & Davies (2006) Thomas D., Davies, R. L., 2006, MNRAS, 366, 510
  * Thomas, Maraston & Johansson (2011) Thomas, D., Maraston, C., Johansson, J., 2011, MNRAS, 412, 2183
  * Tortora et al. (2013) Tortora, C., Romanowsky, A. J., Napolitano, N. R., 2013ApJ, 765, 8
  * Trager et al. (1998) Trager, S. C., Worthey, G., Faber, S. M., Burstein, D., González, J. J., 1998, ApJS, 116, 1
  * Treu et al. (2010) Treu T., Auger M. W., Koopmans L. V. E., Gavazzi R., Marshall P. J., Bolton A. S., 2010, ApJ, 709, 1195
  * Valdes et al. (2004) Valdes, F., Gupta, R., Rose, J. A., Singh, H. P., Bell, D. J., 2004, ApJS, 152, 251
  * van Dokkum & Conroy (2010) van Dokkum P. G., Conroy C., 2010, Nature, 468, 940
  * van Dokkum et al. (2017) van Dokkum, P., Conroy, C., Villaume, A., Brodie, J., Romanowsky, A. J., 2017, ApJ, 841, 68V
  * Vazdekis et al. (1996) Vazdekis, A., Casuso, E., Peletier, R. F., Beckman, J. E., 1996, ApJS,
  * Vazdekis et al. (2003) Vazdekis, A., Cenarro, A. J., Gorgas, J., Cardiel, N., Peletier, R. F., 2003, MNRAS, 340, 1317
  * Vazdekis et al. (2012) Vazdekis, A., Ricciardelli, E., Cenarro, A. J., Rivero-González, J. G., Díaz-Garcá, L. A., Falcón-Barroso, J., 2012, MNRAS, 424, 157
  * Vazdekis et al. (2015) Vazdekis, A., Coelho, P., Cassisi, S., Ricciardelli, E., Falcón-Barroso, J., Sánchez-Blázquez, P., La Barbera, F., Beasley, M. A., Pietrinferni, A., 2015, MNRAS, 449, 1177 (V15)
  * Vazdekis et al. (2016) Vazdekis, A., Koleva, M., Ricciardelli, E., Röck, B., Falcón-Barroso, J., 2016, MNRAS, 463, 3409
  * Wegner et al. (2012) Wegner, G.A., Corsini, E.M., Thomas, J., Saglia, R.P., Bender, R., Pu, S.B., 2012, MNRAS, 144, 78
  * Whitford (1977) Whitford, A. E., 1977, ApJ, 211, 527
  * Widrow, Perrett, Suyu (2003) Widrow, L. M., Perrett, K. M., Suyu, H. S., 2003, ApJ, 588, 311
  * Zieleniewski et al. (2015) Zieleniewski S., Houghton R. C. W., Thatte N., Davies R. L., 2015, MNRAS, 452, 597 (Z15)
  * Zieleniewski et al. (2017) Zieleniewski S., Houghton R. C. W., Thatte N., Davies R. L., Vaughan S. P., 2017, MNRAS, 465, 192

##  Appendix A Emission corrections

Absorption indices analyzed in the present work are corrected for
contamination by (gas) emission lines. Besides the Balmer lines, $\rm H\beta$
and $H\gamma$, contamination affects also the Fe5015 and Mgb absorptions,
because of emission from the [OIII]$\lambda\lambda 4959,5007$ and
[NI]$\lambda\lambda 5198,5200$ doublets, respectively (see e.g. Sarzi et al.
2006 and references therein). Note that, in the present work, we correct for
the effect of emission on absorption features, while we postpone a detailed
analysis of emission lines in the M31 OSIRIS spectra to a forthcoming paper.
The correction is performed as follows:

(i)

$\rm H\beta_{o}$ and $\rm H\gamma_{F}$ are corrected with the same approach as
in LB13 and LB16. For each radially binned spectrum, and each Balmer line, we
perform spectral fitting on a 200 Å-wide region around the line, with a 2SSP
model 141414The fitting allows for a variation of the ages and metallicities,
as well as the relative light-fraction, of the two SSPs., by excluding the
trough of the line, and fitting a Gaussian function to model the emission in
the residual spectrum. The Gaussian fit is removed from the input spectrum to
estimate the emission-corrected line-strength. Uncertainties on emission
correction are estimated by taking into account noise in the input spectrum,
by varying the degree of the polynomial modeling the continuum in the spectral
fitting, as well as by varying the adopted IMF (i.e. the $\rm\Gamma_{b}$) of
the 2SSP models. The corrections turned to be significant (at more than 2.5
sigma level) for all radial bins, amounting, on average, to $\sim 0.3$ Å and
$\sim 0.1$ Å for $\rm H\beta_{o}$ and $\rm H\gamma_{F}$, respectively.

(ii)

To correct Fe5015 and $\rm Mgb5177$, we performed spectral fitting with pPXF
over the spectral region from $\lambda\sim 4800$ to $\sim 5300$ Å, by
excluding regions potentially affected by emission lines, i.e. regions of
$\pm$100 km s-1 around the $\rm H\beta$$\lambda 4861$ line, and around each
line of the [OIII]$\lambda\lambda 4959,5007$ and [NI]$\lambda\lambda
5198,5200$ doublets. We modeled emission lines in the residual spectra with
Gaussian functions, and removed the Gaussian profiles from the input spectra,
to measure the emission-corrected line-strengths. On average, the corrections
amount to $\sim-0.1$ Å, and $+\sim 0.3$ Å, for $\rm Mgb5177$ and Fe5015,
respectively. Note that this method provides also an emission correction to
$\rm H\beta_{o}$ that is similar ($\sim 0.05$ Å lower) to that obtained from
method (i). However, to be consistent with our previous works (e.g. LB19), we
decided to use corrections on $\rm H\beta_{o}$ and $\rm H\gamma_{F}$, from
individual fits around each line, as described above (method i). Our values of
the correction for emission contamination are roughly consistent with those of
$0.29$ Å, $-0.07$ Å and $0.4$ Å, reported by S10 for $\rm H\beta$ (rather than
$\rm H\beta_{o}$), $\rm Mgb5177$, and Fe5015, respectively.

##  Appendix B Correction of $\rm NaD$ for ISM absorption

The NaD stellar absorption is potentially affected by absorption from inter-
stellar matter (ISM), both in the source restframe, as well as, for nearby
systems, in our Galaxy. Although contamination can be severe for late-type
galaxies (Jeong et al., 2013), none of the previous studies has estimated the
effect for the bulge of M31 (e.g. Zieleniewski et al. 2015). Thanks to the
high quality of our OSIRIS spectroscopy, we are able to see the signature of
ISM absorption on the profile of the NaD lines of M31. Fig. 9 illustrates, for
the innermost spectrum of M31, the procedure we adopt to correct the
equivalent width of NaD for ISM absorption. The observed spectrum of M31 is
fitted with pPXF, over a region $\sim 800$ Å–wide around the NaD, as shown by
the black and magenta curves in the top panel of the Figure. The middle and
bottom panels show a zoom-in of the NaD absorption (middle), and the profile
obtained by normalizing the NaD absorption with the pPXF best-fit solution
(bottom). Three narrow “residual absorptions” are detected at $\lambda\sim
5890$, $5895$, and $5902$ Å, roughly corresponding to the positions of the D2
and D1 lines ($\sim 5890$ and $5895$ Å, respectively) in the restframe of M31
(see the vertical thin green segments in the bottom panel), and at the
heliocentric velocity of M31 (about $-300$ km s-1; see the thick green
segments in the bottom panel). Note that in order to minimize the effect of
ISM absorption on the best-fit of the stellar component, we run pPXF by
masking out regions $\pm 100$ km s-1–wide around the three residual
absorptions shown in the middle panel of Fig. 9 (see green lines in the middle
panel). We fit the residual absorptions with three Gaussian functions, and
correct the input observed spectrum by dividing it with the best-fit
Gaussians. The difference of NaD line-strength between the corrected and input
spectra give the NaD correction for ISM absorption, $\rm\delta(NaD)$.

The value of $\rm\delta(NaD)$ is negative, as expected by the fact that ISM
contamination makes the NaD stronger with respect to that of the pure stellar
component. Fig. 10 plots $\rm\delta(NaD)$ as a function of radial distance to
the center of M31. To explore possible systematics in the correction, we
repeated the entire procedure by running pPXF with two different sets of
templates, i.e. MILES stars as well as Na-EMILES SSP models. Also, since the
ISM absorption at $\lambda\sim 5890$ Å is much weaker than the others, and is
barely seen in some spectra (especially for the outer, lower S/N, radial
bins), we repeated the procedure by fitting either three or two Gaussian
functions (in the latter case only the residual absorptions at $\lambda\sim
5895$ and $5902$ Å were fitted). Different methods are shown with different
colours, symbols, and line types in Fig. 10, with average values from
different methods (see the labels) marked with horizontal lines. For all
methods, the correction is significant, and amounts to about $-0.2$ Å. Note
that this value is negligible with respect to the $\rm NaD$ radial gradient as
measured in the present work (see panel g of Fig. 2). As a conservative
approach we averaged out all results shown in Fig. 10, and applied a constant
correction of $-0.215\pm 0.06$Å to all NaD line-strengths analyzed in the
present work 151515In general, the correction procedure is hampered by the
fact that NaD airglow lines do also overlap to the MW ISM absorption lines.
This contribution is negligible in the innermost region, where the signal from
the galaxy is prominent, while it becomes more important in the outermost
bins, hampering a precise estimate of the radial variation of
$\rm\delta(NaD)$. .

Figure 9: Example of the procedure to correct NaD lines for ISM absorption,
for the innermost spectrum of M31. Top panel: observed (black) and pPXF best-
fit (magenta) spectrum, in the wavelength range from $\sim 5600$ to $\sim
6400$ Å. Medium panel: the same as in the top panel but zooming into the
spectral region around the NaD. The green curves mark regions potentially
affected by ISM absorption (see below), that are masked out when running pPXF.
Bottom panel: the black curve plots the observed spectrum normalized to the
pPXF best-fit around the NaD lines. Note the three absorptions at $\lambda\sim
5890$, $5895$, and $5902$ Å, respectively. The thin and thick vertical
segments, labeled as “M31 ISM” and “MW ISM”, mark the expected positions of
the NaD doublet (D2 and D1 lines) in the restframe of M31, and at the
heliocentric velocity of M31 (about $-300$ km s-1), respectively. Regions
potentially affected by ISM absorption (middle panel) are defined by windows
$\pm 100$ km s-1–wide around the positions of the green segments. Red and blue
curves show the multi-Gaussian best-fit to the ISM absorptions and the
normalized residuals, respectively.  Figure 10: Correction to the NaD line-
strength for ISM absorption, $\delta(NaD)$, as a function of galacto-centric
distance, $\rm R$, in units of arcsec. Different symbols and lines correspond
to different methods to compute the correction, as labeled in the upper-right
(see the text for details). Horizontal dot-dashed lines mark the median values
of $\delta(NaD)$ from different methods.

##  Appendix C Age, $\rm[Z/H]$, and $\rm[Mg/Fe]$ grids

Fig. 11 illustrates the estimate of age and metallicity for the M31 spectra
through the $\rm H\beta_{o}$–$\rm[MgFe]^{\prime}$ diagram, using different SSP
models, i.e. base EMILES models with BaSTI and Padova isochrones 161616 Note
that we refer to SSP models based on BaSTI isochrones as “Teramo” models.
(dark- and light-grey grids in the top panel), and $\alpha$–MILES scaled-solar
(green) and $\alpha$-enhanced (red) models in the bottom panel. Since around
solar metallicity, the BaSTI isochrones are hotter than the Padova ones (see
Pietrinferni et al. 2004; V15), Teramo base models are shifted towards higher
$\rm H\beta_{o}$ values than Padova (see the top panel of the Figure). Hence,
Teramo models tend to provide older ages (by $\sim 1$–$2$ Gyr) than Padova
models, as shown in Fig. 4. However, when using $\alpha$–MILES models (see
bottom panel of the Figure), the effect of the isochrones is counteracted by
the fact that $\rm H\beta_{o}$ is expected to decrease with $\rm[Mg/Fe]$ (see
the red and green grids in the bottom panel). Since the M31 spectra are
$\alpha$-enhanced, ages obtained with $\alpha$–MILES models do not differ
significantly from those obtained with Padova base models, as shown in Fig. 4.
We note that the metallicities obtained with Teramo (base) models are also
slightly higher than those for Padova, as the two sets of models differ for
the reference solar metallicity (see V15).

For what concerns the estimates of $\rm[Mg/Fe]$ (bottom panel of Fig. 4), we
find a decrease of $\rm[Mg/Fe]$ in the center of the bulge, which is more
pronounced for $\alpha$–MILES models, and barely seen with base models (where
we use the “proxy” for $\rm[Mg/Fe]$; see LB13 and V15). Fig. 12 plots the
$\rm\langle Fe\rangle$ index as a function of $\rm Mgb5177$. In this diagram
(see e.g. Thomas et al. 2003a), the effect of metallicity and age are almost
parallel, pushing the line-strength of both indices to higher values. On the
other hand, increasing $\rm[Mg/Fe]$ tends to increase $\rm Mgb5177$ and
decrease $\rm\langle Fe\rangle$, so that the effect of $\rm[Mg/Fe]$ is
orthogonal to that of age and metallicity. As seen in the Figure, the radial
gradient for the M31 data-points is along a direction of decreasing
metallicity and increasing $\rm[Mg/Fe]$, with $\rm[Mg/Fe]$ values between
$0.1$ dex for the innermost radial bin, and $\sim 0.25$ dex for the outermost
bins ($r\gtrsim 100$”). Hence, Fig. 12 confirms the presence of a (mild)
positive radial gradient of $\rm[Mg/Fe]$ for the bulge of M31. Note that a
lower $\rm[Mg/Fe]$ in the center of M31 has also been reported by Davidge
(1997).

Figure 11: The age indicator $\rm H\beta_{o}$ is plotted against the
metallicity indicator $\rm[MgFe]^{\prime}$, for EMILES base (top) and
$\alpha$–MILES (bottom) SSP model predictions. In both panels, observed line-
strengths for the M31 bulge are shown with dots and error bars (1-sigma
confidence levels), with colours varying from blue (negative values of the
radial distance $R$) through red (positive values of $R$), as shown in the
bottom-panel inset bar. In each panel, the grids show SSP model predictions
with varying age and metallicity, along the horizontal and vertical
directions, respectively, as shown by the labels. In the top panel, dark- and
light-grey grids correspond to Padova and Teramo models. For $\alpha$–MILES
models (bottom), labels are shown only for the scaled–solar (green) models.
Note that the grids correspond to different values of age and metallicity,
depending on the models’ isochrones (Teramo vs. Padova). The age range from
$\sim 6$ to $\sim 14$ Gyr, while metallicity varies from about $-0.25$ dex
($-0.4$) for Teramo (Padova) models, up to $+0.4$ dex. Predictions for the
highest metallicity values are obtained by extrapolation.  Figure 12: The
$\rm\langle Fe\rangle$ spectral index is plotted as a function of $\rm
Mgb5177$. Observed line-strengths for the M31 bulge are shown with coloured
dots, as in Fig. 11. The dark- and light- grey grids show SSP model
predictions of $\alpha$–MILES models, for an age of $7$ and $14$ Gyr,
respectively, corresponding to the innermost and outer radial bins of M31. The
grids show the effect of varying metallicity, from $-0.2$ to $+0.26$ dex, and
$\rm[Mg/Fe]$, from 0 to $0.4$ dex, as labeled. Note that points for
$\rm[Mg/Fe]$$=0.1$, $0.2$, and $0.3$ dex in the grids, are obtained by linear
interpolation of the models.

##  Appendix D The $\rm TiO1$ and $\rm TiO2$ mismatch at $\rm R=0$

In the innermost radial bins, we are able to match all observed indices of
M31, but $\rm TiO1$ and $\rm TiO2$ (see Fig. 3 and Sec. 4.4). The mismatch
stems from the fact that EMILES models predict both $\rm TiO1$ and $\rm TiO2$
to increase with age and IMF slope, while being independent of $\rm[Z/H]$ (for
$\rm[Z/H]$$\gtrsim 0$). On the other hand, the observed TiO’s show a sharp
rise towards the center of the bulge (within $\sim 40$”), where age is
younger, $\rm\Gamma_{b}$ higher, and metallicity is above solar (see Fig. 4
and 5). Therefore, the best-fit values of both TiO’s underestimate the
observed line-strengths at $\rm R\sim 0$. Note that the mismatch is unlikely
due to the effect of non-solar abundance ratios. In fact, the discrepancy
between observed and model indices does not disappear when including
individual abundance ratios in the fitting procedure (methods I and K).
Moreover, TiO indices are expected to increase with $[\rm\alpha/{\rm Fe}]$,
while $\rm[Mg/Fe]$ abundance is found to decrease in the bulge center (see
Fig. 4).

In the center of M31, stellar metallicity is significantly above solar, with
values as high as $\rm[Z/H]$$=+0.4$ dex. For $\rm[Z/H]$ above solar, EMILES
models are actually computed for $\rm[Z/H]$$=+0.15$, $+0.26$, and $+0.4$,
respectively. However, predictions are safe up to $\rm[Z/H]$$=+0.26$, while
they become unsafe at $\rm[Z/H]$$=+0.4$, due to the dearth of high metallicity
stars in the input stellar libraries (see V15). For this reason, in our
reference approach (e.g. LB17 and LB19), we do not use models with
$\rm[Z/H]$$+0.4$, but rather perform a linear extrapolation of the models
above $+0.26$ (see Sec. 3). To investigate if this extrapolation might be
responsible for the TiO mismatch in the center of M31, we created a set of
EMILES experimental models, with the main aim of improving predictions at
$\rm[Z/H]$$=+0.4$ dex. First, we note that, when computing an SSP model, the
iron metallicity ([Fe/H]), total metallicity ($\rm[Z/H]$), and $\rm[Mg/Fe]$ of
stars to be attached to the isochrone, are related by the equation:

$\rm[Z/H]=[Fe/H]+\beta\cdot[Mg/Fe],$ (2)

where the constant $\beta$ is $\sim 0.75$ for $\alpha$-MILES models (see V15).
Therefore, to compute a given SSP, one needs lower [Fe/H] to achieve the
required $\rm[Z/H]$, if $\rm[Mg/Fe]$$>0$, meaning that one can compute better
quality models if both $\rm[Z/H]$ and $\rm[Mg/Fe]$ are above solar.
Unfortunately, most of MILES stars at [Fe/H]$\gtrsim 0$ have $\rm[Mg/Fe]$
around solar, meaning that, in principle, if one wants to create an
$\alpha$-enhanced model at [Z/H]$\gtrsim 0$, one has to rely on theoretical
differential corrections, with all the uncertainties inherited from
theoretical models of stellar atmospheres (see V15). A way out is that of
using the scatter in $\rm[Mg/Fe]$ of stars at given [Fe/H], to create, by
interpolation, empirical stellar spectra with $\rm[Mg/Fe]$$>0$. Therefore, if
we do not “travel too much above” $\rm[Mg/Fe]$$\sim 0$, we can construct
better-quality $\alpha$-enhanced models at supersolar metallicity, entirely
based on empirical stellar spectra. We decided to follow this approach, by
computing a set of EMILES experimental empirical models with
$\rm[Mg/Fe]$$=+0.1$ 171717These $\alpha$-enhanced models are computed assuming
scaled-solar isochrones, as $\alpha$-enhanced Teramo isochrones are available
only for $\rm[Mg/Fe]$$=0.4$. However, as discussed in V15, the effect of
$\alpha$-enhanced isochrones is secondary compared to that of $\rm[Mg/Fe]$ for
stellar spectra..

Figure 13: SSP model predictions for $\rm TiO1$ (panel a), $\rm TiO2$ (panel
b), $\rm H\beta_{o}$ (panel c), and $\rm[MgFe]^{\prime}$ (panel d), as a
function of total metallicity $\rm[Z/H]$. Dashed and solid lines refer to
models with an age of 7 and 14 Gyr, typical for the innermost and outermost
bins of M31, respectively. Thin and thick lines correspond to a Kroupa-like
($\rm\Gamma_{b}$$=1.3$), and a bottom-heavy ($\rm\Gamma_{b}$$=3$) IMF,
respectively, while black, green, and red curves plot predictions for base,
scaled-solar, and $\alpha$-enhanced EMILES testing models, with
$\rm[Mg/Fe]$$=0.1$, specifically computed for this project (see the text). The
safe range of EMILES is up to $\rm[Z/H]$${}_{\rm safe}=+0.26$ (see V15), with
the $\rm[Mg/Fe]$$=0.1$ testing models having better quality above this
metallicity threshold. Note the trend of $\rm TiO1$ and $\rm TiO2$ to be
constant up to $\rm[Z/H]$safe, and then increase significantly with $\rm[Z/H]$
for $\rm[Z/H]$$>$$\rm[Z/H]$safe.

In Fig. 13 we plot SSP model predictions for both $\rm TiO1$ and $\rm TiO2$,
as well as for the age-sensitive index, $\rm H\beta_{o}$, and for the total
metallicity indicator $\rm[MgFe]^{\prime}$, as a function of $\rm[Z/H]$. We
consider only models computed with Teramo isochrones. Black and green curves
in the Figure show EMILES base and scaled-solar models, respectively, while
red curves correspond to the experimental empirical models with
$\rm[Mg/Fe]$$=+0.1$. Dashed and solid lines in the Figure refer to models with
age of 7 and 14 Gyr, respectively, typical for the innermost and outermost
bins of M31. Thin and thick lines are for a Kroupa-like
($\rm\Gamma_{b}$$=1.3$), and a bottom-heavy ($\rm\Gamma_{b}$$=3$) IMF,
respectively, Fig. 13 shows that both TiO indices do not depend significantly
on metallicity from $\rm[Z/H]\\!\sim\\!0$ to $\rm[Z/H]\\!\sim\\!+0.26$.
Therefore, when extrapolating the models above $\rm[Z/H]$$=+0.26$, as in our
reference approach, one predicts both TiO’s to be constant all the way up to
$\rm[Z/H]$$=+0.4$ dex, and above. On the contrary, base and scaled-solar
models suggest an increase of 0.01–0.015 mag (depending on age and IMF slope)
for $\rm TiO2$, and about 0.005–0.01 mag for $\rm TiO1$, from
$\rm[Z/H]$$=+0.26$ to $\rm[Z/H]$$=+0.4$. Although models at $\rm[Z/H]$$=+0.4$
are unsafe, this trend suggests that our linear extrapolation might be
underestimating the TiO line-strengths at $\rm[Z/H]$$=+0.4$, explaining the
mismatch in the center of M31. The experimental empirical models with
$\rm[Mg/Fe]$$=+0.1$ seem to confirm that this is actually the case, predicting
an even stronger increase of both TiO’s with $\rm[Z/H]$ (see the red curves in
the Figure for $\rm[Z/H]$$>0.26$), with respect to scaled-solar and base
models. However, at the same time, the models with $\rm[Mg/Fe]$$=+0.1$ predict
an increase (decrease) of $\rm H\beta_{o}$ ($\rm[MgFe]^{\prime}$) with
$\rm[Z/H]$ (for $\rm[Z/H]$$>0.26$), that contrasts with the expectation of
$\rm H\beta_{o}$ ($\rm[MgFe]^{\prime}$) being insensitive to metallicity
($\rm[Mg/Fe]$), hence casting some doubts about the reliability of predictions
in Fig. 13. Therefore, while the present analysis seems to support the
conclusion that the mismatch of TiO observed and model line-strengths in the
center of M31 is due to an increase of both $\rm TiO1$ and $\rm TiO2$ with
metallicity, in the highly supersolar regime, current stellar population
models do not allow us to draw firm conclusions.

In order to account for a possible increase of $\rm TiO1$ and $\rm TiO2$ at
$\rm[Z/H]$$>0.26$, we modified our best-fitting method I (see Sec. 4.3). For
$\rm[Z/H]$$>0.26$, we assume a linear parametrization of TiO1 and TiO2 as a
function of $\rm[Z/H]$, with slopes of
$\rm\Delta_{Z}(TiO1)=\delta(TiO1)/\delta([Z/H])$ and
$\rm\Delta_{Z}(TiO2)=\delta(TiO2)/\delta([Z/H])$, respectively. The values of
$\rm\Delta_{Z}(TiO1)$ and $\rm\Delta_{Z}(TiO2)$ are assumed to be independent
of age and IMF slope, and are treated as two extra fitting parameters. From
the data of M31 at $\rm R=0$, we obtain $\rm\Delta_{Z}(TiO1)=0.07\pm 0.014$
mag/dex and $\rm\Delta_{Z}(TiO2)=0.23\pm 0.014$ mag/dex. Note that these
values are remarkably consistent with those one can estimate from the
experimental models in Fig. 13, i.e. $\rm\Delta_{Z}(TiO1)\sim 0.1$ mag/dex and
$\rm\Delta_{Z}(TiO2)\sim 0.25$ mag/dex, respectively. Also, the best-fitting
line-strengths for the modified method I, plotted as green asterisks in Figs.
2 and 3, match all indices of M31 at $\rm R=0$, including the two TiO’s
181818As a further test, we also included the NIR TiO index $\rm TiO0.89$ in
the fitting (see LB17), finding no significant differences in the estimated
values of $\rm\Delta_{Z}(TiO1)$ and $\rm\Delta_{Z}(TiO2)$, as well as in the
best-fit value of $\rm\Gamma_{b}$. . More important for the purpose of the
present paper, when accounting for $\rm\Delta_{Z}(TiO1)$ and
$\rm\Delta_{Z}(TiO2)$, the best-fitting IMF slope, $\rm\Gamma_{b}$, changes
only slightly with respect to method I, decreasing from $\sim 2.4$ to $\sim
2.3$. We conclude that, while better models are certainly needed in the very-
high metallicity regime, our results are robust against the model
uncertainties at high $\rm[Z/H]$.

##  Appendix E Kinematics

To extract the kinematics of the M31 bulge, we have performed an adaptive
binning of the OSIRIS spectra as a function of radial distance, assuming a
minimum $S/N$ ratio of 25 Å-1 (see Sec. 2.3 for details). We verified that
adopting a minimum $\rm S/N$ of 15, rather than 25, does not change
significantly our results. We measured rotation velocity $\rm V_{rot}$,
velocity dispersion, $\rm sigma$, and the higher moments of the line-of-sight
(LOS) velocity distribution, $\rm H_{3}$ and $\rm H_{4}$, by running the
software pPXF (Cappellari & Emsellem, 2004; Cappellari, 2017). We ran pPXF
within four different spectral ranges, i.e. the “blue” ($\lambda\lambda
3900$-$4900$ Å), “Mgb” ($\lambda\lambda 4750$-$5750$ Å), “NaD”
($\lambda\lambda 5700$-$6200$ Å), and “CaT” ($\lambda\lambda 7670$-$8850$ Å)
spectral regions, adopting two sets of templates, i.e. EMILES SSP models with
varying age and metallicity (see Sec. 3), as well as stellar spectra from the
MILES (Sánchez-Blázquez et al., 2006) (for blue, Mgb, and NaD ranges) and CaT
(Cenarro et al., 2001a) (for the CaT range) libraries. For the Mgb range, we
also ran pPXF with $\alpha$–MILES models (including both scaled-solar and
$\alpha$–enhanced models). Fig. 14 compares the kinematics radial profiles
obtained for different methods. Indeed, we find an excellent agreement among
different spectral ranges, and different templates, but for a small offset in
$\sigma$ (amounting to about $-10$ km s-1) for the CaT spectral range (see
pink and light-green curves in panel b of the Figure). Since we obtained
consistent results among methods, we median-combined different results into
final kinematics profiles, shown as black curves in Fig. 14.

Figure 14: The kinematics profiles for the bulge of M31, i.e. rotation
velocity $\rm V_{rot}$(panel a), velocity dispersion $\rm sigma$ (panel b), as
well as $\rm H_{3}$ and $\rm H_{4}$ (panels c and d), are plotted as a
function of galactocentric distance, $\rm R$, in units of arcsec. Lines with
different colours correspond to results obtained by running pPXF on different
spectral ranges, i.e. the “blue”, “Mgb”, “NaD”, and “CaT” regions (see the
text), and with different sets of templates, i.e. EMILES SSP models, MILES and
CaT stars, and $\alpha$–MILES SSP models, as labeled in the lower–left of
panel a. Note that the labels report the fitted spectral range, followed by
the adopted set of templates. Note the small offset of $\sigma$ (about $-10$
km s-1), for results corresponding to the CaT spectral region (see text).
Results from different methods are median-combined into final kinematics
profiles, shown as black curves.

Fig. 15 shows that our kinematics for the M31 bulge (black curves) is in
excellent agreement with that obtained by S10. Rotation is detected at all
radii probed in the present work, and correlates with $\rm H_{3}$. We confirm
that within a region of $\sim 100$”, the galaxy has a velocity dispersion of
$\sim 160$ km s-1, dropping down (by $\sim 10$ km s-1) towards the center (see
Opitsch et al. 2018 and references therein). On the other hand, in the very
central bins, the velocity dispersion increases steeply – a signature of the
central supermassive black hole. Within $\sim 100$”, the $\rm H_{4}$ is
slightly positive, with a typical value of $\sim 0.02$, again in very good
agreement with S10. Note that in our data, we also see a slight trend of $\rm
H_{4}$ decreasing outwards, which is not clearly seen in S10 data.

Figure 15:  Comparison of our kinematics for the bulge of M31 (black circles
and error bars, as in Fig. 14), with that obtained by Saglia et al. 2010,
shown as red stars. Panels a, b, c, and d correspond to the same kinematics
profiles as in Fig. 14. Note the excellent agreement between the two sets of
measurements.

Since the observed line-strengths for CaT are below the predictions of EMILES
models (see, e.g., the red curve in panel f of Fig. 2), the small offset in
$\rm sigma$ for the CaT spectral range (see above) might be due to template
mismatch. In order to further explore this point, we re-ran pPXF in the CaT
region with EMILES models, applying a response for an abundance pattern of
$\rm[Ca/Fe]$$=-0.2$ (as estimated from CvD12a models). In this case, results
for the CaT are in excellent agreement with those obtained for the other
spectral ranges, as shown by the thick magenta line in Fig. 16. We conclude
that template mismatch is the most likely explanation for the (small) offset
of $\rm sigma$ in the CaT region.

Figure 16: The same velocity dispersion profiles as in panel c of Fig. 14 (see
labels in the upper–left), are compared to the $\rm sigma$ profile obtained
with EMILES models, modified for an abundance pattern of $\rm[Ca/Fe]$=$-0.2$
dex, plotted as a thick pink curve (see the text). Note that for models with
$\rm[Ca/Fe]$$=-0.2$ dex, we match the kinematics obtained in the bluer
spectral ranges, removing the $\sigma$ offset seen in panel b of Fig. 14.