Elemental Abundances in M31: Properties of the Inner Stellar Halo^*

We described the CMD position of each star by the parameter X_CMD, which is analogous to photometric metallicity (${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{phot}}$). The advantage of using X_CMD rather than color and magnitude is that we can place all of our stars on the same scale, despite the diversity of the photometric filter sets (Section 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 metal-poor 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}_{\mathrm{CMD}}\lt 0$ are ≳10 times more likely to belong to the MW than M31 (Gilbert et al. 2006). In the sample being evaluated for membership, we classified stars with (${X}_{\mathrm{CMD}}+\delta$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, ${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{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,

$$\begin{eqnarray}&&{\rm{\Sigma }}\mathrm{Ca}=0.5\times {\mathrm{EW}}_{\lambda 8498}+1.0\times {\mathrm{EW}}_{\lambda 8542}+0.6\times {\mathrm{EW}}_{\lambda 8662},\end{eqnarray} \tag{ 2 }$$

following Rutledge et al. (1997). In addition to ΣCa, the calibration to determine [Fe/H]_CaT depends on stellar luminosity,

$$\begin{eqnarray}&&{\left[\mathrm{Fe}/{\rm{H}}\right]}_{\mathrm{CaT}}=-2.66+0.42\times {\rm{\Sigma }}\mathrm{Ca}+0.27\times (V-{V}_{\mathrm{HB}}),\end{eqnarray} \tag{ 3 }$$

where $V-{V}_{\mathrm{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}_{\mathrm{HB}}$) > −5). For fields with CFHT MegaCam ${g}^{{\prime} }$ and ${i}^{{\prime} }$ band photometry (Section 2.3), we approximated $V-{V}_{\mathrm{HB}}$ by ${g}^{{\prime} }-{g}_{\mathrm{HB}}^{{\prime} }$ using an empirical transformation between SDSS and Johnson–Cousins photometry (Jordi et al. 2006),

$$\begin{eqnarray}&&g-V=0.56\times (g-r+0.23)/1.05+0.12.\end{eqnarray} \tag{ 4 }$$

Assuming that g − r = 0.6 and ${g}_{\mathrm{HB}}^{{\prime} }\approx {g}_{\mathrm{HB}}$, we obtained ${g}_{\mathrm{HB}}={V}_{\mathrm{HB}}+(g-V)=25.73$ for stars in M31. For M31's (the MW's) distribution of ΣCa in our model (Section 4.1.3), we found $\langle {\rm{\Sigma }}\mathrm{Ca}\rangle =5.60(3.00)$ Å, with values spanning −4.70 to 21.66 (−4.07 to 10.37) Å and typical errors of $\langle \delta {\rm{\Sigma }}\mathrm{Ca}\rangle \,=1.71(0.88)$ Å.

We measured EW_Na in the sample being evaluated for membership by summing the area under the best-fit Gaussian line profiles fit to the observed spectrum between 8178–8190 and 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 (Section 4.1.3) has $\langle {\mathrm{EW}}_{\mathrm{Na}}\rangle =0.52(3.07)$ Å, with values spanning −10.55 to 9.38 (−3.67 to 8.88) Å and typical errors of $\langle \delta {\mathrm{EW}}_{\mathrm{Na}}\rangle =1.06(0.58)$ Å.

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 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₀ < 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 (Gilbert et al. 2006), ${N}_{{\rm{M}}31}$, to the number of targets with successful radial velocity measurements, ${N}_{\mathrm{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; Dorman et al. 2012, 2015; Gilbert et al. 2012, 2014, 2018; Tollerud et al. 2012). In addition to ${N}_{{\rm{M}}31}$/${N}_{\mathrm{targets}}$ as a function of radius, Figure 3 presents measurements of EW_Na, X_CMD, v_helio, and [Fe/H]_CaT for 1510 secure M31 members and 1794 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 calibration (Rutledge et al. 1997) used in SPLASH (Section 4.1.1).⁹

Zoom In Zoom Out Reset image size

Based on this data, we approximated P(M31) by an exponential distribution,

$$\begin{eqnarray}&&{P}_{{\rm{M}}31}({r}_{j})=\exp (-{r}_{j}/{r}_{p}),\end{eqnarray} \tag{ 5 }$$

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}_{{\rm{M}}31}(r)$ overlaid on the SPLASH survey data.

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, (${{\boldsymbol{x}}}_{j}$, $\delta {{\boldsymbol{x}}}_{j}$), as

$$\begin{eqnarray}&&P({{\boldsymbol{x}}}_{j},\delta {{\boldsymbol{x}}}_{j}| {\rm{M}}31)=\displaystyle \frac{1}{{N}_{i}}\displaystyle \sum _{i=1}^{{N}_{i}}P({{\boldsymbol{x}}}_{j}| {{\boldsymbol{\theta }}}_{i},\delta {{\boldsymbol{\theta }}}_{i}),\end{eqnarray} \tag{ 6 }$$

where (${{\boldsymbol{\theta }}}_{i}$, $\delta {{\boldsymbol{\theta }}}_{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 1510 (1794) for M31 (the MW). The likelihood that a star, j, with unknown membership belongs to the MW, $P({{\boldsymbol{x}}}_{j},\delta {{\boldsymbol{x}}}_{j}| \mathrm{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 nonparametric, N_k-dimensional Gaussian distribution,

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}P({{\boldsymbol{x}}}_{j},\delta {{\boldsymbol{x}}}_{j}| {{\boldsymbol{\theta }}}_{i},\delta {{\boldsymbol{\theta }}}_{i}) & = & -\displaystyle \frac{{N}_{k}}{2}\mathrm{ln}2\pi -\displaystyle \frac{1}{2}\displaystyle \sum _{k=1}^{{N}_{k}}\mathrm{ln}(\delta {x}_{j,k}^{2}+\delta {\theta }_{i,k}^{2})\\ & & -\,\displaystyle \frac{1}{2}\displaystyle \sum _{k=1}^{{N}_{k}}\displaystyle \frac{{\left({x}_{j,k}-{\theta }_{i,k}\right)}^{2}}{\delta {x}_{j,k}^{2}+\delta {\theta }_{i,k}^{2}},\end{array}\end{eqnarray} \tag{ 7 }$$

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 an M31 RGB candidate, as opposed to an MW dwarf candidate, from the odds ratio of the posterior probabilities (Equation (1)),

$$\begin{eqnarray}&&{p}_{{\rm{M}}31,{\rm{j}}}=\displaystyle \frac{P({\rm{M}}31| {{\boldsymbol{x}}}_{j},\delta {{\boldsymbol{x}}}_{j})/P(\mathrm{MW}| {{\boldsymbol{x}}}_{j},\delta {{\boldsymbol{x}}}_{j})}{1+P({\rm{M}}31| {{\boldsymbol{x}}}_{j},\delta {{\boldsymbol{x}}}_{j})/P(\mathrm{MW}| {{\boldsymbol{x}}}_{j},\delta {{\boldsymbol{x}}}_{j})},\end{eqnarray} \tag{ 8 }$$

where the proportionality factor in Equation (1) is equivalent for P(M31$| {{\boldsymbol{x}}}_{j},\delta {{\boldsymbol{x}}}_{j}$) and P(MW$| {{\boldsymbol{x}}}_{j},\delta {{\boldsymbol{x}}}_{j}$).

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}_{{\rm{M}}31,j}$ ≳ 0.75) M31 RGB (${N}_{{\rm{M}}31}=328$) or MW dwarf (N_MW = 92) candidates, excepting six stars with intermediate properties (0.5 < ${p}_{{\rm{M}}31,j}$ ≲ 0.75).¹⁰

Zoom In Zoom Out Reset image size

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}_{\mathrm{helio}}\,\sim \,-130$ km s⁻¹; 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}_{\mathrm{helio}}-\delta {v}_{\mathrm{helio}}$) > −100 km s⁻¹ 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}_{{\rm{M}}31,j}\gt 0.5$) to be M31 members in the following analysis. For stars in common between our data set 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 ($\langle {L}_{\mathrm{splash},j}\rangle \gt 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.

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}_{{\rm{M}}31}\gt 0.5$; Equation (8)). For smooth stellar halo fields, p_sub = 0. For fields with substructure, we computed p_sub using an equation analogous to Equation (8), where the Bayesian odds ratio is substituted with the relative likelihood that a M31 RGB star belongs to substructure versus the halo:

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{sub},j}=\displaystyle \frac{{f}_{\mathrm{sub}}{ \mathcal N }({v}_{j}| {\mu }_{\mathrm{sub}},{\sigma }_{\mathrm{sub}}^{2})}{{f}_{\mathrm{halo}}{ \mathcal N }({v}_{j}| {\mu }_{\mathrm{halo}},{\sigma }_{\mathrm{halo}}^{2})}.\end{eqnarray} \tag{ 9 }$$

