The Chemodynamics of the Stellar Populations in M31 from APOGEE Integrated-light Spectroscopy

Benjamin J. Gibson; Gail Zasowski; Anil Seth; Aishwarya Ashok; Kameron Goold; Tobin Wainer; Sten Hasselquist; Jon Holtzman; Julie Imig; Dmitry Bizyaev; Steven R. Majewski

doi:10.3847/1538-4357/acd9a9

1. Introduction

Roughly half of the stellar mass in the universe is contained in galaxies whose mass is within a factor of 2 of the Milky Way's (e.g., Baldry et al. 2012; Driver et al. 2022). Therefore, by studying the chemical abundances, dynamics, and their joint distributions (chemodynamics) of stars in these galaxies, we can further understand the physical processes that govern the behavior of baryonic matter. We have a detailed understanding of the chemical makeup and kinematics of the stars in the Milky Way (MW), as our perspective from within the galaxy allows us to observe individual stars both photometrically and spectroscopically throughout it. Surveys such as Gaia (Gaia Collaboration et al. 2016), APOGEE (Majewski et al. 2017), LAMOST (Cui et al. 2012), GALAH (De Silva et al. 2015), and several others, have observed billions of stars throughout the MW and provided us with detailed chemodynamical information for each of them. When we break these stars down into groups with similar chemodynamics, we can identify subpopulations that likely were born in different locations at different times (e.g., Helmi et al. 1999, 2018; Belokurov et al. 2006; Koppelman et al. 2019; Myeong et al. 2019). This practice has made the MW the primary basis for understanding how MW-like spiral galaxies form and evolve.

Unfortunately, this perspective from inside the MW has some drawbacks. For example, stars on the farside of the disk and bulge can often be obscured by gas and dust or other foreground stars. Being in the plane of the disk complicates the classification of spiral arms and the bar, leading to different interpretations of galactic structures (e.g., Lee et al. 2015; Ness & Lang 2016). It is also difficult to determine global properties of the MW such as total luminosity, color, mass-to-light ratio, bulge-to-disk mass ratio, etc. as we cannot see all of the MW in the same way we can see external galaxies (Bland-Hawthorn & Gerhard 2016; Fielder et al. 2021).

Luckily, we can work around some of these issues by studying our neighbor, M31, which is believed to be relatively similar to the MW. Despite its distance of ∼785 kpc (e.g., McConnachie et al. 2005) and inclination of 77°, we are able to observe all of M31's disk and bulge. Like the MW, M31 is a large disk galaxy, though it features rings rather than spiral arms. The stellar mass of M31's disk is 5.9 × 10¹⁰ M_⊙ (Yin et al. 2009), about 1.5× that of the MW's at 4 × 10¹⁰ M_⊙ (Sick et al. 2015). M31's disk scale length (5.76 kpc, Dorman et al. 2013) is about 2.2× larger than the MW's (2.6 ± 0.5 kpc, Bland-Hawthorn & Gerhard 2016). These two galaxies are similar in metallicity (e.g., Gregersen et al. 2015) and environment, being in the same galaxy group. M31 has a bar, but it is much weaker and less distinct than the MW's (Athanassoula & Beaton 2006), and M31's bulge is composite: one-third of its stellar mass is in a classical bulge and the other two-thirds is in a boxy/peanut-shaped (b/p) bulge (Blañ et al. 2017). The MW bulge is mostly b/p (Wegg & Gerhard 2013). The bulge-to-disk mass ratio for M31 is 0.43 (Courteau et al. 2011), compared to 0.3 for the MW (Bland-Hawthorn & Gerhard 2016).

The stellar populations of M31's bulge and disk have been the focus of numerous studies over the last 80 yr. Baade (1944) barely resolved stars at the tip of the RGB in M31's bulge and concluded they were more similar to stars seen in the solar neighborhood of the MW, rather than the types of stars seen in globular clusters. In the following decades, several resolved star studies concluded M31's bulge was made of old, metal-rich, Population I stars (e.g., Baum & Schwarzschild 1955; Spinrad & Taylor 1971; Wu et al. 1980; Bohlin et al. 1985; Walterbos & Kennicutt 1987). Bica et al. (1990) measured the metallicity of the bulge and found it to be super-solar ([Z/Z_⊙] ∼ 0.3). The first element abundance measured in the bulge was [Mg/Fe] at 0.3 by Sil'chenko et al. (1998).

Recently, resolved stellar photometry has become available across a vast extent of M31 (e.g., McConnachie et al. 2009; Dalcanton et al. 2012). Olsen et al. (2006) studied resolved near-infrared (NIR) photometry in the bulge and disk, finding that the bulge is old (≥6 Gyr) and has solar metallicity. They also note that the stellar ages are uniform throughout the bulge and inner disk. The Panchromatic Hubble Andromeda Treasury (PHAT; Dalcanton et al. 2012) measured UV, optical, and IR photometry of 117 million individual stars in the bulge and northern half of M31. The star formation history (SFH) of M31 recovered from this survey showed a majority of the stars in the disk are older than 8 Gyr, with a smaller component that has formed in the last 2–4 Gyr (Williams et al. 2017). The 5 kpc inner ring was also found to have ongoing star formation in the last 200 Myr (Lewis et al. 2015). Spectroscopic kinematic measurements of 5800 PHAT stars show a very clear age-velocity dispersion correlation in the disk, where older stars have a higher dispersion (Dorman et al. 2015). A negative metallicity gradient (−0.2 dex kpc⁻¹) was found from 4–20 kpc in the disk, as well as a slight metallicity enhancement along the bar from 3–6 kpc (Gregersen et al. 2015). Analysis of PHAT photometry combined with resolved stellar spectra from the SPLASH survey (Guhathakurta et al. 2005; Gilbert et al. 2006) indicates the disk has a dominant metal-rich population at a median of [Fe/H] = −0.19 and that the disk has been significantly dynamically heated (Escala et al. 2023).

In recent years, several studies have utilized simple stellar population models to analyze integrated-light spectroscopy of M31's bulge and disk. Saglia et al. (2010) examined long-slit spectra of the bulge region and found it to be uniformly old (12 Gyr), of solar metallicity, and to have an alpha abundance of [α/Fe] ∼ 0.2. The very center is younger and slightly more metal-rich. A study of LAMOST (Cui et al. 2012) spectra in the inner 7 kpc agrees that the bulge is uniformly old, and the disk has solar metallicity and a negligible gradient (Zou et al. 2011).

To date, the most detailed integrated-light studies of M31 come from Opitsch et al. (2018, hereafter referred to as O18), who examined optical integral-field unit (IFU) data of the inner bulge and parts of the disk. They constructed detailed kinematic maps of this region and found strong evidence for a bar. Saglia et al. (2018, hereafter referred to as S18) examined the same data to construct maps of the chemistry and age of the stellar populations. They found the bulge to be old (∼13 Gyr), metal-rich (roughly solar abundance) along the bar, and α enhanced (0.2–0.3 dex). Gajda et al. (2021) developed a made-to-measure model of M31 that matched the same data and found that the part of the bulge on either side of the bar was less metal-rich than the bar and disk. They called this pattern a metallicity desert.

In this paper, we analyze high-resolution, NIR, integrated-light spectra from APOGEE (Majewski et al. 2017) of M31 out to a radius of ∼7 kpc⁸ (see Figure 1). Our data cover the entire bar of M31 as well as the inner disk and inner ring at 5 kpc. This data set is unique compared to previous integrated-light spectroscopic studies of M31. No other data set has complete coverage of the inner regions, as most span a few fields in the central regions, or have continuous coverage that does not extend quite so far out. These data are also of higher spectral resolution than previous studies by nearly an order of magnitude, which allows us to measure the average kinematics of the stellar populations with unprecedented accuracy.

APOGEE's NIR coverage contains spectral features from several α elements (Smith et al. 2021), which can be used to trace the SFH and chemical evolution of stellar populations. The NIR is significantly less affected by dust attenuation than the optical and is sensitive to RGB and AGB stars. This allows us to clearly observe and more completely sample the stellar population. Moreover, APOGEE's measurements of ∼650,000 stars in the MW have been leveraged to create a library of single stellar population (SSP) model spectra (Ashok et al. 2021) with the same spectral resolution and wavelength range. This library is unique to our study, affording us an advantage over other studies, as they may have to navigate discrepancies between their data and models, such as differences in spectral resolution, that we do not. Therefore, we can analyze our M31 data with our unique models to enable direct, consistent comparison between the stellar populations of the MW and M31.

In this paper, we present maps of the mean kinematics, abundances, and ages of the stellar population of M31. These maps cover, for the first time, the entirety of the bulge, bar, inner disk, and 5 kpc ring. Section 2 describes the APOGEE observations and data processing, Section 3 details our SSP models, and Section 4 describes the data binning and fitting analysis. The results are presented and discussed in Section 4, and the summary and conclusion are in Section 5.

