The Period-Luminosity Relations of Red Supergiants in M33 and M31

Based on previously selected preliminary samples of Red Supergiants (RSGs) in M33 and M31, the foreground stars and luminous Asymptotic Giant Branch stars (AGBs) are further excluded, which leads to the samples of 717 RSGs in M33 and 420 RSGs in M31. With the time-series data from the iPTF survey spanning nearly 2000 days, the period and amplitude of RSGs are analyzed. According to the lightcurves characteristics, they are classified into four categories in which 84 and 56 objects in M33 and M31 respectively are semi-regular variables. For these semi-regular variables, the pulsation mode is identified by comparing with the theoretical model, which yielded 19 (7) sources in the first overtone mode in M33 (M31), and the other 65 (49) RSGs in M33 (M31) in the fundamental mode. The period-luminosity (P-L) relation is analyzed for the RSGs in the fundamental mode. It is found the P-L relation is tight in the infrared, i.e. the 2MASS $JHK_{\rm S}$ bands and the short-wavelength bands of Spitzer. Meanwhile, the inhomogeneous extinction causes the P-L relation scattering in the $V$ band, and the dust emission causes the less tight P-L relation in the Spitzer/[8.0] and [24] bands. The derived P-L relations in the 2MASS/$K_{\rm S}$ band are in agreement with those of RSGs in SMC, LMC and the Milky Way within the uncertainty range. It is found that the number ratio of RSGs pulsating in the fundamental mode to the first overtone mode increases with metallicity.


INTRODUCTION
Red Supergiants (RSGs) are Population I stars in helium burning stage with a mass range of about 9-27M . They have relatively low effective temperature of ∼ 3000 − 4000 K and the corresponding spectral types of late K-to M-types (Massey et al. 2008). Their radius is large, with the maximum radius being ∼ 1500R (Levesque et al. 2005), so that they have low surface gravity of log g < 1.0, but high luminosity of ∼ 3 500 − 630 000L (Massey et al. 2008;Massey & Evans 2016).
RSG is a phase of very significant mass loss (Stothers & Chin 1996). The Reimers law (Reimers 1975;Kudritzki & Reimers 1978) is an empirical relation between mass loss rate (MLR), stellar luminosity, radius and mass derived from a sample of red giants and RSGs:Ṁ = 5.5 * 10 −13 L R /M . van Loon et al. (2005) also derived the relationship between MLR and stellar parameters, but based on the samples of RSGs and oxygen-rich Asymptotic Giant Branch stars (AGBs) in Large Magellanic Cloud (LMC): log(Ṁ ) = −5.5 + 1.05 log(L/10000L ) − 6.3 log(T eff /3500K). However, Mauron & Josselin (2011) pointed out that the result might be biased toward larger values ofṀ because the samples contained mainly extremely dusty RSG stars and the limitations of the mid-IR data available then. Circumstellar dust forms with stellar mass loss and most RSGs have some amount of circumstellar dust (Verhoelst et al. 2009). These dust have obvious infrared emission, especially in the midinfrared band, which makes the luminosity of such RSGs in infrared band exceed the upper limit of the theoretical luminosity of RSGs, and the contribution can be observed in the Spitzer 8µm and 24µm band.