Here 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 50th 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 recalculated p_sub for such fields (H, S, and D; Figure 1) based on their 50th percentiles, as noted previously. In the following abundance analysis (Section 5), we incorporated p_sub as a weight when determining the chemical properties of the stellar halo.

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 gradient between this work and Gilbert et al. (2014), we measured a ${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{phot}}$ (Section 2.3) gradient from our final spectroscopic sample of 197 halo stars, assuming t = 10 Gyr and [α/Fe] = 0 (Figure 8). We obtained a slope of −0.00 91 ± 0.001 9 dex kpc⁻¹ (−0.00070 ± 0.0016 dex kpc⁻¹) with an intercept of −0.71 ± 0.04 (−0.81 ± 0.03) when including (excluding) substructure. These gradients are inconsistent for r_proj ≳ 20 kpc, where the ${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{phot}}$ gradient including substructure remains flat as the ${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{phot}}$ gradient of the smooth halo declines. Such a difference is not detected from Gilbert et al.'s sample of over 1500 RGB stars spanning 10−90 kpc. However, the ${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{phot}}$ gradients measured in this work provide a more direct comparison to our spectral synthesis–based [Fe/H] gradients.

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 ${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{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, 2007, 2008). If the smooth stellar halo is systematically older than the tidal debris toward 30 kpc, this would steepen the relative ${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{phot}}$ gradient between populations with and without substructure. This agrees with our observation that ${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{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 ${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{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 and 35 kpc. Thus a more likely explanation for the discrepancy in slope between the ${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{phot}}$ and [Fe/H] gradients are uncertainties in the stellar isochrone models at the tip of the RGB.

In general, ${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{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 (Gilbert et al. 2020). 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 (Gilbert et al. 2020). 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).

Alternatively, the discrepancy between the ${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{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 (Section 5.1). We investigated the impact of potential bias from both (1) the exclusion of TiO stars (${b}_{\mathrm{TiO}}$) and (2) S/N limitations (${b}_{{\rm{S}}/{\rm{N}}};$ Section 5.1) on our measured radial [Fe/H] gradients by shifting each [Fe/H] measurement for an M31 RGB star in a given field by ${b}_{\mathrm{TiO}}+{b}_{{\rm{S}}/{\rm{N}}}$. As in Escala et al. (2020), we estimated ${b}_{\mathrm{TiO}}$ from the discrepancy in $\langle [\mathrm{Fe}/{\rm{H}}]{\rangle }_{\mathrm{phot}}$ between all M31 RGB stars in a given field (including TiO stars) and the final sample (excluding TiO stars). We calculated ${b}_{{\rm{S}}/{\rm{N}}}$ from the difference between $\langle [\mathrm{Fe}/{\rm{H}}]\rangle$ for the sample of M31 RGB stars with [Fe/H] measurements in each field (regardless of whether an [α/Fe] measurement was obtained for a given star) and the final sample. These sources of bias do not affect the [α/Fe] distributions.

Incorporating these bias terms yields a radial [Fe/H] gradient of −0.022 ± 0.002 dex kpc⁻¹ (−0.018 ± 0.001 dex kpc⁻¹) 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 and 34 kpc in M31's stellar halo are robust against these two possible sources of bias.

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) to place them on the same scale as M31's halo (Escala et al. 2019, 2020; Gilbert et al. 2019; this work) and other M31 dSph (Kirby et al. 2020) measurements. Kirby et al. (2020) found that the Vargas et al. measurements were systematically offset toward higher [Fe/H] by +0.3 dex compared to their measurements, based on an identical sample of spectra of M31 dwarf-galaxy stars. Kirby et al. did not find evidence of a systematic offset between their [α/Fe] measurements and those of Vargas et al. (2014a). Thus we adjusted [Fe/H] values measured by Vargas et al. (2014a) by the mean systematic offset of −0.3 dex and did not make any changes to their [α/Fe] values. This correction is reflected for the relevant M31 dSphs in Figures 9 and 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 Equation (4) of Escala et al. (2020),¹³

$$\begin{eqnarray}&&{w}_{i,j}=\displaystyle \frac{{M}_{\star ,j}{\rm{\Phi }}({L}_{V,j})/{N}_{j}}{{{\rm{\Sigma }}}_{j=1}^{{N}_{\mathrm{gal}}}{M}_{\star ,j}{\rm{\Phi }}({L}_{V,j})},\end{eqnarray} \tag{ 10 }$$

where ${M}_{\star ,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 distribution of the mock stellar halo. We refer the reader to Escala et al. (2020) for further details.

When resampling the observed inner halo abundance distributions, we assigned each star a probability of being drawn according to its p_halo. Not only did we find that the [α/Fe] distribution of the metal-rich ([Fe/H] > −0.8) and intermediate metallicity (−1.5 < [Fe/H] < −0.8) components of the smooth inner halo are inconsistent with that of present-day M31 dSphs at high significance ($p\ll 1$%), but also that the [α/Fe] distribution of the metal-poor ([Fe/H] < −1.5) component of the smooth halo disagrees with that of low-metallicity ([Fe/H] < −1.5) RGB stars in M31 dSphs ($p\lt 0.3 \%$).

Nevertheless, constructing [α/Fe] distributions according to metallicity bins presents a limited, 1D view of the relationship between [α/Fe] and [Fe/H]. Figure 9 displays [α/Fe] versus [Fe/H] for the inner halo of M31 compared to M31 dSphs with abundance measurements, including NGC 147 (Vargas et al. 2014a) and And X (Kirby et al. 2020). 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.

In order to robustly identify M31 halo stars with abundance patterns similar to that of M31 dSphs (Figure 9), we modeled the observed 2D 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 also 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 Section 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).¹⁴ 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 ($\delta {x}_{i}$, $\delta {y}_{i}$) = (δ[Fe/H]_i, δ[α/Fe]_i)¹⁵ given a model described by means ${\boldsymbol{\mu }}=({\mu }_{x},{\mu }_{y})$ and standard deviations ${\boldsymbol{\sigma }}=({\sigma }_{x},{\sigma }_{y})$ is

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal L }}_{i}({x}_{i},{y}_{i},\delta {x}_{i},\delta {y}_{i}| {\boldsymbol{\mu }},{\boldsymbol{\sigma }})=\displaystyle \frac{1}{2\pi {\sigma }_{{x}_{i}}{\sigma }_{{y}_{i}}\sqrt{1-{r}^{2}}}\\ \ \ \times \ \exp \left(-\displaystyle \frac{1}{2(1-{r}^{2})}\left[\displaystyle \frac{({x}_{i}-{\mu }_{x})}{{\sigma }_{{x}_{i}}^{2}}+\displaystyle \frac{({y}_{i}-{\mu }_{y})}{{\sigma }_{{y}_{i}}^{2}}\right.\right.\\ \ \ \left.\left.-\,\displaystyle \frac{2r({x}_{i}-{\mu }_{x})({y}_{i}-{\mu }_{y})}{{\sigma }_{{x}_{i}}{\sigma }_{{y}_{i}}}\right]\right),\end{array}\end{eqnarray} \tag{ 11 }$$

where i is an index corresponding to an individual RGB star. We incorporated the measurement uncertainties into Equation (11) via the variable ${\sigma }_{{x}_{i}}=\sqrt{{\sigma }_{x}^{2}+{\delta }_{{x}_{i}}^{2}}$, which is defined analogously for ${\sigma }_{{y}_{i}}$. 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 (Equation (11)) using a mixture model,

$$\begin{eqnarray}&&{{ \mathcal L }}_{i}=(1-{p}_{i}^{\mathrm{sub}}){{ \mathcal L }}_{i}^{\mathrm{halo}}+{p}_{i}^{\mathrm{sub}}{{ \mathcal L }}_{i}^{\mathrm{sub}},\end{eqnarray} \tag{ 12 }$$

where ${p}_{i}^{\mathrm{sub}}$ is the kinematically based probability that a star belongs to substructure in M31's stellar halo (Equation (9)). Thus the full stellar halo abundance ratio distribution is represented by an eight parameter model (${{\boldsymbol{\mu }}}_{\mathrm{halo}}$, ${{\boldsymbol{\mu }}}_{\mathrm{sub}}$, ${{\boldsymbol{\sigma }}}_{\mathrm{halo}}$, ${{\boldsymbol{\sigma }}}_{\mathrm{sub}}$), whereas we utilized a four parameter model (${{\boldsymbol{\mu }}}_{\mathrm{dSph}}$, ${{\boldsymbol{\sigma }}}_{\mathrm{dSph}}$) for M31 dSphs. The likelihood of the entire observed data set for a given stellar population (i.e., M31's stellar halo or the dSphs) is therefore the product of the individual likelihoods,

$$\begin{eqnarray}&&{ \mathcal L }=\displaystyle \prod _{i=1}^{N}{{ \mathcal L }}_{i}.\end{eqnarray} \tag{ 13 }$$

Using Bayes's theorem (see also Equation (1)), we evaluated the posterior probability of a particular bivariate model accurately describing a given a set of abundance measurements for a stellar population,

$$\begin{eqnarray}&&\begin{array}{l}P({\boldsymbol{\mu }},{\boldsymbol{\sigma }}{| }_{i=1}^{N}{x}_{i},{y}_{i},\delta {x}_{i},\delta {y}_{i})\\ \ \propto \ P{(}_{i=1}^{N}{x}_{i},{y}_{i},\delta {x}_{i},\delta {y}_{i}| {\boldsymbol{\mu }},{\boldsymbol{\sigma }})P({\boldsymbol{\mu }},{\boldsymbol{\sigma }}),\end{array}\end{eqnarray} \tag{ 14 }$$

where $P{(}_{i=1}^{N}{x}_{i},{y}_{i},\delta {x}_{i},\delta {y}_{i}| {\boldsymbol{\mu }},{\boldsymbol{\sigma }})$ is the likelihood (Equation (13)) and $P({\boldsymbol{\mu }},{\boldsymbol{\sigma }})$ represents our prior knowledge regarding constraints on ${\boldsymbol{\mu }}$ and ${\boldsymbol{\sigma }}$. We implemented noninformative priors over the allowed range for each fitted parameter, where we assumed uniform priors and inverse Gamma priors for nondispersion and dispersion parameters, respectively. We allowed the dispersion parameters in our model to vary between 0.0 < σ < 1.0, the correlation coefficients to vary between −1 < r < 1, and permitted samples of μ to be drawn from $-3.0\lt [\mathrm{Fe}/{\rm{H}}]\,\lt 0.0$ and $-0.8\lt [\alpha /\mathrm{Fe}]\lt 1.2$.

In the case of the variance, ${{\boldsymbol{\sigma }}}^{2}$, the inverse Gamma distribution is described as ${{\rm{\Gamma }}}^{-1}(\alpha ,\beta )$, where α is a shape parameter and β is a scale parameter. This distribution is defined for ${\boldsymbol{\sigma }}\gt 0$ and can account for asymmetry in the posterior distributions for dispersion parameters, which are restricted to be positive. The bootstrap resampled standard deviations (Section 5.3), weighted by the inverse variance of the measurement uncertainty and the probability of belonging to a given kinematical component, are σ([Fe/H]) = 0.53 ± 0.03, σ([α/Fe]) = 0.32 ± 0.02 for the smooth halo, σ([Fe/H]) = 0.51 ± 0.03, σ([α/Fe]) = 0.3 2 ± 0.02 for halo substructure, and σ([Fe/H]) = 0.49 ± 0.02, σ([α/Fe]) = 0.34 ± 0.02 for M31 dSphs. Owing to the similarity in the dispersion for each abundance ratio between the various stellar populations, we fixed the priors on ${{\boldsymbol{\sigma }}}^{2}$ to ${{\rm{\Gamma }}}^{-1}(3,0.5)$ for all ${\sigma }_{[\mathrm{Fe}/{\rm{H}}]}^{2}$ parameters and ${{\rm{\Gamma }}}^{-1}(3,0.2)$ for all ${\sigma }_{[\alpha /\mathrm{Fe}]}^{2}$ parameters. These distributions result in priors that peak at ${\sigma }_{[\mathrm{Fe}/{\rm{H}}]}\sim 0.50$ and ${\sigma }_{[\alpha /\mathrm{Fe}]}\sim 0.30$ with a standard deviation of 0.3 dex.

With the provided formulation, we sampled from the posterior distribution of each model (Equation (14)) using an affine-invariant 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 (Equation (12)), and separately for a model corresponding to M31 dSphs. We evaluated each model using 100 walkers and ran the sampler for 10⁴ steps, retaining the latter 50% of the MCMC chains to form the converged posterior distribution. Table 3 presents the final parameters of each 2D chemical abundance ratio distribution model. We computed the mean values from the 50th percentiles of the marginalized posterior distributions, whereas the uncertainty on each parameter was calculated relative to the 16th and 84th percentiles.

Table 3.  2D Chemical Abundance Ratio Distribution Model Parameters

Model	${\mu }_{[\mathrm{Fe}/{\rm{H}}]}$	${\sigma }_{[\mathrm{Fe}/{\rm{H}}]}$	${\mu }_{[\alpha /\mathrm{Fe}]}$	${\sigma }_{[\alpha /\mathrm{Fe}]}$	r
Smooth Halo	−1.27 ± 0.05	0.54 ± 0.04	0.38 ± 0.03	0.23 ± 0.03	−0.022 ± 0.09
Substructure	−0.88 ± 0.06	0.38 ± 0.05	0.42 ± 0.05	0.18 ± 0.03	$-{0.320}_{-0.14}^{+0.17}$
dSphs	−1.60 ± 0.03	0.45 ± 0.02	0.12 ± 0.02	0.26 ± 0.02	−0.297 ± 0.05

Note. The columns of the table correspond to the model for the observed 2D chemical abundance ratio distributions (Figure 10) for a given stellar population, and the parameters describing the bivariate normal probability density functions (Equation (11)). We fit for the smooth halo and substructure components simultaneously (Equation (12)) using a mixture model. The parameters for the 2D chemical abundance models are the 50th percentiles of the marginalized posterior probability distribution functions (Section 5.4.2), where the uncertainties on each parameter were calculated from the corresponding 16th and 84th percentiles.

Download table as:  ASCIITypeset image

Figure 10 displays the resulting bivariate normal distribution models for the 2D chemical abundance ratio distributions, reflecting the general trends of the observed abundance distributions (Section 5.4.1) for the smooth halo, halo substructure, and dSph populations. The models predict a larger separation between the mean metallicity of the smooth halo and substructure populations, where ${\mu }_{[\mathrm{Fe}/{\rm{H}}]}^{\mathrm{sub}}-{\mu }_{[\mathrm{Fe}/{\rm{H}}]}^{\mathrm{halo}}\,=0.39\pm 0.08$, compared to $\langle [\mathrm{Fe}/{\rm{H}}]{\rangle }_{\mathrm{sub}}-\langle$[Fe/H]${\rangle }_{\mathrm{halo}}=0.09\,\pm 0.06$ from a weighted bootstrap resampling (Section 5.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,

$$\begin{eqnarray}&&{p}_{\mathrm{dSph}}=\left(\displaystyle \frac{{{ \mathcal L }}_{i}^{\mathrm{dwarf}}/{{ \mathcal L }}_{i}^{\mathrm{halo}}}{1+{{ \mathcal L }}_{i}^{\mathrm{dwarf}}/{{ \mathcal L }}_{i}^{\mathrm{halo}}}\right)\times (1-{p}_{\mathrm{sub},{\rm{f}}}),\end{eqnarray} \tag{ 15 }$$

where ${{ \mathcal L }}_{i}$ is the likelihood from Equation (11) and ${p}_{\mathrm{sub},{\rm{f}}}$ is the refined probability that a star belongs to kinematical substructure (see Equation (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}_{\mathrm{sub}}{{ \mathcal L }}_{i}^{\mathrm{sub}}$/(($1-{p}_{\mathrm{sub}}$) ${{ \mathcal L }}_{i}^{\mathrm{halo}}$). Then, we computed ${p}_{\mathrm{sub},{\rm{f}}}$ for each star from the 50th 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}_{\mathrm{sub},{\rm{f}}}\gt {p}_{\mathrm{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 $\langle [\mathrm{Fe}/{\rm{H}}]\rangle$ or $\langle [\alpha /\mathrm{Fe}]\rangle$ 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 higher-dimensional 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⁻¹ < v_helio ≲ −150 km s⁻¹) 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 Sections 6.2 and 6.3.

Potential bias from the omission of red, presumably metal-rich, stars with strong TiO absorption (Section 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 (Section 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 (Kirby et al. 2020), and their abundance measurements are similarly affected by S/N limitations (Section 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.

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 high-resolution spectroscopy in the near-infrared H-band, in contrast with our measurements, which we derive from low- ($R\sim 3000$) and medium- ($R\sim 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 (Holtzman et al. 2015) 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; Escala et al. 2019, 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.

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 buildup 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% and 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. 2013, 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 in situ halo forms mostly from gas stripped from accreted galaxies. For galaxies that have not experienced a recent ($z\lesssim 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] $\gt \ -1$ in the solar neighborhood. Tissera et al. (2012, 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, 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 (Belokurov et al. 2018; Deason et al. 2018; Haywood et al. 2018; Helmi 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; Conroy et al. 2019; di Matteo et al. 2019; Gallart 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, high-eccentricity, α-rich stars on retrograde orbits²² with $[\mathrm{Fe}/{\rm{H}}]\gt -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 (Ibata et al. 2005). Although these stellar structures may reasonably be various in situ 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}_{\mathrm{proj}}\sim 9\,\mathrm{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 models 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. 2016). Indeed, multiple lines of evidence indicate that accretion has played a dominant role in the formation of M31's stellar halo (Section 6.3.2).

Based on this alone, both the dSph-like and non-dSph-like stars in M31's smooth, inner stellar halo (Section 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 progenitor(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 (Section 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 (Section 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 multidimensional 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 general consensus from simulations of MW-like galaxies is that massive dwarf galaxies (${M}_{\star }\sim {10}^{8-10}{M}_{\odot }$) distinct from present-day satellites are the dominant progenitors of the accreted stellar halo (Bullock & Johnston 2005; Robertson et al. 2005; Font et al. 2006; 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. 2006; 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; Fattahi et al. 2020), 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}_{\star }\lt {10}^{8}{M}_{\odot }$ become nonnegligible around ∼100 kpc. The disruption of globular clusters associated with the accretion of low-mass 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 MW-like galaxies with more active accretion histories have higher mass accreted progenitors (${M}_{\star }\sim {10}^{9-10}{M}_{\odot }$) 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) suggest that M31 has likely experienced multiple and continuous contributions to its stellar halo from accreted galaxies (Cooper et al. 2013; Deason et al. 2013; Font et al. 2020). Moreover, the total accreted stellar mass of M31 ($1.5\pm 0.5\times {10}^{10}{M}_{\odot };$ Harmsen et al. 2017) is on the upper end of both observed (Merritt et al. 2016; Harmsen et al. 2017) and predicted (Cooper et al. 2013; 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 (Gilbert et al. 2014, 2020; 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}_{\star }\sim {10}^{9-10}{M}_{\odot }$. 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}_{\star }\sim {10}^{9-10}{M}_{\odot }$ accreted between ∼1 and 4 Gyr ago, depending on the merger scenario (Fardal et al. 2007, 2008; D'Souza & Bell 2018b; Hammer et al. 2018).

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. 2018, 2019). In particular, Leethochawalit et al. (2018, 2019) found that the normalization of the stellar mass–metallicity relation evolves by 0.04 ± 0.01 dex Gyr⁻¹. Assuming that this relation extends to ${M}_{\star }\,\lt {10}^{9}\ {M}_{\odot }$, this means that an SMC-mass galaxy (${M}_{\star }\sim {10}^{8.5}\ {M}_{\odot }$) accreted 10 Gyr ago would have $\langle [\mathrm{Fe}/{\rm{H}}]\rangle =-1.39\pm 0.04$, compared to $\langle [\mathrm{Fe}/{\rm{H}}]\rangle =-0.99\pm 0.01$ at z = 0 (Nidever et al.2020). 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}_{\star }\sim {10}^{5-8}{M}_{\odot }$) contribute the remainder of metal-poor 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 (Section 6.2), owing to their defining similarity to the abundance patterns M31 dwarf galaxies (Section 5.4.2). Given that (1) the chemical abundance patterns of M31's dSph-like population are broadly consistent with the low metallicity ([Fe/H] ≲ −1.5) and early evolution of the MCs (Section 6.2) and (2) we do not detect any kinematical substructure for this population (Section 5.4), the accreted galaxies(s) that contributed to M31's dSph-like population could have been massive (${M}_{\star }\sim {10}^{8-9}{M}_{\odot }$) 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\,\lesssim$ [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 stellar mass of this secondary progenitor should be approximately between the mass of the SMC and LMC (${M}_{\star }\sim {10}^{8.5-9.5}{M}_{\odot };$ Section 6.2) and possibly comparable to Gaia–Enceladus (Section 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 dSph-like 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. 2006; 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 SMC-mass 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)\lesssim \langle [\mathrm{Fe}/{\rm{H}}]\rangle \lesssim -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}_{\star }\sim {10}^{7-8}\ {M}_{\odot }$ may be necessary to explain these low-metallicity stars in an accretion origin scenario.