Evolved Massive Stars at Low-metallicity I. A Source Catalog for the Small Magellanic Cloud

TL;DR

This paper provides a comprehensive multiwavelength catalog of over 45,000 sources in the Small Magellanic Cloud, focusing on evolved massive stars at low metallicity, with detailed classifications and candidate identifications.

Contribution

It introduces a new, well-validated catalog of SMC sources with multi-band data, spectral classifications, and candidate evolved massive stars, facilitating future low-metallicity stellar evolution studies.

Findings

Identified 1,405 RSG, 217 YSG, and 1,369 BSG candidates in the SMC.

Estimated the lower initial mass limit for RSGs as 6-7 solar masses.

Demonstrated the separation of RSGs from AGB stars at faint magnitudes.

Abstract

We present a clean, magnitude-limited (IRAC1 or WISE1 $\leq$ 15.0 mag) multiwavelength source catalog for the SMC with 45,466 targets in total, with the purpose of building an anchor for future studies, especially for the massive star populations at low-metallicity. The catalog contains data in 50 different bands including 21 optical and 29 infrared bands, ranging from the ultraviolet to the far-infrared. Additionally, radial velocities and spectral classifications were collected from the literature, as well as infrared and optical variability statistics were retrieved from different projects. The catalog was essentially built upon a $1''$ crossmatching and a $3''$ deblending between the SEIP source list and Gaia DR2 photometric data. Further constraints on the proper motions and parallaxes from Gaia DR2 allowed us to remove the foreground contamination. We estimated that about 99.5\% of the targets in our catalog were most likely genuine members of the SMC. By using the evolutionary tracks and synthetic photometry from MIST and the theoretical $\rm J-K_S$ color cuts, we identified 1,405 RSG, 217 YSG and 1,369 BSG candidates in the SMC in five different CMDs, where attention should also be paid to the incompleteness of our sample. We ranked the candidates based on the intersection of different CMDs. A comparison between the models and observational data shows that the lower limit of initial mass for the RSGs population may be as low as 7 or even 6 $M_{\odot}$ and the RSG is well separated from the AGB population even at faint magnitude, making RSGs a unique population connecting the evolved massive and intermediate stars, since stars with initial mass around 6 to 8 $M_{\odot}$ are thought to go through a second dredge-up to become AGBs. We encourage the interested reader to further exploit the potential of our catalog.

Tables

Table 1. Table 1: Number of Detected Targets in Each Filter of the SMC Source Catalog

Filter	$Numbe{r}_{detected}$	Absolute percentage	Relative percentagea𝑎aa𝑎a$a$Relative to the filter with the most matches in each dataset.
2MASS_J	45335	99.71%	100%
2MASS_H	45335	99.71%	100%
2MASS_$K_S$	45335	99.71%	100%
IRAC1	44712	98.34%	99.47%
IRAC2	44948	98.86%	100%
IRAC3	38272	84.18%	85.15%
IRAC4	35133	77.23%	18.16%
MIPS24	1261	2.77%	2.81%
WISE1	42460	93.39%	100%
WISE2	42460	93.39%	100%
WISE3	42460	93.39%	100%
WISE4	42460	93.39%	100%
Gaia_G	45466	100%	100%
Gaia_BP	45466	100%	100%
Gaia_RP	45466	100%	100%
VMC_Y	18637	40.99%	99.99%
VMC_J	18639	41.00%	100%
VMC_$K_S$	18638	40.99%	99.99%
IRSF_J	28645	63.00%	100%
IRSF_H	28592	62.89%	99.81%
IRSF_$K_S$	28513	62.71%	99.54%
AKARI_N3	555	1.22%	100%
AKARI_N4	551	1.21%	99.28%
AKARI_S7	274	0.60%	49.37%
AKARI_S11	134	0.29%	24.14%
AKARI_L15	40	0.09%	7.21%
AKARI_L24	18	0.04%	3.24%
HERITAGE_f70	1	—	—
HERITAGE_f100	1	—	—
HERITAGE_f160	2	—	—
HERITAGE_f250	2	—	—
HERITAGE_f350	2	—	—
SkyMapper_$u$	3185	7.01%	7.89%
SkyMapper_$v$	3711	8.16%	9.20%
SkyMapper_$g$	37949	83.47%	94.04%
SkyMapper_$r$	39930	87.82%	98.94%
SkyMapper_$i$	40356	88.76%	100%
SkyMapper_$z$	40171	88.35%	99.54%
NSC_$u$	26425	58.12%	68.24%
NSC_$g$	36700	80.72%	94.77%
NSC_$r$	37025	81.43%	95.61%
NSC_$i$	36297	79.83%	93.73%
NSC_$z$	38724	85.17%	100%
NSC_Y	1250	2.75%	3.22%
M2002_V	11630	25.58%	100%
M2002_B-V (B)	11630	25.58%	100%
M2002_U-B (U)	4000	8.80%	34.39%
M2002_V-R (R)	11563	25.43%	99.42%
GALEX_FUV	95	0.21%	57.93%
GALEX_NUV	164	0.36%	100%

Table 2. Table 2: Observation Epochs of ALLWISE and NEOWISE-R

Beginning	Ending
(MJD-54000)	(MJD-54000)
ALLWISE
	1400
1400	1600
NEOWISE-R
2700	2870
2870	3050
3050	3230
3230	3405
3405	3600
3600	3775
3775	3955
3955	4135
4135	4325
4325

Table 3. Table 3: SMC Source Catalog Contents