2. Observations and Data Processing

Data for this project was taken as part of the Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2017; S. M. Majewski 2023, in preparation). APOGEE, part of the Sloan Digital Sky Survey (SDSS)-III (Eisenstein et al. 2011) and SDSS-IV (Blanton et al. 2017), is a NIR (1.5–1.7 μm), high-resolution (R ∼ 22,500), multi-object spectroscopic survey. APOGEE observed ∼650,000 stars in the MW, as well as individual stars in nearby dwarf galaxies and integrated-light spectra from extragalactic GCs and nearby galaxies (Zasowski et al. 2013, 2017; Beaton et al. 2021; Santana et al. 2021). Data analyzed in this paper were taken using the APOGEE spectrograph (Wilson et al. 2019) on the 2.5 m SDSS telescope at Apache Point Observatory (Gunn et al. 2006). APOGEE also includes data taken from a second APOGEE spectrograph on the du Pont Telescope at Las Campanas Observatory (Bowen & Vaughan 1973).

The standard data processing pipeline for APOGEE is outlined in Nidever et al. (2015), and the latest data release (DR17) is described in Abdurro'uf et al. (2022). The APOGEE Stellar Parameter and Chemical Abundances Pipeline (ASPCAP; García Pérez et al. 2016; Nidever 2021; J. H. Holtzmann et al. 2023, in preparation) determines stellar parameters and ∼20 elemental abundances for stars in the APOGEE data set. Note that ASPCAP was used to develop the models described in Section 3, not for abundance determination of our M31 data.

Throughout this paper, we use fiber position to indicate a location in M31 that has been observed. Each fiber position is observed multiple times, and we will refer to each individual observation as a visit (see Zasowski et al. 2013 Section 2.1 for more details).

2.1. Observations

The M31 observations were performed as part of an ancillary science program of APOGEE (Zasowski et al. 2017). Five sets of targets were designed to observe a total of 1105 locations. These designs were drilled into plug plates for fiber optic cables that gather the light from an object and transmit it to the APOGEE spectrograph. One plate with an M31 design can be seen in Figure 12 in Appendix B. The fiber placement was designed to be as evenly spaced as possible, but was subject to several constraints. APOGEE fibers are 2'' in diameter, and the collision distance between fibers was 72''. The design of the first plate included the very center of M31 and allocated fibers outward from there. The other plates were designed to start 10'' from the center and their fiber positions were no closer than 10'' from the fiber positions in any other plate. Fiber positions avoided bright sources identified with 2MASS (Skrutskie et al. 2006) imagery, as these were assumed to belong to the MW. Sky and telluric calibration fiber positions were determined by standard APOGEE routines (see Sections 5.1 and 5.2 of Zasowski et al. 2013 for more details). They are located outside the disk of M31 and are roughly uniformly dispersed across a circle ∼1.75° in radius from the center of M31. All of the calibration was done using standard APOGEE procedures (see Sections 6.2 and 6.3 of Nidever et al. 2015).

Figure 1 shows the distribution of fiber positions on an optical image of M31. Positions filled with gray are not included in the final analysis (see Section 2.3). For the remainder of this paper, the term bulge will refer to the region inside the black and cyan ellipse in Figure 1. This was chosen to match the region where the bulge makes up >50% of the total luminosity of the galaxy according to Courteau et al. (2011; see Section 5.3.1). With this definition, the bulge dominates the light to ∼1.80 kpc along the disk major axis and ∼5.56 kpc along the disk minor axis, due to the inclination of M31's disk. We note that radii measured in kiloparsec are deprojected onto the disk of M31.

These observations extend to a radius of ∼7 kpc along the bar major axis. Each fiber position was observed 10 times for roughly 66 minutes each time. There are 275 fiber positions in the bulge region; the other 830 are in the inner-disk region, roughly from 3–7 kpc from M31's center. Our fibers cover roughly 150 kpc² of M31's bulge and disk, for an average density of ∼7.4 fibers per square kiloparsec.

2.2. Data Processing

The APOGEE pipeline is optimized for data of individual stars, not integrated light of external galaxies, which often have a low signal-to-noise ratio (S/N) and high velocity dispersion. To account for this, we modified the final step of the standard APOGEE pipeline: the combination of several visits of each fiber position into a single spectrum. The entire process is described in Nidever et al. (2015). We refined this final step as described below to adequately identify and alleviate errors (from skylines, detector persistence, etc.) while preserving as much of our data as possible.

The data gathered for an object observed by APOGEE is condensed into an apStar file (Nidever et al. 2015).⁹ These files contain target flux spectra, sky flux spectra, target flux uncertainty, and pixel bitmasks for each visit to a fiber position, as well as each of these for a single combined spectrum.¹⁰ The apStar files will only include a visit if it is used to create the final combined spectrum, excluding data with very low S/N and poor radial velocity (RV) determination.¹¹ These choices were optimized for stellar spectra, and exclude some of our low S/N and data. For this reason, the standard data release products contain spectra for only 877 fiber positions.¹² We worked with the APOGEE data team to create custom apStar files that removed the S/N cutoff and the RV and barycentric corrections. We received 1082 custom apStar files that contain seven to 10 visits each. These custom apStar files, as well as the final data set outlined in Section 2.3, are available upon request and will be made publicly available as part of a future publication.

2.2.1. Masking Individual Visit Spectra

The APOGEE spectrograph is made of three CCD chips covering different wavelength ranges. There is no data gathered for the wavelengths between the chips. To avoid these chip gaps, we split each spectrum into three separate arrays, one for each chip. Each chip is then masked, combined, and stored individually.

To combine the visit spectra into one final spectrum, we need to identify pixels with unreliable fluxes and significant sky emission. To do this, we first take a 500 pixel running median of the object flux spectrum for each visit to use as a continuum (this same strategy is used by the APOGEE data reduction pipeline). We then calculate ΔFlux = ∣flux − continuum∣, as well as S/O = $\tfrac{{\rm{sky}}\,{\rm{flux}}}{{\rm{object}}\,{\rm{flux}}}$ for each pixel. ΔFlux was chosen as a means to simply remove the pixels that are farthest away from the continuum, which are generally, but not always, due to sky emission. S/O was chosen as a means to also identify pixels with high sky flux. We then calculate the 85th percentile for ΔFlux and S/O and mask pixels that fall above this line for both statistics. Pixels flagged by one statistic are often, but not always, flagged by the other. A percentile cutoff was chosen as it will remove the most affected ∼15% of pixels from a spectrum regardless of the visit's S/N.

We also mask pixels in each visit spectrum according to the visit-level PIXMASK flags.¹³ Specifically, we mask out pixels with bits 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, and 14 set. These bits identify pixels affected by, among other things, medium to high levels of persistence, bad calibration frames, cosmic rays, or saturation.

Upon visual inspection, a number of our spectra contained what appeared to be a prominent emission feature centered near 15930 Å. This feature was identified as an artifact resulting from decreased sensitivity in a small region of the detector; the dome flat shows almost no flux between 15890 and 15960 Å for fibers with IDs between 20 and 30. Since dome flats are removed by dividing, a low amount of flux from the flat will cause a spike in the object flux in the affected region. We mask out this region in affected fibers.

In summary, our final visit-level pixel mask flags pixels affected by severe sky emission, several issues identified by the APOGEE pixel masks, and one region of low detector sensitivity in the green chip. This amounts to masking ∼25% of our available pixels. We note that these issues are specific to our data and not indicative of larger issues with the APOGEE data processing pipeline.

2.2.2. Combining Masked Visit Spectra

The products from the previous section are masked visit spectra of each fiber position. For a given fiber position, each visit's chip spectra are pseudo-continuum normalized using a new 500 pixel running median that ignores the pixels masked in Section 2.2.1. The visits are then combined using an inverse variance weighted average:

$\begin{eqnarray}&&F=\displaystyle \frac{{\rm{\Sigma }}\tfrac{{f}_{i}}{{\sigma }_{i}^{2}}}{{\rm{\Sigma }}\tfrac{1}{{\sigma }_{i}^{2}}};\sigma =\displaystyle \frac{1}{{\rm{\Sigma }}\tfrac{1}{{\sigma }_{i}^{2}}},\end{eqnarray} \tag{ 1 }$

where F and σ are the combined pixel flux and uncertainty, respectively, and f_i and σ_i are the visit pixel continuum-normalized flux and uncertainty, respectively.

If a pixel is masked in at least half of the available visits, then it is masked out in the final combined spectrum. These pixels are assigned a flux value equivalent to the average continuum flux at that pixel and are given a very large uncertainty, we used 10¹⁵ in flux units. This combining process is the same process used in the apStar files¹⁴ (Nidever et al. 2015), with the primary distinctions being that we remove masked pixels from the combined spectrum and leave our final spectra continuum normalized.

