The Hubble Space Telescope Advanced Camera for Surveys Emission Line Survey of Andromeda. I. Classical Be Stars

We observed M31 with the HST in the 12 orbit GO-13857 program (Dalcanton PI). Pure-parallel observations were made with the Advanced Camera for Surveys (ACS)/WFC and WFC3/UVIS camera at three distinct spatial locations that coincided with regions of the galaxy previously observed via the HST/PHAT survey (Dalcanton et al. 2012), yielding six distinct pointings across M31 (Figure 1). We obtained a second epoch of observations of these same six pointings approximately 1 yr later, when HST was able to achieve the same view of these pure-parallel pointings. Observations were made in both cameras using the F658N narrow-band filters and using the F625W broadband filter. A summary of our basic observational design is provided in Table 1. As discussed and explained in Section 3.3, although we present the emission line sample derived from both the ACS/WFC and WFC3/UVIS data, we restricted our Be star statistical analysis to the ACS/WFC data to avoid biasing our results due to the slightly different central wavelengths of the F658N filter between the two cameras.

We utilized the same methodology and pipeline for extracting photometry from our data as used in the PHAT survey, as described in Dalcanton et al. (2012) and Williams et al. (2014). In brief, the images from each epoch of observing (2014 and 2015) for each of the six fields were separately processed through our cloud-based photometry pipeline (Williams et al. 2018), which uses the DOLPHOT photometry package (Dolphin 2000) to force a fit of the precomputed point-spread function on every resolved star in the image stack of greater than 2.5σ significance above the expected background level. The photometry from these stars was then aligned and matched to the full 117 million source catalog from the full PHAT survey as presented in Williams et al. (2014) to allow us to compare the stars' narrow-band properties to their six band spectral energy distributions (SEDs). Overall, the photometric reductions were of comparable quality to those of the PHAT survey, except that we found a systematic uncertainty between epochs. Our analysis showed that our broadband exposures covered a short enough portion of the orbit that the thermal expansion and contraction of the optics caused a significant offset between the orbit-averaged point-spread function (PSF) used by DOLPHOT and that of the F625W images across epochs, which introduced a systematic error that we further discuss in Appendix.

Next, our source catalog was enhanced by adding the output of the Bayesian Extinction And Stellar Tool (BEAST; Gordon et al. 2016) on the PHAT catalog data, which added T_eff, log g, and age information for each source. Gordon et al. (2016) developed the BEAST to model the observed broad SED of an individual star by considering the full observational uncertainties. The BEAST primarily returns the posterior probability distributions for the intrinsic stellar properties of individual stars (e.g., stellar temperature, age, surface gravity) and their line-of-sight dust information (dust column density, grain size distribution, composition). The BEAST has been successfully used for many multiband HST programs such as PHAT (PI: Dalcanton), SMIDGE (PI: Sandstrom), METAL (PI: Roman-Duval), and the Local Volume UV survey (PI: Gilbert). Specifically, Gordon et al. (2016) ran the BEAST for all PHAT stars that were detected in F475W and at least three more bands (among F225W, F336W, F814W, F110W, and F160W) with resolutions of 0.15 in log age, 0.15 in Av, 0.25 in fA, and Z = 0.004, 0.008, 0.019, and 0.03. For this study, we cross-matched this PHAT-BEAST catalog with our classical Be star candidates, and used the BEAST measurements to better categorize and understand Be stars in M31. Finally, we corrected our data for focus-induced offsets, as described in Appendix.

Zoom In Zoom Out Reset image size

The initial reduced data set contained approximately two million stars in each epoch (Figure 2). As this paper focuses solely on the analysis of the classical Be population of the survey, we applied several criteria to isolate the population of main-sequence B-type stars with high quality data from the larger survey database. Our basic process was to utilize the value-added nature of BEAST-derived fundamental stellar parameters to isolate the population of main-sequence B-type stars in our sample, and then use traditional color–magnitude diagrams (CMDs) to confirm the robustness of these inferred stellar parameters.