Column	Name	Description
1	ID	Index
2	R.A.(J2000)	Right Ascension, J2000 (deg)
3	Decl.(J2000)	Declination, J2000 (deg)
4-9	2MASS_J, …, e_2MASS_$K_S$	2MASS photometries with uncertainties (mag)
10-19	IRAC1, …, e_MIPS24	Spitzer photometries with uncertainties (mag)
20-27	WISE1, …, e_WISE4	WISE photometries with uncertainties (mag)
28-33	Gaia_parallax, …, e_Gaia_pmdec	Gaia DR2 astrometric solutions with uncertainties (mas; mas/yr)
34-36	Gaia_G, …, Gaia_RP	Gaia DR2 mean magnitudes (mag)
37-38	Gaia_RV, e_Gaia_RV	Gaia DR2 radial velocity with uncertainty (km/s)
39-44	VMC_Y, …, e_VMC_$K_S$	VMC photometries with uncertainties (mag)
45-50	IRSF_J, …, e_IRSF_$K_S$	IRSF photometries with uncertainties (mag)
51-62	AKARI_N3, …, e_AKARI_L24	AKARI photometries with uncertainties (mag)
63-72	HERITAGE_f70, …, e_HERITAGE_f350	HERITAGE fluxes with uncertainties (mJy)
73-84	SkyMapper_$u$, …, e_SkyMapper_$z$	SkyMapper photometries with uncertainties (mag)
85-96	NSC_$u$, …, e_NSC_Y	NSC photometries with uncertainties (mag)
97-104	M2002_V, …, e_M2002_V_R	M2002 photometries with uncertainties (mag)
105-108	GALEX_FUV, …, e_GALEX_NUV	GALEX photometries with uncertainties (mag)
109	B2011_SAGE_Class	IR color classifications from Boyer et al. (2011)
110	S2013_YSO	YSO classifications from Sewiło et al. (2013)
111	R2015_IRS_Class	MIR spectral classifications from Ruffle et al. (2015)
112	B2010_MK_Class	Optical spectral classifications from Bonanos et al. (2010)
113-115	GF2015_MK_SpVar, …, GF2015_MK_Class	Optical spectral classifications from González-Fernández et al. (2015) (km/s)
116-117	N2018_RV, N2018_MK_Class	Optical spectral classifications from Neugent et al. (2018) (km/s)
118-121	Simbad_RV, …, Simbad_Other_Types	Simbad classifications (km/s)
122-127	Median_WISE1, …, Amp_WISE1	WISE1-band variability statistics (mag)
128-133	Median_WISE2, …, Amp_WISE2	WISE2-band variability statistics (mag)
134-137	Median_IRAC1, …, Amp_IRAC1	IRAC1-band variability statistics (mag)
138-141	Median_IRAC2, …, Amp_IRAC2	IRAC2-band variability statistics (mag)
142-145	Median_VMC_Y, …, Amp_VMC_Y	VMC Y-band variability statistics (mag)
146-149	Median_VMC_J, …, Amp_VMC_J	VMC J-band variability statistics (mag)
150-153	Median_VMC_$K_S$, …, Amp_VMC_$K_S$	VMC $K_S$-band variability statistics (mag)
154	VMC_Type	VMC variable classifications
155-158	Median_IRSF_J, …, Amp_IRSF_J	IRSF J-band variability statistics (mag)
159-162	Median_IRSF_H, …, Amp_IRSF_H	IRSF H-band variability statistics (mag)
163-166	Median_IRSF_$K_S$, …, Amp_IRSF_$K_S$	IRSF $K_S$-band variability statistics (mag)
167-170	Median_Gaia_G, …, Amp_Gaia_G	Gaia G-band variability statistics (mag)
171-174	Median_Gaia_BP, …, Amp_Gaia_BP	Gaia BP-band variability statistics (mag)
175-178	Median_Gaia_RP, …, Amp_Gaia_RP	Gaia RP-band variability statistics (mag)
179-180	Best_Class_Gaia, Best_Class_Score_Gaia	Gaia variable classifications
181-184	Median_NSC_$u$, …, Amp_NSC_$u$	NSC $u$-band variability statistics (mag)
185-188	Median_NSC_$g$, …, Amp_NSC_$g$	NSC $g$-band variability statistics (mag)
189-192	Median_NSC_$r$, …, Amp_NSC_$r$	NSC $r$-band variability statistics (mag)
193-196	Median_NSC_$i$, …, Amp_NSC_$i$	NSC $i$-band variability statistics (mag)
197-200	Median_NSC_$z$, …, Amp_NSC_$z$	NSC $z$-band variability statistics (mag)
201-204	Median_NSC_Y, …, Amp_NSC_Y	NSC Y-band variability statistics (mag)
205-214	OGLE_LPV_Type, …, OGLE_LPV_Iamp3	OGLE long period variables statistics (mag)
215-241	OGLE_CCep_Mode, …, OGLE_CCep_phi31S	OGLE classical Cepheids statistics (mag)
242-252	OGLE_T2Cep_Type, …, OGLE_T2Cep_phi31	OGLE Type II Cepheids statistics (mag)
253-259	OGLE_Ecl_Type, …, OGLE_Ecl_DS	OGLE eclipsing binaries statistics (mag)

Table 4. Table 4: Evolved Massive Star Candidate Selection Criteria

Group	Color Criteria	Magnitude Criteria
BSGGaia		$G\le 2.333\times \left(\mathrm{BP}-\mathrm{RP}\right)+14.617$
		$G\le -3.200\times \left(\mathrm{BP}-\mathrm{RP}\right)+16.830$
YSGGaia		$G\le 14.750$
RSGGaia		$G\le 0.833\times \left(\mathrm{BP}-\mathrm{RP}\right)+13.833$
		$G\le -2.600\times \left(\mathrm{BP}-\mathrm{RP}\right)+19.670$
	$2.2\le \left(\mathrm{BP}-\mathrm{RP}\right)$	$G\le 13.950$
BSGSkyMapper		$i\le 2.933\times \left(r-i\right)+15.797$
		$i\le -6.667\times \left(r-i\right)+15.317$
YSGSkyMapper		$i\le 14.650$
RSGSkyMapper		$i\le 14.650$
		$i\le -5.200\times \left(r-i\right)+18.030$
	$0.9\le \left(r-i\right)$	$i\le 13.350$
BSGNSC		$r\le 2.182\times \left(g-r\right)+15.541$
		$r\le -6.000\times \left(g-r\right)+15.950$
YSGNSC		$r\le 14.750$
RSGNSC		$r\le 0.750\times \left(g-r\right)+14.225$
		$r\le -2.500\times \left(g-r\right)+18.450$
	$1.8\le \left(g-r\right)$	$r\le 13.950$
BSGM2002		$R\le 6.154\times \left(V-R\right)+14.527$
		$R\le -9.000\times \left(V-R\right)+16.800$
YSGM2002		$R\le 14.550$
RSGM2002		$R\le 1.333\times \left(V-R\right)+13.883$
		$R\le 4.000\times \left(V-R\right)+18.150$
	$1.1\le \left(V-R\right)$	$R\le 13.750$
BSG2MASS		$K_S\le 1.200\times \left(J-K_S\right)+14.670$
		$K_S\le -9.333\times \left(J-K_S\right)+16.250$
YSG2MASS		$K_S\le -3.000\times \left(J-K_S\right)+14.350$
RSG2MASS		$K_S\le -1.800\times \left(J-K_S\right)+13.750$
		$K_S\le -6.667\times \left(J-K_S\right)+18.617$
	$1.3\le \left(J-K_S\right)$	$K_S\le 9.950$

Table 5. Table 5: Numbers of Identified Evolved Massive Star Candidates

Dataset	Models	RSGs	YSGs	BSGs
Gaia	MIST	1029	97	952
SkyMapper	MIST	987	74	1035
NSC	MIST	627	34	606
M2002	MIST	836	133	757
2MASS	MIST	1068	93	879
2MASS	CB	1020
Total (cleaned)		1405	217	1369

Table 6. Table 6: Numbers of Ranked Evolved Massive Star Candidates

Ranks	RSGs	YSGs	BSGs
Rank 0	352	6	239
Rank 1	421	17	362
Rank 2	229	44	315
Rank 3	105	51	188
Rank 4	206	99	265
Rank 5	92
Total	1405	217	1369

Full text

1 Introduction

Massive stars are stars born with initial masses $\gtrsim 8\leavevmode\nobreak\ M_{\sun}$. They are relatively rare compared to the large number of low-mass stars due to the initial mass function and their short lifetimes. However, as a result of the intensive interior energy transfer and radiative output, they are responsible for some of the most extreme astrophysics in the Universe, including supernovae (SN), black holes, gravitational waves and long gamma-ray bursts, and are critical for stellar evolution, star formation and chemical evolution throughout the cosmic time (Humphreys & McElroy, 1984; Woosley et al., 2002; Massey, 2003; Meynet et al., 2011; Maeder & Meynet, 2012; Massey, 2013). The first generation of stars in the early Universe is expected to be massive and one of the main contributors of dust content in high-redshift galaxies (Massey et al., 2005; Bromm et al., 2009; Levesque, 2010). Therefore, understanding of the physical properties, evolution and mass loss of massive stars may help to reveal the formation of the primitive cosmic structures (Gall et al., 2011; Smith, 2014; Zhang et al., 2018).

Unfortunately, such stages of the early Universe can not be directly observed. The alternative way is to study the analogs of early Universe in our own cosmic backyard. In recent years, there has been a growing interest in massive stars in the metal-poor environments, due to the advancement of detectors and telescopes, and the needs to extrapolate the astrophysics from the local Universe to the early cosmic epochs. As a consequence, metal-poor star-forming dwarf irregular (dIrr) galaxies in the local Universe serve as an ideal laboratory for investigating the evolution and mass loss of massive stars at low-metallicity, since they mimic the behaviours of galaxies in the early Universe (Kunth & Östlin, 2000; McConnachie, 2012).

