Variability of Red Supergiants in M31 from the Palomar Transient Factory

The evolution of massive stars, namely those with zero age main sequence (ZAMS) spectral types O and B, populating the upper part of the Hertzsprung–Russell (HR) diagram, remains one of the most puzzling fields in stellar physics. Its details are influenced by many uncertain physical processes, which include, among others, mass-loss, convection, and rotation (e.g., Langer & Maeder 1995; Maeder & Meynet 2000; Mauron & Josselin 2011; Georgy 2012; Beasor & Davies 2017; see Martins & Palacios 2013 for a comparison of different models employing varied prescriptions for these physical processes). Furthermore, there are many unknowns surrounding the mapping of these stars to their final fates, with cases of observational results conflicting with theory—e.g., what produces which types of core-collapse supernovae (SNe), and what results in black holes (Heger et al. 2003; Smith et al. 2011; Smartt 2015; Sukhbold et al. 2016).

Despite the incomplete knowledge of these aspects of the evolution of massive stars, they remain from birth to death one of the key players moderating and tracing properties of galaxies, such as the star formation rates, energetics—both radiation and mechanical—and chemical enrichment (see Massey 2013 for a recent review). Red supergiants (RSGs), which are core He-burning (or beyond) stars representing the red-most excursion in the HR diagram during the evolution of massive stars in the ZAMS mass range ∼10–30 M_⊙ (for Population I stars), are the largest and among the optically brightest, and thus easily detected stars in the Local Universe. As a result of direct detections in pre-explosion images, they are recognized as the progenitors of the most abundant type of core-collapse SNe, namely Type II-P (e.g., Van Dyk et al. 2003; Smartt et al. 2009; Groh et al. 2013).

In addition, many RSGs are variable, the properties of which have long been interpreted as pulsations (e.g., Stothers 1969; Guo & Li 2002). The driving mechanism for pulsation in RSGs is not yet fully understood but it is thought to be the κ mechanism in the hydrogen ionization zone, coupled with convection, whose feedback remains unknown (Heger et al. 1997; Yoon & Cantiello 2010). Despite this theoretical uncertainty, the existence of an observed period–luminosity (PL) relation for RSGs entails their potential use as extragalactic distance indicators, which is aided by their high luminosities and their higher prevalence compared to classical Cepheids (e.g., Glass 1979; Feast et al. 1980; Mould 1987).

Most early work on the PL relation focused on the infrared, particularly the K-band, for determining the luminosity, as these objects suffer from strong absorption in the optical, due to, for example, the presence of TiO in their atmospheres (Pierce et al. 2000), and the bolometric correction is small in the K-band. From such a relation, Jurcevic et al. (2000) obtained a distance to the galaxy M101 that was consistent with the Cepheid distance. However, studies investigating the PL relation for RSGs have been just a countable few (e.g., Feast et al. 1980; Wood et al. 1983; Jurcevic et al. 2000; Pierce et al. 2000; Kiss et al. 2006; Yang & Jiang 2011, 2012) in the almost four decades since the pioneering work of Glass (1979).

RSGs have periods of typically several hundreds of days, and a semi-regular and complex variability; consequently, a sufficiently long baseline is needed for their studies. This has been largely made possible either through decades-long collections of the American Association of Variable Star Observers (e.g., Kiss et al. 2006 for the Galactic RSGs) or recent surveys such as ASAS (Pojmanski 2002) and MACHO (Alcock et al. 1997; e.g., Yang & Jiang 2011, 2012 for the RSGs in the Magellanic Clouds). However, until now, there has not been a study of the variability of RSGs in the closest spiral galaxy to us, M31, even though a good number of them have been identified and cataloged (Section 2.1). This was due to a dearth of long-baseline time-domain surveys for M31. However, with recent surveys such as the (intermediate) Palomar Transient Factory (iPTF) survey (Law et al. 2009; Rau et al. 2009; Ofek et al. 2012), the WeCAPP survey (Riffeser et al. 2001) of the bulge of M31 (2000–2003; Fliri et al. 2006), and the PAndromeda project (2010–2012) of Pan-STARRS 1 (Kaiser et al. 2010; Lee et al. 2012), such analyses are now possible. iPTF, in particular, possesses a baseline ≈1.5–2.5 times longer than the other surveys and probes sources down to m_R ≈ 21 mag over the whole galaxy (cf. Section 2.2). In this work, we present the study of the variability properties of RSGs in M31 and a measurement of their PL relation using observations from iPTF.

Furthermore, while the derivation of physical parameters for SN progenitors based on direct detections or upper limits in archival images has become more complex in recent years, such analyses still currently neglect any possible variability of the progenitor star (e.g., Van Dyk et al. 2003; Smartt 2015; Davies & Beasor 2018). Our knowledge of the behavior of the star in the few years leading to the supernova is far from complete, with observational and theoretical work pointing to the possibility of pre-explosion eruptions, (e.g., Arnett et al. 2014; Shiode & Quataert 2014), enhanced mass-loss (e.g., Beasor & Davies 2016, 2017; Morozova et al. 2018), and the growth of pulsations (e.g., Yoon & Cantiello 2010). Even if pulsational variability similar to that observed in known RSG populations were still to be present during the archival imaging of the star, there could be consequence in interpreting its initial mass. In particular, if pulsation amplitude is luminosity-dependent, systematic effects could influence current observations, which point to a deficiency of higher-mass progenitors for SNe II-P (dubbed the "RSG Problem"; e.g., Smartt et al. 2009; see Section 5). As an illustrative example, we also assess quantitatively the possible extent of such an effect in this paper.