2.3. Determining the Final Data Set

For our final data set that we consider in the rest of the paper, we removed fiber positions from our sample for which we were unable to determine reliable kinematics from the final combined spectrum. Specifically, we removed fiber positions that had inconsistent velocities from chip to chip or that had unrealistically high velocity dispersions. The majority of these fiber positions are in the outskirts, and have low S/N spectra that complicate the measurement of their kinematics. We performed the preliminary full-spectrum fitting routine with ppxf (Cappellari 2017) defined in Section 4.2, which matches our data with the model spectral templates described in Section 3; the kinematics are allowed to vary freely from chip to chip in this procedure. From these fits we determined the three templates with the highest weight used in the ppxf fit (weights normalized across each chip). Next, each chip was fit to each of the three templates individually and we determined which combination of templates gave the smallest variation in velocities from chip to chip. We identified 95 fiber positions with RV variations >50 km s⁻¹, then visually inspected each of these spectra and kept seven fiber positions that had otherwise consistent kinematics and RV variations only slightly above 50 km s⁻¹. We also identified 22 fiber positions from the 95 that appeared to have prominent spectral features but nevertheless gave inconsistent kinematics. These 22 spectra were fit again using the Markov Chain Monte Carlo (MCMC) routine described in Section 4.2 and we kept nine whose kinematics converged well.

In the end, 32 fiber positions had no visits with S/N above 1, so these were removed from the sample. We cut 79 fiber positions due to high RV variations between chips. Three other fiber positions had dispersions in one or more chips >220 km s⁻¹, and we do not expect values this high in the disk of M31, so they were cut. Another three fiber positions were cut that had χ² values, which are used to determine the goodness-of-fit, above 5 in the blue and red chips, or above 8 in the green chip (green chip fits are, on average, worse). Lastly, we cut two more fiber positions that were unintentionally centered on GCs. We identified these fiber positions by their low dispersion (<60 km s⁻¹) and verified the GC using imaging from the PHAT survey (Dalcanton et al. 2012). Our final sample consists of 963 fiber positions, or 89% of our initial 1082 fiber position sample (see Figure 1).

2.4. Velocities from Opitsch et al. (2018)

Preliminary testing showed that we could only reliably recover the velocity curve seen in M31 from individual fiber positions with an empirical signal-to-noise ratio (eS/N, see Section 4.2) of ≳20. We also found that an eS/N of ∼60 was needed to determine abundances. Therefore, we bin fiber positions together in the outer disk to increase their eS/N, which requires reliable velocity measurements for each fiber position. We adopt mean stellar radial velocity (RV) data from O18 for shifting and combining binned data (see Section 4.1).

O18 studied the stellar kinematics of M31 using the VIRUS-W optical IFU spectrograph mounted on the 2.7 m telescope at the McDonald Observatory (R ∼ 9000, Fabricius et al. 2012). The O18 data cover M31's bulge as well as azimuthal strips, or spokes, of the disk out to 5.3 kpc along the disk major axis. We take O18's measured line-of-sight velocities and LOESS smooth them using the LOESS Package (Cappellari et al. 2013), which implements the locally weighted regression algorithm of Cleveland (1979) to adaptively smooth over observational errors and unveil the data's underlying structure. We smoothed the data with a degree of 1, fraction of 0.5 (both are parameters in the LOESS package code that help control the smoothness), and included the error bars. This algorithm can extrapolate onto locations that lie beyond the spatial extent of the data; some of our fiber positions lie outside of the bulge region measured by O18 (see Figure 13). The smoothed map was evaluated at our fiber positions to get an RV measurement for each; we call this set of LOESS-smoothed velocity measurements the O18 velocities. These were compared to velocities we measured for spectra we Voronoi binned (in the same manner as Section 4.1) to a target eS/N of 20. Subtracting the former from the latter we found a bias of −0.65 km s⁻¹ and a scatter of ∼19 km s⁻¹. The velocities published in O18 have errors of 3–4 km s⁻¹, which are comparable to our single-fiber position errors. A plot comparing our velocities to those of O18 can be found in Appendix C.

3. SSP Models

3.1. A-LIST

The APOGEE Library of Infrared SSP Templates (A-LIST; Ashok et al. 2021) is a library of SSP template spectra generated from APOGEE stellar spectra. A-LIST spans a wide range of population ages, metallicities, and α abundances. This library of spectral templates was designed specifically to analyze IL spectra from APOGEE.

We only utilize templates with luminosity fraction >0.32, which is a measure of the completeness of the APOGEE sample for a given age, metallicity, and α abundance, as recommended by Section 4.1.1 in Ashok et al. (2021). We also restrict our template sample to have [M/H] ≥ −1.0 and Age ≥ 7 Gyr. The method we use to interpolate between A-LIST templates (described next) can be simplified by restricting our parameter space, since there is an age–metallicity degeneracy in the lowest metallicities. Moreover, we do not expect to find the average population to be more metal-poor or younger in the bulge and inner disk of M31 (e.g., Olsen et al. 2006; Saglia et al. 2018).

3.2. A-LIST Interpolation

We interpolate between the discretely gridded A-LIST templates to use them with the MCMC fitting method described in Section 4.2. We adopt The Cannon (Ness et al. 2015; Casey et al. 2016) for this interpolation. The Cannon is a data-driven method to transfer stellar parameters and abundances (labels) from a training set of spectra to a broader data set. The labels for the training set are determined through independent methods. The Cannon models the flux in each pixel as a linear combination of the labels. In our case, we train the model on our sample of A-LIST spectra using [M/H], [α/M], and log₁₀(age) as labels. This allows us to generate model spectra continuously across our parameter space.

This interpolation routine has the added benefit of fixing a handful of cases where α-element lines do not deepen linearly with increasing α abundance due to the luminosity fraction and ΔT_eff limits of the library. Additionally, The Cannon training procedure allows the user to add weights to pixels. The Cannon procedure ensures that upweighted pixels will be recreated accurately relative to downweighted ones. We upweighted pixels that show the strongest variation of flux with α abundance, $\tfrac{{dF}}{d\alpha },$ ¹⁵ holding age and metallicity constant. In particular, we weight a fraction of the pixels by

$\begin{eqnarray}&&W=-c\displaystyle \frac{{dF}}{d\alpha },\end{eqnarray} \tag{ 2 }$

where W is the weight and c is 10. We apply this to the 10% of pixels with the steepest slopes. The other 90% of pixels are given a weight of 1.

We tested many models, varying the multiplier and percentile, as well as the polynomial order and regularization of The Cannon. The regularization parameter controls the amount to which the flux determined for a given pixel is influenced by the flux of surrounding pixels. Following Ashok et al. (2021), we analyzed APOGEE IL spectra of M31 GCs with [M/H] ≥ −1.0. We found the model that best matched the original A-LIST spectra and recreated the results from Ashok et al. (2021) was an unregularized third-order polynomial with a multiplier of 10 on the slopes of the top 10% steepest pixels.

As can be seen in the left panel of Figure 2, the original A-LIST templates are well reproduced by our interpolation routine. From the right panel, one can see that we are able to reasonably model the flux for abundances between the discrete A-LIST grid points. Figure 3 shows the distribution of the flux in each pixel of all the interpolated models versus the A-LIST models. We show that the interpolated model is a good recreation of A-LIST, with a median absolute deviation of 0.0013 normalized flux units between the two.

**Figure 3.** Comparison of pixel fluxes in the interpolation and the original A-LIST templates. The color bar indicates the number of pixel fluxes at a given point. The small number of outliers to the upper left of (1,1) are due to the low metallicity, low α abundance templates in our sample. There are not many stars with this chemistry in the MW, and as such these templates are more sensitive to fluctuations from star to star. Therefore, these templates are not necessarily as representative of their populations as templates with higher metallicities or α abundances, and they are affected the most by our interpolation routine.
Download figure:
Standard image High-resolution image

4. Binning and Fitting

4.1. Binning

To bin adjacent fiber positions together, we first shift the spectra to a common wavelength array determined by the average of the O18 velocities of the fiber positions in each bin. We then stack the spectra together using the routine outlined in Section 2.2.2, and calculate the eS/N with the first two steps of our general fitting procedure (Section 4.2).

We determine bin members by Voronoi binning in radius and position angle using the code VorBin (Cappellari & Copin 2003). Voronoi binning allows us to increase our eS/N to a target of 70 while maintaining as much spatial resolution as possible. An added desire for our bins is a small spread in the intrinsic velocity dispersions of the spectra in a bin. Velocity dispersion trends closely with radius, so we multiply the radius by 6 so VorBin thinks the fiber positions are further away radially than they actually are. This produces bins as radially narrow as possible. Lastly, we set the zero position angle to the southern disk semiminor axis so that we can have bins that cross over the major axis.