To understand the evolution of massive stars in the low-metallicity environment, both theoretical and observational constraints are needed. Among the many physical parameters of massive stars, one deterministic parameter is the mass loss, which has a profound impact on a star’s lifetime, luminosity (L), effective temperature ($T_{\textrm{eff}}$), radiation field, and the end fate as SN. However, the most important modes of mass loss are the most uncertain (Smith, 2014). It has been known that the mass-loss rates (MLR) adopted in modern stellar evolution codes for standard metallicity-dependent winds of hot main-sequence stars are overestimated by a factor of 2-3, due to the clumped and inhomogeneous stellar winds (Puls et al., 2008). Thus, other factors, e.g., the stellar winds, pulsations, rotation, convection and eruptions of evolved supergiants, as well as binary mass transfer, may have important contributions to the removal of the hydrogen envelope. Or, the overestimation of MLR of the hot stars may just simply raise the mass-limit for what stars can become Wolf-Rayet stars (WRs) through stellar wind mass-loss alone. Meanwhile, it has also been recognized that the unsteady modes of mass loss, like the episodic mass-loss events, are more important than previously thought (Smith & Owocki, 2006; Ofek et al., 2013). Consequently, observation in the infrared wavelengths, especially the relatively longer wavelengths, which are mainly dominated by dust emission, has a great impact on the understanding of MLR of massive stars.

Among the Local Group galaxies, the Large and Small Magellanic Cloud (LMC and SMC) are particularly intriguing due to their low-metallicity environments (about half and one-fifth of the Milky Way; Russell & Dopita 1992; Rolleston et al. 2002; Keller & Wood 2006; Dobbie et al. 2014; D’Onghia & Fox 2016) and close distances, for which individual stars can be resolved, allowing for a detailed analysis of their massive star populations in a variety of ways (Barba et al., 1995; Massey & Olsen, 2003; Evans & Howarth, 2008; Bonanos et al., 2009, 2010; Neugent et al., 2010; Yang & Jiang, 2011, 2012; Bouret et al., 2013; Kourniotis et al., 2014; Hainich et al., 2015; Castro et al., 2018; Yang et al., 2018; Britavskiy et al., 2019; Patrick et al., 2019).

In this paper, we focus on the evolved dusty massive star population in the SMC. We construct a clean catalog of massive stars, including their infrared properties (related to the MLR), astrometric solution regarding their membership to the SMC, time-series data revealing stellar variability, and evolutionary models in relation to their evolutionary stages, aiming to build a comprehensive anchor for future studies. The paper is structured as follows: the multiwavelength source catalog and time-series data are presented in §2 and §3, respectively. The identification of evolved massive star candidates is described in §4. The summary is given in §5.

2 Multiwavelength Source Catalog

Since our goal is to focus on the evolved dusty massive stars in the SMC, we use the Spitzer Enhanced Imaging Products (SEIP) source list as a starting point. SEIP source list contains sources detected with a high signal-to-noise ratio (S/N; 10-sigma level) in at least one channel among 12 near-infrared (NIR) to mid-infrared (MIR) bands of J (1.25 $\mu$m), H (1.65 $\mu$m), KS (2.17 $\mu$m), IRAC1 (3.6 $\mu$m), IRAC2 (4.5 $\mu$m), IRAC3 (5.8 $\mu$m), IRAC4 (8.0 $\mu$m), MIPS24 (24 $\mu$m), WISE1 (3.4 $\mu$m), WISE2 (4.6 $\mu$m), WISE3 (12 $\mu$m) and WISE4 (22 $\mu$m), from Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), Spitzer (Werner et al., 2004) and Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010)111http://irsa.ipac.caltech.edu/data/SPITZER/Enhanced/SEIP/docs/

seip_explanatory_supplement_v3.pdf. We retrieved the initial infrared data from SEIP source list with $3^{\circ}\leq$ R.A. $\leq 25^{\circ}$, $-75.5^{\circ}\leq$ Decl. $\leq-70^{\circ}$ and IRAC1 or WISE1 $\leq$ 15.0 mag, covering almost the whole SMC and also a small part of the Magellanic Bridge (MB). The magnitude cut of IRAC1 or WISE1 $\leq$ 15.0 mag was justified based on a drop-off ($\sim$14.85 mag) in the number counts for 12,748,156 ALLWISE WISE1 single-epoch measurements in the same region as shown in Figure 1. Considering the lower angular resolution of WISE ($\sim 6^{\prime\prime}$) and 2MASS ($\sim 5^{\prime\prime}$) compared to Spitzer ($\sim 2^{\prime\prime}$), and WISE sources within $3^{\prime\prime}$ of a SEIP source were reported, a deblending was performed with a search radius of $3^{\prime\prime}$. Targets having neighbors within $3^{\prime\prime}$ were excluded from the initial dataset, which resulted in 131,233 targets.

Besides the infrared detection, the reliable membership to the SMC is also a crucial factor for our study, for which the astrometric solution from Gaia Data Release 2 (DR2) is vital (Gaia Collaboration et al., 2016, 2018a). Given the very small offset between the Gaia DR2 and the 2MASS (median value of $\sim 0.120^{\prime\prime}\pm 0.157^{\prime\prime}$ for targets with $G\leq 18\leavevmode\nobreak\ mag$ in the same SMC region), as well as the SEIP source list and the 2MASS (median value of $\sim 0.089^{\prime\prime}\pm 0.237^{\prime\prime}$1), we crossmatched the deblended SEIP data and Gaia DR2 with a search radius of $1^{\prime\prime}$ to fix the position, and removed any SEIP targets with multiple counterparts to eliminate the blending. Then, a $3^{\prime\prime}$ crossmatching was performed and SEIP targets with multiple counterparts were removed again, which resulted in 74,237 targets. Since the effective angular resolution of the Gaia DR2 source list has improved to $\sim 0.4^{\prime\prime}$, with incompleteness in close pairs of stars starting below about $2^{\prime\prime}$ (Gaia Collaboration et al., 2018a; Arenou et al., 2018), we may reasonably believe that the vast majority of the blending is removed at the resolutions of both Spitzer and Gaia.

After crossmatching and deblending between SEIP and Gaia data, the Gaia DR2 astrometric solution was used to eliminate the foreground contamination (Lindegren et al., 2018). The selection of SMC members was restricted to targets with errors $<$0.5 mas/yr in proper motion (PM) and $<$0.5 mas in parallax. The first two panels of Figure 2 show the errors versus Gaia PMs in R.A. (left) and Decl. (middle), respectively. We fit a Gaussian profile to PM at each dimension as $0.695\leavevmode\nobreak\ (peak)\pm 0.240\leavevmode\nobreak\ (\sigma)\leavevmode\nobreak\ mas/yr$ in R.A. and $-1.206\pm 0.140\leavevmode\nobreak\ mas/yr$ in Decl., and calculated the limits of $\pm 5\sigma$ shown as the vertical dashed lines. The last panel (right) of Figure 2 shows the errors versus Gaia parallaxes. Similarly, a Gaussian profile fitting was adopted for the parallax as $-0.009\pm 0.066\leavevmode\nobreak\ mas$, while an additional elliptical constraint was also applied with the $5\sigma$ limits of $PM_{R.A.}$ and $PM_{Decl.}$ taken as the primary and secondary radii, respectively. The same criteria of $\pm 5\sigma$ were calculated for the parallax. In brief, the membership of the SMC is constrained by a Gaussian profile of the parallax with an additional elliptical constraint derived from $PM_{R.A.}$ and $PM_{Decl.}$, which results in 45,466 targets. Additionally to the previous astrometric constraints, we also applied another constraint on the radial velocity (RV) for all targets that have Gaia RV measurements, set at RV being larger than $\sim$90 km/s as shown below.