Zoom In Zoom Out Reset image size

We began our sample definition by requiring sources in the new observations to have a signal-to-noise ratio (S/N) of >10 in the F475W, F625W, F658N, and F814W filters, which are the primary filters that we used in our subsequent analysis. Next, we required sources have BEAST-derived 50th percentile effective temperatures that resided within the range expected for B-type stars, 4.00 < Log(T_eff) < 4.51 (Underhill et al. 1979). We also required that sources have a BEAST-derived 50th percentile surface gravity, inclusive of +1σ_err computed from the 84th percentile surface gravity, of Log(g) > 3.75, to isolate main-sequence stars and exclude giants and supergiants (Castelli & Kurucz 2004).

We inspected the distribution of these selected B-star candidates on CMDs to determine whether uncertainties associated with BEAST-derived stellar parameters were introducing large populations of likely false-positive B-type main-sequence stars. We found that adopting a more conservative surface gravity criteria (Log(g)₁₆ > 3.75, where Log(g)₁₆ is the 16th percentile of the Log(g) posterior probability distribution) led to a 4× reduction in the number of sources, but these extra sources that were removed occupied the same region on the CMD as our higher confidence main-sequence B stars. As such, we adopted our less conservative criterion (Log(g)₈₄ > 3.75), which had the net effect of being most effective at excluding only the supergiants from our sample while retaining the main-sequence stars. Based on these inspections, we included additional criteria to reject sources having (F475W−F814W) < 0, which were likely planetary nebulae (see, e.g., Veyette et al. 2014), and also rejected objects having (F475W−F814W) > 1.5, as these sources had BEAST-derived Log(g) values near 3.75 and were clearly isolated from the bulk population with larger Log(g) values, indicating they were likely post-main-sequence giants. Finally, we excluded a small number of sources in our catalog that were assigned the same R.A./decl. coordinates by our matching routine between filters, due to source confusion between different survey tiles. The final data set for pointings 4–6 matching the aforementioned criteria contained 8981 stars in epoch one and 9098 in epoch two, and their CMD positions are plotted in Figure 3.

Zoom In Zoom Out Reset image size

Classical Be stars by definition are stars that exhibit or have exhibited Hα emission, which is interpreted as originating from their circumstellar gas disks (Jaschek et al. 1981). As such, one technique commonly used to identify classical Be stars is the use of photometric two-color diagrams (2-CDs) that utilize a narrow-band filter centered on Hα and an associated filter that samples the nearby continuum region (see, e.g., Grebel et al. 1992; Keller et al. 1999; McSwain & Gies 2005a, 2005b; Wisniewski & Bjorkman 2006). We discuss some of the limitations of this diagnostic in more detail in Section 3.2.1. We adopted a similar diagnostic approach here, using new data in the F625W and F658N HST filters from our GO-13857 program to diagnose the presence of emission or absorption in Hα, along with archival data in the F475W and F814W HST filters from the PHAT survey (Williams et al. 2014).

We first computed synthetic photometry through the known properties of HST's optics and cameras for normal B-type main-sequence stars spanning a range of B subtypes in the F625W and F658N filters. We began by using pysynphot (STScI Development Team 2013) and the Castelli-Kurucz Atlas (Castelli & Kurucz 2004) to create synthetic stellar spectra for five subtypes, B0, B1, B3, B5, and B8, adopting the stellar parameters listed in Table 2, based on average values for each subtype compiled in Underhill et al. (1979). These spectra were created using pysynphot's Icat function, which interpolates between the discrete models in the Castelli-Kurucz Atlas. We adopted a color excess of E(B − V) = 0.4, appropriate for moderately reddened stars in M31 (Clayton et al. 2015), along with R_v = 3.1. This adopted color excess is consistent with the average E(B−V) values computed for our sample of M31 B-type stars via the BEAST (Gordon et al. 2016). We note that the adoption of a color excess appropriate for moderately reddened stars will lead to a more stringent color criterion for the identification of Be stars in our data, such that our Be star rates will be a lower limit. We used the Extinction function to apply reddening to our spectra, yielding photospheric Hα absorption line profiles that are representative of main-sequence B-type stars (or classical Be stars when they are in a "diskless" state). We compile the equivalent widths of these Hα absorption lines in Table 2.