3
RSGs are critical and important as both direct and indirect progenitors of supernovae that spend some time in the RSG phase. Lower mass RSGs are known to be the direct progenitors of corecollapse Type II-P supernovae with clear hydrogen lines in optical spectrum and a distinctive plateau in visual lightcurves. A certain number of Type II-P supernovae have been proved to associate with RSG progenitors by observation (Smartt 2015;Maund 2017). Higher mass RSGs, however, first noticed by Smartt et al. (2009), have not been detected as supernova progenitors. It may indicate that higher mass RSGs are not the direct progenitors of supernova, instead, they evolve back across the H-R diagram, spending a short time as yellow supergiants, blue supergiants or Wolf-Rayet stars before exploding as a supernova (Ekström et al. 2012). But high mass RSGs play an important role in their final fate because the mass-loss processes and dust production in the RSG phase will influence their ending.
Large amount of dust is detected in high redshift galaxies when the low-mass stars have not evolved to the AGB phase to become the dust contributor (Riechers et al. 2013). Therefore, the massive stars are the main if not unique providers of dust then. RSGs and supernovae (SNe) must be the major contributors to the interstellar dust in such distant galaxies.
Most RSGs show some degree of variability in the visual band. The typical amplitude of variation in the V band is about 1 mag, while in the near-infrared the amplitude is smaller, about 0.25 mag in the K band, and even smaller in the mid-infrared band (Levesque et al. 2007;Yang et al. 2018).
According to the characteristics of light variation, the RSG variability is divided into two categories: one is irregular variation which is too complex to be delineated by any period; the other is semiregular variation that can be further divided into two subclasses: one with short period (≤ several hundred days) and the other with Long Secondary Period (LSP) (≥ thousand days). The main difficulty in the study of variability of RSGs is the long time scale of variation. The time-series data need to span over several hundreds of days to cover at least one entire period of variation. If the LSP is present and desired to be determined, an even longer time-baseline is required to last for several years.
Like other types of pulsating variables, RSGs exhibit the period-luminosity (P-L) relation. Both theory and observation prove the existence of the P-L relation in RSGs with semi-regular variation, especially the RSGs with short period variation. Guo & Li (2002) derived the theoretical P-L relations of fundamental (FU), first overtone (1O) and second overtone pulsation models for 15-30M RSGs.
In the Milky Way, based on the light curves collected by the American Association of Variable Star Observers (AAVSO) with a span of 60 years, Kiss et al. (2006) studied 18 red supergiant stars to obtain their P-L relation.
Based mainly on the All Sky Automated Survey (ASAS) (Pojmanski 2002) and MAssive Compact Halo Objects (MACHO) (Alcock et al. 1997) observation, Yang & Jiang (2011, 2012 analyzed the lightcurves of 112 RSGs in Small Magellanic Cloud (SMC) and 169 RSGs in LMC and determined the periods of those semi-regular variables. They derived the P-L relation based on 47 RSGs in LMC and 21 RSGs in SMC with short period semi-regular variation which was consistent with the result of theory (Guo & Li 2002). They also analyzed a sample of 40 RSGs in M33 and obtained their P-L relation in the 2MASS/K S band by using the period from Kinman et al. (1987).
An extension to the environments of various metallicity will examine its effect on the P-L relationship. In the Local Group, the metallicity has a wide range from sub-solar in SMC  Garnett et al. 1997) and then to super-solar in M31 (12 + log(O/H) = 9.00; Zaritsky et al. 1994). It is known that the metallicity affects the ratios of blue-to-red supergiants (B/R) and Wolf-Rayet stars to RSGs (W-R/RSG). The B/R ratio changes by about 7 times (Maeder et al. 1980) and the W-R/RSG ratio changes by about 100 times over 0.9 dex in metallicity (Massey 2002). Unfortunately, differences in W-R/RSG ratio cannot be fully explained by stellar evolutionary models of massive stars. Besides, metallicity has an obvious impact on color indices (or spectral types) and effective temperatures of RSGs. The spectral types of RSGs shift toward earlier at lower metallicities (Elias et al. 1985;Massey & Olsen 2003;Levesque & Massey 2012;Dorda et al. 2016).
The effective temperatures of RSGs shift to warmer at lower metallicities, which was demonstrated 7 of red giants is truly small, on the scale of about one percent. Besides, foreground stars from TCD method are analysed to check the agreement with GAIA distances: 430 of 1677 foreground stars in M33 and 1161 of 3589 foreground stars in M31 are confirmed by the GAIA distances. Removing the foreground red giants from the TCD-selected samples, we are left with the samples of 747 RSG candidates in M33 and 432 RSG candidates in M31 consequently.