Further evaluation of the astrometric excess noise, which measures the disagreement, expressed as an angle, between the observations of a source and the best-fitting standard astrometric model (using five astrometric parameters), shows that about 99.5% of our targets have $astrometric\_excess\_noise\leq 0.5$ mas, indicating good astrometric solutions for the vast majority of our targets (Lindegren et al., 2018). Moreover, Gaia Collaboration et al. (2018b) determined the PMs and and parallaxes for 75 Galactic globular clusters, nine dwarf spheroidal galaxies, one ultra-faint system, and the MCs, where the method was also used by others (see, e.g., Aadland et al. 2018). The basic idea was to first determine the median and robust scatter in PMs and parallaxes by selecting a sample of stars covering a larger field of view (FoV), then further eliminated any sources showing larger scatter in the PMs, and finally constructed a filter based on a covariance matrix of the cleaned sample, which allows one to properly select out likely members. They also provided lists of possible members according to their analysis for all the systems. We crossmatched between our list of 45,466 targets and their list of $\sim$1.4 million sources for the SMC (their data covered larger FoV, which also yielded more identified sources than ours) with a search radius of 1”. The result was extremely good. Out of 45,466 targets in our list, there was only one mismatched target. We also compared their result with our LMC source catalog (in preparation), where 99.98% of our targets were matched. This proved the robustness of our method in selection of extragalactic objects by using Gaia astrometric solution, even without any sophisticated correction or filtering (in general, our method was actually similar to theirs). There were likely four reasons for this. Firstly, we put relatively strict constraints on the Gaia astrometric solution. Secondly, the SMC was close enough to have better quality of the data compared to more distant galaxies. Thirdly, we presumed that, for the current stage, compared to the Galactic targets, extragalactic objects (even as close as in the MCs) were probably only seen as tiny point sources by Gaia, where the fluctuation of photocenter caused by the binarity or convection/pulsation were almost invisible (Chiavassa et al., 2011; Pasquato et al., 2011; Messineo, & Brown, 2019). Fourthly, the crowding issue was largely mitigated since we deblended our sample. However, we would also like to emphasize that, in some rare cases, Gaia PMs and parallaxes might not good enough to distinguish a halo giant (e.g., at a distance of 10-15 kpc) from an SMC supergiant (at a distance of $\sim$60 kpc). Meanwhile, the RV might be not very helpful either, since the SMC’s RV was mostly a reflection of the Sun’s motion around the the Milky Way, and hence a halo star in the same direction might have a similar RV as an SMC member. Even with spectroscopy, it would be hard to distinguish, as a metal-poor giant in the halo might have similar metallicity as an SMC member. Therefore, we used the Besançon models of the Milky Way (Robin et al., 2003) to generate the expected number of Galactic stars as a function of magnitude and color for the same area of the sky as the SMC. The result of simulation indicated that, after applying the same constraints as above, the contamination by the halo stars was less than 0.25% in our sample and could be ignored.

Figure 3 shows $PM_{R.A.}$ versus $PM_{Decl.}$, for which the separation of selected SMC members, NGC104 and NGC362 is clearly shown. Meanwhile, based on the median number density of Figure 3 and the association between the SEIP source list and the 2MASS point source catalog (SEIP sources without valid 2MASS measurements), we estimated the contamination of remaining foreground sources and the possible non-point/background sources for the SMC were around 0.2% ($\sim$98/45,466) and 0.3% ($\sim$131/45,466), respectively, and could be ignored. Background point-like sources, for example, active galactic nuclei, quasars, or blue compact dwarfs, cannot be rejected at this stage, while the non-detection in 2MASS may also be due to the saturation, faintness, or photometric qualities of the targets. The estimation of contamination indicates that about 99.5% of the targets in our source catalog are most likely genuine members of the SMC.

We note that the strict constraints on the astrometric solution and the previous deblending procedure may cause target loss and incompleteness in our catalog to a certain extent, but ensure that we select the true SMC targets. This also can be seen from the histogram of Gaia RVs in Figure 4, where the separation of Milky Way and SMC is clear, and the vast majority of astrometry constrained targets with Gaia RVs larger than $\sim$90 km/s are selected with minimal value of $\sim$95 km/s (targets with both large PMs and RVs are not necessarily members of the SMC, as they could be hypervelocity stars in the Milky Way, runaway stars from the Milky Way or the SMC, or free-floating stars between the Milky Way and the SMC). Figure 5 illustrates the Gaia color-magnitude diagram (CMD) before (gray) and after (red) applying the astrometric constraints, where the large amount of foreground contamination of bright yellow and faint red stars is swept out. In addition, further constraints on the SEIP data, for example, $S/N\geq 3$ for IRAC1 ($i1\_fluxtype=1$), were investigated and resulted only in about 1.66% difference in the total number of sources (45,466 versus 44,712 targets) and were ignored.

Based on this fiducial dataset of 45,466 targets, we retrieved additional optical and infrared data from the following datasets with a search radius of $1^{\prime\prime}$:

• •

18,641 matches (41.00%) from VISTA survey of the Magellanic Clouds system (VMC) DR4. VMC is a NIR YJKS bands multi-epoch survey for the LMC, the SMC and the MB (KS $\lesssim$ 20.3 mag, $\sim$2% photometric and $\sim 0.01^{\prime\prime}$ astrometric precision), using the 4-meter near-infrared optimized Visible and Infrared Survey Telescope for Astronomy (VISTA; Cioni et al. 2011).

• •

28,678 matches (63.08%) from the IRSF Magellanic Clouds point source catalog (MCPS). IRSF MCPS is a NIR JHKS bands photometric catalog for a 40 $deg^{2}$ area in the LMC, an 11 $deg^{2}$ area in the SMC, and a 4 $deg^{2}$ area in the MB (KS $\lesssim$ 16.6 mag, photometric and astrometric accuracies for bright sources are 0.03-0.04 mag and $0.1^{\prime\prime}$, respectively), based on the data from Simultaneous three-color InfraRed Imager for Unbiased Survey (SIRIUS) camera on the InfraRed Survey Facility (IRSF) 1.4-meter telescope (Kato et al., 2007).

• •

625 matches (1.37%) from AKARI SMC bright point source list. The source list represents NIR to MIR imaging and spectroscopic observations of patchy areas in the SMC (N3 $\lesssim$ 16.5 mag, $\sim 0.1\leavevmode\nobreak\ mag$ photometric and $\lesssim 0.8^{\prime\prime}$ astrometric precision), using the Infrared Camera (IRC) aboard AKARI space telescope (Onaka et al., 2007; Murakami et al., 2007; Ita et al., 2010).

• •

4 matches from HERschel Inventory of the Agents of Galaxy Evolution (HERITAGE) band-merged source catalog (units are in flux [mJy] instead of magnitude). HERITAGE catalog is derived based on the data from both Photodetector Array Camera and Spectrometer (PACS; 100 and 160$\mu$m) and Spectral and Photometric Imaging Receiver (SPIRE; 250, 350, and 500$\mu$m) cameras on board the Herschel Space Observatory, to identify dusty objects in the LMC and SMC (the catalog also includes the Spitzer MIPS 70$\mu$m band data, and none of our targets has been detected in the SPIRE 500$\mu$m band; Pilbratt et al. 2010; Meixner et al. 2013; Seale et al. 2014).

• •

40,387 matches (88.83%) from SkyMapper DR1.1 (we apply constraints on the parameters as $flags<4$, $nch\_max=1$, $nimaflags=0$ and $class\_star\geq 0.9$ to retrieve the data; see Wolf et al. 2018 for details). SkyMapper is a southern hemisphere photometric survey in six bands ($u,v,g,r,i,z$; from magnitude 8 to 18, $\sim$1% photometric and $<0.2^{\prime\prime}$ astrometric precision), using the dedicated 1.3-meter SkyMapper telescope (Keller et al., 2007; Bessell et al., 2011; Wolf et al., 2018).

• •

38,759 matches (85.25%) from the NOAO source catalog (NSC) DR1 (we apply constraints on the parameters as $flags<4$ and $class\_star\geq 0.9$ to retrieve the data; see Nidever et al. 2018 for details). NSC is a catalog of sources from most of the public data taken on NOAO’s CTIO-4m+DECam as well as KPNO-4m+Mosaic3 (in $u,g,r,i,z$, Y bands; reach to $\sim$23 mag in most broadband filters with $\sim$1-2% photometric precision, and astrometric accuracy of $\sim$7 mas; Nidever et al. 2018).