The amount of Hα emission that a classical Be star exhibits is to first order representative of the amount of disk material present within the stellocentric distances probed by Hα (Rivinius et al. 2013) at any epoch. We next used the absorption spectra we created in the previous step to create emission line spectra representative of a classical Be star with an arbitrarily strong disk (e.g., having Hα EWs as compiled in Table 2, albeit in emission), and spectra whose photospheric absorption is merely filled in with emission, representative of a classical Be star with a more tenuous disk (e.g., having Hα EWs of 0). We fit a Guassian to the Hα line using curve_fit from scipy.optimize, and then used this expression to create profiles with filled-in photospheric absorption and pure emission signatures. Finally, we computed the synthetic magnitudes in the F625W and F658N filters for the ACS and WFC3 cameras for these spectra for each B subtype to allow us to deduce which of the B-type stars in our sample were plausibly classical Be stars.

Our identification of classical Be stars was made using two-color diagrams (2-CDs) developed for each of our spectral subtypes (e.g., Figures 4–8) set by the 50th percentile T_eff reported by the BEAST. The boundaries of different spectral subtype assignments are simply the midpoint T_eff between subtypes listed in Table 2. Because the T_eff values were inferred from broadband photometric observations by the BEAST, we expect our spectroscopic subtype assignments could differ by 1–2 subtypes. However, because we binned the spectroscopic classifications in Section 3.3, minor misclassifications do not impact our overall results. We also explored effects of including 1σ_err in our 50th percentile T_eff source catalog cuts, and found that this only changed our final fractional Be ratios ±1%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The dotted lines in each 2-CD demarcate the expected values for pure photospheric absorption (black), filled-in absorption (blue), and pure emission (red) profiles from our synthetic spectra. We adopted conservative limits for identifying sources as classical Be stars, requiring their (F625W−F658N) colors to be redder than the photospheric absorption level (e.g., [(F625W–F658N) − $3\sqrt{{({\mathtt{F}}{\mathtt{625}}{{\mathtt{W}}}_{\mathrm{err}})}^{2}+{({\mathtt{F}}{\mathtt{658}}{{\mathtt{N}}}_{\mathrm{err}})}^{2}}$] > photospheric levels). Stated differently, we conservatively defined normal B-type stars without circumstellar disks as either (1) being located below the dotted black line or (2) including the black line within their 3σ error interval in Figures 4–8. Such normal B-type stars are depicted by black circles in these figures.

We further broke down the Be population into two categories: stars whose 3σ error interval extended blueward of the pure emission cutoff (potentially weaker disk systems, depicted by blue triangles) and stars whose 3σ error interval denoted it unambiguously as a strong disk system (depicted by red triangles). We were unable to identify any correlations between the amplitude of (F625W–F658N) excess and (F475-F814W) colors, indicating that the presence of red optical continuum excesses (Rivinius et al. 2013) did not bias us against detecting the strongest disk systems based on our selection criteria.

3.2.1. Limitations of the 2-CD Technique