The AGB stars
Once the foreground stars are removed, the only pollution is the AGB stars in M33 and M31. On the HR diagram, the high luminosity AGB stars have a definite overlap with the low luminosity RSGs (Brunish et al. 1986). We solve this problem by setting the luminosity threshold in multiple bands. According to the mass-luminosity relation of massive stars: where γ ≈ 4 (Stothers & Leung 1971), the RSG luminosity corresponding to the mass range of 9-  Bessell et al. 1998a,b) where the difference between the K and K S bands is very 8 small and ignored in this work. The linear relationship of the apparent magnitude with the apparent bolometric magnitude is derived, which is used to convert the limits of bolometric magnitudes to that in the corresponding band. The results of linear fitting between bolometric magnitude and that in the λ band are shown in Fig.3 and 4, e.g. the upper and lower limits in the J band is derived to be 13.54 mag and 18.73 mag which correspond to the limits of apparent bolometric magnitudes of m max bol = 15.09 and m min bol = 19.86 for M33. The case of M31 is presented in Table 2 and Fig.4.
It is worthy to note that the linear fitting has relatively large dispersion in the Spitzer 8µm and 24µm bands. This can be understood by the contribution of circumstellar dust to the flux at 8µm and 24µm in addition to the stellar photosphere. This scattering leads to relatively uncertain limits of brightness in the 8µm and 24µm bands.
Once the range of the brightness in some bands, mainly the lower limit of the brightness, is set up, the locations of the RSG stars can be delineated in the CMD with the help of the color index and the AGB stars can be separated. Several near-and mid-infrared CMDs are shown in Fig.5 and 2. The minimum luminosity of all sources is about one magnitude above the theoretical limit of RSGs corresponding to 9M . The reason may be the selection effect of observation that the fainter RSGs are not detected due to the sensitivity limit.
3. In the CMD of K S vs J − K S , the red limit is at J − K S = 1.60, which is the limit of the carbonrich and oxygen-rich stars from Hughes & Wood (1990). The blue limit is at J − K S = 0.50, which is the observational limit of the RSGs in the Galactic halo obtained by Josselin et al. (2000). of the luminosity. Once the red limit is set, the blue limit is set to ensure the probability that lies in the range between the limits of blue and red colors is 95% based on the distribution. to the impact of circumstellar dust, we only mark by experience the range of these three color index.
For M31, 22 sources are bluer than the blue limit of J − K S or J − [3.6], however, 11 of them have been identified as K-type or M-type supergiant stars by spectroscopy (Massey & Evans 2016) so that they are kept, and 1 of these sources is also fainter than the lower limit of Spitzer/IRAC [3.6] band magnitude. Besides, 1 source is fainter than the lower limit of the Spitzer/MIPS [24] band magnitude and removed. In total, 12 sources are removed.
We are left with the samples of 717 RSGs in M33 and 420 RSGs in M31 at last.

The iPTF data
The time-series data is taken from the iPTF survey, a wide-field survey for an exploration of optical transients using the 1.2-meter Samuel Oschin Telescope at Palomar Observatory (Law et al. 2009;Rau et al. 2009). This survey uses a large field camera covering 7.8 square degrees with 11 2048 × 4096 CCDs at a resolution of 1.01 arcsec/pixel. In the iPTF survey, 2 and 6 frames covered the M33 and M31 sky area respectively from 2009 Aug to 2015 Jan (about 2000 days in total) with a typical seeing of 2 arcseconds. Although the images were taken in both R and g band, the R band observations are dominant and make use of >80% time, reaching 20.5 mag at a 5-sigma level.
Therefore we choose the R band iPTF data to analyze the variation of RSGs. The RSGs are selected in the iPTF images within three arcseconds from the position of the star. Fig.7 and Fig.8 show the iPTF images that include 658 and 377 RSG stars in M33 and M31 respectively, whose light variations are further analyzed.

Photometry
The CCD images were processed in a standard way. We measured all the iPTF images available (796 images of M33, 1980 images of M31) using the Sextractor code (Bertin & Arnouts 1996) with fixed coordinates. For a precise photometry, a stable reference star should be chosen for differential photometry, and its brightness and color should be similar to the target stars in order to cancel 11 the influence of the instruments and atmosphere as much as possible. For this purpose, we made an artificial reference star whose brightness is the average of all the red supergiant stars in each image, i.e. the average of approximately 700 RSGs in M33 and 100 RSGs in M31 within one image.
Naturally, this artificial star has the very typical color of RSGs. The large scale of the sample should smooth out the variation of all red supergiant stars, which makes the artificial reference star stable.
The stability of the brightness is checked by comparing with a star from the Tycho reference catalogue (Hog et al. 1998) with a brightness of V T = 12.608 mag and a color of V T − BT = 1.704 similar to RSGs. As shown in Fig.9, the standard error of differential photometry between the artificial reference star and the Tycho reference star is about 0.04 mag, which implies the photometric accuracy to be better than 0.04 mag. Considering this Tycho reference star is much brighter than the RSGs, the real accuracy should be even better. In order to achieve a good photometry quality, the size of the aperture is adjusted according to the full width at half-maximum (FWHM) of the PSF in each image by R aperture = 1.5 + 1.2 × FWHM.