We made some manual edits to the final binning results to even out the eS/Ns, and in the end we have 164 bins with an eS/N ≳ 60. The final bins are shown as the colored polygons in Figure 5; the fiber positions included in each bin are denoted by the black and white dots inside each polygon. In Section 5.5, we discuss the robustness of this method to the choices involved.

4.2. General Fitting Procedure

We utilize full-spectrum fitting via Penalized Pixel-Fitting (ppxf; Cappellari 2017) to determine the kinematics, abundances, and ages of the stellar populations in M31. ppxf is a reduced χ² minimization software that shifts and broadens template spectra (in our case, the interpolated A-LIST model) to fit data. We assume that the stellar population making up our IL spectrum has the line-of-sight velocity and velocity dispersion found by the ppxf fitting procedure. The light-weighted mean age, metallicity, and α abundance of the population at that fiber position is the age, metallicity, and α abundance of the best-fitting interpolated SSP template. That template is determined by an MCMC routine that searches the entire parameter space for the set with the maximum likelihood, L = −χ²/2, using the emcee software package (Foreman-Mackey et al. 2013).

The fitting procedure is as follows:

1.
Mask the following pixels in each spectrum (in addition to Section 2.2.1):
- (a)
  The 250 pixels closest to the red end of each chip, and the 100 pixels closest to the blue end of each chip. This is to avoid the chip gaps present in our templates that get blueshifted into the fit by ppxf.
- (b)
  Pixels in each spectrum within 2.5 Å of skylines with an intensity larger than 0.5 from Rousselot et al. (2000).¹⁶
- (c)
  Pixels in the following ranges that were repeatedly found to be poorly fit by the A-LIST templates: 15360–15395, 16030–16070, and 16210–16240 Å (see Section 5.5.5).
2.
Calculate eS/N: Using ppxf, fit each chip individually with an mdegree of 40, which uses a multiplicative polynomial of order 40 to represent the continuum of the spectrum. This value was chosen as it was found to produce fits with no structure in the residuals. Fit the data with templates interpolated at the discrete A-LIST grid points. Use the χ² from this fit to calculate the rank-based estimate of the standard deviation of the residuals over the uncertainties, q:
$\begin{eqnarray}&&{\sigma }_{G}=0.7413({q}_{75}-{q}_{25}),\end{eqnarray} \tag{ 3 }$
where q₇₅ and q₂₅ are the upper and lower quartiles of q. Then, multiply the uncertainty arrays by σ_G to make the χ² value for this fit equal 1. Use these scaled uncertainties for the final fit in step 3 and to calculate the eS/N of the spectrum.
3.
Perform the fit: Initialize emcee routine bounds and other parameters. Set the walker bounds as follows: Velocity Dispersion between 60 and 200 km s⁻¹, [M/H] between −1.0 and 1.0, [α/M] between 0 and 0.5, and age between 7 and 12 Gyr. Initialize walkers randomly between these bounds. Set the mdegree for the ppxf fit to 2. Run the emcee routine using StretchMove¹⁷ with the desired number of kinematic and stellar population parameters; this can be four or five (see next paragraph). Use 30 walkers for 200 iterations and a burn-in period of 125 iterations.¹⁸ Take the final results for each parameter as the median values in the ensemble chain after the burn-in period.

For bins comprised of more than one fiber position, we fit the spectra with a 4D MCMC routine to determine the mean dispersion, metallicity, α abundance, and age of the stellar population. The RV is fixed to be the eS/N weighted mean of the O18 velocities of the fiber positions in the bin.

For bins that contain just one fiber position, we fit the spectra with a 5D MCMC routine to determine the velocity in addition to the four parameters above. For these fits, the velocity walkers are initialized in a uniform distribution ±15 km s⁻¹ around the O18 velocity for that fiber position.

A fit to a bin near the Voronoi target eS/N of 70 is shown in Figure 4; just the blue chip is shown. Bin 48 consists of a single fiber due west of the center of M31. The χ² of this fit was 1.99. The measured RV is −367 km s⁻¹ and the dispersion is 126 km s⁻¹.

4.3. Measurement Error Determination

The errors quoted in Section 5 are determined through a jackknife resampling procedure where data points are systematically removed from the sample and the rest are reanalyzed to determine the variations between data points (Tukey 1958). For our data the final spectrum for single-fiber-position bins is made by stacking visit spectra, whereas the final spectrum for multi-fiber-position bins is made by shifting and stacking the spectra for fiber positions in that bin. We vary these constituent spectra in our jackknife process to determine the errors. We calculate the jackknife-determined errors and an error scaling term, which is the median of the ratio of jackknife-determined errors to MCMC-determined errors, for bins with five or more constituent spectra. For all other bins, we multiplied the MCMC-determined errors by the scaling term to get the quoted error on each measurement. In general, our jackknife-determined errors are larger by a factor of 1–5 than the MCMC-determined errors, with no strong correlation between eS/N and error size. Jackknife-determined errors measure the intrinsic variation in the sample, rather than the MCMC-determined errors which are statistical, hence the former being larger. The median scaling used for each parameter was 2.13 for velocity, 2.43 for dispersion, 2.47 for metallicity, and 1.93 for α abundance.

5. Results and Discussion

5.1. Mean Trends in the Bulge and Disk

The primary results of this work are presented in Figure 5. This figure shows the light-weighted average kinematics, chemistry, and ages of the stellar populations in the bulge and inner disk. The bottom-right panel of Figure 5 shows the eS/N of each bin. Median errors for each parameter are shown as black lines in the color bars in Figure 5.

In Sections 5.1.1–5.1.5, we discuss the mean trends for each parameter in Figure 5. In Section 5.2, we discuss the central five fiber positions that likely host significant nuclear populations. Section 5.3 describes our results for abundance gradients in the bulge. We compare our results to the literature in Section 5.4 and examine other analysis choices we investigated in Section 5.5.

5.1.1. Radial Velocity

The top left panel of Figure 5 shows our RV results for M31, which have typical errors of 3.8 km s⁻¹. Our results are consistent with a systemic velocity of −300 ± 4 km s⁻¹ for M31 (de Vaucouleurs et al. 1991; McConnachie et al. 2005). In Figure 6 we LOESS smooth these results in the bulge. The raw data is shown as colored circles with black outlines, and the smoothing is shown as the red, white, and blue backgrounds. The lime line is a twist in the RV profile along the direction of the bar. This is a well-known bar effect (e.g., Stark et al. 2018) and has previously been seen by O18. See Section 5.4.3 for further discussion.

**Figure 6.** LOESS-smoothed RV in the bulge. Fiber positions used in the smoothing are indicated by black-outlined circles, colored by the value of the parameter used in the smoothing. For velocity, we smooth using a fraction of 0.5 and a degree of 1. The lime line on the left-side panel indicates where the smoothed velocity equals −300 km s⁻¹, showing the "s" shaped velocity profile that is an effect of M31's bar. As in Figure 5, the solid and dashed black lines indicate the disk and bar major axes, respectively.
Download figure:
Standard image High-resolution image

5.1.2. Velocity Dispersion

Results for velocity dispersion are shown in the top-right panel of Figure 5. We show that the bulge has a dispersion >130 km s⁻¹ that decreases roughly linearly (∼−10 km s⁻¹ kpc⁻¹) along all position angles. The inner disk of M31 has a velocity dispersion ≤100 km s⁻¹. Our RV dispersion map peaks at ∼160 km s⁻¹ off-center to the southeast. The dynamics of the very center are complex (e.g., Lockhart et al.2018), and our results here are discussed further in Section 5.2.

5.1.3. Metallicity

Our metallicity results are presented in the left middle panel of Figure 5. The mean metallicity of the bulge is generally subsolar ([M/H] = −0.149 dex); the inner disk shows near-solar metallicity ([M/H] = −0.006). The nucleus is enhanced relative to the rest of the bulge ([M/H] = −0.055). We emphasize that these metallicities are the H-band-light-weighted mean metallicities of what is likely a range of population metallicities in these regions of the galaxy.

We find a clear pattern of metal deficiency in the bulge perpendicular to the bar's major axis, which is shown in Figure 7. This structure was uncovered by Gajda et al. (2021) as so-called metallicity deserts. S18 determined that the bar of M31 is metal-enhanced relative to the off-bar regions, but we do not find as clear of a signature. We attribute this to a lack of spatial resolution sufficient to entirely negate the effects of stochasticity (see Section 5.5.1 for more). We examine the metallicity gradients in Section 5.3.

5.1.4. α Abundance

The right middle panel of Figure 5 shows our map of the α abundance in M31. The bulge of M31 is enhanced in α abundance with a median [α/M] = 0.281. The results in this region exhibit statistically significant variations in the α abundance of 2.5–3 times our jackknife-determined errors. These errors accurately quantify observational uncertainties, but do not account for stochastic variation in the stellar population from point to point. These statistical fluctuations do not appear to trace a coherent spatial pattern. This is consistent with the bulge having no α abundance substructure, but rather having a spread in α abundance that is evident from stochastic variations in the fraction of light from stars of different α abundance that are in each fiber.