Our methodology identified 552 Be stars out of 8981 (B + Be) stars in epoch 1, and 542 Be stars out of 9098 (B + Be) stars in epoch 2. The resultant Be fraction (# Be/(# normal B + # Be)) for each pointing is shown in Table 3, and is broken down as a function of spectral subtype in Table 4. We note that pointings 1–3, obtained from the WFC3 camera (see Section 2.1), are excluded from these statistical analyses for the remainder of the paper, as the filter utilized (F658N) only sampled the wing of Hα, therefore biasing us against robustly detecting weaker disk signals. We provide comprehensive information about each normal B-type star and classical Be star in the catalog in Table 8.

The overall average # Be/(# normal B + # Be) stars (in pointings 4–6) is ∼6%, with no statistically significant difference between epoch 1 (6.15% ± 0.26%) and epoch 2 (5.96% ± 0.25%). These overall Be fractions are strongly impacted by the larger number of later spectral type sources we detected, as seen in Table 4 and Figure 9. In particular, Figure 9 illustrates the clear trend of the Be fraction as a function of spectral type. The Be fraction is largest for B0-type stars and progressively decreases for later subtypes.

Zoom In Zoom Out Reset image size

We caution that our detection rates as a function of spectral type in Table 4 and Figure 9 are likely subject to several biases. The total main-sequence B-type star population (i.e., column 4 in Table 4) exhibits incompleteness based on expectations from an initial mass function (IMF). Assuming no incompleteness at the B0 subclass and adopting a Kroupa IMF (Kroupa 2001), our B-type star population is incomplete by a factor of 1.2× (B1), 1.6× (B3), 2.1× (B5), and 6.3× (B8). Of course if this incompleteness affects our recovery of Be stars at the same rate it does normal B-type stars, the fractional Be content we nominally derived in Table 4 and Figure 9 would remain unchanged. We also note that the larger (F625W−F658N) photometric errors present for later spectral subtypes acts as another factor that biases us against detecting weaker disk systems in later spectral subtypes. Since we do not a priori know the intrinsic prevalence of strong versus weaker disk systems as a function of spectral subtype, we do not quantify an upper limit to our quoted Be fractions owing to this systematic.

The wide (year-long) time baseline of our two epochs of observations allows us to assess the statistical prevalence of disk-loss and disk-renewal episodes, albeit the poor temporal sampling limits our ability to constrain the detailed timescale of these events. Because of the circumstellar disk's Hα-emitting properties, large changes in the size of these disks will cause sources to significantly move in 2-CDs. We classified a Be star as exhibiting a disk-loss or disk-renewal event if it was determined to be a Be star in one epoch and its (F625W−F658N) color was at or bluer than the expected value for a normal B-type star for the other, and vice-versa. We show the evolution of such disk-loss and disk-renewal events in 2-CD space for B0V-type stars in Figure 10. Following McSwain et al. (2009), we defined the percentage, p, of transient Be stars, as the ratio of Be stars that exhibit either disk-loss or disk-renewal events, n_Betrans, relative to the total number of Be stars, n_Betotal:

$$\begin{eqnarray}&&p=\displaystyle \frac{{n}_{\mathrm{Betrans}}}{{n}_{\mathrm{Betotal}}}.\end{eqnarray} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

We also defined the error in this percentage to be:

$$\begin{eqnarray}&&{\sigma }_{p}=\sqrt{\displaystyle \frac{p(1-p)}{{n}_{\mathrm{Betotal}}}}.\end{eqnarray} \tag{ 2 }$$

Approximately 13% of our sample exhibited a disk-loss event while ∼10% of our sample exhibited a disk-renewal event (Table 5). Roughly 77% of our sample maintained their status of having a disk at some level through the 1 yr baseline of our data. As seen in Figure 10 for our B0-type stars, such nontransient systems still exhibited evidence of variability in their Hα emission. Note that the total number of Be stars reported in both epochs in Table 5 is lower than the number reported in epoch # 1 or # 2 in Table 4 as some sources had too low S/N to meet our source catalog requirements in one of the two epochs.

Ascertaining whether the stars in our catalog belong to distinct clusters or a general field population is of interest as the uniform age of cluster members can allow one to assess how the Be phenomenon evolves with a star's fractional main-sequence lifetime. We matched our source catalog against clusters identified by the citizen science–based Andromeda Project and compiled in Johnson et al. (2015). Specifically, we required our sources to be within the quoted aperture radius for each cluster, R_ap, as compiled in Johnson et al. (2015).

Of the 8981 B + Be stars we detected in epoch 1, we classified 420 cluster stars in 85 unique clusters and 8561 field stars. We computed a corresponding fractional Be content of 16.19% ± 1.80% in cluster environments and a much smaller fractional Be content of 5.65% ± 0.25% in field environments (Table 6). Our results from epoch 2 are consistent with those found in epoch 1 to within errors. In epoch 2, of the 9098 B + Be stars, we classified 417 as cluster stars in 86 unique clusters and 8681 as field stars, and computed a fractional Be content of 14.39% ± 1.72% in cluster environments and a much smaller fractional Be content of 5.55% ± 0.24% in field environments (Table 6). Be stars are located in 23 unique clusters in epoch 1 and 27 unique clusters in epoch 2. We caution that our photometry is likely incomplete or of worse quality in the densest cluster environments, as source contamination from blends could be more detrimental.

Table 6.  Prevalence of Be Stars in Cluster vs. Field Environments

Epoch	Location	# Be Stars	# Normal B Stars	Be Fraction $\pm 1\sigma$
1	Field	484	8077	5.65%  ± 0.25%
1	Cluster	68	352	16.19%  ± 1.80%
2	Field	482	8199	5.55% ± 0.24%
2	Cluster	60	357	14.39% ± 1.72%

Note. The fractional Be content (# Be/(# normal B + # Be)) for objects identified as likely belonging to cluster and field environments is tabulated for both epochs of our observations.

Download table as:  ASCIITypeset image

The distribution of ages for our Be stars located in clusters, adopted from the cluster ages determined and tabulated by Johnson et al. (2016), is shown in Figure 11. We also computed the fractional main-sequence lifetime for each cluster star, defined as the ratio of a star's age to its main-sequence lifetime. This was done by assigning a mass to each inferred spectral subtype, taken from Table 1 of Silaj et al. (2014), and adopting main-sequence lifetimes for solar metallicity stars, assuming Ω/${{\rm{\Omega }}}_{\mathrm{crit},\mathrm{initial}}=0.8$, from interpolating Table 2 of Georgy et al. (2013). The resultant distribution of our cluster sample as a function of fractional main-sequence lifetime is shown in Figure 11. We remark that this methodology initially assigned 12 of our sources with fractional main-sequence lifetimes exceeding 1, likely due to improper main-sequence lifetime assignment related to the 1–2 spectral subtype uncertainty of our sample (Section 3.2), and were reassigned fractional main-sequence lifetimes of 1. We further discuss and interpret trends in these data in Section 4.2.

Zoom In Zoom Out Reset image size

A remarkable aspect of the Be phenomenon is dramatic episodes of complete disk-loss and disk-regeneration. Long-term photometric, spectroscopic, and polarimetric studies have suggested that complete dissipation of some disks occurs over timescales as long as one to several years (e.g., Wisniewski et al. 2010; Carciofi et al. 2012), whereas partial disk-growth episodes imprint a variety of clear observational signatures over timescales of days to months to years (Rivinius et al. 1998; Grundstrom et al. 2011; Draper et al. 2014; Labadie-Bartz et al. 2017, 2018). Both observational and theoretical studies have suggested that the disk dissipation phase can last for ∼2× longer than the typical timescale for disk build-up phases (Haubois et al. 2012; Vieira et al. 2017; Labadie-Bartz et al. 2018). Since each observational bandpass/diagnostic probes different physical regions within Be disks (Haubois et al. 2012; Rivinius et al. 2013), it is important to test this different disk growth versus dissipation timescale prediction of the viscous decretion disk model (Haubois et al. 2012) with a variety of observational techniques.

These events plausibly represent epochs in which the amplitude or efficacy of the disk driving mechanisms change. Diagnosing both the timescale and frequency of such events therefore seems likely to elucidate the fundamental mechanism(s) driving the phenomenon. Early studies of the frequency of transient events were limited, based only on small samples of Galactic clusters (e.g., 45 Be stars distributed across seven clusters; McSwain et al. 2008, 2009), revealing a wide dispersion in the observed transient percentage (0% ± 18% to 75% ± 22%, with a composite mean of 51% ± 7%; McSwain et al. 2009). Granada et al. (2018) characterized the IR color variability and location in IR CMD and 2-CD space of a sample of Be stars in open clusters to infer that ∼9%–15% of their sample of Be stars had a dissipating or small disk. Larger volume surveys have quantified that long-term variations (LTVs) in the R-band occur over many years in 37% (81/217) of Be stars (Labadie-Bartz et al. 2017). However, LTVs are not classified on the strict basis of complete dissipation of an observational signature of a disk or the regeneration of a disk from a diskless state, and include partial dissipation and build-up events.

Our survey provides a significant statistical boost to previous studies of the frequency of transient events, i.e., a 10× larger sample size than McSwain et al. (2009). Our sample includes both cluster and field stars, which differs from earlier studies that focused solely on cluster populations (e.g., McSwain et al. 2009). We found 22% ± 2% (100/452 cluster and field Be stars; Table 5) of Be stars experienced a disk-loss or disk-renewal episode within 1 yr, which is lower than the composite mean of 51% ± 7% (23/45 Be stars) reported by the smaller McSwain et al. (2009) study, although these transient rates do technically agree within 3σ of the quoted errors. The rate of disk-loss episodes that we observed across our entire sample (∼13%; Table 5) is similar to the fraction of dissipating or small disks inferred from the IR photometric study of Granada et al. (2018). We do caution that Be stars have been found to vary over periods of many years, and our one-year survey may not have captured the full transition period for many of these stars. Indeed, Figure 10 illustrates that many of our nontransient sources exhibit variability, some of which could be revealed to be the onset of transient events given higher photometric accuracy and/or additional epochs of observations.

We note that the compilation of transient Be fractions from McSwain et al. (2008, 2009) was computed from four years (NGC 3766) and two years (Collinder 272, Jogg 16, IC 2581, NGC 3293, NGC 4755, NGC 6231, and NGC 6664) of observations, so it is not clear how much the different durations of these surveys affect these statistical comparisons. If we make the simple assumption that, to first order, the Be transient fraction over several years scales linearly by a rate set by the duration of observations, the resultant yearly transient Be fractions are remarkably consistent. Namely, the 22% ± 2% yr⁻¹ we determined from our present study is consistent with that derived from NGC 3766 (17.3% ± 3% yr⁻¹; McSwain et al. 2008) and the remaining six clusters in McSwain et al. (2009); (20.5% ± 4.5%). Since the rate of disk-outburst events has been noted to be more prevalent in early-type stars compared to midtype Be stars (Labadie-Bartz et al. 2018), it could be interesting to further build up statistics on the rate of the transient Be fraction as a function of spectral type, to see if a similar trend emerges.

Finally, we remark that our detection of modestly more disk-loss events (57; Table 5) compared to disk-renewal events (43; Table 5) is unexpected given the duration of our sample (1 yr) and previous reports that the disk dissipation phase can last for ∼2× the typical timescale for disk build-up phases (Vieira et al. 2017; Labadie-Bartz et al. 2018). If disk dissipation is a slow process, we would expect to preferentially detect the more rapid disk renewal events, which is the opposite of what we observe. We caution that our adoption of conservative 2-CD criteria for the presence of a disk, requiring a 3σ elevation above Hα photospheric levels, could be introducing an observational bias in our data that misclassifies disks decreasing in mass as disk-loss events, and conversely underestimates small-scale disk-renewal events.

We explored how Be stars identified in cluster environments depended on cluster age and fractional main-sequence lifetimes in Section 3.5. Interpreting trends in these data first requires a discussion of the known limitations of our data set. Fundamentally, our analysis is enabled by characterization of the fundamental stellar parameters of our data set by use of the BEAST Bayesian toolset. However, as noted in Section 3.2, we do expect our spectral subtyping to have uncertainties of 1–2 subtypes. We see likely byproducts of these uncertainties in Figure 11. For example, the presence of early subtype Be stars in clusters >30 Myr old could imply these clusters have experienced multiple episodes of star formation and/or be a reflection that the spectral subtypes assigned to this population are incorrect by 1–2 subtypes. Similarly, we noted in Section 3.5 that we reassigned 12 cluster stars to have a fractional main-sequence lifetime of 1 (from >1), which was a likely byproduct of small mismatches in our spectral subtyping.

With this caveat in mind, we do see several interesting trends appear in Figure 11. First, we observe a clear population of Be stars at early fractional main-sequence lifetimes. This supports the idea that a subset of classical Be stars must emerge on the ZAMS as rapid rotators, as suggested by Wisniewski et al. (2007). Although the enhancement in the number of early-type Be stars observed at late fractional main-sequence lifetimes could support the idea that the Be phenomenon is enhanced with evolutionary age (Fabregat & Torrejón 2000; Martayan et al. 2010), we caution that we do not have sufficient statistics with our existing M31 data set to uniformly sample how the frequency of early-type Be stars changes as a function of fractional main-sequence lifetime.

Many studies of Galactic, LMC, and SMC Be stars have found that the fractional Be content of an environment is related inversely to its metallicity (Maeder et al. 1999; Martayan et al. 2006, 2010; Wisniewski & Bjorkman 2006; Iqbal & Keller 2013). The reported amplitude of this trend varies substantially in the literature, likely due to differences in completeness, filter systems, disparate stellar ages, and nonuniform ranges of spectral types among the data sets. This trend has been suggested to arise from differences in the average rotation rates of stars in these environments, with low metallicity environments having a larger population of rapid rotators (e.g., Maeder & Meynet 2000, 2001; Martayan et al. 2007b). Since the present day metallicity of M31 is approximately 1.5× solar (Clayton et al. 2015), our data set allows us to explore this trend on an extended metallicity range.

The spectral type dependence of the fractional Be content has been extensively studied in the Galaxy, LMC, and SMC (see, e.g., Zorec & Briot 1997; McSwain & Gies 2005b; Wisniewski & Bjorkman 2006; Mathew et al. 2008; Martayan et al. 2010; Rivinius et al. 2013). Although the decrease in percentage of Be stars with spectral type that we observed in M31 (Figure 9) is generally consistent with that reported in the literature in other environments, significant attention must be paid to factors such as sample size, completeness, and the methodology used to quantify Be stars in each study.

The number of single epoch Be stars we identified for statistical analysis in our M31 sample (∼550 Be stars; Table 4), is 1.6× larger than the full Galactic sample used by McSwain & Gies (2005b), 2.0× larger than the full Galactic sample used by Mathew et al. (2008), 1.7× larger than the slitless SMC sample of Martayan et al. (2010), 1.2× larger than the photometric SMC sample of Wisniewski & Bjorkman (2006), and 1.2× larger than the photometric LMC sample of Wisniewski & Bjorkman (2006). However, despite our large sample size, the total number of main-sequence B-type stars (# normal B + # Be) that we found as a function of spectral subtype (Table 4) is flat through the mid- to late-B-type stars. This deviation from that expected based on a nominal IMF distribution indicates our M31 sample is incomplete at mid- to late- spectral subtypes, similar to that reported for all SMC and LMC studies (see, e.g., Wisniewski & Bjorkman 2006; Martayan et al. 2010).

If we limit the detailed comparison of our M31 results with the literature to our earliest spectral subtypes, where completeness is less of an issue, we do find some interesting results. Figure 12 depicts the fractional Be content for early-type (B0–B3) stars as a function of metallicity, using Galactic Be star data from McSwain & Gies (2005b), LMC and SMC Be star data from Wisniewski & Bjorkman (2006), Maeder et al. (1999), and Keller et al. (1999), and M31 Be star data from the present study. The fractional Be content we observe in M31 is lower than that observed in the SMC or LMC, in agreement with metallicity trends previously noted in the literature (Maeder et al. 1999; Martayan et al. 2006, 2010; Wisniewski & Bjorkman 2006; Iqbal & Keller 2013). While Be stars are more prevalent in the SMC than in M31 for the earliest subtypes by a factor of 2.8×, this ratio is still lower than the 3–5× enhancement between the SMC and the Milky Way reported by Martayan et al. (2010). By contrast, our data indicate no clear change in the frequency of early-type Be stars in M31 and the Milky Way, despite their different metallicity. Since the current star formation rate of M31 is lower than that of the Milky Way (Yin et al. 2009; Robitaille & Whitney 2010; Lewis et al. 2015), it is possible that our results are a reflection that the Be phenomenon is enhanced by evolutionary age as stars progress through their main-sequence evolutionary lifetimes.

Zoom In Zoom Out Reset image size

The SPLASH survey (Dorman et al. 2012) conducted a Keck II/DEIMOS optical spectroscopic survey of ∼5000 isolated stars in M31, including regions observed by the PHAT survey. The SPLASH survey took into account the amount of crowding present during targeting to reduce the role of nearby neighbor contamination, and targeted sources brighter than F814W = 22 or F475W = 24 (Prichard et al. 2017). The survey identified a total of 146 main-sequence Be stars that overlapped with the PHAT footprint and determined a net fraction Be content (# Be/(# normal B + # Be)) of ∼12% (Prichard et al. 2017).

We identified a total of 14 emission line stars in the SPLASH Hα emission survey that overlapped with our HST footprint. These include 10 stars classified as Be stars, 2 classified as RHeB, 1 classified as AGB, and 1 classified as T-MS (Figure 13). SPLASH utilized 08 wide slitlets on Keck/DEIMOS. To investigate the potential role of source confusion and contamination in these spectra, we computed the number of sources in our HST footprint that fell within the Keck/DEIMOS aperture size for each of the 10 Be stars identified by Prichard et al. (2017) that fell within our FOV. We found 24 HST-detected sources in total were present within the Keck aperture for these 10 Be stars, and provide a summary of each of the HST-detected sources that fell within these Keck/DEIMOS apertures in Table 7. Most of these HST-detected sources had BEAST-derived stellar parameters that were inconsistent with them being main-sequence B-type stars: 7 HST-detected sources had T_eff more consistent with O-type stars and 14 HST-detected sources had Log(g) values more consistent with post-main-sequence giants. As seen in Table 7, each SPLASH Be star was matched with one source having a bright F475W magnitude, along with several less luminous stars. It is therefore possible that the Hα EWs reported by Prichard et al. (2017) have slightly elevated (or suppressed) signal from these fainter neighbors that could lead to spurious detections of Be stars in the survey. We note that none of the stars that we identified as Be stars in the analysis of our HST data were matched with Be stars identified by the SPLASH survey.

Zoom In Zoom Out Reset image size

Given the clear dependence of the Be phenomenon as a function of spectral subtype that we see in our M31 sample (Figure 9 and Table 4) and the shallower depth of source brightnesses probed by the SPLASH survey, we caution that the global (spanning all spectral subtypes) M31 Be fraction of ∼12% cited by Prichard et al. (2017) is strongly biased by overweighting contributions from the earliest B subtypes.

In a broader context, we noted in Section 3.5 that our HST survey also likely suffered from source confusion, especially in the densest cluster regions of M31. If the effects of source confusion serve to produce a greater frequency of false-positive detections of classical Be stars, as seems possible from our overlapping HST coverage of the SPLASH Hα footprint, this would imply that the Be fractions we derived in our study could be similarly affected by false-positive detections, especially in the densest stellar environments.