The paper is organized as follows. In Section 2.1, we introduce the sample of RSGs in M31. We describe the iPTF data in Section 2.2 and show the RSG light curves constructed from these data in Section 2.3 and the method to extract their periods in Section 2.4. We then derive the PL relation in Section 3 and compare it with theoretical MESA models in Section 4, and touch upon the implication of the variability of RSGs for SN progenitors in Section 5. We end with our conclusions in Section 6.

2.1. The Sample of RSGs and Their Physical Parameters

2.2. The iPTF Data

The iPTF survey, carried out with the 1.2 m Samuel Oschin telescope at Palomar Observatory, continued imaging M31 (including the outskirts of this galaxy) from its predecessor, the PTF survey (Law et al. 2009; Rau et al. 2009; Ofek et al. 2012), with the same wide-field detector covering more than 7 deg² with 11 active CCDs, and a typical spatial sampling of 2'' full-width at half maximum (FWHM). The baseline of the data for M31 analyzed here extends from 2012 May until 2017 February. The imaging was done in two passbands, R and g, with the largest fraction (>80%) in R at a limiting magnitude reaching ≈21 mag and an average cadence of 1 day (typically taking 2 observations per night).

We thus use only the R-band PTF/iPTF (hereafter simply PTF) data set of about 10,000 images. This set comprises over 1500 observation epochs of M31 covering 1.8 × 2.4 deg² with 6 CCDs (as shown in Figure 1), over the approximately 5 years of baseline. Difference imaging was performed on all these images via the difference imaging pipeline of the PTF collaboration, the IPAC/iPTF Discovery Engine (PTFIDE; Masci et al. 2018), and this pipeline also carried out detections of varying sources on the difference images, which we term here as "raw detections."

Zoom In Zoom Out Reset image size

2.3. Optical Light Curves of RSGs and Their Variability

The PTFIDE raw detections catalog, in principle, provides us with a means to determine whether a source in M31 is variable or not. In fact, we find from the results of cross-matching (with a search radius of 2'', consistent with the typical FWHM) that all RSGs lying within the PTF footprint (253) are flagged as variable by the PTFIDE pipeline. However, as already established by Masci et al. (2017) and Soraisam et al. (2018), these "raw detections" are largely dominated by artifacts of image differencing such as, for example, dipoles from imperfect PSF-matching, edges of masked image regions due to saturated stars, CCD defects, etc.

For our science goal, we do not want to compromise completeness by resorting to thresholding on output parameters of image differencing. Rather, we implement forced photometry to construct the RSG light curves, using the difference images themselves, and analyze their variability. Difference imaging significantly alleviates the problem of crowding in the M31 fields for measuring their fluxes. Subtraction artifacts in individual difference images will largely be reflected as rogue points in the light curves, which we effectively deal with (e.g., by masking those points) without discarding the source.

In the same manner as Soraisam et al. (2017), we perform aperture photometry at the positions of the M31 RSGs and apply a curve-of-growth correction to the measured fluxes. The latter is a different factor for each star applied to correct for the flux missed due to the limited size of the aperture (taken to be the FWHM). The subtracted fluxes of the sources are then added using the template/reference image PSF-fit photometry catalog and calibrated using the relevant zero-point in the science images (see Masci et al. 2017 for details). We drop 9 RSGs with bad photometry in a large section of their light curves, due to, for example, their location close to masked parts of images or in regions close to the bulge where the quality of image differencing is poor. Example calibrated light curves are shown in Figure 2.¹¹

Zoom In Zoom Out Reset image size