The nucleus is more α-poor than the rest of the bulge, but still super-solar (median [α/M] = 0.179), and the central fiber position has [α/M] = 0.006. The disk is also enhanced (median [α/M] = 0.274), but shows more clear spatial substructure. The LOESS-smoothed map of the bulge is shown in Figure 8. Among the features visible are lower α abundances along the minor axis and asymmetry along the bar major axis; a clear signature of this α-poor region to the northwest is also seen by S18. It could be that this is the end of the bar protruding out from the bulge into the disk; however, an equivalent feature is not seen as strongly on the opposite side of the bulge. Analysis of the data with symmetric binning across the disk major axis was unable to produce a similar deficiency on both sides of the bar. Dust lanes in the region could obscure the stellar populations to the southwest and increase the apparent enhancement, but further analysis of this asymmetry is beyond the scope of this work.

**Figure 8.** LOESS-smoothed α abundances in the bulge, as in Figure 6. The bulge is uniformly enhanced relative to the disk by roughly 0.1 dex. The central fiber position has an abundance of 0.063 ± 0.008 dex, and the other four nuclear fiber positions are also deficient. The region to the northeast is also α-poor and is possibly the end of the bar. The parameters used for smoothing the α abundance are fraction = 0.5 and degree = 2. As in Figure 5, the solid and dashed black lines indicate the disk and bar major axes, respectively.
Download figure:
Standard image High-resolution image

5.1.5. Stellar Age

The bottom-left panel of Figure 5 shows our results for the age of the stellar populations in M31. We find the bulge to be uniformly 12 Gyr old. Our model grid extends only to 12 Gyr, so our results are a lower bound. In the disk, we find most bins to be old as well, with a few exceptions. We note that the A-LIST models are not particularly sensitive to age, as ages for spectra used to create A-LIST are from astroNN (Mackereth et al. 2019), which trains a neural network on asteroseismic data to determine ages from spectra that are accurate to ∼30%–35%. Additionally, two models with mono-abundance and ages only off by 1 Gyr are very similar (median flux differences below 0.001). Therefore, we trust the bins with ages >11 Gyr; the median error on these ages is 0.31 Gyr. Bins with ages <11 Gyr have elevated errors (median 4.35 Gyr) and higher χ²s, as the MCMC routine cannot locate a maximum likelihood age for these few bins. Analysis of these bins with age fixed to 12 Gyr generally decreased metallicity and α abundance, but these results remain in line with expectations (see Section 5.5.4 for more details). We identified three bins with ages <8 Gyr and χ² > 3.5 whose fixed-age χ² was greatly improved and will utilize those results for our analysis.

5.2. Nuclear Regions

The nuclear star cluster (NSC) of M31 dominates the light out to 5'' (≃0.02 kpc; Kormendy & Bender 1999; Peng 2002), so our central fiber position is sensitive to only NSC stars. We have four other fiber positions from 11'' to 12'' (≃0.045 kpc) that are likely also sensitive to nuclear stars. Thus, we call these five fiber positions the nuclear fiber positions. The stellar populations in this region are distinct from the bulge and inner-disk populations. We find these fiber positions to be metal-enhanced ([M/H] ≃ −0.05) and α deficient ([α/M] ≃ 0.18) relative to the rest of the bulge.

The central 4 pc of M31 was investigated by Lockhart et al. (2018) who show that the dynamics of this region are complex. They find a quickly rotating eccentric disk here with asymmetric velocity magnitudes on either side. They find the dispersion to be peaked at 381 ± 55 km s⁻¹ 0 farcs 13 from the center, though this high-dispersion area is only ∼0 farcs 2 in diameter. APOGEE fibers have an on-sky diameter of 2'', implying our result for this fiber position should be the average value of their data. They find that a larger area is redshifted from −300 km s⁻¹ than is blueshifted. Additionally, the majority of their dispersion plot is near 150 km s⁻¹. This lines up with our result for the central fiber position, which has an RV of $-{225.56}_{-0.459}^{+0.486}$ km s⁻¹ and a velocity dispersion of ${151.91}_{-0.307}^{+0.364}$ km s⁻¹. This indicates the NSC of M31 is toward the Sun slower than the rest of M31. O18 found a ring of increased dispersion (up to ∼180 km s⁻¹) a few hundred arcseconds (∼1 kpc) from the nuclear regions, which have a lower dispersion (around 150 km s⁻¹). Saglia et al. (2010) show that the velocity dispersion reaches a minimum around 5'' from the center before rising again further out. This is consistent with our analysis, which indicates the peak dispersion being off-center by a few hundred arcseconds. We would likely see a ring of increased dispersion too with better spatial sampling. Both our results and those of O18 are subject to stochasticity. S18 show an increase in metallicity to [Z/H] = 0.35 in the very center of M31 but no significant decrease in α abundance. The central 5'' hosts a mostly old stellar population (7–13 Gyr Saglia et al. 2010).

We find the NSC of M31 to be old (≃12 Gyr) and to have an RV of −225.56 ± 3.23 km s⁻¹, velocity dispersion of 156.91 ± 2.45 km s⁻¹, [M/H] of −0.055 ± 0.013, and [α/M] of 0.063 ± 0.008 (the most α-poor fiber position in our sample). Any further interpretation of these results is beyond the scope of this paper.

5.3. Abundance Gradients

We derived gradients for the mean metallicity and α abundance along three axes in the bulge: the disk major and minor axes, as well as the bar major axes (Figure 10). Gradients are calculated by assigning each fiber position the abundance calculated for its bin. We isolate the fiber positions within 0.01° on-sky of each axis (Figure 9) and calculate the average radius for these fiber positions. This radius is used to derive the gradient by fitting a straight line to the plot of abundance versus radius. We only fit to fiber positions falling within the bulge region. The central fiber positions (red open circles in Figure 10) are located in the nuclear regions of M31 (see Section 5.2) which is made up of different stellar populations than the bulge and inner disk. These five fiber positions are excluded from the gradient calculations.

**Figure 9.** This figure demonstrates how the gradients in Figure 10 are calculated. Fiber positions along each axis are assigned the abundance for their bin, then the radius for each fiber position is calculated. The radius for the bin is the mean radius for fiber positions in that bin along that axis. Unused fiber positions are denoted as gray dots.
Download figure:
Standard image High-resolution image

**Figure 10.** Metallicity (top) and α abundance (bottom) gradients for the fiber positions shown in Figure 9. The calculated gradients (blue lines) are the slopes of a straight line fit to the data taking abundance errors into account and are quoted in the title of each panel. Red open circles denote the five nuclear fiber positions, which are not included in gradient calculations. Red circles denote Voronoi bins with only one fiber position, black circles denote Voronoi bins with more than one fiber position. The radius error bars correspond to the minimum and maximum radius of fiber positions along that axis in each bin, and the abundance errors were determined by the jackknife procedure in Section 4.3. The shaded blue regions show the errors on our gradient fits. The dashed lines show the region along each axis where the bulge makes up >50% of the total light, dotted lines indicate this region for >90%.
Download figure:
Standard image High-resolution image

We find strongly positive metallicity gradients along the disk and bar major axes (0.043 and 0.027 dex kpc⁻¹, respectively) with a flattening beyond ∼3 kpc around solar metallicity. The spread in abundance for red points (corresponding to single-fiber-position bins) is likely stochastic in nature (see Section 5.5.1). The gradient along the disk minor axis is noisy but still generally positive (0.037 dex kpc⁻¹). There is a notable decrease in metallicity of ∼0.1 dex from the nuclear fiber positions to bins between 1 and 3 kpc, corresponding to the metallicity desert perpendicular to the bar.

Gradients for the α abundance are much more scattered than those for metallicity. They are negative, at ∼−0.017 dex kpc⁻¹ along the disk major and minor axes and −0.0085 dex kpc⁻¹ along the bar. The single-fiber-position bins (generally associated with the bulge) show higher maximum enhancement than the multi-fiber-position bins in the disk, indicating overall enhancement of the bulge relative to the disk. The effects of stochasticity are still apparent. The bottom-right panel of Figure 10 clearly shows the asymmetry along the bar, which affects the gradient calculated here, decreasing its magnitude.

Our abundance gradients in the bulge region are likely the result of the superposition of a metal-rich, α-poor disk and a less metal-rich, α-rich disk. Light in the central regions will be dominated by the bulge, but the relative fraction of this component's contribution to the total light will decrease steadily with distance. As a result, the mean metallicity gradients are steeply positive, but the intrinsic gradients for each component (classical bulge, b/p bulge, inner disk, etc.) may not be.

5.3.1. Bulge–Disk Decomposition