PERIOD DETERMINATION
Before we start to determine the period of light variation, the original photometric data from differential photometry are processed further. First the outlying points are removed by using the locally weight linear regression (LWLR), with its most common method known as locally estimated scatterplot smoothing (LOESS). The weight function used for LOESS is the Gaussian function: ω(x, z) = exp(−(x − z) 2 /2κ 2 ) where we adopted κ = 4.5 with eye-check after trying the value in the range of [0.01,10]. The points lying more than 3σ away from the locally fitting line are removed.
Then the moving-average filter was used to smooth the light curve, in which a window of 9 days is taken. Fig.10 (Stellingwerf 1978) method from PyAstronomy for the period determination. The parameters used in the methods are shown in Table 3 in details. The final selection of the period depends on the goodness of fitting the lightcurve. A period is assigned to an object with semi-regular variation only when it is detected by at least two methods whose difference in period is usually (> 70%) less than 40 days while up to 50-90 days in a few cases. The period is searched continuously in the power spectrum of the residual data after subtracting the light variation with the most prominent period until the local signal to noise ratio of the power spectrum is less than 4. Fig.11 shows an example of period determination, where both the PDM and GLS methods find a period at ∼ 347 days, and the WWZ method also confirms that this period is present in most of the observational time. Based on the period and the regularity of light variation, these RSGs are classified into five categories: 1. Semi-regular variables, marked as S, which have definitely a period and are used for the followup analysis of the P-L relations. There are 84 (56)

PERIOD-LUMINOSITY RELATION
The P-L relation is derived for the bands with photometry data available, from the optical V band through 2MASS near-infrared JHK S bands to the Spitzer/IRAC and Spitzer/MIPS bands. As mentioned earlier, we only adopted the RSGs that show semi-regular variation and whose period are found by at least two methods.
The observational results are compared with the model by Guo & Li (2002) in Fig.14 and Fig.15.
The absolute magnitude in the K S band of the model is obtained from the absolute bolometric luminosity with the help of the bolometric correction in the K S band by Josselin et al. (2000).
Although the models are not perfectly matched with the observation, two sequences are clear in that most of the stars follow the sequence of the fundamental-mode pulsation and a few stars are closer to the first overtone pulsation mode. The three regular variables all locate close to the first overtone model. Because the pulsations of the fundamental and first overtone modes follow different P-L relations, they must be separated to study this relation. We identify 19 stars in M33 and 7 stars in M31 pulsating in the first overtone mode, and their P-L relation is not analyzed due to too few stars. Nevertheless, there is no clear indication of a P-L relation for them. The P-L relation of the 65 RSGs in M33 and 49 RSGs in M31 that pulsate in the fundamental mode is then analyzed.
In principle, the brightness should be the mean magnitude when the P-L relation is dealt with.
However, the brightness in the V band is taken from the single measurement by LGGS. Since the amplitude of variation in the V band is ∼ 1 mag, this can cause some dispersion in the P-L relation in the V band. Although the brightness in the infrared bands is also taken from the single measurement, the small amplitude of variation makes the single measurement approximate to the mean brightness and do not bring significant scatter in the P-L relation.
Following the convention, a linear function is taken to fit the relation between the absolute magnitude and log P. The results are displayed in Fig.16, Fig.18 and Table 4 for the M33 sample, and   15 in Fig.17, Fig.19 and Table 5 for the M31 sample. It is apparent that there is a tight relation in the near infrared and the mid-infrared band. Quantitatively the correlation coefficients are around 0.8 (see Table 4 and   et al. 2014;Clayton et al. 2015), and it is not well measured in M33. We take the value of 3.1 for simple and reasonably good approximation. Fig.18 and 19 illustrate the relation of W BV with log P.
The relation is much improved that the correlation coefficient is increased to around 0.4 (see the last row of Table 4 and Table 5), which verifies that the lack of P-L relation in the V band is partly caused by interstellar extinction. For a correct P-L relation in the V band, an accurate correction of interstellar extinction should be made to each star individually and accurately. In summary, the P-L relation exists between the photosphere brightness and the period of light variation for RSGs with semi-regular variation in M33 and M31. There is a tight P-L relation for apparent brightness in the 2MASS/JHK S and Spitzer/IRAC1-3 bands. Meanwhile, the P-L relation should exist and can show up in the V band when the interstellar extinction is corrected, and it will become tighter in the  (1987): M K = (−3.59 ± 0.41) × logP + (−0.88 ± 0.69). The P-L relation is in agreement with this work within a reasonable error range as shown in Fig.20. However, we derived the P-L relation based on a larger sample of 84 sources with semi-regular variation and longer time baseline. Moreover, the P-L relation is derived in many more bands, including the visual V band, other 2MASS bands and the Spitzer bands.