• •

11,630 matches (25.58%) from a UBVR CCD survey of the Magellanic Clouds by Massey (2002) (M2002; note that due to the existence of some extreme large errors in the catalog, we replace them, e.g., error $>$ 1.0 mag, with NULL values). M2002 is a survey of 14.5 $deg^{2}$ region in the LMC and 7.2 $deg^{2}$ region in the SMC (V $\lesssim$ 18.0 mag, $<3\%$ photometric and $\sim 0.3^{\prime\prime}$ astrometric precision), using 0.61-meter Curtis Schmidt telescope at CTIO.

• •

164 matches (0.36%) from revised GALEX source catalog for the All-Sky Imaging Survey (GUVcat_AIS). GUVcat_AIS is a science-enhanced, “clean” catalog of GALEX ultraviolet (UV) sources with typical depth of 20.8 and 19.9 mag, position accuracy of $0.32^{\prime\prime}$ and $0.34^{\prime\prime}$ in NUV and FUV bands, respectively (Morrissey et al., 2007; Bianchi et al., 2017).

All the datasets are also deblended with a search radius of $3^{\prime\prime}$, prior to crossmatching with the SEIP-Gaia dataset. Figure 6 shows the normalized transmission curves of all filters used in our study. In total, we have 50 filters including 21 optical (including the two UV filters) and 29 infrared filters. The spatial distributions of the additional optical (left) and infrared (right) datasets are shown in Figure 7. The GALEX and HERITAGE data are not shown in the diagram due to the paucity of matches.

Besides the photometric data, additional classifications were also retrieved from the literature with a search radius of $1^{\prime\prime}$, including:

• •

37,375 matches from Boyer et al. (2011) (IR color classifications). It was an investigation of the IR properties of cool, evolved stars in the SMC, including the red giant branch (RGB) stars and the dust-producing red supergiant (RSG) and asymptotic giant branch (AGB) stars using observations from the Spitzer Space Telescope Legacy Program entitled “Surveying the Agents of Galaxy Evolution in the Tidally Stripped, Low Metallicity SMC (SAGE-SMC)”.

• •

30 matches from Sewiło et al. (2013) (IR color and SED classifications). They used CMDs based on the multiwavelength photometric data and SED fitting to identify a population of $\sim$1,000 intermediate- to high-mass young stellar objects (YSOs) in the SMC.

• •

43 matches from Ruffle et al. (2015) (MIR spectral classifications). They classified 209 point sources observed by Spitzer Infrared Spectrograph (IRS; Houck et al. 2004) using a decision tree method, based on infrared spectral features, continuum and spectral energy distribution shape, bolometric luminosity, cluster membership and variability statistics (all the targets from Kraemer et al. 2017 were also included).

• •

695 matches from Bonanos et al. (2010) (optical spectral classifications). It was a catalog of 5,324 massive stars from the literature with accurate spectral types, and a multiwavelength photometric catalog for a subset of 3,654 of these stars in the SMC, intending to study their infrared properties.

• •

198 matches from González-Fernández et al. (2015) (spectral variability flag, radial velocities and optical spectral classifications; for convenience, we only kept the first spectral classification for targets with multiple measurements). They studied physical properties of about 500 RSGs in the LMC and SMC by using NIR/MIR photometry and optical spectroscopy, aiming at exploring the fainter end of RSGs and extrapolating their behavior to other environments by building a more representative sample.

• •

113 matches from Neugent et al. (2018) (data originated from Neugent et al. 2010; radial velocities and optical spectral classifications). They spectroscopically observed 176 near-certain (Category 1) SMC yellow supergiant stars (YSGs) among $\sim$500 candidates to test against the evolutionary model in the low-metallicity environment.

• •

39,295 matches from Simbad (Wenger et al., 2000). The radial velocities, optical spectral classifications, main object types, and auxiliary object types were retrieved.

The unmatched targets are likely due to larger PMs, blends, or the quality cuts in SEIP or Gaia catalog.

This multiwavelength source catalog with 45,466 targets serves as the backbone of our study. The sample consists of “bona-fide” and “dusty” SMC targets determined by both astrometric measurement and infrared detection. Table 1 shows the absolute and relative percentage (relative to the filter with the most matches in each dataset) of detected targets in each filter. Figure 8 shows the histograms of magnitude distribution for each dataset (for convenience, the HERITAGE data is not shown here). The bin size is 0.1 mag, except for the GALEX (0.5 mag) and AKARI (0.25 mag) data. For WISE3 and WISE4 bands, due to the fact that the majority of the targets ($\sim$75% in WISE3 and $\sim$95% in WISE4 bands) have low S/N ($<2$) and are derived with a 95% confidence brightness upper limit, the histograms only show targets with $S/N\geq 2$.

3 Multiwavelength Time-Series Data

Following Yang et al. (2018), the MIR time-series data of WISE1 (3.4 $\mu$m) and WISE2 (4.6 $\mu$m) bands for all the 45,466 targets were collected from both ALLWISE (Cutri & et al., 2013) and Near-Earth Object WISE Reactivation mission (NEOWISE-R; Mainzer et al. 2014), with a search radius of $1^{\prime\prime}$ and the following parameter constraints, $qi\_fact>0$, $saa\_sep>0$, $moon\_masked=0$, $qual\_frame>0$, $det\_bit=3$, $\chi^{2}\leq 10$, $S/N\geq 3$ (see Yang et al. 2018 for details), which resulted in about 2.28 million measurements from ALLWISE and 8.78 million measurements from NEOWISE-R. With the NEOWISE 2019 data release, the total frame coverage is about twelve major epochs spanning $\sim$3,200 days ($\sim$8.8 years) with two epochs from ALLWISE and ten epochs from NEOWISE-R separated by an approximately three-year gap. The beginning and ending of each epoch set by us are given in Table 2. We binned the data within each epoch by using the median values of the date and magnitudes. For each epoch, we required at least five valid points to calculate the median value. Figure 9 shows examples of the original lightcurves overlapped with the binned lightcurves. The median absolute deviation (MAD) and standard deviation (SD) were used to calculate the long-term (full lightcurve) and short-term (within single epoch) variability of each target with at least five valid points, where the former was more resistant to outliers than the latter (Rousseeuw et al., 2009). Although it is possible to derive periods based on the current data, the period search may be highly contaminated by the strong alias structures due to the very low sampling of WISE data as shown in Figure 9, that we only have twelve epochs spanning $\sim$3200 days and each epoch only covers about 5 to 10 days (we refer interested readers to Chen et al. 2018 for the WISE catalog of periodic variable stars). In addition to MAD and SD, the full amplitude ($\rm Amp=max_{mag}-min_{mag}$) was also calculated for each target. More details about the WISE time-series data reduction can be found in Yang et al. (2018). In total, there are 39,694 (87.30%) targets having variability statistics in both WISE1 and WISE2 bands covering almost the whole area of our sample. The lack of variability statistics for some targets is likely due to either faint magnitudes or the quality cuts we adopted.

Besides the WISE data, we also collected other sets of variability data from different projects, including IR data from SAGE-Var program (Riebel et al., 2015), VMC, IRSF, and optical data from Gaia, NSC and OGLE. SAGE-Var is a follow-up to the Spitzer legacy program SAGE (Meixner et al., 2006; Gordon et al., 2011), for which six total epochs of photometric observations at IRAC1 (3.6 $\mu$m) and IRAC2 (4.5 $\mu$m) bands are obtained covering the bar of the LMC and the central region of the SMC with 15 different timescales ranging from $\sim$20 days to $\sim$5 yr. We collected the SMC data and calculated the median magnitudes, MADs, SDs and Amps for targets with all six epochs available. In total, there are 7,160 (15.75%) targets in IRAC1 band and 5,894 (12.96%) targets in IRAC2 band matched with our source catalog within $1^{\prime\prime}$. The left panel of Figure 10 shows the spatial distribution of those targets overlapped on our source catalog.