If we assume that the metallicity gradient in the disk is the result of the superposition of a subsolar metallicity bulge and solar metallicity disk, then we can use the photometric bulge–disk decomposition from (Courteau et al. 2011) to determine the uniform metallicity of the bulge and inner disk. We use their Equations (2) and (4) to determine the fractional intensity of light coming from the bulge as a function of radius. We fit their Sérsic bulge and exponential disk model to our data and find our data are well represented by a bulge with [M/H] = −0.219 ± 0.008 and a disk with [M/H] = 0.0584 ± 0.025. This relationship along the disk major axis is shown in Figure 11.

**Figure 11.** Decomposition of the bulge metallicity gradient along the disk major axis. Points were selected in the same as the top left panel of Figure 10. Nuclear fiber positions are not included in the fit. All other points are colored by the fraction of their light coming from the bulge. The blue line shows how our modeled metallicity changes as a function of radius.
Download figure:
Standard image High-resolution image

5.4. Comparison to Literature

5.4.1. Chemodynamics

The most similar analysis to ours found in the literature is from O18 (kinematics) and S18 (abundances and ages). As mentioned in Section 2.4, these data are from the VIRUS-W optical IFU. These IFU data are much more densely packed in the bulge than our APOGEE data, but their results are still subject to stochastic effects. We calculated the difference between abundance in neighboring bins for all of the bins in their central region, excluding the spokes. The standard deviation of this distribution was σ[Z/H] = 0.136 and σ[α/Fe] = 0.108. Despite this, their spatial sampling was dense enough to uncover detailed spatial structure. O18 performed full-spectrum fitting using ppxf with 41 kinematic standard stars and constructed detailed kinematic maps of this region and found strong evidence for a bar from cuts of the kinematics through the disk major axis. S18 examined the same data using Lick indices to determine the chemistry of M31. They found the inner 2.3 kpc (where the light is dominated by stars in the classical bulge) to be old (12.9 Gyr), metal-rich ([Z/H] up to 0.35, averaging to be slightly super-solar) and α enhanced ([α/Fe] of 0.28). The b/p bulge is the same age at 13.1 Gyr, but is slightly subsolar with [Z/H] = −0.04 and less α abundant at [α/Fe] = 0.25. The bar is slightly younger, but has similar abundance to the classical bulge, and is metal-rich.

In general, our results agree with this analysis. Our kinematics show strong evidence of bar effects in both the "s"-shaped velocity curve and the velocity versus position angle plot. The spatial structure we find in our abundances is also correlated with those found in the above. They see enhancement in metallicity along the spurs of their data that extend into the disk along the disk major axis, and depression for the spurs along the disk minor axis, much as we do. We also find the metallicity deserts on either side of the bar, though we hesitate to identify this structure as an enhancement along the bar. Our data are too sparsely sampled spatially, and as such are too subject to stochastic effects, to say this for sure (see Section 5.5 for further discussion). We find our mean metallicities are lower than those of S18 in the inner ∼2 kpc by 0.1–0.15 dex, but that this pattern flips beyond 2 kpc, and instead our results are 0.05–0.1 dex higher than S18's. S18 notes near-uniform enhancement in α abundance in the bulge with the exception of the nuclear regions, much as we do (excepting the statistical fluctuations), and see a reduction in abundance along the northeast side of the bar. They also find the bulge to be older than 12 Gyr.

5.4.2. Gradients

Many studies have analyzed the metallicity gradient in the disk of M31. Studies of H ii regions in the disk of M31 (∼4–20 kpc) have indicated a slightly negative metallicity gradient of −0.023 ± 0.002 dex kpc⁻¹ (Zurita & Bresolin 2012). Analysis of planetary nebulae in M31's outer disk shows a more shallow metallicity gradient around −0.011 dex kpc⁻¹ from 18–43 kpc (Kwitter et al. 2012). The PHAT survey shows a negative metallicity gradient from 4–20 kpc that matches that of the H ii analysis at −0.020 ± 0.004 dex kpc⁻¹ (Gregersen et al. 2015). The gradients found by these studies are consistent with simulations of the merger history of M31 that indicate it had a recent major interaction, potentially with M32, 2–4 Gyr ago (e.g., D'Souza & Bell 2018; Hammer et al. 2018).

These studies have analyzed and drawn conclusions from metallicity gradients in the disk >4 kpc from the center of M31. Gradients inside this radius may result from different physical processes, as they are inside a different structure. We do not interpret our positive gradients in the context of M31's disk formation or merger history. Instead, as was shown in Section 5.3.1, they can be modeled as an effect of the changing light contribution from the bulge and disk structural components.

The gradient in the bulge was studied by S18, who calculate metallicity and α abundance gradients in the central 500'' (roughly 1.9 kpc). They call the inner ∼78'' (∼0.3 kpc) the classical bulge and calculate steeply negative and shallow positive gradients for [Z/H] and [α/Fe], respectively. S18 find a shallow positive metallicity gradient out to 500'' for the bar and a shallow negative metallicity gradient off the bar. Both on- and off-bar regions show a shallow negative α abundance gradient. Gregersen et al. (2015) quote their gradient for the outer disk, but note that their results within 4 kpc show a positive metallicity gradient closer to our results (see their Figure 4, note that their results in this region are subject to crowding and completeness effects). Outside of 4 kpc, our results show a roughly flat gradient out to 7 kpc (see the top-right panel of Figure 10 for the best example), much more in line with previous studies.

5.4.3. Bar and Bulge Properties

The twist in the velocity field along the disk minor axis seen in Figure 6 is often called nonaxisymmetric motion or a noncircular velocity field. It is seen in perhaps 20%–50% of disk galaxies, most often in the gas motions (Haynes et al. 1998; Swaters et al. 1999; Kornreich et al. 2000; Andersen & Bershady 2013; Bloom et al. 2017), and identified as early as 1978 (e.g., Peterson et al. 1978; Sanders & Tubbs 1980). The strongest deviations from circular gas motions in CALIFA (Sánchez et al. 2012) galaxies are associated with the presence of a bar (García-Lorenzo et al. 2015). The most detailed analysis of asymmetric stellar motions was done by Stark et al. (2018). They investigate ∼2800 MaNGA galaxies and find this feature (which they call inner bends) in ∼20%–25% of their sample. They find it is most closely associated with having a strong bar, however roughly 45% of the strongly barred galaxies in their sample show no distortion. They note these numbers may be skewed by the difficulty of visually identifying bars.

Our analysis finds inconclusive evidence of the metal-rich bar seen in S18. However, our results may line up well with the stellar chemodynamics of the MW's bar and bulge. Rojas-Arriagada et al. (2020) found that the MW's bulge can be split into three spatially overlapping populations at [Fe/H] = +0.32, −0.17, and −0.66 dex in order of decreasing population concentration. They confirm the interpretation that metal-rich stars in the disk secularly evolve into the boxy/peanut bulge and argue that the intermediate metallicity population could have formed in situ at high redshift. The metal-poor stars could be associated with the early thick disk. Similarly, Wegg et al. (2019) find that bulge stars on the most bar-like orbits have the highest metallicities ([Fe/H] = 0.3) and inner-disk stars have lower metallicity ([Fe/H] = 0.03). These results, as well as those of Queiroz et al. (2021), line up well with ours, but the nature of the MW bulge is still up for debate: Bovy et al. (2019) define the bar by its kinematics and find it to be more metal-poor and α-rich on average compared to the inner disk.

An N-body simulation of the bar's evolution from Fragkoudi et al. (2018) indicates the bar is made up of metal-rich ([Fe/H] > 0.0) and metal-poor ([Fe/H] < 0.0) populations whose fractional contribution to the bar varies spatially. The metal-rich population makes a longer, more narrow bar, whereas the metal-poor population is shorter and wider. This results in a nonuniformly enhanced bar that is metal-poor in the center and metal-rich at the ansae. The regions perpendicular to the bar are also metal-poor. This final result may line up with ours best, as we find the nonnuclear fiber positions in the center to be generally the most metal-poor, as well as a positive radial metallicity gradient along the bar.

5.5. Examining Robustness to Analysis Choices

Above, we have presented an analysis whose details depend on the methodological choices outlined in Section 3, which we will call our default method. To test our robustness to these methodological choices, we reanalyzed our data in several different ways. In this section, we discuss why we chose our default method and the robustness of the features seen in our results to these methodological choices.

5.5.1. Criteria for Our Presented Results

Our primary focuses when choosing which analysis variation to present were to ensure that results were consistent between neighboring bins and that gradients were reasonably tight.

To quantify bin-to-bin consistency we calculated the difference between abundances in neighboring bins (as in Section 5.4.1 for S18) and determined the standard deviation of the distribution. For our default method, we find σ[M/H] = 0.08, and σ[α/M] = 0.05. This was similar to or less than our results for the other analysis variations described in the following subsections, which varied from 0.10 to 0.16 in metallicity and 0.03 to 0.06 in α abundance.

