Red Supergiants in M31 and M33. I. The Complete Sample

The aim of this paper is to establish a complete sample of red supergiants (RSGs) in M31 and M33. The member stars of the two galaxies are selected from the near-infrared (NIR) point sources after removing the foreground dwarfs from their obvious branch in the J − H/H − K diagram with the archival photometric data taken by the UKIRT/WFCAM. This separation by NIR colors of dwarfs from giants is confirmed by the optical/infrared color–color diagrams (r − z/z − H and B − V/V − R) and the Gaia measurement of parallax and proper motion. The RSGs are then identified by their outstanding location in the members’ J − K/K diagram due to high luminosity and low effective temperature. The resultant sample has 5498 and 3055 RSGs in M31 and M33 respectively, which should be complete because the lower limiting K magnitude of RSGs in both cases is brighter than the complete magnitude of the UKIRT photometry. Analysis of the control fields finds that the pollution rate in the RSG sample is less than 1%. The by-product is the complete sample of oxygen-rich asymptotic giant branch stars (AGBs), carbon-rich AGBs, thermally pulsing AGBs, and extreme AGBs. In addition, the tip-RGB is determined together with its implication on the distance modulus to M31 and M33.


Introduction
The red supergiants (RSGs) are Population I massive stars in the core-helium burning stage. It is generally believed that the initial mass of an RSG is at least ∼8M e . But the lower limit of the initial mass for the RSG population may be as low as 7M e or even 6M e . The radius of the RSGs can reach ∼1500R e (Levesque et al. 2005), and they have low surface gravity and high luminosity of 3500−630,000L e (Massey et al. 2008;Massey & Evans 2016). A complete catalog of RSGs is the basis for studying the properties of RSGs more accurately, such as to examine massive-star evolution as a function of metallicity (Maeder et al. 1980;Massey 2002Massey , 2013, to estimate the total contribution of dust by RSGs to interstellar dust (Reimers 1975;Kudritzki & Reimers 1978;Gordon et al. 2016), and to calibrate the periodluminosity (P-L) relations of RSGs (Kiss et al. 2006;Yang & Jiang 2011Soraisam et al. 2018;Chatys et al. 2019;Ren et al. 2019) and the scaling relations between granulation and stellar parameters (Ren & Jiang 2020).
The Small Magellanic Cloud (SMC), Large Magellanic Cloud (LMC), Triangulum Galaxy (M33), and Andromeda galaxy (M31) are all near enough that RSGs can be detected and resolved individually, providing important cases to learn the statistical properties of RSGs in a galaxy. There have been some collections of RSG samples in these galaxies. The sample of RSGs was on the scale of a few tens of objects in the early studies of the LMC and SMC (Feast et al. 1980;Catchpole & Feast 1981;Wood et al. 1983;Pierce et al. 2000), and increased to a couple of hundred later (Massey 2002;Massey & Olsen 2003;Yang & Jiang 2011Neugent et al. 2012;González-Fernández et al. 2015). Recently, Yang et al. (2019Yang et al. ( , 2020a identified 1405 and 2974 RSGs in the SMC and LMC, respectively, which is a drastic increase from previous studies and estimated to be about 90% complete. This revolutionary progress comes from both much more collections of data and the method to remove the foreground stars. Specifically, they combine a variety of color-magnitude diagrams (CMDs) to identify RSGs and remove the foreground contamination by Gaia's proper motion and parallax.
M31 and M33 are much more distant than the MCs with distance moduli larger by about 5 mag, which calls for alternative methods to identify RSGs. Previously, Massey et al. (2006Massey et al. ( , 2007Massey et al. ( , 2009, Drout et al. (2012) selected the initial RSG sample by V<20 and V−R0.85 for M31 and V−R0.6 for M33 through the Local Group Galaxies Survey (LGGS) observation. The V-band criterion was set to ensure sufficient brightness to avoid confusion with asymptotic giant branch stars (AGBs), and the color index was set to limit the stars to K and later types. The foreground dwarfs were further removed from this initial sample by the B−V/V−R diagram (Massey et al. 2009;Drout et al. 2012). Consequently, they identified 437 and 776 RSG candidates in M31 and M33, respectively. Massey & Evans (2016) measured radial velocities and determined spectral types for 255 (about 60%) of these stars and confirmed they are truly RSGs after comparing their radial velocities with the expected values of the Population I objects in M31 (Massey et al. 2009).
The RSG samples in M31 and M33 identified by Massey et al. (2009) are far from complete. Ren et al. (2019) found that the minimum luminosity of RSGs in these samples is about 1 mag above the theoretical limit of RSGs for a 9 M e star, which indicates that the samples missed the faint RSGs. In addition, Messineo & Brown (2019) discovered 889 RSG candidates in the Milky Way galaxy from Gaia Data Release 2 (DR2). Considering that both M31 and M33 are spiral galaxies like the Milky Way, the number of RSGs should be comparable. Moreover, the sample of 889 RSGs in the Milky Way cannot be complete. In comparison with the newly found large sample (a couple of thousands) of RSGs in the SMC and LMC, M31 and M33 with a much larger geometrical size are expected to host many more RSGs. As mentioned in Massey & Evans (2016), their completeness limit was set by V20, which implies that the RSG population of M31 was only complete down to ∼15M e , corresponding to L L log 4.7 according to the mass-luminosity relation of massive stars Stothers & Leung 1971). If we set out to complete the identification of RSGs in M31 and M33 down to ∼7M e , i.e., the bolometric magnitudes M bol of ∼−3.71 mag, the V-band magnitudes will be 22.63 and 22.89 for M31 and M33, respectively, by taking the distance modulus of 24.40 (Perina et al. 2009) and 24.66 (Orosz et al. 2007), the bolometric correction BC V of −0.94 for a 4000K RSG (Massey et al. 2009), and the extinction A V =1 (Massey et al. 2009). With the increase of the photometric error to the faint objects, it becomes more and more difficult to distinguish the foreground dwarfs from the member RSGs by using the B−V/V−R diagram. Thus, it is hard to yield a complete and pure identification of faint, cool RSGs in M31 and M33 by using optical data like LGGS.
In this work, we try to establish a complete sample of RSGs in M31 and M33 using near-infrared data in a new way. The paper is organized as follows: Section 2 for the data, Section 3 for the method to remove the foreground stars, and Section 4 for how to identify RSGs.

Data and Reduction
The JHK brightness comes from images taken with the Wide Field Camera (WFCAM) from mid-2005 to 2008 on the 3.8 m United Kingdom Infra-Red Telescope (UKIRT) located in Hawaii (Irwin 2013). WFCAM consists of four Rockwell Hawaii-II (HgCdTe 2048 × 2048) detectors, each covering ¢ 13 .65 on sky with 0 4/pixel. For some exposures with microstepping, which is used to recover some of the lost resolution when observing conditions are undersampled, the microstepped frames are interwoven to give an effective sampling of 0 2/pixel in the 2×2 microstep mode or 0 133/pixel in the 3×3 microstep mode. 4 For the images we used in this work, the average seeing on all frames varied between ∼0 7−1 2. The images were processed by the Cambridge Astronomical Survey Unit (CASU) and made available via the WFCAM Science Archive 5 (deprecated frames are purged in this work). The data products include the calibrated stacked images and the corresponding source catalogs. We further convert the FITS-format source tables into ASCII-format catalogs, applying all necessary corrections using the program provided by CASU. 6 The resultant ASCIIformat catalogs contain R.A., decl., magnitude, magnitude error, stellar classification flag, etc. We cross-match the results between the J,H, and K bands with a radius of 1″. The sources with the stellar classification flag of −1 (stellar), −2 (probably stellar), and −7 (source with bad pixels) are kept and regarded to be point sources. For the cases where the flags disagree in the three bands, at least two of the JHK bands must meet the above conditions to be selected. We add an "N_Flag" index to the JHK catalogs to indicate the number of bands in which the source is identified as a stellar source, i.e., 3 means all three bands are labeled "stellar." Finally, there are 1,245,930 and 203,486 point sources in M31 and M33, respectively. Different from other works, the source with a flag of −7 is kept in order to guarantee the completeness of the sample. For M31, the percentage of "−7" sources is about 5% in the J band, while it is up to about 40% in the H and K bands in the initial sample. The large fraction of "−7" sources in the HK band and difference with the J band are caused by the H-and K-band data used being all taken in the 3×3 microstep mode, while all of the J-band data are taken in the non-microstep mode. As a result, the percentage of "−7" sources increases several times in the H and K bands because the data are flagged as having bad pixels in the core radius. The case of M33 is similar. We checked the original source catalogs and found that most of the sources with a fractional number of bad pixels is due to soft-edged aperture. So, the photometry of "−7" sources is also considered reliable.
In order to examine the photometry accuracy of UKIRT, the J and K magnitudes from UKIRT are compared with those of the Two Micron All-Sky Survey (2MASS; Skrutskie et al. 2006). Because the UKIRT astrometry and photometry are calibrated with the 2MASS point-source catalog, the UKIRT photometric results agree very well with those of 2MASS as expected. In the case of M31, the 2MASS 6x point-source catalog can reach J∼20 mag and K∼18.5 mag, while the sources flagged "AAA" 7 reach J∼17.5 mag and K∼16 mag. The 1σ differences between 2MASS and UKIRT are ∼0.2 mag at J∼17.5 mag and K∼16.5 mag. In the case of M33, the main 2MASS point-source catalog can reach J∼17.5 mag and K∼16.5 mag, while the "AAA" sources to J∼16 mag and K∼15 mag with the 1σ difference better than 0.1 mag. These differences are well within the claimed uncertainty.
The distributions of the K-band magnitude and corresponding error are shown in Figure 1. In general, the photometry is better for M33 than M31 in both depth and accuracy because the star field in M33 is less crowded. If the magnitude that is 0.5 mag brighter than the drop-off point is considered to be complete, then sources brighter than 17.94 and 18.22 are complete in M31 and M33, respectively. The spatial locations of the selected stars are shown in Figure 2, where M32 and M110 are labeled as well by their location and size. It can be seen that the samples cover almost all the regions of M31 and M33, with a very small part of M33 missed. There are some additional fields beyond M31 and M33 which will be used as control fields to estimate the pollution rate of RSGs. The observed CMDs of the initial sample are shown in Figure 3, where the recently identified member stars are decoded by color.

Removing the Foreground Stars
Although M31 and M33 are not located in the Galactic plane, Massey et al. (2007) showed that the contamination by foreground stars is serious toward these sight lines. There are a few ways to remove foreground stars. Yang et al. (2019Yang et al. ( , 2020aYang et al. ( , 2020b separate efficiently the SMC and LMC members from foreground stars by using astrometric solution from Gaia/DR2 because the MC members concentrate on the values expected from the motion of MCs relative to the Galaxy. Unfortunately, this method cannot be applied to M31 and M33 effectively because they are so distant that the proper motions and parallaxes are too small to be measurable by Gaia. In the sample, only 6.7% (∼83,260/1,245,930) and 9.3% (∼18,901/203,486) of sources have Gaia parallaxes and proper motions measurements, many of which are unreliable. For M31 and M33, Massey & Evans (2016) removed foreground stars by radial velocities and spectral type from optical spectroscopy, which can identify RSGs and membership correctly but only for very bright sources. As mentioned earlier, the optical B−V/V−R diagram is used to remove foreground dwarfs. This two-color method is deeper than spectroscopy and works very well to distinguish foreground dwarfs from RSGs. But this method is useful only for highaccuracy photometry and limited to bright sources (see Massey et al. 2009;Drout et al. 2012;Massey & Evans 2016 and our later analysis).
Instead of optical observation, we rely mainly on the nearinfrared photometric data. Because the effective temperature of RSGs is in the range 3000-5000 K (Massey et al. 2008;Neugent et al. 2010;Yang et al. 2019), their major radiation goes around near-infrared so that RSGs will be most easily detectable in the J, H, and K bands. For an RSG with color index V−K∼4.0, its K magnitude would be about 18 at V=22. The infrared band also has much less extinction than the optical (Wang & Chen 2019). We take advantage of the near-infrared bands to identify RSGs by using the UKIRT observation of M31 and M33. Besides, the Gaia DR2 data are used to remove the foreground giant stars, though practically no object is removed in this way. We remove the foreground dwarfs using the near-infrared color-color diagram. The study of Bessell & Brett (1988) found that the intrinsic color indexes have clear bifurcations on the J−H/H−K two-color diagram for giant and dwarfs. Dwarfs have higher surface gravity than giants or supergiants; the collision rate between atoms is higher, and molecules are easier to form (Allard & Hauschildt 1995). This makes molecules form at relatively high temperatures in dwarfs, causing absorption in the H band and darkening the H-band brightness and eventually leading to smaller J−H and bigger H−K than giants.
The borderline between dwarfs and giants are redefined intentionally and specifically, though Bessell & Brett (1988) already obtained the intrinsic color indexes. For this purpose, the high-accuracy photometric data are chosen with the error of J-, H-, and K-band photometry less than 0.05 mag and N_Flag=3, which means the object is identified as "stellar" in all three bands. The J−H/H−K diagrams for these accurate photometries are shown in Figure 4 for M31 and M33. It can be seen that there is a very clear boundary between giants and dwarfs. With the increase of H−K, J−H increases to about 0.7 and then turns down for dwarfs; meanwhile, the J−H of red giants and supergiants starts from about 0.7 and increases monotonically. The trends and values coincide very well with the result of Bessell & Brett (1988).
The dividing line is defined quantitatively. In a step of 0.01 in the range of H−K from 0.05 to 0.3, the point of the maximum surface density is calculated on the dwarf branch, and then the piecewise function is used to fit those points. When H−K0.13, the function is linear; for H−K>0.13, the function is quadratic. This piecewise function represents the relation of the J−H and H−K of dwarfs. Because the foreground extinction is pretty small (an average A V ∼0.17; Schlegel et al. 1998;Schlafly & Finkbeiner 2011;Cordiner et al. 2011), these colors basically represent the intrinsic color indexes of dwarfs. We shift the piecewise functions up by 0.12 and 0.09 mag by eye-check and take them as the dividing lines between giants and dwarfs in the M31 and M33 fields, where the difference in the shift in the vertical axis is caused by the different foreground extinction to the two galaxies, i.e., A V ∼ 0.17 mag and 0.11 mag for M31 and M33, respectively, according to SFD98 (Schlegel et al. 1998). The function forms and coefficients of the adopted dividing lines are listed in the first rows of Table 1.
With the dividing line determined, the criterion is applied to the entire initial sample to remove foreground dwarfs. In addition, we remove sources with H−K<0.1 that is apparently bluer than RSGs. This action certainly also removes the blue and yellow supergiants in the galaxy, which are absent in the final catalog of the members in M31 and M33. Finally, 414,490 and 77,091 foreground dwarfs are removed in the M31 and M33 fields, i.e., 33.3% and 37.9% of the initial samples.
The above method is expected to remove all of the foreground dwarfs in the sample, but the uncertainty of color indexes would move some sources around the borderline. The completeness and pollution rate of the selected member stars are estimated by Monte Carlo simulation. First, the sources with "N_Flag=3" and σ J,H,K <0.05 mag are taken as a no-error "perfect" sample so that their locations in the J−H/H−K diagram can absolutely decide being a dwarf or a giant. Then, we perform 5000 simulations for a random error with a two-dimensional Gaussian distribution for each source whose width is four times the error of J−H and H−K. Finally, the UKIRT/NIR criterion is used to divide dwarfs and giants to compute the completeness and pollution rate of giants when the JHK photometric errors are limited to less than 0.2 mag. For M31 and M33, the simulation results show that the completeness of the selected stars is about 93%, and the pollution rate is about 9%.
Actually, the pollution rate calculated by this method will be overestimated, and the completeness will be underestimated. On one hand, we multiply the errors of the sources whose JHK photometric errors are less than 0.05 mag by 4 to simulate the case where the JHK photometric errors are limited to be less than 0.2 mag, which means increasing the photometric errors systemically. On other hand, the distribution of sources with JHK photometric errors less than 0.05 mag is already scattering  due to the photometric errors. Therefore, the true pollution rate of the selected giants is smaller than the simulation value, and the completeness is larger than the simulation value.

Double-check by Optical/Infrared Color-Color Diagrams
Although the J−H/H−K diagram works very efficiently and is applicable to all the sample stars, the identification deserves to be checked by other methods. One purpose is to confirm the identification, and the other is to further remove some foreground stars that are very close to the borderline in the NIR color-color diagrams (CCD) but may be significantly distinguishable in optical bands in particular for some relatively blue stars.
1. The r−z/z−H diagram. By convoluting the spectrum of the MARCS model with the transmission functions of the r,z, andH filters, Yang et al. (2020c) found that the r−z/z−H diagram can distinguish dwarfs from giants well. We introduce the r−z/z−H diagram only as an auxiliary method to remove foreground dwarfs, because the information used to distinguish dwarfs from giants here, H-band photometry, has already been used in the J−H/H−K CCD (see Section 3.1.1), and only those foreground dwarfs with good photometric quality will be taken into account. The PS1/DR2 data are used, where the forced mean PSF magnitude is taken for its consideration of both photometric depth and accuracy. The data flags of PS1/ DR2 are complex and only those with good measurements are used. The cross-match between UKIRT and PS1 by a radius of 1″ resulted in 184,750 and 37,304 stars in M31 and M33 with good photometry. Similar to the method in Section 3.1.1, we select the "good-measurement" sources with r-, z-, and H-band photometric errors less than 0.05 mag to determine the borderline between dwarfs and giants in the r−z/z−H diagram. The positions of the maximum surface density of the dwarf   branch on r−z/z−H diagram are calculated in a step of 0.01 of r−z from 0 to 2.5, and then a piecewise function is used to fit those points. When r−z1.1, the sigmoid function, which is also known as the logistic function, is used for fitting; otherwise, the function is a quadratic curve. As shown in Figure 5, the piecewise functions are shifted up by 0.12 and 0.10 mag for M31 and M33 to become the dividing lines between giants and dwarfs. It should be noted that a slight difference in the intrinsic color indexes r−z or z−H appears between our results and the MARC model though they are in general agreement. The functions used are listed in the second rows of Table 1. The criterion is applied to the PS1 sources with "good measurement," and the photometric error in the r, z, and H bands is less than 0.1 mag. The sources below the dividing line or with r−z<0.3 are removed. This removes 102,174 and 17,648 dwarfs in the M31 and M33 fields, respectively, in which 93,938 (91.9%) and 15,568 (88.2%) are also removed by the NIR CCD. In other words, an additional 8,236 and 2,080 stars are removed, which are mostly around the borderline in the NIR CCD or in the area close to the center of the galaxy where the photometry is of poor quality. 2. The B−V/V−R diagram. Massey (1998)  The logistic function is taken to fit these points and shifted up by 0.12 and 0.10 mag for M31 and M33, respectively, respectively, as the dividing lines shown in Figure 6. The function forms and coefficients of dividing lines are listed in the third rows of Table 1. The LGGS criteria work very well for bright sources when Massey et al. (2009) and Massey & Evans (2016) limit the sources by V brighter than 20 mag, i.e., photometric error less than 0.01. But for faint sources with a slightly large photometric error, to distinguish dwarfs from giants is very hard. As shown in the left panel of Figure 7, a large portion of the foreground dwarfs selected by the NIR CCD are located in the B−V/V−R region of dwarfs, confirming the consistency between optical and NIR criteria. But when the photometric error increases to 0.05 mag, many of the NIR-selected foreground dwarfs fall into the giants region of the LGGS diagram as shown in the right panel of Figure 7, needless to say that many of the objects have photometric uncertainty of ∼0.1 mag. This can be understood that the difference between giants and dwarfs in the B−V/V−R diagram is too small to tolerate any significant photometric error. Comparing Figure 6 with Figure 4 and Figure 5 demonstrates clearly that the difference of colors in the optical (mostly ∼0.1 mag) is much less significant than in near-infrared (mostly ∼0.2 mag). The cross-match of the LGGS catalog with UKIRT by a radius of 1″ results in 92,441 (7.4% of the initial NIR sample) and 45,279 (22.3%) associations. These associations are a small part of the sample, and taking the risk of mixing dwarfs and giants together into account, the LGGS criteria are not used to remove dwarfs in this work.

Removing the Foreground Giants
It is worth noting that our UKIRT and PS1 criteria can be used only to remove foreground dwarfs, but foreground giant stars cannot be removed in this way. Instead, we make use of the Gaia astrometric information to remove them because the foreground giants should present measurable motions.
Together with some foreground dwarfs, the foreground giants are searched for using parallaxes and proper motions from Gaia/ DR2 (Gaia Collaboration et al. 2018). Stars are considered to be foreground objects if they satisfy either parallax or proper motion constraints. Specifically, we remove the sources whose distances are less than the Milky Way scale (i.e., 50 kpc; Liu et al. 2017) with astrometric solution relative error less than 20% (i.e., |σ ω /ω|, | | s m m a a * * and | | s m m d d are both smaller than 20%). The distances here are calculated with the Smith-Eichhorn correction method from Gaia-measured parallax and its error (Smith & Eichhorn 1996). Besides, a source is also removed if the measured proper motion is greater than that expected for a star with a velocity of  where μ α* and μ δ are proper motions in R.A. and decl.. The only sources removed by Gaia criteria are expected to be foreground giants. Among the 83,260 (M31) and 18,901 (M33) Gaia-UKIRT cross-identified stars, 14,795 and 2,433 stars are removed by the Gaia criteria as well as by UKIRT. There are 15 and 5 sources identified as foreground stars only by Gaia. However, the stars removed only by Gaia are faint in the K band, mostly fainter than 14 mag. A typical red giant star is as bright as -4 mag in K, and 14 mag means a distance of 40 kpc. Thus, these sources are very unlikely to be foreground red giants. We tend to believe that the astrometric information of Gaia for most of these sources is uncertain. As these stars actually have color-magnitude similar to RSGs in the host galaxies, this uncertainty may be caused by the variation of the photocenter due to the large-scale convection of RSGs (Chiavassa et al. 2011). Therefore, the Gaia-only sources are not removed from the sample, which means no foreground red giants are removed from the sample. This result confirms the argument that the contamination of foreground red giants is very small toward the sight lines of M31 and M33 by Massey & Evans (2016).

Comparison with the Besançon Model
The Besançon Milky Way stellar population synthesis model (Robin et al. 2003) is introduced to examine whether the foreground stars are removed correctly. Within the distance range of 0-50 kpc, the expected foreground stars at the direction of M31 and M33 with a sky area of 12 deg 2 and 3 deg 2 , approximately the coverage of the UKIRT data, are computed by the Besançon model. Their distribution in the J−K/K diagram is displayed in Figures 8 and 9, where the red giants (red dots) are diagnosed by < g log 2.5 (from G5 later on) because the blue (super)giants with large g log are few and might be removed by our criteria. For the M31 area, the number of foreground giants and foreground dwarfs with Kband magnitude brighter than 20 mag and fainter than 12 mag is 55 and 214,157 respectively. For the M33 area, the numbers are 12 and 26,046, respectively. In comparison, our criterion removed 422,741 and 79,176 dwarfs in the sky area of M31 and M33, which is basically consistent with but still more than that predicted by the Besançon model. On the contrary, we removed no foreground red giants, which can be understood as the foreground red giants being so bright that they should be saturated in the UKIRT observation with a stellar classification flag of −9 (possibly saturated objects) and have been excluded from the initial sample.
As shown in Figures 8 and 9, the foreground dwarfs removed by the UKIRT and PS1 criteria are consistent with the Besançon model, but those faint sources in the K band have larger dispersions than the "no-error" Besançon model due to the effect of photometric error. Meanwhile, some stars removed appear on the red side with J−K>1.0, which is not present in the Besançon model. Re-examination found that these sources are removed by the r−z/z−H criterion, and they are mostly located close to the center of the galaxy whose photometry suffers relatively large uncertainty. Nevertheless, they have little effect on the RSG sample because they are too red at given brightness for RSGs, and they may be AGB or RGB stars.
The number of foreground stars removed by the three methods is summarized in Table 2. Apparently, the NIR colorcolor diagram works much more effectively than the other two methods because the distinction of giants from dwarfs is significant and RSGs are bright in the near-infrared.    Although this work focuses on the catalogs of RSGs in M31 and M33, other evolved populations can be identified using the CMDs of the member stars after removing the foreground stars, which helps to identify RSGs. For this purpose, we choose the sources with "N_Flag=3" and JHK photometric errors less than 0.1 mag to define the regions of various populations in the near-infrared CMDs shown in Figure 10, and then apply the criterion to classify other less accurately measured sources. In order to ensure that the objects are point sources, we drop those with "N_Flag=2" (i.e., only two of the JHK bands are marked as point sources) but marked as extended source by PS1. On this basis, we define the sources with "N_Flag=2" as "Rank 2" sources and sources with "N_Flag=3" as "Rank 1" sources.
The division between various populations in the CMD is mainly guided by previous studies of evolved populations in the SMC and LMC by Yang et al. (2019Yang et al. ( , 2020aYang et al. ( , 2020b. Additionally, the MIST (MESA Isochrones and Stellar Tracks) model (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015Dotter 2016;Choi et al. 2016) of a 5M e and 7M e star is referenced for M31 and M33 with [Fe/H]=0.3 and [Fe/H]=0.1, respectively. Above these empirical and theoretical results, the division is decisively determined by the density contour of the stars in the J−K/K diagram. First, the position of the Tip-RGB (TRGB) is determined, which will be discussed in Section 4.4, and sources fainter than TRGB are RGB stars. The sources brighter than TRGB are divided into two major groups: relatively bluer RSGs and relatively redder AGB stars. The former is the leading group to be discussed in detail in the following. AGB stars are subdivided into oxygen-rich AGBs (O-rich AGBs), carbon-rich AGBs (C-rich AGBs), extreme AGBs (X-AGBs), and thermally pulsing AGBs (TP-AGBs). The locations of O-rich, C-rich, and X-AGB stars in the J−K/K diagram have been previously identified by several works, such as Cioni et al. (2008) and Yang & Jiang (2011), and our results agree with theirs. But the location of TP-AGB stars is suggested for the first time. It can be seen from Figure 10 that the TP-AGB branch is almost parallel to and very close to but slightly redder than the RSG branch, which led to previously classifying them as RSGs. However, there is a clear gap from RSGs, and they coincide very well with the 5M e MIST TP-AGB model. We expect that spectroscopy and light variation of these objects (see Ren & Jiang 2020) would further confirm the nature of these sources.
The total numbers of different stellar populations are listed in Table 3. Also listed are the numbers of "Rank 1" sources in brackets, whose type flags are all "stellar" in the JHK bands, implying the lower limit of the number of sources of each type. As will be shown in Section 4.4, the TRGB of M31 and M33 is 17.62 mag and 18.11 mag in K band, which is brighter than the completeness magnitude (17.94 and 18.22). The samples of AGB and RSG stars should be more or less complete, while the samples of RGB stars that are incomplete for the fainter ones are not all detectable. M31 is similar to our Milky Way galaxy in type and metallicity and even size; these numbers should be a valuable reference for studying the evolved stellar populations in our Galaxy as a whole.

RSGs
The RSG branch is very obvious in the J−K/K diagram after removing foreground stars shown in Figure 3. The magnitude and color criteria can then be used to define the RSG region in this diagram. As pointed out in previous section, RSGs stand out in the J−K/K diagram above the TRGB. We take the K-band magnitude of TRGBs (17.62 mag for M31 and 18.11 mag for M33) as the lower limit of RSGs. Although this limit may be questioned, we have two reasons. One is that the configuration in the J−K/K diagram is continuous for the RSG branch until the TRGB where the break occurs at the lowest density. The other is that the core-helium-burning stage of the 7M e star in the MIST track starts from about this position, which indicates that the lower mass limit of RSG is about 7M e .
Here we take two straight lines as the red and blue boundaries enclosing the RSG branch. The quantitative forms of the blue and red boundaries are 20.00 33.00, 3 8.00 25.00. 4 Previous samples of RSGs in M31 and M33 were obtained from optical observations based on the B−V/V−R diagram (Massey et al. 2009;Drout et al. 2012). These photometrically classified RSGs were further checked by radial-velocity information for membership determination. Massey & Evans (2016, ME16 for short) and Drout et al. (2012, D+12 for short) identified 255 RSGs in M31 and 204 RSGs in M33 with 189 Rank 1 highly likely supergiants and 15 Rank 2 possible supergiants. Among them, 240 in M31 and 201 in M33 radialvelocity-confirmed RSGs are in our initial sample, of which 180 and 154 sources are considered as member stars, and finally 180 and 147 sources are also identified as RSGs in our work. These stars are labeled "D" in the column "LGGSType" of Table 4 and 5.
These previously identified RSGs in M31 (240) and M33 (201) are compared with the whole sample in Figure 10. It can be seen that previous samples miss the red and faint RSGs and include some objects too blue or too red to be an RSG. Note. a The numbers in the brackets include only the sources labeled as "stellar" in all of the UKIRT/JHK bands.

Completeness and Pureness
As mentioned in Section 2, sources brighter than 17.94 and 18.22 in the K band are complete. With the lower limits of RSGs in M31 and M33 being K=17.62 and 18.11, the samples are considered to be complete, i.e., there are 5498 and 3055 RSGs in M31 and M33, respectively. RSGs in M31 and M33 are listed in Table 4 and 5 with R.A. and decl. coordinates, magnitudes, magnitude errors, astrometric information, etc. One thing to note is that our preliminary selection of stars is not very strict in that the stars labeled "−7" for bad pixel are included and that only two of the JHK bands labeled "stellar" are required. If all three JHK bands are required to be "stellar," then the number of RSGs is reduced to 3268 and 2804 for M31 and M33, respectively. The number of reduction for M31 is significant and very minor for M33, again due to the much more crowded field in M31. We think these numerals must have underestimated the sample. Moreover, if the sources with two of the three bands labeled "−7" are removed, the number of RSGs is further reduced to 3154 and 2635, respectively.
Comparing the location of various types of evolved stars in the J−K/K diagram with the limiting magnitude of 2MASS, it can be inferred that some RSGs and AGBs are detectable by 2MASS, but not all of them. If the quality flag is further restricted to "AAA," then only bright AGBs and RSGs in M31 with K<16 can be recognized by the 2MASS photometry. For M33, the situation is reduced because the main 2MASS point-source catalog has a brighter limiting magnitude of K∼15 for the "AAA" sources.
The pollution rate of the RSG sample is estimated from the control field. First, the UKIRT fields are divided into geometric regions of M31, M32, M110, and M33 and the control fields. The region of a galaxy is demarcated by an ellipse whose major and minor axes are determined by the = B 25 mag arcsec 2 isophotes in the B band shown in Figure 11. The fields beyond the geometric regions of the galaxies are classified into control fields. The major axes, minor axes, and position angles of the galaxies are listed in Table 6. The sky area of galaxies is calculated by S=πab, where a and b are the half-major and half-minor axes. The sky area of the control fields is the total sky area minus the sky area of the galaxies. The sky area of different regions and the number of RSGs are listed in Table 6. If all of the RSGs in the control fields are foreground stars, or fake RSGs, the pollution rate can be estimated. First, the fake RSGs per sky area are calculated by where N galaxy is the number of RSGs within the galaxy's area. For M31, N galaxy =N M31 +N M32 +N M110 . The pollution rate turns out to be 1.30% and 0.49% for the catalog of RSGs in M31 and M33, respectively. Indeed, the fake RSGs in the control fields located around the rim of the galaxies as shown in Figure 11, which indicates some of them are actually the members of M31 or M33. In another word, M31 and M33 extend to a larger area than the labeled ellipse defined by 25 mag arcsec 2 isophotes. The true pollution rate should be smaller than the above derived values.

Spatial Distribution
The spatial distribution of the selected RSGs is shown in Figure 11 with the GALEX ultraviolet image as the background to display the massive-star regions. The locations coincide very well with the spiral arms, which are expected for RSGs as massive stars. This structure supports our identification of RSGs.

Density of RSGs as a Function of Metallicity
The number of RSGs is 5225 within 2.567 deg 2 of M31 and 3001 within 0.644 deg 2 of M33. Apparently, the surface density of RSGs is not the same. It is well known that metallicity influences the time a massive star spends in the RSG stage. In order to characterize the density and the massive-star formation rate, the number of RSGs is normalized to the stellar mass of the galaxy. In addition to M31 and M33, we supplement the other galaxies (SMC, LMC, and MW) whose RSG samples are systematically studied.  Gehrz (1989) predicted at least 5000 RSGs. Both values are displayed in Figure 12, and the former one  Figure 12. Variation of the number of RSGs per stellar mass with metallicity for five galaxies. For the Milky Way, the lower limit comes from the already identified number of RSGs by Messineo & Brown (2019,) while the higher value comes from the predicted number by Gehrz (1989). indicates the lower limit. The adopted stellar masses are 3.1×10 8 M e , 1.5×10 9 M e , 5.2×10 10 M e , 2.6×10 9 M e , and 1.0×10 11 M e for the SMC, LMC, MW, M31, and M33, respectively (Besla 2015;Fattahi et al. 2016 (Garnett et al. 1997), and 9.00 (Zaritsky et al. 1994) for the SMC, LMC, MW, M33, and M31, respectively.
The RSG density per 10 8 M e is presented in Figure 12. With increasing metallicity, the RSG density per stellar mass decreases rapidly. When the metallicity ( ) + 12 log O H increases by 0.9 dex, the number of RSGs per stellar mass decreases by about 60 times. This can be understood as metallicity affecting the time stars spent in different evolutionary stages. Previous studies have shown that metallicity affects the ratio of blue-to-red supergiants (B/R) and Wolf-Rayet stars to RSGs (W-R/RSGs). When metallicity increases by 0.9 dex, the B/R ratio and the W-R/RSG ratio increase by about 7 times (Maeder et al. 1980) and 100 times (Massey 2002) respectively. Our conclusion that the density of RSGs decreases with increasing metallicity is consistent with this scenario. On the other hand, metallicity cannot be the only factor to influence the density of RSGs. M33 is much more metal-rich than LMC, but holds very similar RSG density. In fact, M33 is the only SAcd-type galaxy in them; the strong starforming activity may account for the higher massive-star formation rate. Among these galaxies, M31 should be the most similar to our galaxy in both type and metallicity. Taking M31 as the reference, the number of RSGs in our galaxy should be ∼2800 after scaling by stellar mass. But this is only the lower limit of the number of RSGs in our galaxy. Because the metallicity of the Milky Way is lower than that of M31, the number of RSGs per stellar mass should be higher than that of M31. As shown in Figure 12, the Gehrz (1989) prediction of the number of 5000 RSGs in our galaxy agrees well with the overall trend of the number of RSGs per stellar mass with metallicity, but is still lower than the overall trend. A detailed study of this problem will be presented in our future work.

Tip-RGB of M31 and M33
The photometry depth of UKIRT covers the TRGB. In Figure 10, the density of stars decreases gradually and then increases in the J−K/K diagram from RGB to AGB, and the location of TRGB is at the lowest density. In mathematics, the saddle point needs to satisfy two conditions: one is that this  The apparent magnitude and color index of saddle points obtained by this algorithm are J−K=1.20 and K = 17.62 for M31 and, J−K=1.09 and K = 18.11 for M33 shown in Figure 13. In principle, the apparent magnitude depends on the distance and the metallicity of the galaxy, and the color index depends only on the metallicity. This provides a route to derive the metallicity and distance of M31 and M33. With the increase of metallicity, the J−K of the TRGB increases and the K-band absolute magnitude of the TRGB becomes brighter (Bellazzini et al. 2004;Górski et al. 2018). Using the relation between the K-band absolute magnitude and the ( ) -J K 0 of the TRGB from Górski et al. (2018), the absolute magnitude can be calculated as

Summary
The archival photometric data taken by UKIRT/WFCAM from mid-2005 to 2008 are used to select RSGs in M31 and M33, which is supplemented by the PS1, LGGS, and Gaia photometry and astrometric information. The foreground dwarfs are removed mainly by their obvious branch in the J−H/H−K diagram due to the significant darkening in H at the higher effective temperature of dwarfs than giants. This identification of dwarfs in the NIR color-color diagram is examined further by optical/infrared colors, specifically in the r−z/z−H and B−V/V−R diagram, and also supported by the Gaia measurement of parallax and proper motion. The depth of photometry, complete to about K=18 mag, combined with the criteria limited within the NIR colors, brings about a complete sample of RSGs in M31 and M33. The RSGs are identified in the members' J−K/K diagram from their outstanding locations caused by high luminosity and low effective temperature. The final sample includes 5498 and 3055 RSGs in M31 and M33, respectively. The control fields are used to estimate the pollution rate of the RSGs in our sample, which is found to be less than 1%. By comparing with the LMC, the SMC, and the MW galaxy, it is found that the number of RSGs per stellar mass decreases with metallicity, which can be understood as the metallicity's effect on the duration of the RSG phase for a star. In addition, the type of galaxy may also play a role in that an Sc type hosts more RSGs than an Sb galaxy.