Elemental Abundances in M31: Properties of the Inner Stellar Halo

We present measurements of [Fe/H] and [$\alpha$/Fe] for 128 individual red giant branch stars (RGB) in the stellar halo of M31, including its Giant Stellar Stream (GSS), obtained using spectral synthesis of low- and medium-resolution Keck/DEIMOS spectroscopy ($R \sim 3000$ and 6000, respectively). We observed four fields in M31's stellar halo (at projected radii of 9, 18, 23, and 31 kpc), as well as two fields in the GSS (at 33 kpc). In combination with existing literature measurements, we have increased the sample size of [Fe/H] and [$\alpha$/Fe] measurements from 101 to a total of 229 individual M31 RGB stars. From this sample, we investigate the chemical abundance properties of M31's inner halo, finding $\langle$[Fe/H]$\rangle$ = $-$1.08 $\pm$ 0.04 and $\langle$[$\alpha$/Fe]$\rangle$ = 0.40 $\pm$ 0.03. Between 8--34 kpc, the inner halo has a steep [Fe/H] gradient ($-$0.025 $\pm$ 0.002 dex kpc$^{-1}$) and negligible [$\alpha$/Fe] gradient, where substructure in the inner halo is systematically more metal-rich than the smooth component of the halo at a given projected distance. Although the chemical abundances of the inner stellar halo are largely inconsistent with that of present-day dwarf spheroidal (dSph) satellite galaxies of M31, we identified 22 RGB stars kinematically associated with the smooth component of the stellar halo that have chemical abundance patterns similar to M31 dSphs. We discuss formation scenarios for M31's halo, concluding that these dSph-like stars may have been accreted from galaxies of similar stellar mass and star formation history, or of higher stellar mass and similar star formation efficiency.