To determine the tightness of gradients we calculated the residuals from our gradient fits along the disk major axis for each analysis variation. The standard deviation of our residuals for our default method are σ[M/H] = 0.05 and σ[α/M] = 0.03. These numbers for the other variations ranged from 0.04–0.17 in metallicity and 0.03–0.05 in α abundance. Our presented results have the best combination of these statistics.

These statistics capture the stochasticity in our results, which is largely a combined effect of our irregular spatial sampling and APOGEE's small fiber diameter (2'' on-sky). Since APOGEE's fibers are so small, we are sensitive to fewer stars in each position than those of O18 and S18.

5.5.2. Binning Effects

The binning scheme used in our final results was motivated to obtain as many radially narrow bins as possible at or above an eS/N of 60 (see Section 4.1). What we lose by using this routine, however, is the ability to identify more detailed azimuthal substructure in the disk, as well as symmetry along the disk minor axis. We created a different binning scheme that maintained symmetry across both the disk major and minor axes, but sacrificed radial narrowness and eS/N. Analysis of this binning strategy produced similar maps to Figure 5 and gradients to Figure 10, but introduced more stochasticity into our measurements. Bins that were next to each other azimuthally often did not have similar abundance. This strategy also did not show an equivalent to the α deficiency seen on the northeast side of the bar on the southwest side.

5.5.3. Continuum Treatment

The handling of the continuum shape of our spectra has an effect on the scaling of our abundance results. There are two main ways to handle the continuum using ppxf. The first one, which we use in our default method, is to manually normalize the spectrum using a running median. The second is to use the mdegree keyword in ppxf, which applies a multiplicative Legendre polynomial to the model spectrum to mimic the continuum shape of the data. ppxf has full control over the coefficients of this polynomial; mdegree simply tells ppxf what order polynomial to apply.

We tested the effects of continuum and mdegree on our results using our models and real data. To test the models, we took the templates with metallicity between −0.6 and +0.4 and α abundance of 0.1 or 0.3, and added noise and broadening to simulate our observed spectra. This was done by Gaussian broadening the templates to a dispersion consistent with our results for M31, and adding noise corresponding to the uncertainty array of a random binned spectrum. We call these templates artificially noisy templates, or ANTs. For one test we did no modification to the continuum of the ANTs; for the other test, we added in the continuum from the binned spectrum we used to add noise (for reference, our continuum calculations are described in Section 2.2.1). The ANTs were analyzed using values of mdegree between 0 and 10. We found that the best-fitting value of dispersion, metallicity, or α abundance across all values tested did not change with changing mdegree, and got back our input results within errors. For the ANTs with continuum added, the dispersion increased by roughly 20 km s⁻¹ from mdegree of 0–10, and the metallicity increased by 0.2–0.5 dex, depending on the input metallicity.

We performed a similar set of tests on our data, handling continuum three different ways: continuum normalized with a second-order multiplicative polynomial (mdegree = 2), continuum normalized with mdegree = 4, and noncontinuum normalized with an mdegree of 6, 8, or 10, depending on the number of pixels in the chip.

In general, the structures seen in our results were consistent with each method, though the stochasticity in our individual measurements increased with increasing mdegree. For instance, the statistics calculated in Section 5.5.1 were generally higher when using an mdegree of 4. The neighboring bin inconsistency for metallicity increased to 0.0614, and the gradient residuals increased to 0.118. We also found that the scaling of our dispersion and metallicity was correlated with mdegree. Using an mdegree of 4 resulted in mean metallicities that were 0.151 ± 0.066 dex higher than using an mdegree of 2. Dispersions were ∼5 km s⁻¹ higher, and α abundance was minimally changed. However, using the higher values of mdegree without continuum normalizing resulted in dispersions ∼50 km s⁻¹ and metallicities ∼0.5 dex higher than our default method. These increased metallicities fell outside of the A-LIST grid, and our interpolated model had to extrapolate to create these templates.

Structurally, we found that these higher mdegree fits resulted in a stronger, but still rather weak signature of bar metallicity enhancement and a more uniformly α-enhanced bulge. However these metallicity results were also more stochastic everywhere, as were the ages. Ultimately, we chose to use an mdegree of 2 to reduce model complexity, provide more robust results, and bring our results more in line with the A-LIST parameter space and previously published kinematics and abundances.

5.5.4. Fixing Age to 12 Gyr

As stated previously, our models are not particularly sensitive to stellar age differences of a few gigayears. We find that most of our bins end up with ages of 12 Gyr (the edge of our model grid), with a few exceptions. We performed our primary analysis method with the restriction that every bin must be fixed to 12 Gyr and found minimal changes to our results. There are 10 bins with ages below 11 Gyr. Fixing the age in these bins caused a decrease in metallicity and α abundance in every bin. The average change in abundance was Δ[M/H] = −0.107 and Δ[α/M] = −0.039. Additionally, the majority of bins with low ages have elevated errors in metallicity, α abundance, and age relative to the 12 Gyr bins.

5.5.5. Effects of Masking

Our 85th percentile cutoff for masking out skylines as described in Section 2.2.1 had the potential to mask out the deepest sections of absorption features in our spectra. To test the effects of our masking, we ran two tests on the ANTs from Section 5.5.3. For the first tests, we performed the fits using the mask arrays from the binned spectra used to make the ANTs. For the second test, we masked out no pixels other than the misfit regions and chip gaps. In both cases we found minimal bias (−0.01 and −0.02 dex, respectively) of our output metallicities from the input ones, indicating our results are not affected by our masking routine.

We also found several regions in our spectra that were consistently poorly fit across the majority of our bins. These regions fell into two categories: known skylines and poorly modeled lines.

We were not able to determine every skyline using the strategy outlined in Section 2.2.1, and sometimes for skylines that we did detect, we did not mask out all of the affected pixels. To account for this, we mask out pixels within 2.5 Å of known skylines prior to fitting. We utilize the skylines from Rousselot et al. (2000) with intensities >0.5.

The three larger regions that are masked are due to poor data-model matching. We were consistently unable to fit these regions across a majority of our spectra. The region from 15360–15395 Å is a deep, broad line whose depth was not well fit. The region from 16030–16070 Å is another deep line, parts of which were sometimes masked by our skyline identification and other times not. Lastly, the region from 16210–16240 Å is a series of smaller lines with inconsistent fits from spectrum to spectrum.

6. Summary and Conclusions

In this paper we analyzed 963 high-resolution, NIR, integrated-light APOGEE spectra of the bulge and inner disk of M31. We spatially binned these spectra into 164 bins with an eS/N ≳ 60. These binned spectra were fit with interpolated SSP model spectral templates from A-LIST (Ashok et al. 2021). The results include maps of the mean RV, velocity dispersion, metallicity, α abundance, and age of the stellar populations in the inner ∼7 kpc of M31 (Figure 5). We also computed mean metallicity and α abundance gradients along the disk major and minor axes and the bar major axis. A summary of our primary conclusions is below.

1.
Abundance gradients: We observe a positive mean metallicity gradient in the bulge region (inner ∼4 kpc). We model this gradient two ways: with a linear model (Section 5.3) and as a superposition of two components with uniform metallicity (Section 5.3.1). First, the linear gradient is steepest along the disk major axis at Δ[M/H] = 0.043 ± 0.021 dex kpc⁻¹. The α abundance gradients in this region are negative, though shallower, and noisy. The linear gradient along the disk major axis is Δ[α/M] = −0.018 ± 0.016 dex kpc⁻¹. Second, the superposition gradient models the increasing metallicity as the changing light contributions from a uniformly metal-rich inner disk ([M/H] = 0.0584 ± 0.025) and less metal-rich bulge ([M/H] = −0.219 ± 0.008).
2.
Abundance features: The mean metallicity of the bulge of M31 is subsolar ([M/H] = $-{0.149}_{-0.081}^{+0.067}$ ), and the disk is, on average, nearly solar ([M/H] = $-{0.023}_{-0.052}^{+0.050}$ ). The metallicity on either side of the bar is subsolar ([M/H] ∼ −0.2), which is less than the mean metallicity of the entire bulge, and a potential indicator of a metal-enriched bar. The bar of M31 has been found to be enhanced in metallicity relative to the rest of the bulge (S18), but we are unable to clearly identify this pattern. It should be noted that our results in the bar are stochastic, making it difficult to clearly identify fine substructure in the bulge (see Section 5.5.1).The bulge of M31 is enhanced in α elements ([α/M] = ${0.281}_{-0.038}^{+0.035}$ ). Statistical fluctuations are noted in results for individual fiber positions in the bulge, but these fluctuations do not appear to trace a coherent spatial pattern. The disk is similarly enhanced on average ([α/M] = ${0.274}_{-0.025}^{+0.020}$ ) but has a more defined structure. We find a region of decreased α abundance to the northeast of the bulge along the bar major axis. This could be the end of the bar sticking out of the bulge above the disk, but further analysis is needed to interpret this feature.
3.
Kinematics: The line where the RV is equal to the systemic velocity of M31 (−300 km s⁻¹) is twisted in the direction of the bar major axis in the bulge. This feature is strongly associated with the presence of a bar (e.g., Stark et al. 2018). The velocity dispersion of M31 does not peak in the very center. Our results find a peak ∼1 kpc from the center to the southeast. This is consistent with the findings of O18, who identify a ring of increased dispersion this far from M31's center. The dispersion decreases with distance from M31's center at a rate of ∼10 km s⁻¹ kpc⁻¹.
4.
Stellar age: Both the bulge and disk are uniformly old at ≥12 Gyr, with the exception of a small number of disparate bins. Fixing the age of these disparate bins to 12 Gyr decreased both the metallicity and α abundance results for these bins to be more in line with adjacent bins.
5.
Nuclear region: M31's nucleus contains a stellar population distinct from the bulge and the disk. The nucleus is metal-enhanced ([M/H] ≃ −0.05) and α deficient ([α/M] ≃ 0.18) relative to the rest of the bulge. We find the very center of M31 is the most α-poor of any region (0.006 dex), and is moving away from the Sun relative to the rest of M31, with a velocity of $-{225.56}_{-0.459}^{+0.486}$ km s⁻¹.