M31
The P-L relation of RSGs in M31 was obtained previously in the K band by Soraisam et al. (2018) who also used the iPTF data. Their data spanned from MJD 56000 to 58000. Here we used the iPTF data from MJD 55050 to 57050 by the limit to data access. Soraisam et al. (2018) compared with the theoretical P-L relation from Modules for Experiments in Stellar Astrophysics (MESA) (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018. Consistent with our results, the majority of sources pulsate in the fundamental mode. After ignoring 10 RSGs in the first overtone pulsation mode and 2 RSGs with periods < 100 days, they obtained the P-L relation: M K = (−3.38 ± 0.27) × logP + (−1.32 ± 0.75).
This P-L relation is also in agreement with our work within the error range as shown in Fig.20 Fig.21. It can be seen that the P-L relation in M33 is very similar to SMC. The P-L relations seem to be similar in these galaxies except for the Milky Way. The P-L relation of RSGs in the Milky Way appears a little different that may be caused by the dispersion of estimated distance and interstellar extinction. It may be argued that the P-L relations in K band is universal in these nearby galaxies. It should be mentioned that the optical bands would be better to examine the dependence of the P-L relation on metallicity because the P-L relation is less dependent on metallicity in near-infrared bands than optical bands (Bhardwaj et al. 2017). Nevertheless, interstellar extinction is serious to make the determination of P-L relation more uncertain in optical bands.
is taken into account, the P-L relation shows up although it is not as tight as in the infrared bands.

19
The results are compared with the P-L relations of RSGs in other galaxies, i.e. SMC, LMC and the Milky Way. It is found that the P-L relation in K band in these galaxies are similar, showing no apparent dependence on the metallicity.
From the P-L relation, the pulsation mode of RSGs in the M33 and M31 is compared with the theoretical model and found to be very possibly in the fundamental pulsation mode. When compared with the SMC, LMC and Galaxy, it seems that the pulsation mode shows dependence on the metallicity. This phenomenon may be related to convection and its dependence on metallicity and deserves further investigation.
Due to the limitations of the iPTF data available at now, the LSP of RSGs cannot be determined with reasonable accuracy. We may obtain long baseline data from the combination of the iPTF data and the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) DR2 data to determine the LSP in future.
We thank Prof. Yan Li for helpful discussions. Special thanks go to the anonymous referee for his/her very helpful suggestions, which improved the paper significantly. This work is supported by           (right). The grey shaded area represents the instability strip (Tammann et al. 2003). Red dots are the regular variables, and black dots with ID from Table 6 and 7 are the RSGs in the instability strip.     The symbols are the same as in Fig.16 39 Figure 19. The same as Fig. 18, but for M31.   .          Table 6 continued on next page          Table 6 continued on next page                                            Table 6 is published in its entirety in the machine-readable format online.          Table 7 continued on next page  Table 7 continued on next page  Table 7 continued on next page    Table 7 continued on next page    Table 7 continued on next page             Table 7 is published in its entirety in the machine-readable format online.