We currently rely on the available information from VMC DR4 regarding the variability statistics, as the final release (which may also include the time-series data) with global photometric and astrometric calibration will be made upon completion of the survey222http://horus.roe.ac.uk/vsa/. Since the observational cadences for different targets and filters were irregular and the median values of cadences varied from few hours to hundreds of days, for each filter, we constrained all the values of median magnitudes, MADs, SDs, and Amps to be within the range of 0 to 99 with at least five good measurements ($*nGoodObs\geq 5$) to avoid any unphysical values. In total, there are 11,197 (24.63%) targets matched with our source catalog within $1^{\prime\prime}$ as shown in the middle panel of Figure 10. There are 11,197 (24.63%) targets in Y band, 11,125 (24.47%) targets in J band, and 7,574 (16.66%) targets in $K_{S}$ band. Additionally, there are 204 targets that have been classified as Cepheids by the VMC team matched with our source catalog within $1^{\prime\prime}$.

For IRSF data, we retrieved the time-series data of $\sim$1,000 targets for each filter from Ita et al. (2018), where a very long-term (2000-2017) near-infrared variable star survey towards an area of $3\leavevmode\nobreak\ deg^{2}$ along the bar in the LMC and an area of $1\leavevmode\nobreak\ deg^{2}$ in the central part of the SMC was carried out with more than one hundred times repeated observations for each area. The median magnitudes, MADs, SDs, Amps were calculated for each targets. There are 160 targets in J band, 161 targets in H band, and 154 targets in KS band matched with our source catalog within $1^{\prime\prime}$ as shown in the right panel of Figure 10.

Gaia DR2 provides classifications for more than 550,000 variable sources consisting of different types of variables. However, only a subset of the variable stars classified as a certain type are characterized in detail and a fraction of the classifications may well be wrong (Gaia Collaboration et al., 2019; Mowlavi et al., 2018). Since the time-series data would be provided in the future release, we retrieved the variability statistics (including MADs, SDs and Amps) and classifications from Gaia Archive333https://gea.esac.esa.int/archive/ with typically $\sim$30 measurements spanning $\sim$620 days. However, we note that the MAD is scaled by 1.4826, so that the expectation of the scaled MAD at large number of measurements equals the standard deviation of a normal distribution444https://gea.esac.esa.int/archive/documentation/GDR1/datamodel

/Ch2/phot_variable_time_series_gfov_statistical_parameters.html. For consistency of the dataset, we reversed the scaled MAD by dividing it with 1.4826. There are 1,379 (3.03%; blue) targets with variability statistics which are matched with our source catalog within $1^{\prime\prime}$ as shown in the left panel of Figure 11, while 1,277 (2.81%; orange) of them are classified, and 868 (1.91%; red) of them having $best\_class\_score\geq 0.5$ ($best\_class\_score$ is a quantity between 0 and 1 provided by the classification pipeline to estimate the confidence of the classifier in the identification of various variability types), for which the classification may be acceptable.

The NSC time-series data of the SMC are mainly from the Survey of the Magellanic Stellar History (SMASH; Nidever et al. 2017). We retrieved the data from NOAO Data Lab555https://datalab.noao.edu/ with the same sky coverage and constraint parameters ($flags<4$ and $class\_star\geq 0.9$). By crossmatching with our source catalog using a search radius of $1^{\prime\prime}$, we retrieved 132,342 measurements in $u$-band, 355,823 measurements in $g$-band, 497,443 measurements in $r$-band, 210,496 measurements in $i$-band, 404,087 measurements in $z$-band, and 16,916 measurements in Y-band. However, due to the irregular sampling in different filters, the observational timescales vary from tenths of a day to more than one thousand days with several to hundreds measurements. Thus, we required at least five valid measurements for each target in individual filter in order to calculate the MAD, SD and Amp. We need to emphasize that as a result of the irregular sampling, and particularly, a strong systematic effect we discovered during the data processing (instead of a relatively uniform distribution, different SMC regions show variable levels of average variabilities), therefore, the calculated variability may not fully represent the true variability of the target and the user should be cautious when using these values. The middle panel of Figure 11 shows the spatial distribution of targets with variability statistics in different filters, including 17,548 (38.60%) targets in $u$-band, 25,897 (57.00%) targets in $g$-band, 29,572 (65.04%) targets in $r$-band, 19,632 (43.18%) targets in $i$-band, 35,850 (78.85%) targets in $z$-band, and 1,212 (2.67%) targets in Y-band.

Finally, we obtained additional data from the Optical Gravitational Lensing Experiment (OGLE; Udalski et al. 1992; Szymanski 2005; Udalski et al. 2008, 2015) by using a search radius of $1^{\prime\prime}$ from both OGLE-III Catalog of Variable Stars (O3CVS) and OGLE-IV Collection of Variable Stars (O4CVS), which resulted in 8,956 (19.70%) Long Period Variables (LPVs) from O3CVS, and 482 (1.06%) Classical Cepheids (CCeps), 11 (0.02%) Type II Cepheids (T2Ceps), and 87 (0.19%) Eclipsing Binaries (Ecls) from O4CVS (Soszyński et al., 2011, 2015; Pawlak et al., 2016; Soszyński et al., 2018) as shown in the right panel of Figure 11. However, since the OGLE data are calculated by using Fourier analysis, which is different from our dataset, only the classifications are used in the further analysis. Figure 12 shows Gaia color-magnitude diagrams with variable classifications from Gaia (left; only shows targets with $best\_class\_score\geq 0.5$), OGLE (middle), and VMC (right).

We also checked the time-series data from SkyMapper DR1.1. However, this first data release provides only few measurements ($\sim$1-2) covering $\sim$300 days, which is not suitable for robust variability calculation.

Our SMC source catalog of 45,466 targets is available in its entirety in CDS. Table 3 documents the columns found in the catalog. Targets without errors indicate either a 95% confidence upper limit, or the errors are simply too large to be reliable (e.g., $>$1.0 mag).

4 Identifying Evolved Massive Star Candidates on the Color-Magnitude Diagrams