These results are robust to a number of different methodological choices. Our findings agree well with those of previous resolved and integrated-light studies of M31 and are consistent with our understanding of the MW bulge and inner disk. The data used to produce the maps in Figure 5 are published with this work as formatted in Table 1 in Appendix A.

Future work will reanalyze these spectra using an expanded methodology that will fit two spectral templates to the data in order to characterize multiple kinematically distinct co-spatial stellar populations in M31. Such populations could include the separate classical and boxy/peanut-shaped components in the bulge and the thick and thin disks, with a particular emphasis on identifying a trend like the so-called α-bimodality that is seen in the MW (e.g., Haywood et al. 2013; Hayden et al. 2015).

Acknowledgments

A very heartfelt thanks is extended to the APO plate pluggers who figured out how to plug our atypical plate designs (see Figure 12)—their work was vital to this project. We thank the referee whose insightful comments and suggestions have improved the clarity of the paper. The authors thank Nicholas Boardman and Jianhui Lian for their helpful discussions on this work. B.J.G., G.Z., and A.S. acknowledge support for this project was provided by NSF grant No. AST-1911129.

Funding for SDSS-IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS website is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

The Digitized Sky Survey was produced at the Space Telescope Science Institute under U.S. Government grant No. NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. The plates were processed into the present compressed digital form with the permission of these institutions.

Appendix A: Data Tables

The results of our analysis will be published in electronic data tables taking the same form as Table 1. These tables can be used to create maps similar to those in Figure 5.

Table 1. Description of Bins and Chemodynamics Results from Our Analysis

Bin	# of Fibers	R.A.		Decl.		Radius		Position Angle		χ²	eS/N	Velocity		Dispersion		[M/H]		[α/M]		Age
BIN	NFib	MinRA	MaxRA	MinDec	MaxDec	MinRad	MaxRad	MinPos	MaxPos			v	e_v	d	e_d	m	e_m	a	e_a	A	e_A
		(deg)	(deg)	(deg)	(deg)	(kpc)	(kpc)	(deg)	(deg)			(km s⁻¹)	(km s⁻¹)	(km s⁻¹)	(km s⁻¹)	(dex)	(dex)	(dex)	(dex)	(Gyr)	(Gyr)
1	1	10.685	10.685	41.269	41.269	0.0	0.0	0.0	0.0	12.251	419.807	−225.559	3.235	151.913	2.449	−0.055	0.013	0.063	0.008	11.988	0.006
11	1	10.645	10.645	41.239	41.239	0.832	0.832	186.791	186.791	1.908	71.253	−368.24	10.809	122.906	5.11	−0.028	0.152	0.235	0.011	11.926	0.024
15	1	10.65	10.65	41.218	41.218	1.056	1.056	169.237	169.237	1.553	50.449	−399.717	8.577	133.01	4.683	−0.158	0.069	0.352	0.025	7.425	4.283
37	1	10.639	10.639	41.264	41.264	1.954	1.954	224.39	224.39	2.466	103.435	−337.403	2.498	139.653	4.221	−0.23	0.049	0.305	0.011	11.924	0.024
57	101	10.427	10.675	41.03	41.27	3.578	5.218	145.944	232.709	3.471	74.998	−443.396	nan	83.272	5.09	0.013	0.049	0.26	0.016	11.971	0.059
69	14	10.603	10.65	41.303	41.356	4.781	5.56	263.392	305.275	2.449	79.582	−280.668	nan	92.891	4.041	−0.058	0.042	0.304	0.018	11.962	0.024
72	48	10.386	10.568	41.073	41.29	6.101	8.106	189.877	245.273	2.566	58.017	−410.188	nan	82.715	5.449	0.136	0.132	0.274	0.013	12.0	nan
78	2	10.656	10.66	41.317	41.317	3.232	3.385	297.8	301.01	1.874	72.098	−271.675	nan	111.532	1.765	−0.086	0.023	0.268	0.009	11.954	0.009
96	1	10.675	10.675	41.289	41.289	1.329	1.329	302.063	302.063	2.218	118.54	−268.25	2.928	131.968	2.13	−0.238	0.013	0.291	0.005	11.976	0.011
123	1	10.736	10.736	41.269	41.269	2.396	2.396	52.166	52.166	1.843	81.905	−274.832	5.836	131.064	1.118	−0.138	0.092	0.224	0.031	11.851	4.779
139	1	10.73	10.73	41.313	41.313	0.88	0.88	359.736	359.736	1.58	51.68	−242.241	6.168	124.011	5.113	−0.087	0.042	0.26	0.011	11.884	0.15
150	1	10.706	10.706	41.294	41.294	0.473	0.473	354.708	354.708	1.985	99.435	−248.027	5.946	132.642	4.73	−0.138	0.014	0.293	0.008	11.967	0.018
156	1	10.693	10.693	41.25	41.25	1.217	1.217	123.823	123.823	1.92	104.577	−318.641	3.931	153.07	3.635	−0.35	0.045	0.343	0.014	11.911	0.05
164	6	10.76	10.776	41.233	41.253	4.854	5.27	65.274	84.65	2.227	76.315	−296.462	nan	100.632	4.095	−0.077	0.018	0.259	0.03	11.975	0.025

Notes. The top row of the header describes the type of quantity given in the columns. The middle row gives the column name in the machine-readable version. Columns with a "e_" prefix correspond to the error bounds on our measurements. Units for each quantity are given in the bottom row of the header. Bins 57, 69, 72, 78, and 164 have nans for e_v because they are multi-fiber-position bins, and as such their velocities were determined from the average O18 velocity of all the fiber positions in each bin (see Section 2.4). Bin 72 has nan for e_A because its age was fixed to 12 Gyr (see Sections 5.1.5 and 5.5.4).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

Appendix B: Plug Plate

One of the five plates designed to take observations of M31 can be seen in Figure 12. This plate is unusual compared to the thousands of other APOGEE plug plates that are typically designed to observe individual stars across the whole plate, rather than a galaxy in the middle. The M31 fiber positions are the densely packed positions just left of the center of this plate. Other positions on this plate were used for sky and telluric flux measurements, taking spectra of M31 halo objects such as luminous blue variables or GCs (Sakari et al. 2016), and measurements of galaxies for MaNGA (Bundy et al. 2015).

Appendix C: Comparison of O18 Velocities to Measured Velocities

We compare the O18 velocities from Section 2.4 to ones we determined independently from our data to ensure the validity of our usage of the O18 velocities in this work. We Voronoi binned our fiber positions in the same way as Section 4.1, but for a target eS/N of 20, as we cannot determine reliable kinematics below this limit. To determine kinematics, we perform the general fitting routine from Section 4.2, but run the MCMC in just two dimensions, one for velocity and the other for dispersion. For each step in the MCMC, we fit the data with the three templates with the highest weighting in the preliminary fit, much the same as in Section 2.3. The result is a tight correlation between our results and the O18 velocities with a bias of −0.653 km s⁻¹ and a scatter of ∼19 km s⁻¹. Figure 13 shows our results for this exercise on the left and the velocities from O18 on the right. Note the data on the right are NOT the LOESS-smoothed O18 velocities from Section 2.4, but rather the raw data from the paper. This figure also shows the differences in our data coverage of M31's bulge and disk.

**Figure 13.** Comparison of parameters we determined (left) to those published in O18 (right). Note the differences in data coverage between the two data sets. Our data spans more of the inner disk continuously opposed to the spokes in O18's data set. The O18 data much more densely samples the bulge than ours.
Download figure:
Standard image High-resolution image