To examine the variability of these RSGs over the full baseline, we compute the root-mean-square (rms) deviations from the mean for all of the extracted light curves, and compare them to those of static stars. We define a static star as a star in the catalog of the reference image used in the image subtraction having no detection from the PTFIDE pipeline in all of the input images used in the subtraction (also termed science images by Masci et al. 2017). We extract the light curves for these static stars in the same manner used for the RSGs. Since the surface brightness of M31 is not uniform over the face of the whole galaxy, there is a gradient in the local background noise of a resulting difference image, which deteriorates severely, particularly toward the bulge. The RSGs, on the other hand, are largely distributed in the disk of the galaxy. To obtain a sample of static stars with a similar background noise distribution as the population of RSGs, we select the static stars that are located in proximity (within about 1') to any of the RSGs.

In Figure 3 we plot the distribution of the rms deviations found for these static stars as a function of their average R-band magnitudes over the baseline of the survey (blue points). Rms values increase toward fainter magnitudes, and the distribution we find for stars in M31 agrees well with that for M-dwarf stars in PTF derived by Law et al. (2012). Also shown in Figure 3 are the M31 RSGs (red circles), and we see that the variabilities (as determined by the rms deviations) of some of the RSGs, around the mean $\langle {m}_{R}\rangle$ range of 18.0–19.0, are comparable to those of static stars, i.e., their variability is consistent with noise. The filled circles indicate RSGs with measured periods (Section 2.4).

Zoom In Zoom Out Reset image size

We visually separate the RSGs with significant variability from the RSGs whose rms deviations are consistent with noise. To this end, we draw a straight line in the plane of $\langle {m}_{R}\rangle -{\rm{\Delta }}{m}_{R}$, with the slope being similar to the linear fit through the static stars, and try different intercepts to most efficiently select significantly variable RGSs. This demarcation is shown with the orange line in Figure 3. This method is not very stringent, but we verified that the main results of this paper remain the same, with small shifts in this line. Adopting this demarcation, we find that 167 RSGs, constituting ≈70% of our sample, have measurable variability from PTF. These are examined in the following section to determine what periodicity, if any, is present in their multi-year light curves.

In Figure 4, we plot the rms deviations for our sample of M31 RSGs versus their M_K magnitudes (left panel) and luminosities (right panel), whenever available from ME16. RSGs with significant variability are plotted in red, while those with PTF upper limits on variability (those below the orange line in Figure 3) are plotted in black.

Zoom In Zoom Out Reset image size

As can be seen, at low luminosities there is an overlap between the RSGs with detected variability and those below our detection threshold, despite the fact that they fully separate in the $\langle {m}_{R}\rangle -{\rm{\Delta }}{m}_{R}$ plane. This can be attributed to the varying degrees of extinction in the R-band suffered by individual RSGs, which subsequently modulate their signal-to-noise ratios (S/Ns). However, despite this overlap, it is notable that all of the "non-varying" RSGs are characterized by low intrinsic luminosity. All RSGs in our sample that are brighter than M_K ≈ −10 (and log L/L_⊙ ≈ 4.8) show variability, and many of them show variability at the level Δm_R > 0.1. We assume this variability of RSGs to be associated with pulsations (e.g., Heger et al. 1997; Guo & Li 2002). Furthermore, our result points to variability/pulsations already occurring for log L/L_⊙ ≳ 4.6. Theoretical calculations of Yoon & Cantiello (2010) could only find pulsations above log L/L_⊙ > 4.95 (e.g., their Figure 1). This is, however, due to the fact that the implicit scheme adopted in stellar evolution calculations introduces substantial damping and thus they could only place a lower limit on the existence of pulsations in RSGs.¹² Our result, on the other hand, shows that stars with masses lower than the 16–17 M_⊙ limit found by Yoon & Cantiello (2010) can also pulsate. In addition, we find a hint of a positive correlation between variability amplitude and brightness for the RSGs with detected variability. In Section 5, we will explore the influence that this observed variability can have on progenitor mass estimates of core-collapse supernovae, in particular the trend of larger amplitude variations for more luminous RSGs.

2.4. Period Determination

The light curves of RSGs are known to be semi-regular, and some of them are completely irregular (see Figure 2). The periodicity of unevenly sampled regular time-series/light curves has been efficiently handled by conventional standard algorithms such as, for example, the Lomb–Scargle periodogram (LS; Lomb 1976; Scargle 1982), phase dispersion minimization (PDM; Stellingwerf 1978), analysis of variance (Schwarzenberg-Czerny 1989), and the hybrid algorithm of Saha & Vivas (2017). However, for the RSGs, these methods tend to fall short, as these objects are characterized by one of the most complex light curve morphologies. This includes a strong red-noise component speculated to arise from the convective cells on the surface and multi-periodicity. For example, Kiss et al. (2006) found two distinct periods for a significant fraction of their Galactic RSG sample over a baseline extending to more than 5 decades: one on the order of a few hundred days, which is typically associated with RSG pulsations, and the other greater than 1000 days, similar to the long secondary periods in Miras with an unknown nature (see also Percy & Abachi 2013). Compared to this, our baseline will not be sensitive enough to pick up these longer (>1000 days) periods.

Bound by the limitations of the popular standard methods for our application, we turn to a relatively recent approach based on Gaussian Process (GP) modeling (a non-parametric statistical model that assumes the distribution over a set of random function values as a multi-variate Gaussian, typically with a mean of 0; see Rasmussen & Williams 2005). In principle, this Bayesian technique provides a statistically robust and powerful way to model any phenomenology, particularly ones for which prior knowledge may not be available. For a GP, determining the covariance matrix is the crux of the modeling, which in turn is determined by the power spectrum of the signal. In almost all of the existing off-the-shelf GP-based methods used by the stellar community, a functional form is assumed for generating the covariance matrix elements, assuming stationarity for the light curve (i.e., its statistical properties are assumed constant) or even periodicity with a single frequency (e.g., Wang et al. 2012). Given the complexities of the RSG light curves as mentioned above, these algorithms are not accurate for their analysis.

For our case, the GP modeling is performed using the so-called critical filter algorithm of signal reconstruction as implemented by Oppermann et al. (2013) using the NIFTy package of Selig et al. (2013) aimed toward cosmological signals, and applied to period analysis by N. Oppermann et al. (2018, in preparation). Simply put, the critical filter provides a methodology to reconstruct the underlying signal, whose power spectrum is unknown, by simultaneously also reconstructing this spectrum from the observed data. The computation is iterative. We start with an initial guess of the power spectrum that defines the prior distribution, and reconstruct the signal, thus obtaining a posterior distribution conditioned on the observed data and the given power spectrum, which is again a Gaussian. The latter distribution gives an improved estimate of the power spectrum, which is used to update the prior distribution, and then this is used to reconstruct the signal and so on. The iteration is stopped when the new estimate of the power spectrum converges. Since Kiss et al. (2006) have illustrated the presence of a 1/f component in the power spectra of these RSGs, this provides a natural choice for our initial guess, but we also impose the condition that there must have been at least five iterations before convergence. We obtain a similar reconstruction even when starting with a flat spectrum.

Examples of the resulting power spectra of the light curves are shown in the right panels of Figure 5, with the red-noise component evident. To extract the peaks in a spectrum, we fit a power law to the background red-noise and use that to set an S/N threshold of 7 above the background. We obtain the slopes for these RSGs in the range −1.51 to −0.86. We then select the peaks above this threshold with values greater than 10⁻³ day⁻¹ (constrained by our baseline). For cases with more than one peak, we choose the one with maximum power (i.e., the dominant mode) and for a few sources with blended peaks, we take their excess power-weighted average. We adopt these frequencies to define the periods of the light curves used for deriving the PL relation (Section 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The methodology we have used reconstructs the light curve along the complete time axis, fully taking into account the observational sampling. We have verified by shuffling the magnitude measurements while keeping the same irregular time sampling and calculating the power spectrum of the resulting light curve, that no artificial periods corresponding to the observational pattern emerge as power peaks.

A detailed analysis of how the results from other period-finding algorithms fare against that of the GP modeling is beyond the scope of this paper (it will be presented by N. Oppermann et al., 2018, in preparation). Nevertheless, for a simple comparison, we use available python implementations of common period-finding algorithms, specifically LS, supersmoother¹³ and PDM, to estimate the periods of some example light curves of M31 RSGs. We consider the RSGs with IDs 227 and 189 shown in Figure 5, with quite complicated shapes. The different algorithms give a range of periods: 139 days (LS), 1113 days (supersmoother), 472 days (PDM), 505 days (this work) for RSG_227; 263 days (LS), 1329 days (supersmoother), 786 days (PDM), 331 days (this work) for RSG_189. As can be seen from the power spectra of these two RSGs in Figure 5, the results from the other methods, though associated with power peaks in most cases, do not pick up the dominant peaks barring PDM for RSG_227.

We successfully obtain the pulsation periods (less than ≈1000 days) for 63 RSGs (the remaining 104 do not have a significant peak), and these are shown in Table 1. We mark these RSGs in Figure 4 with the filled red circles, and as is evident, some bright, highly variable sources do not have detectable periods, an example of which is shown in the last panel of Figure 5. It is quite likely that with an extended baseline (for example with the Zwicky Transient Facility, Bellm 2014), a period >1000 days could be found for such cases (see also Section 3). Another example of a bright but less variable, non-periodic source is also shown in Figure 5. A detailed investigation of why some of the sources in our sample show detectable periods and others, though variable, show none, is beyond the scope of this paper. As described above, the mechanism of pulsation in RSGs, and its interaction with convection, are not well understood.

Table 1.  Properties of the RSGs in M31^a

IDb	Starc	MKc	Teff (K)c	log(L/L⊙)c	$\langle {m}_{R}\rangle$	ΔmRd	Period (days)
1	J003950.86+405332.0	−10.03	3850.0	4.86	18.63	0.27	391
2	J003950.98+405422.5	−10.28	3650.0	4.89	18.46	0.11	⋯
3	J003957.00+410114.6	−9.4	3650.0	4.54	18.80	0.12	⋯
4	J004015.18+405947.7	−10.22	3700.0	4.89	18.62	0.23	391
5	J004015.86+405514.1	−8.86	3950.0	4.42	18.74	(0.08)	⋯
6	J004019.15+404150.8	−9.87	3750.0	4.76	18.72	0.10	⋯

Notes.

^aMeasurements that are not available are indicated by a dash. ^bOur ID follows the order in which the RSGs appear in the ME16 catalog. ^cThese values are obtained from ME16. ^dThe bracketed value is for an RSG whose variability is below our threshold (see Figure 3).

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

Download table as:  DataTypeset image

Figure 6 shows the periods of the RSGs obtained in the previous section and their absolute K-band magnitudes from ME16 (Section 2.1). Of the 63 RSGs with identified periods, 51 appear as a coherent band between log(P/days) ≈ 2.4 at M_K ≈ −9.5 and log(P/days) ≈ 3.0 at M_K ≈ −11.5. However, 12 RSGs follow a separate sequence (copper circles, leftmost portion of the plot), with 10 of them appearing at similar M_K < −9.5 as the bulk of the RSGs in the adjacent band, but with an offset of log(P/days) ≈ 0.3. The nature of these stars will be examined in Section 4.

Zoom In Zoom Out Reset image size

The remaining two stars of the 12 have periods <100 days (IDs 52 and 166, cf. Table 1). Short periods of around 100–150 days have been found for RSGs (e.g., Kiss et al. 2006), but no RSGs in the literature have pulsation periods as short as these stars, prompting us to examine them in detail. Both stars possess extremely regular light curves as compared to the other RSGs (see Figure 5), relatively high T_eff values (4300 K on the temperature scale of ME16), and relatively low luminosities (M_K ≥ −9.5 mag). These physical properties overlap with those observed for ultra-long-period Cepheids (Fiorentino et al. 2012). Indeed, star RSG 166 was classified as a Cepheid by Riess et al. (2012). As the RSG sample of ME16 was simply selected to have a temperature cooler than around 4300 K and be a member of M31 (as opposed to strictly stars on the RSG branch), it is not wholly unexpected that a small number of stars on the luminous and cool edge of the Cepheid instability strip would be identified.

Ignoring the 12 RSGs from above, we obtain a simple linear fit through the rest as

$$\begin{eqnarray}&&{M}_{K}=(-3.38\pm 0.27)\times \mathrm{log}P+(-1.32\pm 0.75),\end{eqnarray} \tag{ 1 }$$

shown by the blue line in Figure 6. The RSG P–M_K relation we have derived here is the first for M31. We see a dispersion around this relation of 0.29 mag. A part of it is possibly due to the errors in the adopted K-band magnitudes, which could reach a tenth of a magnitude (for example, the 2MASS values given by Massey et al. 2009 and used in ME16), and the fact that these K-band magnitudes are from single (nightly) observation epochs, while the amplitude in this band could be as high as 0.25 mag (Wood et al. 1983) (thus, a contribution of <0.18 mag to the dispersion). Also, circumstellar reddening may contribute to this dispersion; ME16, however, applied an average correction for M31 RSGs obtained by Massey et al. (2009) to the M_K values of these stars that should reduce this effect. But a possible part of the cause could also be that the dispersion is intrinsic for these objects—an apparent effect of additional parameters (cf. Section 4).

Given the broad trend of increasing variability with brightness (Figure 4) discussed in Section 2.3, and the positive correlation of period with brightness (Figure 6), the same broad trend can be expected between the periods and the rms deviations. This is indicated in Figure 6 by the gradient in the face-color of the points. In light of such a trend, it is feasible that a period longer than the range probed in this study is present in the highly variable sources without detectable periods (Figure 5, Section 2.4). Furthermore, the theoretical result of increasing pulsation amplitude with initial masses or luminosities (e.g., Yoon & Cantiello 2010) appears to be consistent with our observational result.

Figure 7 shows a comparison of our results for the P–M_K relation with those of others for different galaxies. Our slope and intercept are consistent with those of Kiss et al. (2006) for the Galactic RSG sample with 13 sources and a dispersion of 0.46 mag. Furthermore, they are in good agreement with those of the LMC, SMC, and M33 derived by Yang & Jiang (2012) using 47, 21, and 40 sources, respectively, with dispersion values of ≈0.3 mag (they claim their relation for the Galactic sample is less accurate due to uncertainties in distances).¹⁴ Thus, the data from several galaxies characterized by different metallicities (from sub-solar by a factor of ≈4 in the Small Magellanic Cloud to super-solar by a factor of ≈2 in M31, e.g., Massey 2003), are all consistent with a universal PL relation for RSGs, although the uncertainties still allow a wide range of parameters. Such universality may not be unexpected given current ideas on what drives the pulsation in RSGs (κ mechanism in hydrogen, Section 4), and this potentially could have important implications for the life of very low-metallicity/Population III stars (e.g., Moriya & Langer 2015).

Zoom In Zoom Out Reset image size

We investigate the pulsation mode represented by the periods we have found for the RSGs with theoretical models constructed using Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015, 2018), a standard tool for stellar evolutionary calculations. One of the latest updates presented by Paxton et al. (2015) includes fully incorporating existing methodology for asteroseismology computations (GYRE software of Townsend & Teitler 2013), and thus facilitates comparison of our observational results with theory.

To produce theoretical predictions for RSG pulsation periods, we use version 9793 of the MESA software. We consider stellar models with ZAMS masses of 12–24 M_⊙, which are evolved from ZAMS to depletion of Carbon in the core. Our MESA settings follow the recent work of Chun et al. (2018), including a calibrated mixing length parameter α_MLT = 2.7 to match inferred RSG temperatures in M31 at Z = 0.04. We also include step overshoot with parameter f_ov = 0.15 and mass-loss according to the "Dutch" wind prescription (de Jager et al. 1988; Vink et al. 2001). More details on MESA settings can be found in Chun et al. (2018). Our models do not include any rotation.

After the models leave the main sequence, we include adiabatic GYRE calculations in each step for T_eff < 4500 K. These calculations find all radial modes with periods in the range 10–10,000 days, which is more than adequate to identify the fundamental mode and first overtone in these RSG models. The left panel of Figure 8 shows the resulting period tracks for several different mass models, along with points representing 20 equally spaced time intervals along each of the tracks for the fundamental mode periods. When comparing to the observed RSG data in the right panel of the same figure, we show the same tracks only after the onset of core helium burning, where the models spend more than 95% of their RSG lifetime.

Zoom In Zoom Out Reset image size

Our modeling makes no attempt to discern which modes should be excited. The work of Heger et al. (1997) suggests that RSG pulsations are excited by the κ mechanism in the hydrogen ionization zone. Heger et al. (1997) also find that pulsation periods roughly scale as P ∝ L/M, in agreement with the analytic predictions of Gough et al. (1965). Our MESA models roughly agree with this scaling.

The right panel of Figure 8 shows that the bulk of the pulsating RSGs in our sample are consistent with the fundamental mode pulsation of the theoretical models. Furthermore, the sequence of 10 stars offset to shorter periods for a given M_K is broadly consistent with the first-overtone mode.

The HR diagram for cool supergiants is shown in Figure 9. Colored lines represent the same MESA models used to produce Figure 8, and the RSGs with measured pulsation periods are shown in gray (other M31 RSGs are plotted as blue circles for reference). The observed RSGs lie in proximity to MESA model tracks for broadly the same initial masses, whether period or T_eff is used as a parameter (right panel of Figures 8 and 9). The 10 stars with measured periods offset from the bulk of the sample in the PL plane that are not Cepheids are shown by the star symbols in the HR diagram. These stars with shorter periods overlap with the bulk of the RSG population with measured periods in their temperatures and luminosities. This strengthens our association of the variability observed in these stars with the first-overtone mode of stars with the same properties as the bulk population.

Zoom In Zoom Out Reset image size

As is already known, there is no single-valued mass–luminosity relation for the RSGs. Consequently, tracks of different masses at the same luminosity can smear out the period for the RSGs, and may be a reason for the scatter seen in Figure 6. Thus, the mass of the star, and more importantly its evolutionary stage, naturally appear as additional parameters contributing to some of the dispersion in the P–M_K relation. As a result, the precision of any extragalactic distance estimate based on this relation can never be better than the limit imposed by the intrinsic dispersion.

RSGs as progenitors of SNe II-P/L (e.g., Ekström et al. 2012) have been established based in part on the direct detection of some progenitor stars in the archival images of nearby galaxies hosting these SNe (e.g., Van Dyk et al. 2003; Smartt 2009). A recent comprehensive review of such observations spanning a 15-year period, and their implications for the progenitor masses of hydrogen-rich core-collapse supernovae, was presented by Smartt (2015; hereafter SS15). The general approach to determining progenitor masses from these pre-explosion observations is to convert the photometric measurements to bolometric luminosities and then compare to the initial mass–final luminosity relation from stellar models. One of the important results from such an analysis is that the maximum-mass for these progenitors is claimed to be at a much lower value (≈18 M_⊙; Smartt et al. 2009; SS15) than the high-mass limit of RSGs (25–30 M_⊙) predicted by stellar evolution theory.

However, various systematic errors, often resulting from the uncertain final stages in RSG evolution, can influence this result. In particular, when converting observed fluxes in a limited (often single) photometric band to bolometric luminosities, treatments of extinction and bolometric corrections are critical. For example, Davies & Beasor (2018) showed that the spectral types of RSGs evolve to later types as they approach core-collapse, resulting in a larger bolometric correction than typically assumed for the imaged progenitor stars. They thus pointed out that failing to account for the spectral evolution of RSGs can introduce a systematic error in the adopted bolometric correction.

Most of the archival detections of supernova progenitors are from single-epoch observations, and we have demonstrated in Section 2.3 (see Figure 4) that while lower luminosity RSGs have observed variability of ≲(5–10)%, a large fraction of the RSGs with luminosities above log(L/L_⊙) ≈ 4.8 show a significant variability. Thus, if our observed sample of M31 RSGs is representative, the intrinsic variability of RSGs will preferentially affect the initial mass estimates of high luminosity (and therefore high initial mass) progenitor stars. It is therefore pertinent to examine how big such an effect can be.

Recently, Johnson et al. (2017; see also Kochanek et al. 2017) examined the variability of four SN II-P progenitor stars, from approximately four to eight years prior to the supernovae. They derived upper limits on their variability at (5–10)% of their R-band luminosities, corresponding to Δm_R < 0.05–0.10. Three of these SN progenitors had sufficient observations to derive bolometric luminosities, all of which were in the range log(L/L_⊙) = 4.5–4.9 (Davies & Beasor 2018; Kochanek et al. 2017; van Dyk et al. 2017). Thus, as is evident from Figure 4 (right panel), the upper limits on their variability are also consistent with our results for the RSGs in M31.

In this section, we make a quantitative assessment of whether intrinsic variability of RSGs can have any significant effect on the SN progenitor mass estimates.

In converting a photometric measurement in some filter (m_λ) to the bolometric luminosity (${ \mathcal L }=\mathrm{log}(L/{L}_{\odot })$), a correction (dm_λ to obtain the mean value $\langle {m}_{\lambda }\rangle$) for possible variability of the progenitor can be applied as

$$\begin{eqnarray}&&2.5\cdot \langle { \mathcal L }\rangle ={M}_{\odot }-({m}_{\lambda }+{{dm}}_{\lambda }-{A}_{\lambda }+{\mathrm{BC}}_{\lambda }),\end{eqnarray} \tag{ 2 }$$

where the term inside the brackets on the right side is the bolometric magnitude of the star, A_λ and BC_λ are the extinction and bolometric corrections, respectively, for the star in the given filter, and M_⊙ is the bolometric magnitude of the Sun. In the absence of color information to constrain the BC_λ of the star, typically the value for an M0 supergiant is assumed, with its 1σ error as 0.3 mag, corresponding to the standard deviation of the bolometric corrections spanning the spectral types from late K to late M (see Smartt et al. 2009, SS15).

The detections and upper limits of the progenitor stars have been acquired mostly in the I- and R-bands. Since we do not have information on the intrinsic variability of our RSG sample in the I-band, we assume a similar level of variability as in the neighboring R-band (Figure 4). Then, given a random-phase measurement m_λ from the pre-explosion image of the star, and hence, a random-phase luminosity ${ \mathcal L }$ computed from it, we are interested in the mean luminosity $\langle { \mathcal L }\rangle$ for the star. To this end, we use Bayes' Theorem to derive the posterior probability distribution for the mean luminosity $\langle { \mathcal L }\rangle$, given the observed random-phase luminosity ${ \mathcal L }$, according to

$$\begin{eqnarray}&&P(\langle { \mathcal L }\rangle | { \mathcal L })=\displaystyle \frac{P({ \mathcal L }| \langle { \mathcal L }\rangle )\cdot P(\langle { \mathcal L }\rangle )}{P({ \mathcal L })}.\end{eqnarray} \tag{ 3 }$$

The luminosity distribution of our RSG sample is taken as the prior $P(\langle { \mathcal L }\rangle )$. We obtain a similar result when using a flat prior. To construct the likelihood $P({ \mathcal L }| \langle { \mathcal L }\rangle )$, we use our full sample of M31 RSGs with luminosity available from ME16, i.e., all the RSGs in Figure 3 both below and above the PTF variability sensitivity limit. We bin the stars in log L with bin-width ≈0.1 (Figure 4), and compute the distribution of ${dm}=m-\langle m\rangle$ magnitudes using the light curves of all the stars in each bin. Our construction of relative fluxes in the PTF R-band assumed these relative fluxes are similar in the other filters where most progenitor images were taken. We scale the resulting distribution in the corresponding bin containing $\langle { \mathcal L }\rangle$ by a factor 2.5 and shift it by $\langle { \mathcal L }\rangle$ (cf. Equation (2)) to obtain the likelihood. The denominator in Equation (3) is just a normalization constant. With all these ingredients, we compute the posterior distribution $P(\langle { \mathcal L }\rangle | { \mathcal L })$ and then obtain the mean and standard deviation of this distribution.

The result is shown in Figure 10 by the red points. The diagonal dotted line in the figure shows the 1:1 relation between the random-phase ${ \mathcal L }$ and the mean $\langle { \mathcal L }\rangle$, and thus represents the typical assumption made when computing supernova progenitor masses from pre-explosion images.

Zoom In Zoom Out Reset image size

From this figure, we see that $\langle { \mathcal L }\rangle$ versus ${ \mathcal L }$ for our RSG sample is consistent with the 1:1 relation. The vertical bars of the red points represent standard deviations of the $\langle { \mathcal L }\rangle$ distributions, and as can be seen, they become larger with increasing luminosity, reflecting the fact that the observed variability in our sample is larger for higher-luminosity stars (Figure 4). We also plot the measured luminosities along with the measurement errors of the SN II-P progenitor stars given by SS15 with blue points based on the 1:1 relation. As is evident, the measurement errors (from the photometry, extinction, and bolometric corrections, cf. Equation (2)) of the random-phase luminosities of the progenitors are much larger than the standard deviations of the $\langle { \mathcal L }\rangle$ distributions. The overall shift in the masses of the progenitors accounting for the intrinsic variability (a mere 1–2 M_⊙) is not significant, given the rather large uncertainty on the applied corrections (ΔA_λ and ΔBC_λ) and photometric measurements (Δm_λ). The uncertainty due to variability (dm), however, is an inevitable physical effect that cannot be overcome by taking more precise measurements and thus will eventually become the dominant effect.

The analysis presented above rests on the assumption that our full sample of RSGs in M31 is representative of the variability that RSGs exhibit in the final years before core-collapse. If these properties change as an RSG enters the final nuclear burning stages, our results and conclusions could be influenced. Furthermore, we have assumed that the variability of RSGs in the PTF R-band is similar in the filters (largely I) in which the SN II-P progenitors were imaged. If RSGs were markedly more variable in the I-band, then our results could be affected. Nevertheless, there is no obvious argument for variability amplitude to be larger in I than in R.

While intrinsic variability of RSGs with respect to the luminosity estimates of SN II-P progenitors is shown here to not be significant, it could, however, play an important role in other physical aspects of RSGs of consequence to their initial mass estimates, in particular a possible enhancement of the mass-loss (Heger et al. 1997; Yoon & Cantiello 2010). A similar effect is known to take place in Miras (e.g., Bowen 1988) and in Cepheids (Neilson & Lester 2008), wherein variability due to radial pulsations and the associated lifting and lowering of extended atmospheric layers is regarded to be an important driver for mass-loss. Multi-epoch monitoring, in multiple filters, for the variability of potential SN II-P progenitors, like the efforts of Gerke et al. (2015), Kochanek et al. (2017), and Johnson et al. (2017), will continue to be beneficial and have implications beyond the photometric uncertainty of the progenitor stars.

The Palomar Transient Factory's nearly 2000 days of monitoring of M31 has allowed for the first measurements of the prevalence of RSG variability in this important nearby galaxy. In this first study, we confined our efforts to the 255 bright RSGs that were cataloged and typed by the LGGS group (ME16). As was known for our galaxy, the Magellanic Clouds, and M33, these evolved massive stars can be highly variable, and often periodic. Based on the rms deviations of the PTF R-band light curves, 167 of the 255 RSGs exhibit variability above the noise limit of PTF data (Figure 3). For these sources, we see a broad positive correlation between the variability amplitudes, represented by the rms values, and their brightness.

We found that the light curves of these variables are characterized by semi-regular (and irregular for some cases) long-period pulsations over a strong red-noise (Kiss et al. 2006). Using a robust, probabilistically motivated methodology for reconstruction of their power spectra, we determined significant pulsational periods for 63 of the variable RSGs. Using the absolute magnitudes measured in the K-band by ME16, we derived the P–M_K relation for the RSGs in M31 and showed that it is consistent with that of the Galaxy, the Magellanic Clouds, and the M33 RSG samples, which represent a wide range of metallicities. This points toward a universality of this relation, at least for the Local Universe. It will be interesting to explore this relation with a large sample of RSGs from a larger number of distant resolved galaxies with the upcoming large telescopes like the Large Synoptic Survey Telescope (LSST; Ivezic et al. 2008) and the European Extremely Large Telescope (Greggio et al. 2012).

Comparison of the measured period–luminosity relation with the theoretical one based on the MESA models shows that the pulsations are consistent with the fundamental radial mode. At the same time, a clear sequence of outliers is seen in the distribution with our observed sample, which is found to be consistent with the first-overtone mode based on the models—another novel feat of this work.

Finally, we explored whether the measured RSG variability could be large enough to confound the inference of stellar mass from the observed supernova progenitor luminosity. We found that this effect is small compared to current observational uncertainties.

Future insights on RSG variability in M31 with the existing PTF data could be gained from large-scale variability studies of bright stars not necessarily in the ME16 catalog. In addition, with the recent commissioning of the ZTF instrument, continued monitoring of all RSGs will be possible and hopefully reveal longer period pulsations than possible with the original data set. Moreover, variability observations of different RSG populations spanning near and far host galaxies will be an asset for conducting many important studies of these stars. In the near future, LSST, with its superb photometry and expected cadence of three days, has the potential to generate such a product from monitoring a select group of galaxies.

This work was supported by the National Science Foundation through grants PHY 11-25915, PHY 17-148958, and ACI 16-63688, and funded in part by the Gordon and Betty Moore Foundation through Grant GBMF5076. M.S. is grateful to Niels Oppermann for the detailed discussions on statistical approaches and methods of proper signal reconstruction. M.S. also thanks Mathew Graham and Ashish Mahabal for helpful discussions on the period-finding algorithms, and Matteo Cantiello for the helpful suggestions and comments on the manuscript. E.B. thanks Sang-Hyun Chun for providing MESA inlists necessary to reproduce their RSG models, and Matteo Cantiello for discussion of RSG evolution in MESA. PyFITS is a product of the Space Telescope Science Institute, which is operated by AURA for NASA. Support for this work was provided to M.R.D. by NASA through Hubble Fellowship grant NSG-HF2-51373 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. M.R.D. acknowledges support from the Dunlap Institute at the University of Toronto. M.G. acknowledges the hospitality of the Kazan Federal University (KFU) and support by the Russian Government Program of Competitive Growth of KFU.

Software: numpy (van der Walt et al. 2011), scipy (Jones et al. 2001), astropy (Astropy Collaboration et al. 2013), pyfits, pandas (McKinney 2010), mpi4py (Dalcín et al. 2005), matplotlib (Hunter 2007), DAOPHOT (Stetson 1987), DAOGROW (Stetson 1990), NIFTy (Selig et al. 2013).