INTRODUCTION
In the ΛCDM cosmological paradigm, L galaxies like the Milky Way (MW) and Andromeda (M31) form through hierarchical assembly (e.g., White & Rees 1978). Debris from mergers across cosmic time are deposited within the extended stellar halo, where they remain observationally identifiable because the phasemixing timescales in the outskirts of galaxies are long compared to the age of the universe (e.g., Helmi et al. 1999;Bullock & Johnston 2005). Simulations of stel-lar halo formation in MW-like galaxies have shown that the mass and accretion time distributions of progenitor dwarf galaxies can imprint strong chemical signatures in a galaxy's stellar population (e.g., Robertson et al. 2005;Font et al. 2006b;Johnston et al. 2008;Zolotov et al. 2010;Tissera et al. 2012), particularly in terms of [Fe/H] and [α/Fe]. Measurements of α-element abundance (O, Ne, Mg, Si, S, Ar, Ca, and Ti) and iron (Fe) abundance encode information concerning the relative timescales of Type Ia and core-collapse supernovae (e.g., Gilmore & Wyse 1998), such that galactic systems with different evolutionary histories will have distinct chemical abundance patterns. In this way, stellar halos serve as fossil records of a galaxy's accretion history. This theory has been extensively put into practice in the MW, where the differing patterns of [α/Fe] and [Fe/H] between its stellar halo and satellite dwarf galaxies have revealed their fundamentally incompatible enrichment histories (Shetrone et al. 2001;Venn et al. 2004).
Studies of the kinematics and chemical composition of individual stars in the MW have provided a detailed window into the formation of its stellar halo (e.g., Carollo et al. 2007Carollo et al. , 2010Nissen & Schuster 2010;Ishigaki et al. 2012;Haywood et al. 2018;Helmi et al. 2018;Belokurov et al. 2018Belokurov et al. , 2020. However, the MW is a single example of an L galaxy. Observations of MW-like stellar halos in the Local Volume have revealed a wide diversity in their properties, such as stellar halo fraction, mean photometric metallicity, and satellite galaxy demographics, which likely results from halo-to-halo variations in merger history (Merritt et al. 2016;Monachesi et al. 2016a;Harmsen et al. 2017;Geha et al. 2017;Smercina et al. 2019). In these studies, both the MW and M31 have emerged as outliers at the quiescent and active ends, respectively, of the spectrum of accretion histories for nearby L galaxies.
Owing to its proximity (785 kpc; McConnachie et al. 2005), M31 is currently the only L galaxy that we can study at a level of detail approaching what is possible in the MW. M31's nearly edge-on orientation (i = 77 • ; de Vaucouleurs 1958) provides an exquisite view of its extended, highly structured stellar halo (e.g., Ferguson et al. 2002;Guhathakurta et al. 2005;Kalirai et al. 2006b;Gilbert et al. 2007Gilbert et al. , 2009bIbata et al. 2007Ibata et al. , 2014McConnachie et al. 2018). Most notably, M31's stellar halo contains a prominent tidal feature known as the Giant Stellar Stream (GSS; Ibata et al. 2001), where the debris from this event litters the inner halo (Brown et al. 2006;Gilbert et al. 2007). Since the discovery of M31's halo (Guhathakurta et al. 2005;Irwin et al. 2005;Gilbert et al. 2006), its global metallicity, and kinematical properties have been thoroughly char-acterized from photometric and shallow (∼1 hr) spectroscopic surveys (Kalirai et al. 2006a;Ibata et al. 2007;Koch et al. 2008;McConnachie et al. 2009;Gilbert et al. 2012Ibata et al. 2014). However, it is only recently that Vargas et al. (2014a,b) made the first chemical abundance measurements beyond metallicity estimates in M31's halo and dwarf galaxies.
We have undertaken a deep ( 6 hour) spectroscopic survey using Keck/DEIMOS to probe the formation history of M31 from the largest sample of [Fe/H] and [α/Fe] measurements in M31 to date. This has resulted in the first [α/Fe] measurements in the GSS , the inner halo , and the outer disk , in addition to an expanded sample of [α/Fe] and [Fe/H] measurements in M31 satellite galaxies  for individual stars; Wojno et al. 2020 for coadded groups of spectra) and the outer halo . Some of our key results include (1) evidence for a high efficiency of star formation in the GSS progenitor Escala et al. 2020) and the outer disk ), (2) the distinct chemical abundance patterns of the inner halo compared to M31 satellite galaxies Kirby et al. 2020), (3) support for chemical differences between the inner and outer halo Gilbert et al. 2020), and (4) chemical similarity between MW and M31 satellite galaxies Wojno et al. 2020). In this contribution, we analyze the global chemical abundance properties of the kinematically smooth component of M31's inner stellar halo.
This paper is organized as follows. In § 2, we present our recently observed M31 fields. We provide a brief overview of our chemical abundance analysis in § 3 and develop a statistical model to determine M31 membership in § 4. We present our [α/Fe] and [Fe/H] measurements for 128 M31 RGB stars and analyze the combined sample of inner halo abundance measurements (including those in the literature) in § 5. Finally, we compare M31 to the MW and the Magellanic Clouds, and place our results in the context of stellar halo formation models in § 6.  Table 1 presents previously unpublished deep ( 5 hr) observations of six spectroscopic fields in M31. Fields f109_1, f123_1, f130_1, a0_1, a3_1, and a3_2 were observed in total for 7.0, 6.25, 6.74, 6.79, 6.44, and 6.60 hours, respectively. For five of these fields, we utilized the Keck/DEIMOS (Faber et al. 2003) 600 line mm −1 (600ZD) grating with the GG455 order blocking filter, a central wavelength of 7200 Å, and 0.8" slitwidths. We observed a single field (f123_1) with the 1200 line mm −1 (1200G) grating with the OG550 order blocking filter, a central wavelength of 7800 Å, and 0.8" slitwidths. We observed each spectroscopic field using two separate slitmasks that are identical in design, excepting slit position angles. This difference minimizes flux losses owing to differential atmospheric refraction at blue wavelengths via tracking changes in parallatic angle. The spectral resolution of the 600ZD (1200G) grating is approximately 2.8 (1.3) Å FWHM, or R∼3000 (6500) at the Ca II triplet region (λ ∼ 8500 Å). Similarly deep observations of DEIMOS fields in M31, which we further analyze in this work, were previously published by Escala et al. (2020Escala et al. ( , 2019 (600ZD) and Gilbert et al. (2019) (1200G).
One-dimensional spectra were extracted from the raw, two-dimensional DEIMOS data using the spec2d pipeline (Cooper et al. 2012;Newman et al. 2013), including modifications for bright, unresolved stellar sources (Simon, & Geha 2007). Kirby et al. (2020) provides a comprehensive description of the data reduction process. We included additional alterations to the data reduction pipeline to correct for the effects of differential atmospheric refraction, which preferentially affects bluer optical wavelengths ). This correction is equally applied to both 1200G and 600ZD spectra, although it is most significant for 600ZD spectra with λ 5000 Å.

Field Properties
Fields f109_1, f123_1, f130_1, a0_1, and a3 are located at 9, 18, 23, 31, and 33 kpc away from the galactic center of M31 in projected distance. We assumed that M31's galactocenter is located at a right ascension of 0.71 hours and a declination of 41.3 degrees. Table 1 provides the mask center and mask position angle of each field on the sky. Figure 1 illustrates the locations of our newly observed DEIMOS fields relative to the center of M31, including previously observed fields Gilbert et al. 2019) further analyzed in this work, overlaid on a RGB star count map from the Pan-Andromeda Archaeological Survey (PAndAS; McConnachie et al. 2018). We also show the locations of pencil-beam HST/ACS fields (Brown et al. 2009), some of which overlap with DEIMOS fields , from which Brown et al. (2006Brown et al. ( , 2007Brown et al. ( , 2008 derived stellar age distributions. The deep fields presented in this work were previously observed using shallow (∼1 hr) DEIMOS spectroscopy with the 1200G grating to obtain kinematical information for each field (Gilbert et al. 2007(Gilbert et al. , 2009a. Based x 4') and orientation of the DEIMOS fields, whereas the dashed magenta line corresponds to 50 projected kpc. We include fields observed with both the low-resolution (600ZD; magenta) and medium-resolution (1200G; green) gratings on DEIMOS. Nearby HST/ACS fields (202"x 202"; Brown et al. 2009) are shown as gold stars. Observations first presented in this work are outlined in white (Table 1; § 2). Our additional deep fields include 9 kpc, 18 kpc, 23 kpc, and 31 kpc in the stellar halo and 33 kpc in the Giant Stellar Stream.
on the shallow 1200G observations, Gilbert et al. (2007Gilbert et al. ( , 2009a argued that fields f109_1, f123_1, f130_1, and a0_1 probe the properties of the stellar halo of M31, whereas a3 probes the GSS. However, multiple kinematical components can be present in a given field. For example, f123_1 contains substructure known as the Southeast shelf, which is likely associated with the GSS progenitor (Gilbert et al. 2007;Fardal et al. 2007;Escala et al. 2020). Unlike fields along the GSS located closer to M31's center, a3 does not show evidence for the secondary kinematically cold component (KCC) of unknown origin (Kalirai et al. 2006a;Gilbert et al. 2009aGilbert et al. , 2019. We expect that f130_1 is associated with the "smooth", relatively metal-poor halo of M31, based on the known properties of the overlapping DEIMOS field f130_2 ( Figure 1; Escala et al. 2019). Similarly, the velocity distributions of a0_1 and f109_1 are fully consistent with that of M31's stellar halo (Gilbert et al. 2007). Despite being the innermost M31 field in our sample (r proj = 9 kpc, or r disk = 38 kpc assuming i = 77 • ),  Escala et al. (2019Escala et al. ( , 2020 and Gilbert et al. (2019). a Slitmasks labeled as "a", "b", etc., are identical, except that the slits are tilted according to the parallactic angle at the approximate time of observation for the slitmask.
f109_1 does not show any kinematical evidence for significant contamination ( 10%) by M31's extended disk (Gilbert et al. 2007). This also applies to the Southeast shelf, which was predicted to extend to the location of f109_1 with a large velocity dispersion (∼350 km s −1 ; Fardal et al. 2007;Gilbert et al. 2007), such that any hypothetical SE shelf stars in this field could not be kinematically separated from the halo population.

Photometry
The photometry for the majority of fields published in this work were obtained from MegaCam images in the g , i filters using the 3.6 m Canada-France-Hawaii Telescope (CFHT). The MegaCam images were obtained by Kalirai et al. (2006a) and reduced with the CFHT MegaPipe pipeline (Gwyn 2008). For 9 targets in field f109_1 absent from our primary catalogs, we sourced g and i band photometry from the Pan-Andromeda Archaeological Survey (PAndAS) point source catalog (McConnachie et al. 2018). For field f130_1, the g , i magnitudes were transformed to Johnson-Cousins V, I using observations of Landolt photometric standard stars (Kalirai et al. 2006a). In the case of fields a0_1 and a3, the original photometry was obtained in the Washington M, T 2 filters by Ostheimer (2003) using the Mosaic camera on the 4 m Kitt Peak National Observatory (KPNO) telescope and subsequently transformed to the Johnson-Cousins V, I bands using the relations of Majewski et al. (2000).  (Table 1). For reference, we overplot PARSEC isochrones (black lines; Marigo et al. 2017), assuming 9 Gyr ages and [α/Fe] = 0, with [Fe/H] = −2.2, −1.0, −0.5, and 0 (from left to right). Stars without successful radial velocity measurements ( § 3) are represented as black points. Stars that are likely to be MW foreground dwarfs are indicated as grey diamonds, whereas likely M31 RGB stars ( § 4.1) are represented as blue circles. Figure 2 presents the extinction-corrected colormagnitude diagrams (CMDs) for each field in the relevant photometric filters used to derive quantities based on the photometry such as the photometric effective temperature (T eff,phot ), surface gravity (log g), and metallicity ([Fe/H] phot ). We show all stars in a given field for which we extracted 1D spectra (M31 RGB stars, MW foreground dwarf stars, and stars for which we were unable to evaluate M31 membership owing to failed radial velocity measurements; § 4).
For fields f109_1, f130_1, a0_1, and a3, for which stellar spectra were obtained using the 600ZD grating ( § 2.1), we calculated T eff,phot , log g, and [Fe/H] phot following the procedure described by Escala et al. (2020). In summary, the color and magnitude of a star are compared to a grid of theoretical stellar isochrones to derive the above quantities. We utilized the PARSEC isochrones (Marigo et al. 2017), which include molecular TiO in their stellar evolutionary modeling, and assumed a distance modulus to M31 of m − M = 24.63 ± 0.20 (Clementini et al. 2011). For the 600ZD fields, we assumed 9 Gyr isochrones based on the mean stellar ages of the stellar halo and GSS, as inferred from HST CMDs (9.7, 11.0, 10.5 Gyr for the 11, 21, and 35 kpc ACS fields, and 8.8 Gyr for the ACS stream field; Brown et al. 2006Brown et al. , 2007Brown et al. , 2008. The photometric quantities (T eff,phot , log g) for the single 1200G-based field, f123_1, were derived following the procedure outlined by Kirby et al. (2008), assuming an identical distance modulus and 14 Gyr isochrones from a combination of model sets (Girardi et al. 2002;Demarque et al. 2004;VandenBerg et al. 2006). As described in detail by Escala et al. (2019) and summarized in § 5, T eff,phot and log g are used as constraints in measuring T eff , [Fe/H], and [α/Fe] from spectra of individual stars, where these measurements are insensitive to the employed isochrone models and assumed stellar age.

CHEMICAL ABUNDANCE ANALYSIS
We use spectral synthesis of low-and mediumresolution stellar spectroscopy to measure stellar parameters (T eff ) and abundances ([Fe/H] and [α/Fe]) from our deep ( 5 hr) observations of M31 RGB stars. For a detailed description of the low-and medium-resolution spectral synthesis methods, see Escala et al. (2019Escala et al. ( , 2020 and Kirby et al. (2008), respectively. The lowand medium-resolution spectral synthesis procedures are nearly identical in principle, excepting differences in the continuum normalization given the differing wavelength coverage between the low-and medium-resolution spectra (∼4500−9100 vs. 6300−9100 Å). For 1200G spectra, the continuum is determined using "continuum regions" defined by Kirby et al. (2008), whereas such regions would be unilaterally defined for 600ZD spectra owing to the high density of absorption features toward the blue wavelengths. Chemical abundances for individual stars measured using each technique are generally consistent within the uncertainties .
Prior to the chemical abundance analysis, the radial velocity of each star is measured via cross-correlation with empirical templates (Simon, & Geha 2007;Kirby et al. 2015;Escala et al. 2020) observed in the relevant science configuration ( § 2.1). Systematic radial velocity errors of 5.6 km s −1 (Collins et al. 2011) and 1.49 km s −1 (Kirby et al. 2015) from repeat measurements of identical stars are added in quadrature to the random component of the error for observations taken with the 600ZD and 1200G gratings, respectively. The spectral resolution is empirically determined as a function of wavelength using the width of sky lines (Kirby et al. 2008), and in the case of the bluer 600ZD spectra, arc lines from calibration lamps McKinnon et al., in preparation). The observed spectrum is then corrected for telluric absorption using a template of a hot star observed in the relevant science configuration (Kirby et al. 2008;Escala et al. 2019), shifted into the rest frame based on the measured radial velocity, and an initial continuum normalization is performed.
We measured the spectroscopic effective temperature, T eff , informed by photometric constraints, and fixed the surface gravity, log g, to the photometric value ( § 2.3). We simultaneously measured [Fe/H] and [α/Fe] from regions of the spectrum sensitive to Fe and αelements (Mg, Si, Ca -with the addition of Ti for medium-resolution spectra), respectively, by comparing to a grid of synthetic spectra degraded to the resolution of the applicable DEIMOS grating (600ZD or 1200G) using Levenberg-Marquardt χ 2 minimization. The grids of synthetic spectra utilized were generated for 4100−6300 Å and 6300−9100 Å, respectively, by Escala et al. (2019) and Kirby et al. (2008)  ).

MEMBERSHIP DETERMINATION
Separating M31 RGB stars from the intervening foreground of MW dwarf stars has served as one of the primary challenges for spectroscopic studies of individuals stars in M31. The colors and heliocentric radial velocity distributions of MW and M31 stars exhibit significant overlap (e.g., Gilbert et al. 2006), thus the difficulty in disentangling the two populations when the distances to such faint stars are unknown. Early spectroscopic studies of M31 employed simple radial velocity cuts to exclude MW stars (e.g., Reitzel, & Guhathakurta 2002;Ibata et al. 2005;Chapman et al. 2006), resulting in kinematically biased populations of M31 RGB stars and relatively uncontaminated, albeit incomplete samples of M31 stars. Gilbert et al. (2006) performed the first rigorous, probabilistic membership determination in M31 using various diagnostics, including (1) heliocentric radial velocity, (2) the strength of the surface-gravity sensitive Na I λλ8190 doublet, (3) CMD position, and (4) the discrepancy between photometric and calcium triplet based metallicity estimates as a distance indicator.
Given the variety of photometric filters utilized, Gilbert et al.'s method cannot be uniformly applied to all spectroscopic fields analyzed in this work. For inner halo fields (r proj 30 kpc), Escala et al. (2020) illustrated that a binary membership determination using the aforementioned diagnostics is sufficient to recover the majority of stars classified as likely M31 RGB stars by the more sophisticated method of Gilbert et al. (2006) with minimal MW contamination. However, this binary determination excludes all stars with v helio > −150 km s −1 , where some of these stars may be M31 members at the positive tail of the stellar halo velocity distribution. It also does not allow us to assign a degree of certainty to our membership determination for each star. Thus, we used Bayesian inference to assign a membership probability to each observed star with a successful radial velocity determination ( § 3).

Membership Probability Model
We evaluated the probability of M31 membership for all stars with successful velocity measurements based on (1) a parameterization of color (X CMD ), (2) the strength of the Na I absorption line doublet at λλ8190 (EW Na ), (3) heliocentric radial velocity ( § 3; v helio ), and (4) an estimate of the spectroscopic metallicity based on the strength of the calcium triplet ([Fe/H] CaT ).
Thus, according to Bayes' theorem, the posterior probability that a star labeled by an index j is an M31 member-given independent measurements x j = (EW Na,j , X CMD,j , v helio,j , [Fe/H] CaT,j ) with uncertainties δ x j -is proportional to, where P (M31) is the prior probability that a star observed on a given DEIMOS slitmask is a member of M31, and P ( x j |M31) is the likelihood of measuring a given set of membership diagnostics, x j , assuming the star is a M31 member. Analogously, we can also construct P (MW| x j ), the posterior probability that a star belongs to the MW foreground population given a set of diagnostic measurements, where the sum of P (MW| x j ) and P (M31| x j ) is unity.

Measuring Membership Diagnostics
We described the CMD position of each star by the parameter X CMD , which is analogous to photometric metallicity ([Fe/H] phot ). The advantage of using X CMD rather than color and magntiude is that we can place all of our stars on the same scale, despite the diversity of the photometric filter sets ( § 2.3). Assuming 12.6 Gyr isochrones and a distance modulus of 24.47, we defined X CMD = 0 as the color of the most metalpoor PARSEC (Marigo et al. 2017) isochrone ([Fe/H] = −2.2) in the relevant photometric filter and X CMD = 1 as the most metal-rich PARSEC isochrone ([Fe/H] = +0.5). Then, we used linear interpolation to map the color and magnitude of a star to a value of X CMD . The uncertainty on X CMD was derived from the photometric errors by using a Monte Carlo procedure. This normalization provides the advantage of easily identifying stars that are bluer than the most metal-poor isochrone at a fixed stellar age by negative values of X CMD , where stars with X CMD < 0 are 10 times more likely to belong to the MW than M31 ). In the sample being evaluated for membership, we classified stars with (X CMD + δX CMD ) < 0 as MW dwarf stars.
By including [Fe/H] CaT in combination with X CMD as an independent diagnostic in our model, we can additionally distinguish between foreground stars at unknown distances and distant giant stars. Given that our assumed distance modulus is appropriate only for stars at the distance of M31, [Fe/H] phot , and analogously, X CMD , measurements are fundamentally incorrect for MW stars. Therefore, the two stellar populations will appear distinct in X CMD versus [Fe/H] CaT space. In order to compute [Fe/H] CaT in the sample being evaluated for membership, we fit Gaussian profiles to 15 Å wide windows centered on Ca II absorption lines at 8498, 8542, and 8662 Å. Then, we calculated a total equivalent width for the calcium triplet from a linear combination of the individual equivalent widths, (2) following Rutledge et al. (1997a). In addition to ΣCa, the calibration to determine [Fe/H] CaT depends on stellar luminosity, where V −V HB is the Johnson-Cousins V -band apparent magnitude above the horizontal branch, assuming that V HB = 25.17 for M31 (Holland et al. 1996) and that a star is within the RGB ((V −V HB ) > −5). For fields with CFHT MegaCam g and i band photometry ( § 2.3), we approximated V − V HB by g − g HB using an empirical transformation between SDSS and Johnson-Cousins photometry (Jordi et al. 2006), Assuming that g − r = 0.6 and g HB ≈ g HB , we obtained g HB = V HB + (g − V ) = 25.73 for stars in M31. We measured EW Na in the sample being evaluated for membership by summing the area under the bestfit Gaussian line profiles fit to the observed spectrum between 8178−8190, 8189−8200 Å with central wavelengths of 8183, 8195 Å using least-squares minimization. This parameter is sensitive to temperature and surface gravity (Schiavon et al. 1997), thus functioning as a discriminant between MW dwarf stars and M31 RGB stars. The distribution of EW Na for M31 (the MW) in our model ( § 4.1.3) has EW Na = 0.52 (3.07) Å, with values spanning −10.55-9.38 (−3.67-8.88) Å and typical errors of δEW Na = 1.06 (0.58) Å.

Prior Probability of M31 Membership
The probability of a given star observed on a DEIMOS slitmask belonging to M31, P (M31), increases with decreasing projected radial distance from the center of M31, owing to the increase in M31's stellar surface density. The trend of increasing probability of M31 membership with decreasing projected radius is augmented compared to expectations from M31's surface brightness profile (Courteau et al. 2011;Gilbert et al. 2012) by our photometric pre-selection of DEIMOS spectroscopic Figure 3. Properties of spectroscopically confirmed M31 (grey circles, purple contours) and MW (black crosses, blue contours) member stars from the SPLASH (Guhathakurta et al. 2005;Gilbert et al. 2006) survey. The contour levels correspond to 16 th , 50 th , and 84 th percentiles of the distributions. Histograms of star counts are also shown for both M31 and MW members. For clarity, we omit showing uncertainties in each parameter. The middle right panel illustrates the number of spectroscopically confirmed M31 members  relative to the number of targets with successful radial velocity measurements in a given field as a function of projected radius, as observed by the SPLASH survey, where we have fit the relationship with an exponential distribution, PM31(r) ( § 4.1.2). We used heliocentric radial velocity (v helio ), a parameterization of color (XCMD), Na I λλ8190 equivalent width (EWNa), and calcium triplet metallicity estimates ([Fe/H]CaT) as diagnostics to determine membership ( § 4.1.1, § 4.1.3) for our sample of stars with successful radial velocity measurements.
targets. These selection criteria are designed to include M31 members and exclude MW foreground stars. The photometry of the targets spans the magnitude range characteristic of M31 RGB stars (20 < I 0 < 22.5). We also considered narrow-band, Washington DDO51 photometry in the fields a0_1 and a3. The DDO51 filter isolates the Mg b triplet, which acts as discriminant between MW dwarf and M31 giant stars due to its sensitivity to surface gravity. We expect that an enhancement in the probability of observing an M31 RGB star owing to magnitude cuts is fairly uniform at a factor of 2 within r proj 30 kpc (Gilbert et al. 2012), where the majority of our targets are located. For fields toward the outer halo, such as a0_1 and a3, DDO51-selection bias increases the likelihood of observing an M31 RGB star by a factor of 3−4 (Gilbert et al. 2012).
Thus, we parameterized P(M31) empirically using the ratio of the number of secure M31 RGB stars , N M31 , to the number of targets with successful radial velocity measurements, N targets , from the Spectroscopic and Phomtometric Landscape of Andromeda's Stellar Halo (SPLASH; Guhathakurta et al. 2005;Gilbert et al. 2006) survey. The SPLASH survey consists of tens of thousands of shallow (∼1 hr exposures) DEIMOS spectra of stars along the line of sight toward M31's halo, disk, and satellite dwarf galaxies (e.g., Kalirai et al. 2010;Tollerud et al. 2012;Dorman et al. 2012Dorman et al. , 2015Gilbert et al. 2012Gilbert et al. , 2014Gilbert et al. , 2018. In addition to N M31 /N targets as a function of radius, Figure 3 presents measurements of EW Na , X CMD , v helio , and [Fe/H] CaT for 1,510 secure M31 members and 1,794 secure MW members across 29 spectroscopic fields in M31's stellar halo. We controlled for differences in methodology between this work and SPLASH by re-determining X CMD homogeneously for the SPLASH data from the original Johnson-Cousins photometry and measuring [Fe/H] CaT for our sample based on the same   (Table 1) with successful radial velocity measurements ( § 3). We also include previously observed fields H, S, D, and f130_2 ( Figure 1; Escala et al. 2020). For clarity, we omit showing measurement uncertainties. Each point is color-coded by its probability of belonging to M31 (Eq. 8; § 4.1.3). The SPLASH distributions utilized in our probability model ( § 4.1) are underlaid as contours. For field D, we used a distinct velocity criterion owing to the presence of M31's disk at ∼ −130 km s −1 . We identified 328 (210) secure M31 RGB stars and 92 (75) secure MW dwarf stars across the spectroscopic fields presented in this (prior) work.
calibration (Rutledge et al. 1997a) used in SPLASH ( § 4.1.1). 1 Based on this data, we approximated P(M31) by an exponential distribution, where r j is the projected radius of a star from M31's galactic center and r p = 45.5 kpc. Figure 3 includes P M31 (r) overlaid on the SPLASH survey data.

Likelihood of M31 Membership
1 Our equivalent width measurement procedure differs from that utilized in the SPLASH survey. As opposed to summing the flux decrement in a window centered on a given absorption feature, we performed Gaussian fits. However, we do not expect this difference in metholodgy to significantly affect the usability of the SPLASH data to construct the likelihood ( § 4.1.3), given that our measured EW Na and [Fe/H] CaT distributions are consistent with SPLASH ( Figure 4).
Given the wealth of existing information on the properties of M31 RGB stars (and the MW foreground dwarf stars characteristic of our selection function) from the SPLASH survey, we assigned membership likelihoods to individual stars informed by this extensive data set. We described the likelihood that a star, j, with unknown membership belongs to M31 given its diagnostic measurements and uncertainties, ( x j , δ x j ), as, where ( θ i , δ θ i ) is a set of four diagnostic measurements and uncertainties for a star, i, from the SPLASH survey that is a secure M31 member. The total number of secure SPLASH member stars, N i , equals 1,510 (1,794) for M31 (the MW). The likelihood that a star, j, with unknown membership belongs to the MW, P ( x j , δ x j |MW), is defined analogously. Assuming normally distributed uncertainties, the log likelihood that a star, j, is a member of either M31 or the MW given a single set of SPLASH measurements, i, is a non-parametric, N kdimensional Gaussian distribution, where k corresponds to a given membership diagnostic (EW Na , X CMD , v helio , or [Fe/H] CaT ) and N k is the number of diagnostics utilized for a given star. For some stars in our sample, we were unable to measure EW Na and/or [Fe/H] CaT as a consequence of factors such as weak absorption, low S/N, or convergence failure in the Gaussian fit. For such stars, we excluded EW Na and/or [Fe/H] CaT as a diagnostic, such that N k = 2−3. In four dimensions, the separation between the M31 red giants and MW foreground dwarfs ( Figure 3) is analogous to ∼2.8σ (∼7.4σ) in units of the covariance matrix of M31's (the MW's) distribution. Finally, we computed the probability that a star is a M31 RGB candidate, as opposed to a MW dwarf candidate, from the odds ratio of the posterior probabilities (Eq. 1), where the proportionality factor in Eq. 1 is equivalent for P (M31| x j , δ x j ) and P (MW| x j , δ x j ).

Results of Membership Determination
Figure 4 summarizes our membership determination for 426 total stars with successful radial velocity measurements across the six spectroscopic fields first presented in this work. The probability distribution is strongly bimodal, where the majority of stars are either secure (p M31,j 0.75) M31 RGB (N M31 = 328) or MW dwarf (N MW = 92) candidates, excepting 6 stars with intermediate properties (0.5 < p M31,j 0.75). 2 Figure 4 also includes a homogeneously re-evaluated membership determination for fields H, S, D, and f130_2 ( Figure 1; Escala et al. 2020), which we further analyze in this work. Across these four fields, 346 stars have successful radial velocity measurements, 210 (75) of which are classified as secure M31 (MW) stars. In total, 61 stars in these fields have intermediate properties, most of which (60) are located in field D. Owing to the presence of M31's northeastern disk at MW-like line-of-sight velocity (v helio ∼ −130 km s −1 ; Escala et al. 2020), we calculated M31 membership probabilities in D without the use of radial velocity as a diagnostic. Then, we classified all stars with (v helio − δv helio ) > −100 km s −1 as MW contaminants, as in Escala et al. (2020).
In order to maximize our sample size across all spectroscopic fields, we considered stars that are more likely to belong to M31 than the MW (p M31,j > 0.5) to be M31 members in the following analysis. For stars in common between our dataset and SPLASH, we recovered 91.1% of stars classified as M31 RGB stars by Gilbert et al. (2006), where we used an equivalent definition of membership ( L splash,j > 0). The excess MW contamination is 0.30%, in addition to the expected 2-5% from Gilbert et al.'s method. The 8.6% discrepancy results from stars at MW-like heliocentric velocities that we conservatively classified as MW stars, where these stars are considered M31 RGB stars in SPLASH. Figure 5 illustrates the heliocentric radial velocity distributions for all stars with successful measurements ( § 5), including both M31 RGB stars and MW dwarf stars ( § 4.1), across the spectroscopic fields. We also show the adopted Gaussian mixture models ) describing the velocity distribution for each field, which were computed using over 5,000 spectroscopically confirmed M31 RGB stars across 50 fields in M31's stellar halo. Gilbert et al. omitted radial velocity as a membership diagnostic ( § 4) in their analysis and simultaneously fit for contributions from M31 and MW components to obtain kinematically unbiased models for each field. Table 2 presents the parameters characterizing the velocity model for each field, where we assumed 50 th percentile values of Gilbert et al.'s marginalized posterior probability distribution functions. For the stellar halo components, we transformed the mean velocity from the Galactocentric to heliocentric frame using the median right ascension and declination of all stars in a given field .

Kinematics of M31 RGB Stars
We confirmed that the observed velocity distribution for each field is consistent with its velocity model using a two-sided Kolmogorov-Smirnov test. As discussed in § 2.2, f109_1, f130_1, and a0_1 probe the "smooth" stellar halo of M31 with no detected substructure, whereas fields f123_1 and a3 show clear evidence The Southeast shelf, a tidal feature potentially originating from the GSS progenitor Gilbert et al. 2007), is evident in f123_1, whereas the substructure in fields a3_1 and a3_2 correspond to the GSS.
of substructure known to be associated with the Southeast shelf Gilbert et al. 2007Gilbert et al. , 2019Escala et al. 2020) and the GSS. Owing to the spatial proximity ( Figure 1) and kinematical similarity (Gilbert et al. 2007(Gilbert et al. , 2009a between fields f130_1 (this work) and f130_2  and fields a3 (this work), we consider them together in our subsequent abundance analysis ( § 5).

Substructure Probability
Based on the velocity model for each field, we assigned a probability of belonging to substructure in M31's stellar halo, p sub , to every star identified as a likely RGB candidate (p M31 > 0.5; Eq. 8). For smooth stellar halo fields, p sub = 0. For fields with substructure, we computed p sub using an equation analogous to Eq. 8, where the Bayesian odds ratio is substituted with the relative likelihood that a M31 RGB star belongs to substructure versus the halo, where v j is the heliocentric velocity of an individual star and µ, σ, and f are the mean, standard deviation, and fractional contribution of the substructure or halo component in the Gaussian mixture describing the velocity distribution for a given field. In contrast, Escala et al. (2020) calculated p sub based on the full posterior distributions, as opposed to 50 th percentiles alone, of their velocity models for fields with kinematical substructure. Given that we included abundance measurements of M31 RGB stars from Escala et al. (2020) in this work, we re-calculated p sub for such fields (H, S, and D; Figure 1) based on their 50 th percentiles, as above. In the following abundance analysis ( § 5), we incorporated p sub as a weight when determining the chemical properties of the stellar halo. cluding previously published inner halo fields with abundance measurements Gilbert et al. 2019). We calculated each average chemical property from a bootstrap resampling of 10 4 draws of the final sample for each field, including weighting by the inverse square of the measurement uncertainty and the probability that a star belongs to a given kinematical component ( § 4.2.1).
Hereafter, we predominantly restricted our abundance analysis to RGB stars ( § 4.1) from all six fields that are likely to be dynamically associated with the kinematically hot stellar halo of M31. Additionally, we incorporated likely M31 halo stars from inner halo fields with existing abundance measurements Gilbert et al. 2019) into our final sample. We will analyze the chemical composition of RGB stars likely belonging to substructure, such as the Southeast shelf (f123_1) and outer GSS (a3), in a companion study (I. Escala et al., in preparation). A catalog of stellar parameters and abundance measurements for M31 RGB stars in the six spectroscopic fields is presented in Appendix A.
This section is organized as follows. § 5.1 defines our final spectroscopic sample and discusses potential observational biases, and § 5.2 presents a brief overview of our results pertaining to each newly published spectroscopic field. In § 5.3, we measure spectral synthesis based radial abundance gradients for M31's inner stellar halo and compare to metallicity gradient predictions from photometry. Lastly, § 5.4 compares our enlarged sample of abundance measurements in M31's stellar halo to its dSphs.

Sample Selection and Potential Biases
We vetted our final sample to consist only of reliable abundance measurements ( § 3) for M31 RGB stars. Similar to Escala et al. (2019Escala et al. ( , 2020, we restrict our analysis to M31 RGB stars with δ(T eff ) < 200 K, δ([Fe/H]) < 0.5, δ([α/Fe]) < 0.5, and well-constrained parameter estimates based on the 5σ χ 2 contours for all fitted parameters (T eff , [Fe/H], and [α/Fe]). We also require that convergence is achieved in each of the measured parameters ( § 3). Unreliable abundance measurements often result from an insufficient signal-to-noise (S/N) ratio, translating to an effective S/N threshold of 8 Å −1 for robust measurements of [Fe/H] and [α/Fe]. Such S/N limitations result in a bias in our final sample against metal-poor stars with low S/N spectra, but do not affect the [α/Fe] distributions. This bias is negligible for our innermost halo fields (f109_1, f123_1), whereas it is on the order of 0.10−0.15 dex for our remaining fields (f130_1, a0_1, a3).
We also manually inspected spectra to exclude stars with clear signatures of strong molecular TiO absorption in the wavelength range 7055−7245 Å from our final sample. It is unclear whether abundances measured from TiO stars are accurate owing to the lack of (1) the inclusion of TiO in our synthetic spectral modeling (Kirby et al. 2008;Escala et al. 2019) and (2)  appropriate calibration sample. For stars with successful radial velocity measurements, 31.5%, 16.7%, 34.5%, 32.5%, 30.4%, and 31.3% of stars in fields f109_1, f123_1, f130_1, a0_1, a3_1, and a3_2 have clear evidence of TiO in their spectra. To be conservative, we excluded 54 M31 RGB stars that showed TiO but otherwise would have made the final sample. In total, 128 M31 RGB stars across the six spectroscopic fields pass the above selection criteria, thereby constituting our final sample. 3 Figure 7 presents a CMD of M31 RGB stars in each field, color-coded by the probability that a given star belongs to kinematically identified substructure in the stellar halo ( § 4.2.1). We have also highlighted our final sample of stars with reliable abundance measurements and indicated which stars were excluded as a consequence of TiO absorption strength.
As discussed in detail by Escala et al. (2020), removing stars on the basis of TiO results in a bias against stars with red colors, which translates to a bias in photometric metallicity (

Abundance Distributions for Individual Fields
In this section, we provide a brief overview of the chemical abundance properties of the six spectroscopic fields first presented in this work (Table 1; Figure 1). We discuss the global properties of the inner stellar halo in § 5.3 and § 5.4.
• f109_1 (9 kpc halo field): Field f109_1 is the innermost region of the stellar halo of M31 yet including weighting by the inverse variance of the measurement uncertainty and the probability that a star belongs to a given kinematical component. a The parameters of the velocity model are the 50 th percentiles of the marginalized posterior probability distribution functions. These were computed by Gilbert et al. (2018) for all fields except H, S, and D. We have transformed the 50 th percentile values for the stellar halo components from the Galactocentric to heliocentric frame, based on the median right ascension and declination of all stars in a given field. For fields H, S, and D, Escala et al. (2020) fixed the stellar halo component to the parameters derived by Gilbert et al. (2018) to independently compute the posterior distributions. b Chemical abundances for fields f109_1, f123_1, f130_1, a0_1, and a3 are first presented in this work. We have included chemical properties of previously published M31 fields H, S, D, and f130_2  and f207_1  for reference. We further analyze the halo populations of H, S, f130_2, and f207_1 in this work. c We combined the chemical abundance samples for fields f130_1 (this work) and f130_2  given their proximity ( Figure 1) and the consistency of their velocity distributions (Gilbert et al. 2007. The same is true for fields a3_1 and a3_2 (this work), where Gilbert et al. (2009aGilbert et al. ( , 2018 illustrated the similarity in kinematics between these fields. probed with chemical abundances (Gilbert et al. 2007. It does not contain any detected kinematical substructure ( § 2.2, 4.2). Excepting fields along the GSS, which have an insufficient number of smooth halo stars to constrain the abundances of the stellar halo component, f109_1 is more metal-rich on average ( [Fe/H] = −0.93 +0.08 −0.09 ) than the majority of halo components in fields at larger projected distance (Table 2; see also § 5.3). Based on the abundances for this field, the inner halo of M31 may be potentially less α-enhanced on average ( [α/Fe] = +0.32 ± 0.08) than the stellar halo at larger projected distances, although we did not measure a statistically significant [α/Fe] gradient in M31's inner stellar halo ( § 5.3). Additionally, we did not find evidence for a correlation between [Fe/H] and [α/Fe] for this field, where we computed a distribution of correlation coefficients from 10 5 draws of the measured abundances perturbed by their (Gaussian) uncertainties.
• f123_1 (18 kpc halo field): Field f123_1 is dominated by the smooth stellar halo, but it also has a clear detection of substructure ( § 4.2; Table 2) known as the Southeast shelf ( § 2.2; Fardal et al. 2007;Gilbert et al. 2007). We defer further analysis of this component to a companion paper (I. Escala et al., in preparation)  • f130_1 (23 kpc halo field): Similar to f109_1, f130_1 does not possess detectable substructure. We combined the chemical abundance samples for this field with f130_2  due to their proximity ( Figure 1) and the consistency of their velocity distributions (Gilbert et al. 2007 Table 2) agree within the uncertainties with previous determinations from an 11 star sample by Escala et al. (2019Escala et al. ( , 2020. The lack of a significant trend between [α/Fe] and [Fe/H] is also maintained by the larger sample. The inner halo as probed by field f130 at 23 kpc appears to be more metal-poor and α-enhanced than field f109_1 at 9 kpc, possibly representing a population that assembled rapidly at early times ). The star formation history inferred for this region of M31's halo indicates that the majority of the stellar population is over 8 Gyr old (Brown et al. 2007), suggesting that an accretion origin would require the progenitor galaxies to have quenched their star formation at least 8 Gyr ago.
• a0_1 (31 kpc halo field): The stellar population at the location of a0_1 is solely associated with the kinematically hot component of the stellar halo ( § 4.2). The 10 RGB star sample in this field suggests a positive trend between [α/Fe] and [Fe/H] (r = 0.24 +0.12 −0.24 ) that is likely a consequence of small sample size. Despite its similar projected distance from the center of M31 and kinematical profile, a0_1 appears to be more metal-rich ( [Fe/H] = −1.35 ± 0.10) than f130, although it may be comparably α-enhanced (Table 2). Nearby HST/ACS fields at 35 kpc ( Figure 1) imply a mean stellar age of 10.5 Gyr (Brown et al. 2008) for the vicinity, compared to a mean stellar age of 11.0 Gyr at 21 kpc (near f130). The full age distributions sug-gest that the stellar populations are in fact distinct between the 21 and 35 kpc ACS fields: the star formation history of the 35 kpc field is weighted toward more dominant old stellar populations and is inconsistent with the star formation history at 21 kpc at more than 3σ significance (Brown et al. 2008). If applicable to a0_1, this suggests that its stellar population is both younger and more metal-rich than that at 23 kpc in field f130.
• a3 (33 kpc GSS fields): Field a3 is dominated by GSS substructure ( § 2.2, § 4.2), such that it may not provide meaningful constraints on the smooth stellar halo in this region. However, the stellar halo at 33 kpc along the GSS in field a3 may be more metal-poor ( [Fe/H] = −1.48 +0.12 −0.13 ) than at fields f207_1 and S, which are located at 17 and 22 kpc along the GSS ( Table 2). The abundances in this field clearly show a declining pattern of [α/Fe] with [Fe/H], which is characteristic of dwarf galaxies. A comparison between [Fe/H] for the GSS between fields f207_1 and S indicates that the GSS likely has a metallicity gradient, as found from photometric based metallicity estimates (Ibata et al. 2007;Gilbert et al. 2009a;Conn et al. 2016;Cohen et al. 2018). We will quantify spectral synthesis based abundance gradients in the GSS and make connections to the properties of the progenitor in future work (I. Escala et al., in preparation).
In agreement with previous findings Gilbert et al. 2019), the M31 fields are αenhanced with a significant spread in metallicity. This implies that stars in the inner stellar halo formed rapidly, such that the timescale for star formation was less than the typical delay time for Type Ia supernovae. For an accreted halo, a spread in metallicity coupled with high α-enhancement can indicate contributions from multiple progenitor galaxies or a dominant, massive progenitor galaxy with high star formation efficiency (Robertson et al. 2005;Font et al. 2006b;Johnston et al. 2008;Font et al. 2008). We further discuss formation scenarios for M31's stellar halo in § 6.3.

The Inner vs. Outer Halo
The existence of a steep, global metallicity gradient in M31's stellar halo is well-established from both Ca triplet based (Koch et al. 2008)  plored. From a sample of 70 M31 RGB stars, including an additional 21 RGB stars from Gilbert et al. (2019), with spectral synthesis based abundance measurements, Escala et al. (2020) found tentative evidence that the inner stellar halo (r proj 26 kpc) had higher [α/Fe] than a sample of four outer halo stars (Vargas et al. 2014b) drawn from ∼70−140 kpc. Gilbert et al. (2020) increased the number of stars in the outer halo (43−165 kpc) with abundance measurements from four to nine. In combination with existing literature measurements by Vargas et al. (2014b), Gilbert et (Vargas et al. 2014a;Kirby et al. 2020) and the MW halo (e.g., Ishigaki et al. 2012;Hayes et al. 2018). In comparison, [α/Fe] is higher at fixed [Fe/H] in M31's inner halo than in its dwarf galaxies   (Vargas et al. 2014b;Escala et al. 2019Escala et al. , 2020Gilbert et al. 2019Gilbert et al. , 2020, we re-determined the average values for the inner halo, obtaining [Fe/H] = −1.08 ± 0.04 (−1.17 ± 0.04) and [α/Fe] = 0.40 ± 0.03 (0.39 ± 0.04) for M31's inner halo when excluding (including) substructure. From this enlarged sample of inner halo abundance measurements, we assessed the presence of spectral synthesis based radial gradients in [Fe/H] and [α/Fe] over the radial range spanned by our spectroscopic sample (8-34 kpc).
Taking stars in the kinematically hot stellar halo to have p sub < 0.5 (Eq. 9), we identified 122 (75) stars that are likely associated with the "smooth" stellar halo (kinematically cold substructure) within a projected distance of r proj 35 kpc of M31. We emphasize that these numbers are simply to provide an idea of the relative contribution of each component to the stellar halo-in the subsequent analysis, we do not employ any cuts on p sub , but rather incorporate M31 RGB stars, regardless of halo component association, by using p sub as a weight. We have excluded all M31 RGB stars in field D (Figure 1) from our analysis sample owing to the presence of M31's northeastern disk. To determine the radial abundance gradients of the smooth inner halo, we fit a line using an MCMC ensemble sampler (Foreman-Mackey et al. 2013), weighting each star by p halo = 1 − p sub and the measurement uncertainty. We also determined the radial gradients in the case of including substructure by removing p sub as a weight.
In addition to showing our measured radial gradients, we also include the photometric metallicity gradient of Gilbert et al. (2014) in Figure 8. Using a sample of over 1500 spectroscopically confirmed M31 RGB stars, Gilbert et al. (2014) measured a radial [Fe/H] phot gradient of −0.011 ± 0.001 dex kpc −1 between 10−90 projected kpc. The radial photometric metallicity gradients measured by Gilbert et al. with and without kinematical substructure were found to be consistent within the uncertainties. In contrast to Gilbert et al. (2014), our results suggest that the radial [Fe/H] gradient of the smooth halo is inconsistent with that of the halo including substructure over the probed radial range. The substructures in our inner halo fields are likely GSS progenitor debris Gilbert et al. 2007Gilbert et al. , 2009a, thus the change in slope at its inclusion may reflect a convolution with the distinct metallicity gradient of the GSS progenitor (I. Escala et al., in preparation).

Spectroscopic vs. Photometric Metallicity Gradients
In order to control for differences in sample size, target selection, number and locations of spectroscopic fields utilized, and the radial extent of the measured gradi- A possible explanation for the difference in trends with substructure between spectral synthesis and CMD-based gradients is the necessary assumption of uniform stellar age and α-enhancement to determine [Fe/H] phot . Although we did not measure a significant radial [α/Fe] gradient (Figure 8), M31's inner stellar halo has a range of stellar ages present at a given location based on HST CMDs extending down to the main-sequence turn-off (Brown et al. 2006(Brown et al. , 2007(Brown et al. , 2008. If the smooth stellar halo is systematically older than the tidal debris toward 30 kpc, this would steepen the relative [Fe/H] phot gradient between populations with and without substructure. This agrees with our observation that [Fe/H] phot −[Fe/H] is increasingly positive on average toward larger projected distances, where the discrepancy is greater for the smooth halo than in the case of including substructure. Additionally, a comparatively young stellar population at 10 kpc compared to 30 kpc could result in a steeper [Fe/H] phot gradient in better agreement with our measured [Fe/H] gradient. The difference between mean stellar age at 35 kpc and 10 kpc is 0.8 Gyr (Brown et al. 2008), though the mean stellar age of M31's halo does not appear to increase monotonically with projected radius. Assuming constant [α/Fe], this mean age difference translates to a negligible gradient between 10−35 kpc. Thus, a more likely explanation for the discrepancy in slope between the [Fe/H] phot and [Fe/H] gradients are uncertainties in the stellar isochrone models at the tip of the RGB.
In general, [Fe/H] phot is offset toward higher metallicity, where the adopted metallicity measurement methodology can result in substantial discrepancies for a given sample (e.g., Lianou et al. 2011). For example, the discrepancy between the CMD-based gradient of Gilbert et al. (2014) and literature measurements from spectral synthesis (Vargas et al. 2014b;Gilbert et al. 2019;Escala et al. 2020) decreases when assuming an α-enhancement ([α/Fe] = +0.3) in better agreement with the inner and outer halo . Despite differences in the magnitude of the slope, the behavior with substructure, and intercept of the [Fe/H] gradient from different metallicity measurement techniques, our enlarged sample of spectral synthesis based [Fe/H] measurements provides further support for the existence of a large-scale metallicity gradient in the inner stellar halo, where this gradient extends out to at least 100 projected kpc in the outer halo . For a thorough consideration of the implications of steep, large-scale negative radial metallicity gradients in M31's stellar halo, we refer the reader to the discussions of Gilbert et al. (2014) and Escala et al. (2020).

Metallicity Gradients
Alternatively, the discrepancy between the [Fe/H] phot and [Fe/H] gradients could be partially driven by the bias against red, presumably metal-rich, stars incurred by the omission stars with strong TiO absorption from our final sample ( § 5.1). We investigated the impact of potential bias from both (1)  Incorporating these bias terms yields a radial [Fe/H] gradient of −0.022 ± 0.002 dex kpc −1 (−0.018 ± 0.001 dex kpc −1 ) and an intercept of −0.60 ± 0.03 (−0.58 ± 0.03) without (with) substructure. The primary effect of including bias estimates is a shift toward higher metallicity in the overall normalization of the gradient. The gradient slopes calculated including bias estimates are consistent with our previous measurements, which did not account for potential sources of bias. We can therefore conclude that the slopes of radial [Fe/H] gradients between 8-34 kpc in M31's stellar halo are robust against these two possible sources of bias.

Comparing the Stellar Halo to M31 dSphs
In ΛCDM cosmology, simulations of stellar halo formation for M31-like galaxies predict that the chemical abundance distributions of an accreted component should be distinct from the present-day satellite population of the host galaxy (Robertson et al. 2005;Font et al. 2006b;Tissera et al. 2012), as is observed for the MW (Unavane et al. 1996;Shetrone et al. 2001Shetrone et al. , 2003Tolstoy et al. 2003;Venn et al. 2004). This chemical distinction is driven by the early assembly of the stellar halo, where its progenitors were accreted ∼8-9 Gyr ago, as opposed to ∼4-5 Gyr ago for surviving satellite galaxies (Bullock & Johnston 2005;Font et al. 2006b;Fattahi et al. 2020). Accordingly, Escala et al. (2020)  . This could be a consequence of an insufficient sample size at low metallicity, where they identified 29 RGB stars likely belonging to the smooth stellar halo (p halo > 0.5; Eq. 9). Another possibility is that such low-metallicity stars, with lower average αenhancement ( [α/Fe] ∼ 0.25; Escala et al. 2020), may represent a stellar population in the stellar halo more similar to present-day M31 dSphs.

1-D Comparisons
To investigate whether the [α/Fe] distributions of the stellar halo and dSphs at low-metallicity are in fact statistically indistinguishable, we repeated the analysis of Escala et al. (2020) with our expanded sample of 197 inner halo RGB stars (excluding field D; Figure 1) with abundance measurements.
In contrast to Escala et al. (2020), we applied corrections to the abundances measured by Vargas et al. (2014a)  [α/Fe] values. This correction is reflected for the relevant M31 dSphs in Figures 9 and 10.
We constructed [α/Fe] distributions for a mock stellar halo built by destroyed dSph-like progenitor galaxies, weighting the contribution of each M31 dwarf galaxy by (Eq. 4 of Escala et al. 2020), 5 where M ,j is the stellar mass and L V,j is the V-band luminosity of the host dSph, Φ is the V-band luminosity function of present-day M31 dSphs (McConnachie 2012), N j is the number of RGB stars with abundance measurements in galaxy j, and N gal is the total number of M31 dSphs contributing to the abundance distribu-  (Vargas et al. 2014a) and And X . Evidently, the majority of inner halo RGB stars are inconsistent with the stellar populations of M31 dSphs, as previously found using the combined 91 RGB star sample of Escala et al. (2020) and Gilbert et al. (2019), and reinforced by this work. In general, [α/Fe] tends to be higher for M31's inner halo at fixed metallicity than for the dSphs. Many of the high-metallicity, α-enhanced stars in the inner halo have a high probability of being associated with kinematical substructure, such as the GSS.
In addition to these trends, Figure 9 suggests the existence of a stellar population preferentially associated with the dynamically hot halo that also has chemical abundance patterns similar to M31 dSphs. These stars are not contained within a well-defined metallicity bin (such as [Fe/H] < −1.5, the metal-poor bin utilized by Escala et al. 2020), but rather span a region of [α/Fe] versus [Fe/H] space coincident with the mean trend of the dSphs.

Modeling 2-D Chemical Abundance Distributions
In order to robustly identify M31 halo stars with abundance patterns similar to that of M31 dSphs (Figure 9), we modeled the observed 2-D chemical abundance ratio distributions, as advocated by Lee et al. (2015). Appendix B presents a proof of concept of the statistical method utilized in the following analysis, which is based on the chemical abundance distributions of the MW's stellar halo. Appendix B illustrates that differences in populations can be meaningfully distinguished given our typical uncertainties.
We considered the smooth component of the inner stellar halo, substructure in the inner halo, and M31 dSphs as distinct populations in our modeling. We expanded upon the M31 dSph sample analyzed in § 5.4.1 by incorporating abundance measurements from galaxies that were previously excluded on the basis of their small sample sizes: And X (N stars = 9 ; Kirby et al. 2020) and NGC 147 (N stars = 7 ;Vargas et al. 2014a). 6 Taking δ([Fe/H]) < 0.5, δ([α/Fe]) < 0.5, and [Fe/H] < −0.5, these two additional galaxies result in a sample size of 293 abundance measurements for M31 dSphs.
Assuming the form of a bivariate normal distribution, the likelihood of observing a given set of abundance measurements (x i , y i ) = ([Fe/H] i , [α/Fe] i ) and uncertainties (δx i , δy i ) = (δ[Fe/H] i , δ[α/Fe] i ) 7 given a model described by means µ = (µ x , µ y ) and standard deviations σ = (σ x , σ y ) is, 6 We did not include M32 in our M31 dSph abundance sample. Vargas et al. (2014a) measured abundances for 3 stars in M32 with [Fe/H] < −0.5 and concluded that they were not representative of the galaxy. 7 We assumed independent errors (δx i , δy i ) in our model. In actuality, we expect that some amount of covariance, δxy, exists between our measurement uncertainties. The net result of such dependent uncertainties could be a perceived change in the correlation coefficient, r (Eq. 11). However, we found that accounting for δxy does not significantly alter our error distribution, so we anticipate that the effect of this covariance on r is minimal.
where i is an index corresponding to an individual RGB star. We incorporated the measurement uncertainties into Eq. 11 via the variable σ xi = σ 2 x + δ 2 xi , which is defined analogously for σ yi . The correlation coefficient, r, is an additional model parameter that accounts for covariance between x and y.
For the full stellar halo sample, we modeled the smooth and substructure components simultaneously by combining the respective likelihood functions (Eq. 11) using a mixture model, where p sub i is the kinematically-based probability that a star belongs to substructure in M31's stellar halo (Eq. 9). Thus, the full stellar halo abundance ratio distribution is represented by an eight parameter model ( µ halo , µ sub , σ halo , σ sub ), whereas we utilized a four parameter model ( µ dSph , σ dSph ) for M31 dSphs. The likelihood of the entire observed data set for a given stel-  Figure 10) for a given stellar population, and the parameters describing the bivariate normal probability density functions (Eq. 11). We fit for the smooth halo and substructure components simultaneously (Eq. 12) using a mixture model. The parameters for the 2-D chemical abundance models are the 50 th percentiles of the marginalized posterior probability distribution functions ( § 5.4.2), where the uncertainies on each parameter were calculated from the corresponding 16 th and 84 th percentiles.
lar population (i.e., M31's stellar halo or the dSphs) is therefore the product of the individual likelihoods, Using Bayes' theorem (see also Eq. 1), we evaluated the posterior probability of a particular bivariate model accurately describing a given a set of abundance measurements for a stellar population, where P ( N i=1 x i , y i , δx i , δy i | µ, σ) is the likelihood (Eq. 13) and P ( µ, σ) represents our prior knowledge regarding constraints on µ and σ. We implemented noninformative priors over the allowed range for each fitted parameter, where we assumed uniform priors and inverse Gamma priors for non-dispersion and dispersion parameters, respectively. We allowed the dispersion parameters in our model to vary between 0.0 < σ < 1.0, the correlation coefficients between −1 < r < 1, and permitted samples of µ to be drawn from −3.0 < [Fe/H] < 0.0 and −0.8 < [α/Fe] < 1.2.
In the case of the variance, σ 2 , the inverse Gamma distribution is described as Γ −1 (α, β), where α is a shape parameter and β is a scale parameter. This distribution is defined for σ > 0 and can account for asymmetry in the posterior distributions for dispersion parameters, which are restricted to be positive. The bootstrap resampled standard deviations ( § 5. With the above formulation, we sampled from the posterior distribution of each model (Eq. 14) using an affineinvariant Markov Chain Monte Carlo (MCMC) ensemble sampler (Foreman-Mackey et al. 2013). That is, we solved for the parameters describing the two-component stellar halo abundance distribution (Eq. 12), and separately for a model corresponding to M31 dSphs. We evaluated each model using 100 walkers and ran the sampler for 10 4 steps, retaining the latter 50% of the MCMC chains to form the converged posterior distribution. Table 3 presents the final parameters of each 2-D chemical abundance ratio distribution model. We computed the mean values from the 50 th percentiles of the marginalized posterior distributions, whereas the uncertainty on each parameter was calculated relative to the 16 th and 84 th percentiles. Figure 10 displays the resulting bivariate normal distribution models for the 2-D chemical abundance ratio distributions, reflecting the general trends of the observed abundance distributions ( § 5.4 .3). However, the models predict similar mean α-enhancements between the two halo populations.
Based on these models, we calculated the probability that a star is both kinematically associated with the smooth component of the stellar halo and has abundance patterns similar to M31 dSphs, where L i is the likelihood from Eq. 11 and p sub,f is the refined probability that a star belongs to kinematical substructure (see Eq. 9). For every sampling of the posterior distribution for each star, we calculated the odds ratio of the Bayes factor, K, for the substructure versus halo models, where K = p sub L sub i /( (1 − p sub ) L halo i ). Then, we computed p sub,f for each star from the 50 th percentile of these distributions. The outcome of this procedure is an updated determination of the substructure probability for each star that incorporates both kinematical and chemical information. For stars with nonzero p sub , the net effect is p sub,f > p sub for metal-rich stars ([Fe/H] −1.0) and p sub,f < p sub for metal-poor stars. This change occurs because the substructure model (Table 3) has a higher mean [Fe/H] than the smooth halo model, such that metal-rich stars are more likely to belong to substructure. Figure 10 illustrates the selection of smooth halo stars with M31 dSph-like chemical abundances, where we have identified 22 RGB stars that securely fall within this category (p dSph > 0.75). These stars are at least three times more likely to have abundances consistent with M31 dSphs than the bulk of the smooth halo population. Hereafter, we refer to this group of stars as a dSph-like population. We refrain from characterizing [Fe/H] or [α/Fe] for these dSph-like stars given that our sample is likely incomplete. On the high metallicity end, we cannot attribute stars to the dSph-like population because the M31 dSph PDF utilized in its detection cuts off at [Fe/H] −1.0 (Figure 9) owing to the mass range and metallicity distribution functions of the dSphs. In part due to this incompleteness, we emphasize that the identification of a subset of M31 RGB stars as dSph-like is not akin to a measurement of M31's accreted halo fraction. This is further made the case by (1) our limited sample size, (2) the lack of available higherdimensional phase space information, and (3) the fact that we have defined the dSph-like sample by de facto excluding stars with a high probability of belonging to kinematically identifiable halo substructure.
By definition, dSph-like stars belong to the dynamically hot component of the stellar halo in which we do not detect any kinematical substructure. These stars span the same range of parameter space in line-of-sight velocity (−700 km s −1 < v helio −150 km s −1 ) and projected distance (8 kpc < r proj < 34 kpc) as the full sample of RGB stars in M31's inner halo. Unsurprisingly, dSph-like stars are more likely to be found in spectroscopic fields along the minor axis dominated by the smooth halo, such as f130 and a0_1, as opposed to fields dominated by tidal debris along the high surface-brightness core of the GSS. We investigate possible origins for the dSph-like stars in § 6.2 and § 6.3.

The Effect of Potential Sources of Bias on Population Detections
Potential bias from the omission of red, presumably metal-rich, stars with strong TiO absorption ( § 5.1) does not affect the robust identification of smooth halo stars with chemical abundance patterns similar to M31 dSphs. Introducing maximal bias estimates for [Fe/H] based on this source alone ( § 5.3.2) shifts the mean metallicity of the halo and substructure models (Table 3) to −1.05 ± 0.05 and −0.64 ± 0.07, respectively. TiO stars are not a significant source of bias in M31 dSphs , and their abundance measurements are similarly affected by S/N limitations ( § 5.1) compared to M31 RGB stars. Thus, the net effect of the omission of TiO stars would be an increased separation between the dSph and halo PDFs. This would result in the classification of 36 secure dSph-like stars, thereby reinforcing the detection of this population.  Hayes et al. (2018) illustrated that their data set is equivalent to the two distinct populations of Nissen & Schuster (2010, and by extension contains Gaia-Enceladus stars (Haywood et al. 2018). Gaia-Enceladus was likely comparable to the Small Magellanic Cloud (SMC; M ∼ 10 8.5 M ) in stellar mass at infall, with M ∼10 8−9 M and was accreted ∼10 Gyr ago (e.g., Helmi et al. 2018;Belokurov et al. 2018;Gallart et al. 2019;Mackereth et al. 2019;Fattahi et al. 2019). 8 In contrast to earlier work (using DR12) on metalpoor MW field stars (Hawkins et al. 2015) that relied on kinematical selection, Hayes et al. (2018) used a larger sample of stars in combination with a data-driven approach to identify stellar populations that were distinct in multi-dimensional chemical abundance space. This  In order to perform a more direct comparison between M31's and the MW's chemical abundance distributions, we perturbed the MW halo abundances by uncertainties from 10 4 draws of the empirical error distribution of our measurements. 9 During each draw, each MW halo star was perturbed by random values sourced from a normal distribution defined by the uncertainties on a single randomly selected M31 RGB star and the MW halo star. 10 The M31 RGB star is randomly selected from a metallicity bin within 0.2 dex of the MW halo star to preserve the metallicity-dependence of the M31 We show an example of a single perturbation of the MW halo stars in the upper right panel of Figure 11. The primary effect of this perturbation is that the MW halo stars now span a similar range in [α/Fe] as the M31 halo stars. If the MW was observed at the distance of M31 (i.e., if the MW abundances had an M31-like error distribution), we would measure  Figure 12). In this figure, we also qualitatively compared M31 RGB stars that likely belong to M31's dSph-like population (p dSph > 0.75; Eq. 15) to the MW PDFs. Given that these M31 stars chemically resemble dSphs, it is interesting to compare to the accreted, LMg population in the MW's stellar halo. A subset of M31's dSph-like stars appear similar to the LMg population, but it is unclear if all of M31's dSph-like stars could originate from a LMg-like population. In particular, M31's halo may have stars that are more metal-poor ([Fe/H] −1.8) than those observed in Hayes et al.'s MW halo sample. 11 We note that the APOGEE abundances attributed to Gaia-Enceladus from Helmi et al. (2018) extend down to [Fe/H] ∼ −2.5 and up to [α/Fe] ∼ 0.5, potentially encapsulating the low-metallicity RGB stars in M31's dSph-like population. We explore origin scenarios for M31's dSph-like stars in § 6.2 and § 6.3.

Potential Sources of Systematic Offsets
We refrained from drawing detailed conclusions based on these comparisons owing to potential systematic offsets in the abundance distributions between the MW and M31 data sets. As previously discussed, the APOGEE chemical abundances are based on highresolution spectroscopy in the near-infrared H-band, in 11 We do not expect any systematic biases in APOGEE elemental abundances to result in the dearth of data at [Fe/H] −1.8 and above the ASPCAP grid edge at [Fe/H] = −2.5. APOGEE is biased against metal-rich stars in fields that are either distant (i.e., dominated by cool giants) or enriched to supersolar metallicity (Hayden et al. 2014(Hayden et al. , 2015, but there is no quantified bias for metal-poor stars in APOGEE with [Fe/H] −2.5. contrast with our measurements, which we derive from low-(R ∼ 3000) and medium-(R ∼ 6000) resolution spectroscopy at optical wavelengths with comparatively low S/N. Additionally, the APOGEE abundances are internally calibrated against observations of metal-rich open clusters  and externally calibrated to correct for systematic offsets. Although our spectral synthesis method relies on MW globular clusters to determine the systematic uncertainty on our abundance measurements (Kirby et al. 2008;, 2020), we do not calibrate the stellar parameters and elemental abundances outputted by our abundance pipelines to external observations. The selection functions also differ between Hayes et al. (2018) and this work, where we have avoided selection criteria based on chemical abundances. Although we acknowledge these caveats, an exploration of the impact of these possible systematic offsets and selection effects between the MW and M31 data sets is beyond the scope of this paper.   Kirby et al. 2020). Because the chemical abundance distributions of M31's dSphs are consistent with that of the MW's dSphs at fixed stellar mass , these dSph-like stars in M31 are also similar to the MW's dSphs. This implies that these dSph-like stars may have been accreted onto M31's stellar halo from galaxies of similarly low stellar mass and star formation history (LM-SFH) as the dSphs, or of higher stellar mass and similar star formation efficiency (HM-SFE). Although the identified dSph-like stars are statistically likely to have been accreted, they do not represent a complete sample of such stars ( § 5.4.2), and are therefore not to be confused with the hypothetical true accreted fraction of M31's stellar halo ( § 6.3.1) Given that the dSph-like stars could in large part originate from massive dwarf galaxies ( § 6.3.2), we compared the abundances of M31's dSph-like stars to the Magellanic Clouds (MCs), which have both high stellar masses (2.7 × 10 9 M and 3.1 × 10 8 M for the LMC and SMC, respectively; van der Marel et al. 2002;Stanimirović et al. 2004) and pronounced metal-poor populations. We used the MC abundance sample of Nidever et al. (2019) from APOGEE DR16 (Ahumada et al. 2019) obtained as part of SDSS-IV (Blanton et al. 2017) through the installation of a second APOGEE spectrograph in the Southern hemisphere (Wilson et al. 2019). 12 Figure 13 shows [α/Fe]  Chemical evolution models for the LMC (and by extension, the SMC, owing to their similar chemical abundance patterns) predict the dominance of Type Ia supernovae over core-collapse supernovae during this quiescent epoch (Bekki & Tsujimoto 2012;Nidever et al. 2019), resulting in the observed decline in [α/Fe] with [Fe/H] in the metal-poor tails. These tails are comparably ancient to the dSph-like stars: the predicted agemetallicity relation from Nidever et al. (2019) has the LMC reaching [Fe/H]∼ −1.5 at 6 Gyr into its evolution (Nidever et al. 2019), suggesting that the majority of its metal-poor tail was in place by ∼7 Gyr ago. 13 Based on their chemical abundances, Nidever et al. (2019) concluded that the early evolution of the MCs was characterized by low star formation efficiency (gas mass converted into stellar mass). This low efficiency may result from the MCs having evolved in isolation, given that they are likely on first infall into the MW's potential well (Besla et al. 2007(Besla et al. , 2012. The relatively quiescent star formation history of the MCs until ∼4 13 The models of Bekki & Tsujimoto (2012) reach [Fe/H] ∼ −1.5 by ∼ 9−11 Gyr ago, but were not well-constrained at low metallicity by the available data for the LMC at the time of their study.
Gyr ago may also be a consequence of such isolation (e.g., Smecker-Hane et al. 2002;Harris & Zaritsky 2009;Weisz et al. 2013). Thus, we can infer that the galaxie(s) that were accreted onto M31's halo to form its dSph-like population, which is characterized by low star formation efficiency, experienced chemical evolution similar to the early evolution of the isolated MCs. We evaluate whether the scenario of multiple low-mass progenitors (LM-SFH) or a dominant high-mass progenitor (HM-SFE) is more likely to explain the origin of M31's dSph-like stars in § 6.3.2.

Formation Scenarios for M31's Inner Stellar Halo
In this subsection, we explore the implications of our abundance measurements in M31's stellar halo for insitu and accreted formation channels ( § 6.3.1). We also investigate possible origins of M31's dSph-like population in the context of the stellar mass function of accreted galaxies ( § 6.3.2).

The In-Situ vs. Accreted Halo
Simulations predict the existence of an in-situ component of the stellar halo, which formed in the main progenitor of the host galaxy, in addition to an accreted component formed from external galaxies (Zolotov et al. 2009;Font et al. 2011;Tissera et al. 2013;Cooper et al. 2015). In addition to a component formed via dissipative collapse, the in-situ halo includes contributions from kinematically heated stars originating in the disk (Purcell et al. 2010;McCarthy et al. 2012;Tissera et al. 2013;Cooper et al. 2015). The relative contributions of these two formation channels to stellar halo build-up depend on the stochastic accretion history of the host galaxy, where more active histories correspond to lower in-situ stellar halo mass fractions (Zolotov et al. 2009). However, the in-situ fraction also depends on the numerical details of a given simulation (Zolotov et al. 2009;Cooper et al. 2015). The more recent simulations of Cooper et al. (2015) have found that the in-situ fraction typically ranges between ∼30-40%. Despite variations in the predictions of this fraction, multiple studies have found that the in-situ component tends to dominate within the inner ∼20-30 kpc of the stellar halo (Zolotov et al. 2009;Font et al. 2011;Tissera et al. 2013;Cooper et al. 2015).
Predicted chemical signatures for the in-situ and accreted stellar halo can also vary depending on the simulation. In some studies (Zolotov et al. 2009;Font et al. 2011;Tissera et al. 2013Tissera et al. , 2014, in-situ star formation produces more metal-rich stellar halos than accretion alone. In contrast, Cooper et al. (2015) argued that the in-situ and accreted halo may be indistinguishable based on metallicity alone, assuming that the insitu halo forms mostly from gas stripped from accreted galaxies. For galaxies that have not experienced a recent (z 1) major merger, Zolotov et al. (2010) found that the in-situ halo has higher [α/Fe] at fixed [Fe/H] than the accreted halo for stars with [Fe/H] > −1 in the solar neighborhood. Tissera et al. (2012Tissera et al. ( , 2013 found the opposite trend, where accreted stars tend to be more αenhanced (but also more metal-poor) than in-situ stars.
On the observational front, a bimodality in the metallicity of the MW's stellar halo with radius was interpreted as evidence in favor of a two-component halo defined by in-situ and accreted populations (Carollo et al. 2007(Carollo et al. , 2010. After Gaia DR2, it has become apparent that the MW's presumably in-situ inner stellar halo is in actuality the remnants of Gaia-Enceladus Haywood et al. 2018;Belokurov et al. 2018;Deason et al. 2018). A veritable in-situ component of the inner stellar halo has recently emerged via its distinct chemical composition and kinematics in relation to Gaia-Enceladus tidal debris (Hayes et al. 2018;Haywood et al. 2018;Di Matteo et al. 2019;Gallart et al. 2019;Conroy et al. 2019, with earlier suggestions from Gaia DR1 by Bonaca et al. 2017). In particular, Belokurov et al. (2020) identified this component from old, higheccentricity, α-rich stars on retrograde orbits 14 where with [Fe/H] > −0.7 in the solar neighborhood, associating it with formation from the MW's proto-disk following its last significant merger.
At this time, the presence of a significant in-situ component in the inner stellar halo of M31 is less clear than in the case of the MW. Along the major axis of M31's northeastern disk, there is compelling kinematical evidence for both a rotating inner spheroid (Dorman et al. 2012) and dynamically heated stars originating from the disk (Dorman et al. 2013). An extended disk-like structure possibly formed from M31's last significant merger has also been detected within the inner ∼40 kpc of M31's disk plane . Although these stellar structures may reasonably be various insitu components, they have not been detected within the radial range spanned by our data (∼8-34 kpc in M31's southeastern quadrant, along its minor axis; Figure 1). Even at the distance of our innermost spectroscopic field (r proj ∼ 9 kpc, or r disk ∼ 38 kpc assuming i = 77 • ), the contribution from the extended disk is expected to be 10% (Guhathakurta et al. 2005;Gilbert et al. 2007). Interestingly, M31's global stellar halo properties appear to agree with predictions from accretion-only mod-els for stellar halo formation (Harmsen et al. 2017), along with other MW-like, edge-on galaxies in the Local Volume from the GHOSTS survey (Radburn-Smith et al. 2011;Monachesi et al. 2016a). Indeed, multiple lines of evidence indicate that accretion has played a dominant role in the formation of M31's stellar halo ( § 6.3.2).
Based on this alone, both the dSph-like and non-dSphlike stars in M31's smooth, inner stellar halo ( § 5.4.2) could be accreted populations resulting from chemically distinct groups of progenitors. For example, many of the non-dSph-like stars have high substructure probabilities, and are therefore likely associated with GSS tidal debris. The remainder of the non-dSph-like stars, which are associated with the smooth component of the stellar halo, could plausibly originate from separate massive progentior(s) with high star formation efficiency in a pure accretion scenario. As for the dSph-like smooth halo stars, the similarity of their chemical abundance patterns to M31's dSphs provides strong evidence in favor of an accretion origin in either the LM-SFH or HM-SFE scenarios ( § 6.2). However, we acknowledge the possibility that the few low-metallicity ([Fe/H] −2) dSph-like stars with high α-enhancement ([α/Fe] 0.3) could have formed in-situ. In this regime, dSph and stellar halo chemical abundance patterns are generally similar (e.g., Venn et al. 2004), such that [α/Fe] is degenerate with star formation history (Lee et al. 2015).
Conversely to the pure accretion origin scenario, the overlap between M31's broader stellar halo population and the MW's in-situ HMg population in abundance space ( § 6.1; Figure 11) suggests that M31 could reasonably possess a significant, ancient in-situ population, instead of having a predominately accretion origin in distinct progenitor(s). In an in-situ formation scenario, the bulk of M31's metal-rich stars in the smooth halo may have originated from a proto-stellar disk, similar to the MW's in-situ halo population (Belokurov et al. 2020). Given that we are lacking multi-dimensional kinematical information or higher-dimensional chemical information, it is difficult to distinguish between accreted and in-situ formation scenarios for the bulk of M31's stellar population without performing a detailed comparison to simulations, which is beyond the scope of this paper.

The Mass Function of Destroyed Galaxies
The general consensus from simulations of MWlike galaxies is that massive dwarf galaxies (M ∼ 10 8−10 M ) distinct from present-day satellites are the dominant progenitors of the accreted stellar halo (Bullock & Johnston 2005;Robertson et al. 2005;Font et al. 2006b;Cooper et al. 2010;Deason et al. 2016;D'Souza & Bell 2018a;Fattahi et al. 2020). These galaxies have a median accretion time of 9 Gyr ago (Bullock & Johnston 2005;Font et al. 2006b;Fattahi et al. 2020). Furthermore, a few massive progenitors are predicted to form the bulk of the accreted halo (Cooper et al. 2010;Deason et al. 2016;D'Souza & Bell 2018a;, where the secondary progenitor is on average half as massive as the primary progenitor (D'Souza & Bell 2018a). In particular, Fattahi et al. (2020) found that such massive progenitors dominate the stellar mass budget within the inner 50 kpc of the stellar halo, where contributions from progenitors with M < 10 8 M become non-negligible around ∼100 kpc. The disruption of globular clusters associated with the accretion of lowmass dwarf galaxies and/or ancient systems that formed in-situ likely plays a subdominant role in stellar halo formation, with contributions to the stellar mass budget of ∼2-5% based on recent observational (Koch et al. 2019;Naidu et al. 2020) and theoretical (Reina-Campos et al. 2020) studies.
The details of the stellar mass function of destroyed constituent galaxies are dictated by the accretion history of the host galaxy. Deason et al. (2016) found that MWlike galaxies with more active accretion histories have higher mass accreted progenitors (M ∼ 10 9−10 M ) and larger accreted halo mass fractions than their quiescent counterparts. The wealth of substructure visible in M31's stellar halo (e.g., Ferguson et al. 2002;Ibata et al. 2007;McConnachie et al. 2018), including the GSS (Ibata et al. 2001), and the lack of an apparent break in its density profile (Guhathakurta et al. 2005;Irwin et al. 2005;Courteau et al. 2011;Gilbert et al. 2012) suggests that M31 has likely experienced multiple and continuous contributions to its stellar halo from accreted galaxies (Deason et al. 2013;Cooper et al. 2013;Font et al. 2020). Moreover, the total accreted stellar mass of M31 (1.5 ± 0.5 × 10 10 M ; Harmsen et al. 2017) is on the upper end of both observed (Merritt et al. 2016;Harmsen et al. 2017) and predicted Deason et al. 2016;Monachesi et al. 2019) accreted stellar halo mass fractions for MW-like galaxies.
The large-scale negative metallicity gradient of M31's stellar halo  this work) may also lend support to it being dominated by a few massive progenitors, although this trend could also result from a significant in-situ component in the inner regions of the halo (Font et al. 2011;Tissera et al. 2014;Monachesi et al. 2019). Thus, M31's accreted mass function is likely weighted toward progenitor galaxies with M ∼ 10 9−10 M . This agrees with the hypothesis that the main contributor to M31's stellar halo (including substructure) is the GSS progenitor, a massive galaxy with M ∼ 10 9−10 M accreted between ∼1-4 Gyr ago, depending on the merger scenario (Fardal et al. , 2008Hammer et al. 2018;D'Souza & Bell 2018b).
The mass spectrum of progenitor galaxies dictates the metallicity of the accreted stellar halo through the stellar mass-metallicity relation for galaxies (Gallazzi et al. 2005;Kirby et al. 2013) and its redshift evolution (Gallazzi et al. 2014;Ma et al. 2016;Leethochawalit et al. 2018Leethochawalit et al. , 2019. In particular, Leethochawalit et al. (2018Leethochawalit et al. ( , 2019 found that the normalization of the stellar-mass metallicity relation evolves by 0.04 ± 0.01 dex Gyr −1 . Assuming that this relation extends to M < 10 9 M , this means that an SMC-mass galaxy (M ∼ 10 8 Observations of the MW-like GHOSTS galaxies (Harmsen et al. 2017) have firmly established such a relationship between stellar halo mass and metallicity, as first suggested by Mouhcine et al. (2005). The physical driving force for this relation is the fact that the most massive progenitors dominate halo assembly (D'Souza & Bell 2018a; Monachesi et al. 2019).
In addition to setting the stellar halo metallicity, the most massive progenitors largely determine the full metallicity distribution function of the stellar halo (Deason et al. 2016;Fattahi et al. 2020). For a typical progenitor mass of 10 9−10 M , these galaxies contribute 90% of metal-poor stars ([Fe/H] < −1) and ∼20-60% of stars with [Fe/H] < −2 (Deason et al. 2016). Deason et al. also found that classical dwarf galaxies (M ∼ 10 5−8 M ) contribute the remainder of metalpoor stars, with negligible contributions from ultra-faint dwarf galaxies.
Indeed, M31's dSph-like stars most likely have an accretion origin in either the LM-SFH or HM-SFE scenarios ( § 6.2) owing to their defining similarity to the abundance patterns M31 dwarf galaxies ( § 5.4.2). Given that (1) the chemical abundance patterns of M31's dSphlike population are broadly consistent with the lowmetallicity ([Fe/H] −1.5), early evolution of the MCs ( § 6.2) and (2) we do not detect any kinematical substructure for this population ( § 5.4), the accreted galaxies(s) that contributed to M31's dSph-like population could have been massive (M ∼ 10 8−9 M ) if they evolved in an isolated, low star formation efficiency environment and were accreted 8 Gyr ago.
Given M31's halo properties and inferred accretion history (based on halo formation models), a possible scenario for the formation of the dSph-like population (−2.5 [Fe/H] −1.0) is the accretion of a massive, secondary progenitor with low star formation efficiency (e.g., the HM-SFE scenario). Assuming that the GSS progenitor is the dominant progenitor, the stel-lar mass of this secondary progenitor should be approximately between the mass of the SMC and LMC (M ∼ 10 8.5−9.5 M ; § 6.2) and possibly comparable to Gaia-Enceladus ( § 6.1). An early, massive accretion event such as this would also deposit its debris closer to the center of the host potential in the inner halo owing to dynamical friction (Cooper et al. 2010;Tissera et al. 2013;Fattahi et al. 2020). Alternatively, M31's dSphlike population could have formed exclusively from the accretion of multiple progenitors similar in stellar mass and star formation history to classical dwarf galaxies (e.g., the LM-SFH scenario). Although this hypothesis cannot be rejected, it is less favored by the predictions of stellar halo formation in a cosmological context (Bullock & Johnston 2005;Font et al. 2006b;Deason et al. 2016;D'Souza & Bell 2018a).
To attempt to discern between the LM-SFH and HM-SFE hypotheses, we consider the accretion of an SMCmass galaxy onto M31 ∼10 Gyr ago. Assuming a significant metallicity dispersion as observed in LG galaxies (∼0.4 dex; e.g., Kirby et al. 2013;Ho et al. 2015;Kirby et al. 2020), such a galaxy would span −1.9 (−2.2) [Fe/H] −0.9 (−0.6) within 1σ (2σ). It is unclear if such a galaxy could account for the most metal-poor stars observed in M31's dSph-like population ([Fe/H] −2.2). Additional progenitors with M ∼ 10 7−8 M may be necessary to explain these low-metallicity stars in an accretion origin scenario.

SUMMARY
We have presented measurements of [α/Fe] and [Fe/H] for 128 individual M31 RGB stars from low-(R ∼ 3000) and medium-(R ∼ 6000) resolution Keck/DEIMOS spectroscopy. With a combined sample of 197 M31 RGB stars with abundance measurements in inner halo fields Gilbert et al. 2019; this work), we have undertaken an analysis of the chemical abundance properties of M31's kinematically smooth stellar halo between 8-34 kpc. Our primary results are the following: 7. We concluded that M31's dSph-like population was most likely accreted ( § 6.3) onto the stellar halo from progenitors either similar in stellar mass (M ∼ 10 6−8 M ) and star formation history to M31 dSphs (the LM-SFH scenario), or progenitors of higher mass (M 10 8−9 M ) and similarly low star formation efficiency (the HM-SFE scenario). The remaining, dominant population of stars in the smooth component of M31's halo may result from progenitor(s) distinct from those of the dSphlike population, and/or represent an in-situ stellar halo component.

ACKNOWLEDGMENTS
We thank the referee for their comments, which improved this paper. The authors would additionally like to thank Miles Cranmer and Erik Tollerud for insightful discussions regarding statistics and David Nidever for providing a catalog of APOGEE data for the Magellanic Clouds. We also thank Stephen Gwyn for reducing the photometry for slitmasks f109_1 and f123_1 and Jason Kalirai for the reduction of f130_1.
I.E. acknowledges support from a National Science Foundation (NSF) Graduate Research Fellowship un-der Grant No. DGE-1745301, in addition to a Carnegie-Princeton Fellowship through the Carnegie Observatories. This material is based upon work supported by the NSF under Grants No. AST-1614081 (E.N.K., I.E.) AST-1614569 (K.M.G, J.W.), and AST-1412648 (P.G.). E.N.K gratefully acknowledges support from a Cottrell Scholar award administered by the Research Corporation for Science Advancement, as well as funding from generous donors to the California Institute of Technology. E.C.C is supported by a Flatiron Research Fellowship at the Flatiron Institute. The Flatiron Institute is supported by the Simons Foundation. The analysis pipeline used to reduce the DEIMOS data was developed at UC Berkeley with support from NSF grant AST-0071048.
We are grateful to the many people who have worked to make the Keck Telescope and its instruments a reality and to operate and maintain the Keck Observatory. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.    (Kirby et al. 2008;Escala et al. 2019Escala et al. , 2020 Note-(This table is available in its entirety in machine-readable form.).

B. IDENTIFYING PROBABLE ACCRETED STARS IN THE STELLAR HALO OF THE MW
In § 5.4.2, we modeled the 2-D chemical abundance distributions of M31's stellar halo and its dSphs in or-der to statistically identify M31 RGB stars with abundances similar to M31 dSphs. The suggestion that such stars may exist in M31's stellar halo was first made apparent by Figure 9 through a comparison to abundance The unperturbed MW halo abundance distribution color-coded according to a star's probability of being SMC-like (defined analogously to Eq. 15), as evaluated from the perturbed abundance distributions. Magenta stars are SMC-like, whereas yellow stars are not. measurements in M31 dSphs obtained using the same techniques. In this Appendix, we use APOGEE chemical abundances of the MW's stellar halo (Hayes et al. 2018) and the ancient, metal-poor tail of the SMC ( § 6.2; Nidever et al. 2019), which we perturb by an M31-like error distribution ( § 6.1), to illustrate the efficacy of this method for successfully recovering dwarf-galaxy-like stars in a stellar halo population. Figure 14 outlines the steps in this proof of concept. First, we perturbed the MW's halo abundance distribution by M31-like errors, as described in § 6.1. Then, we fit a 2D PDF, following § 5.4.2, to the perturbed MW halo abundances without separating it into the LMg or HMg populations of Hayes et al. 2018. The left panel of Figure 14 illustrates the resulting PDF (yellow contours). Next, we isolated the metal-poor tail of the SMC abundance distribution ([Fe/H] −1.4). As discussed in § 6.2, chemical evolution models for the Large Magellanic Cloud (Bekki & Tsujimoto 2012;Nidever et al. 2019) indicate that its metal poor tail was in place by ∼7-9 Gyr ago at the latest. Given the similarly low star formation efficiency of the SMC (Nidever et al. 2019), we can assume that its metal-poor tail is similarly ancient. The more ancient, metal-poor portion of the SMC abundance distribution is the relevant portion to compare to the predominately old stellar halo of the MW (e.g., Gallart et al. 2019;Bonaca et al. 2020). We perturbed the SMC metal-poor tail by M31-like errors and fit a 2D PDF to resulting abundance distribution (left panel of Figure 14, magenta contours).
The middle panel of Figure 14 shows the original unperturbed MW halo and SMC abundance distributions compared to the 2D PDFs calculated from the perturbed distributions. We used the MW halo PDF and the SMC metal-poor (MP) PDF to calculate a probability that each MW halo star has abundances similar to the SMC (analogously to Eq. 15). To calculate this probability, we used the perturbed values of the abundances for a given star. The right panel of Figure 14 shows the unperturbed MW halo abundances color-coded according to this probability calculated using the perturbed equivalent.
We compared the stars classified as "SMC-like" (p SMC > 0.75) by the perturbation analysis to the stars that are truly SMC-like (e.g., the LMg stars). The LMg population corresponds to stars that were accreted during the Gaia-Enceladus merger event ( § 6.1; Helmi et al. 2018;Belokurov et al. 2018;Haywood et al. 2018) of a galaxy at least as massive as the SMC. Assuming that stars with p SMC > 0.75 (0.5) are SMC-like, this results in the correct classification of 99% (98%) of SMC-like stars as belonging to the LMg population. If we assume that stars with p SMC < 0.25 (0.5) are not SMC-like, this results in the correct classification of 38% (32%) of not-SMC-like stars as belonging to the HMg population. This is because the classification of SMC-like stars is incomplete at high metallicity beyond the bounds of the SMC MP PDF. However, the relevant point for our analysis is that the classification of stars as SMC-like is highly accurate. Repeating this exercise for different perturbations of the MW halo and SMC abundance distributions does not alter these percentages by more than 0.5%.  Figure 15. Stars identified as SMC-like (pSMC > 0.75; pink circles) compared to MW halo stars that were not classified as SMC-like separated according their membership in the LMg or HMg population (black and brown points, respectively; Hayes et al. 2018). Some stars classified as SMC-like reside close to the HMg region of the MW halo abundance distribution because they were classified based on their perturbed abundances and uncertainties. 99% of SMC-like stars are correctly classified as originating from the accreted LMg population. (Right panel) 2D PDFs fit to abundance distributions of the perturbed MW halo (solid yellow contours) and SMC MP tail (solid pink contours), compared to 2D PDFs fit to the unperturbed equivalents (dashed gold and purple contours, respectively). The true, underlying abundance distribution is constrained by both the perturbed and unperturbed abundance distributions. Figure 15 further supports this analysis. The left panel shows stars that we have identified as SMC-like overlaid on the full MW halo abundance distribution, where it is evident that the overwhelming majority of SMC-like stars correspond to the LMg population. The right panel shows PDFs fit to the perturbed and unperturbed abundance distributions for the MW halo and SMC MP tail. The primary difference between the perturbed and unperturbed PDFs is that the unperturbed PDFs are narrower in the [α/Fe] direction, whereas the shapes, angles, and metallicity spreads between each case are similar. Both the perturbed and unperturbed distributions represent approximations of the true, underlying parent distribution. The fact that they show significant overlap demonstrates that we can in fact reasonably constrain the underlying PDFs describing the abundance distributions in spite of our uncertainties in [α/Fe].
In this example, the statistical methodology developed in § 5.4.2 is successful at identifying MW stars that are similar in their 2D [α/Fe] vs. [Fe/H] abundances to stars in the metal-poor tail of the SMC abundance distribution, using abundances perturbed by the typical uncertainties in our M31 sample. This gives us confidence that we can statistically identify M31 halo stars that are similar in their 2D [α/Fe] vs. [Fe/H] abundances to stars in M31's dSph galaxies, and thus may have formed in similar environments and been accreted onto M31's halo. This is true despite the fact that we did not assume any knowledge about the true underlying abundance distributions, and despite the size of the uncertainties [α/Fe]. Thus, this statistical method can successfully identify dwarf-galaxy-like stars not only in the MW's halo, but also in the smooth component of M31's stellar halo.