As we focus on the evolved dusty massive stars, the primary task is to identify them. We utilize the evolutionary tracks and synthetic photometry from Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015, 2018) Isochrones & Stellar Tracks (MIST666http://waps.cfa.harvard.edu/MIST/; Choi et al. 2016; Dotter 2016), which covers a wide range of ages, masses and metallicities by using solar-scaled abundance under a single computational framework, to identify evolved massive star candidates on the CMDs of our multiwavelength source catalog.

We used a canonical value of 18.95 as the distance modulus of the SMC (Graczyk et al., 2014; Scowcroft et al., 2016). Since the metallicity of SMC is about 10% to 20% solar (Russell & Dopita, 1992; Dobbie et al., 2014; D’Onghia & Fox, 2016), we adopted the chemical composition of -1.0 to -0.7 dex for [Fe/H]. The non-rotation and rotation ($V/V_{crit}=0.40$) models of 7 to 40 $M_{\sun}$ were computed with no extinction and extinction of $A_{V}=1.0$ mag, respectively (Cardelli et al., 1989; Zaritsky et al., 2002; Haschke et al., 2011; Gao et al., 2013). We chose the color-magnitude combinations based on the available synthetic photometry in MIST with relative percentage $>90\%$ (see Table 1), at longer wavelengths (some dusty massive stars might not be identified in the shorter wavelengths due to higher extinction and reddening). In those models were also clearly shown the Yellow Void between blue supergiant stars (BSGs) and RSGs in order to identify YSGs.

Figure 13 shows multiple optical CMDs of Gaia, SkyMapper, NSC and M2002 datasets, overlapped with the evolutionary tracks of 7, 9, 12, 15, 20, 25, 32, and 40 $M_{\sun}$ generated from MIST synthetic photometry. The tracks are color coded based on the equivalent evolutionary phases (EEPs) from core helium burning to carbon burning, and $T_{\rm eff}$ of $7500\leavevmode\nobreak\ K<T_{\rm eff}$ (blue; BSGs), $5000<T_{\rm eff}\leq 7000\leavevmode\nobreak\ K$ (yellow; YSGs) and $T_{\rm eff}\leq 5000\leavevmode\nobreak\ K$ (red; RSGs) (Neugent et al., 2010). The regions of each type of evolved massive stars are outlined by the dashed lines with color and magnitude criteria listed in Table 4. Targets without errors are not shown in the CMDs, while the average photometric uncertainties are indicated when available. Due to the relatively large MLR during the RSGs phase, the star could be heavily obscured by the surrounding dust envelope (Smith et al., 2001; Massey et al., 2005; Levesque et al., 2006; Yang et al., 2018; Ren et al., 2019). Thus, we also empirically extended the RSGs region from the reddest and faintest points of the models to the even redder but not fainter area, in order to avoid the contamination from extreme AGBs (x-AGBs, stars are most likely to be experiencing a “superwind”, where the mass-loss rate can increase by a factor of 10, and a thick dust envelope obscures the star at optical wavelengths; van Loon et al. 2006; Boyer et al. 2011) as shown in the diagrams. It may be that some super-AGBs (stars in the mass range of $\sim 7-10\leavevmode\nobreak\ M_{\sun}$ that represent a transition to the more massive supergiant stars and are characterised by degenerate off-centre carbon ignition analogous to the earlier helium flash; Herwig 2005; Siess 2006, 2007, 2010; Groenewegen et al. 2009; Doherty et al. 2017) are also selected by the extension. However, we expect that they can be rejected by using several methods as shown in Yang et al. (2018). Moreover, inevitably, there will be contamination of the main sequence massive stars at the blue end, which cannot be easily disentangled. The diagrams show a clear bimodal distribution of the BSGs and RSGs candidates with few YSGs candidates lying between them. We notice that the conventional lower limit of the stellar mass for “massive star” is usually defined as eight solar mass. However, our initial tests show that the observational data are fitted better with the evolutionary track of seven solar mass (or even lower) as shown in the diagrams. This may be due to the treatment of parameters in the model such as convective overshooting, rotation, mixing, metallicity, and MLR, or the uncertainties of extinction and/or bolometric correction, or it may be true that “massive star” do have lower limit of stellar mass, which are complicated and need further investigation.

It can be seen from the Figure 13 that in the optical bands, generally, massive stars evolve horizontally across the upper part of CMD, while it is slightly different in the NIR bands. The left panel of Figure 14 shows the 2MASS CMD overlapped with MIST tracks, for which the tracks extend from the fainter and bluer region to the brighter and redder area due to the displacement of $T_{\rm eff}$ and intrinsic emission peaks. The right panel of Figure 14 shows a different method to classify RSGs originated from Cioni et al. (2006a) and Boyer et al. (2011) (hereafter CB method), where Oxygen-rich AGBs (O-AGBs) are defined by KS $<$ KS-band TRGB (KS-TRGB$\approx$12.7 mag) and K1 $<\rm J-K_{S}<$ K2, Carbon-rich AGBs (C-AGBs) are defined by KS $<$ KS-TRGB and K2 $<\rm J-K_{S}<$ 2.1 mag, x-AGBs are defined by $\rm J-K_{S}>2.1$ mag (all AGB stars are brighter than the K0 line, except x-AGBs), and RSGs are defined by $\rm\Delta(J-K_{S})=0.25$ mag from the O-AGBs shown as the dashed line in the diagram (the distance and 0.05 mag for the metallicity between LMC and SMC are corrected; Cioni et al. 2006b). It occurs to us that for the CB method, it covers almost the whole magnitude range of RSGs population down to the KS-TRGB, where the MIST tracks are unable to cover. However, MIST tracks are more broadened with a part of the tracks also extending to the O-AGBs region. This discrepancy between MIST tracks and CB method will be addressed in our next paper. A simple calculation by using a constant bolometric correction (BC) of $\rm BC_{K_{S}}=2.69$ (Davies et al., 2013) and $\rm A_{K_{S}}=0.1$ mag shows that the RSGs candidates close to KS-TRGB (e.g., $\rm K_{S}=12.6$ mag) only correspond to $\sim 10^{3.4}$ ($\sim$2,500) solar luminosity ($L_{\sun}$).

We also notice one important observational piece of evidence that, for the RSGs population, there is a distinct branch stretching continuously from the top of the luminous cool region towards the relatively faint warm area, reaching approximately to the tip of the red giant branch (TRGB) without blending into the AGBs population. It is obviously beyond the limit of 7 $M_{\sun}$ track and most probably reaching down to 6 $M_{\sun}$. False detection is ruled out due to two reasons. Firstly, all the Gaia, SkyMapper, NSC and 2MASS data show the same tendency except the M2002 data, which is likely due to the photometric sensitivity. Secondly, our sample is strictly constrained by the Gaia astrometric solution and verified by comparing with the result of Gaia Collaboration et al. (2018b), which indicate that the luminosities and colors of our targets are almost certain and the photometric errors are also very small at the magnitude range around the TRGB. Figure 15 shows the zoom in regions of $\rm 0.8<BP-RP<2.5$ and $\rm 14.0<G<18.0$ for Gaia (upper left), $\rm 0<r-i<1.2$ and $\rm 13.5<i<17.5$ for SkyMapper (upper right), $\rm 0.3<g-r<2.0$ and $\rm 14.0<r<18.0$ for NSC (bottom left), and $\rm 0.3<J-K_{S}<1.5$ and $\rm 10.5<K_{S}<15.5$ for 2MASS (bottom right). From the diagram and the observational point of view, there is no doubt that the separation between the RSG and AGB populations is quite clear, even at faint magnitude. Additionally, the distribution of a large sample of spectroscopic SMC RSGs (more than 300 and brighter than $\rm K_{S}\approx 11.0\leavevmode\nobreak\ mag$) on the 2MASS CMD follows almost exactly the MIST tracks, indicating the excellence of the MIST model prediction (details will be presented in our next paper). Thus, even though there may still be inevitable contamination of AGBs at the faint magnitude end of the RSG sample, we do believe that it should be small, and will be easily resolved by the next generation of large-scale spectroscopic surveys (e.g., 4MOST, MOONS; Cirasuolo et al. 2012; de Jong et al. 2012). Despite of the small contamination, the overlapping of lower (initial) masses of RSGs and the upper limits of AGBs, and the clear separation between them at faint magnitude range are also very interesting and puzzling. It means that, at given magnitude, some stars evolve to RSGs and the others evolve to AGBs. Currently, we have no certain explanation for this. It may be due to the rotation, binarity, chemical composition, MLR, or even more sophisticated evolution. Nevertheless, this may indicate the uniqueness of RSGs population, which connects the evolved massive and intermediate stars, since stars with initial mass around 6 to 8 $M_{\sun}$ are thought to go through a second dredge-up to become AGBs (Eldridge et al., 2007). The low mass RSGs may be also related to the intermediate luminosity optical transients (ILOTs; Prieto et al. 2008; Bond et al. 2009; Berger et al. 2009). Still, more investigations are needed to confirm the true nature of RSGs.

We combined the candidates of each type of evolved massive stars from different datasets and removed the duplications, which resulted in 1,405 RSG, 217 YSG and 1,369 BSG candidates listed in Table 5. Since the candidates were mostly selected based on the MIST model prediction in different datasets with variety of filters, sky coverage, photometric sensitivities and qualities, it was difficult to judge how reliable the candidacies were. Therefore, we ranked (Rank 0 to 5) the candidates based on the intersection between different CMDs, where Rank 0 indicated that a target was identified as the same type of evolved massive star in all five datasets (Gaia, SkyMapper, NSC, M2002 and 2MASS) by the MIST models and so on, and Rank 5 indicated the additional RSG candidates identified by the CB method but not recovered by the MIST models. Table 6 shows the numbers of ranked candidates for each type of evolved massive stars. Figure 16 illustrates two CMDs (Gaia and 2MASS) for all the candidates with ranks, where RSG, YSG, and BSG candidates are color coded in red, yellow, and blue colors ranging from dark (Rank 0) to light (Rank 5). Detailed information about each type of evolved massive star candidates are presented in separate tables available in CDS with the similar formats as Table 3, except the Rank has been added as an additional column in the end of each table. It has to be emphasized that some candidates may have different classifications at the same time, which is likely due to either the inevitable slightly overlapping of adjacent types of massive stars, or the larger uncertainties of photometry at the fainter magnitudes. Moreover, there are very few stray candidates way into other populations deserved to be investigated. For example, a couple of BSG candidates appearing in the RSGs region in the 2MASS CMD are likely due to their thick circumstellar dust envelope, since they are also very bright in the MIR bands. Similarly, we notice that there are few RSG candidates scattered in the much fainter and redder regions in the Gaia CMD. However, simultaneous inspection of the optical and MIR CMDs shows that almost all the scattered candidates in the bottom right region show infrared excess and/or high MIR luminosity, which may indicate that the dimming in the optical band is caused by the circumstellar dust envelope. Finally, Figure 17 shows the spatial distribution of evolved massive star candidates. It is clearly shown that due to the interaction between LMC and SMC, a bunch of candidates are stretched towards the MB. These populations follow the distribution of star formation in the SMC which is predominant along the bar of the SMC and extend to the LMC through the MB. This is depicted from the star cluster (especially for those with ages $<$100 Myr) distribution of the SMC (Bitsakis et al., 2018). Further detailed analysis of identified massive star populations will be presented in our following papers.

Finally, we would like to emphasize again one important thing that, our purpose was to study the evolved dusty massive stars in the SMC by identifying BSG, YSG, and RSG populations primarily based on the infrared detections, and constrained by the astrometric solutions and the evolutionary models. In that sense, our sample is NOT complete due to the strict constraints on the astromety, deblending, limitations of the models, and more important, the (presumably) serious differences in the luminosity completeness limits as a function of temperature. As we adopted the magnitude cut at the MIR wavelength (IRAC1 or WISE1 $\leq$ 15.0 mag), hotter stars would suffer more from incompleteness due to the weaker radiation at the far end of the Rayleigh-Jeans tail than the cooler stars. Even for the hot stars themselves, the degeneracy in the infrared detection would occur for hotter stars with larger MLR (less contribution from the stellar emission and more contribution from the dust emission) and cooler stars with smaller MLR (vice versa), since they might have similar radiative intensities at the same MIR wavelengths, which also happens in all magnitude ranges (not to mention the effect of binarity). One should NOT simply conclude that the BSGs to RSGs ratio (B/R ratio) is about 1:1 in the SMC, since both the observation and theoretical prediction indicates that there will be more BSGs than the RSGs (e.g., B/R ratio$\sim$4 or more) in the SMC (Meylan & Maeder, 1982; Humphreys & McElroy, 1984; Guo, & Li, 2002).

5 Summary

We present a clean, magnitude-limited (IRAC1 or WISE1 $\leq$ 15.0 mag) multiwavelength source catalog for the SMC with 45,466 targets in total. We intend to build our catalog as a comprehensive dataset serving as an anchor for future studies, especially for massive star populations at low-metallicity. The catalog contains data in 50 different bands including 21 optical and 29 infrared bands, retrieved from SEIP, VMC, IRSF, AKARI, HERITAGE, Gaia, SkyMapper, NSC, M2002, and GALEX datasets, ranging from ultraviolet to far-infrared. Additionally, radial velocities and spectral classifications were collected from the literature, as well as the infrared variability statistics, including MAD, SD, and Amp, derived from WISE, SAGE-Var, VMC and IRSF, and the optical variability statistics derived from Gaia, NSC and OGLE.

The catalog was essentially built upon a $1^{\prime\prime}$ crossmatching and a $3^{\prime\prime}$ deblending between SEIP source list and Gaia photometric data. We further constrained the PMs and parallaxes from Gaia DR2 to remove the foreground contamination, by applying a Gaussian profile in parallax with an additional elliptical constraint derived from $PM_{R.A.}$ and $PM_{Decl.}$. We estimated that about 99.5% of the targets in our catalog were most likely to be genuine members of the SMC.

By using the evolutionary tracks and synthetic photometry from MIST and also the theoretical $\rm J-K_{S}$ color cuts from CB method, we identified three evolved massive star populations in the SMC, namely the BSGs, YSGs and RSGs, in five different CMDs. There are 1,405 RSG, 217 YSG and 1,369 BSG candidates, respectively. We emphasized that our sample was NOT complete due to several reasons. We ranked the candidates based on the intersection of different CMDs, where the source with the most intersections was given the highest rank. A comparison between the models and observational data shows that, the lower limit of initial mass for the RSGs population may reach down to 7 or even 6 $M_{\sun}$ and the RSG is well separated from the AGB population even at faint magnitude, making RSGs a unique population connecting the evolved massive and intermediate stars, since stars with initial mass around 6 to 8 $M_{\sun}$ are thought to go through a second dredge-up to become AGBs.

We encourage the interested reader to further exploit the potential of our catalog, including, but not limited to, massive stars, supernova progenitors, star formation history, stellar population, stellar kinematics, chemical evolution, individual/integrated spectral energy distribution, time-domain astronomy, and so on. A further detailed analysis of the RSG population in the SMC will be presented in our next paper.

6 Acknowledgments

We would like to thank the anonymous referee for many constructive comments and suggestions. We acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement number 772086), and from Hubble Catalog of Variables project funded by the European Space Agency (ESA) under contract No.4000112940. B.W.J. and J.G. gratefully acknowledge support from the National Natural Science Foundation of China (Grant No.11533002 and U1631104). We thank Man I Lam and Stephen A. S. de Wit for helpful comments and suggestions.

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology. It is funded by the National Aeronautics and Space Administration. This publication makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

This research uses services or data provided by the NOAO Data Lab. NOAO is operated by the Association of Universities for Research in Astronomy (AURA), Inc. under a cooperative agreement with the National Science Foundation.

The national facility capability for SkyMapper has been funded through ARC LIEF grant LE130100104 from the Australian Research Council, awarded to the University of Sydney, the Australian National University, Swinburne University of Technology, the University of Queensland, the University of Western Australia, the University of Melbourne, Curtin University of Technology, Monash University and the Australian Astronomical Observatory. SkyMapper is owned and operated by The Australian National University’s Research School of Astronomy and Astrophysics. The survey data were processed and provided by the SkyMapper Team at ANU. The SkyMapper node of the All-Sky Virtual Observatory (ASVO) is hosted at the National Computational Infrastructure (NCI). Development and support the SkyMapper node of the ASVO has been funded in part by Astronomy Australia Limited (AAL) and the Australian Government through the Commonwealth’s Education Investment Fund (EIF) and National Collaborative Research Infrastructure Strategy (NCRIS), particularly the National eResearch Collaboration Tools and Resources (NeCTAR) and the Australian National Data Service Projects (ANDS).

This research has made use of the SIMBAD database and VizieR catalog access tool, operated at CDS, Strasbourg, France, and the Tool for OPerations on Catalogues And Tables (TOPCAT; Taylor 2005).

Based on data products from observations made with ESO Telescopes at the La Silla or Paranal Observatories under ESO programme ID 179.B-2003.