Testing Evolutionary Models with Red Supergiant and Wolf–Rayet Populations

We often refer to massive stars as "cosmic engines" as they are responsible for manufacturing many of the heavy elements in the universe (Johnson 2019). In addition, their strong stellar winds, and eventual disruption as core-collapse supernovae, provide most of the mechanical energy input to the interstellar medium, triggering new generations of star formation (see, e.g., Abbott 1982). Further, their UV radiation creates the H ii regions that populate late-type galaxies. Massive stars are also the source of the most energetic phenomenon yet found, emitting gamma-ray bursts (GRBs) as they collapse into black holes (Woosley 1993; Bloom et al. 2002). Mergers of their stellar remnants (black holes and neutron stars) disturb the curvature of spacetime, propagating as gravitational waves. Understanding massive star evolution, and getting the evolutionary models right (and knowing their limitations), is important not only for understanding the evolution of massive stars per se, but also for interpreting the light from distant galaxies (e.g., Stanway 2017) and a host of other astrophysically interesting phenomena.

Massive star evolution is not a solved problem. The physics of massive stars poses many challenges to modeling their evolution; chief among these is mass loss. The first stellar UV spectroscopic observations showed that mass loss is ubiquitous among massive OB stars (Morton 1967; Snow & Jenkins 1977). The incorporation of mass loss in early stellar models showed the dramatic effect this had on their evolution (see, e.g., Chiosi & Nasi 1974; De Loore et al. 1977, 1978; Chiosi et al. 1979; De Loore 1980; Lamers et al. 1980; Brunish & Truran 1982a, 1982b), broadening the main sequence and causing stars to evolve at higher luminosities than they would without mass loss. It also provides a mechanism to produce Wolf–Rayet (WR) stars through single-star evolution (Conti 1975; Chiosi et al. 1979; Chiosi & Maeder 1986). Aside from mass loss, rotation is another important ingredient in the models. Although we often think of stellar evolution as being determined by the initial mass and chemical composition (i.e., the Russell-Vogt theorem), the initial angular momentum has also been shown to be a third important parameter due primarily to the effects of rotationally induced mixing, but its inclusion in evolution is complicated. (See, e.g., Maeder & Meynet 2000 and Maeder 2009.) Finally, mass transfer and tidal interaction due to a binary companion will also have profound effects on the subsequent evolution of massive stars, as reviewed by Langer (2012) and De Marco & Izzard (2017). The inclusion of interacting binaries in massive star evolution models strongly affects the resulting stellar populations; see, e.g., Eldridge & Stanway (2020).

Despite these complications, the current generation of evolutionary models has enjoyed many successes in matching the observed distributions and properties of massive stars. Single-star evolutionary models have been shown to match the luminosity distribution of the very short-lived yellow supergiant phase in the LMC and in M33 (Drout et al. 2012; Neugent et al. 2012b). Both single-star and binary evolution models match the relative number of WN- and WC-type WR stars in nearby galaxies (Neugent et al. 2012a; Eldridge et al. 2017; Dorn-Wallenstein & Levesque 2018; Stanway et al. 2020), as well as the shape of the luminosity function of red supergiants (RSGs) in M31 (Neugent et al. 2020b). It should be emphasized that such tests involving the numbers of evolved massive stars are very exacting, as they turn a "magnifying lens" on stellar evolution calculations (Kippenhahn & Weigert 1990; Meynet & Maeder 2005).

Possibly the most vexing failure has been the inability of these models to correctly predict the relative number of RSGs and WRs. In the standard single-star evolutionary models (Ekström et al. 2012; Eldridge et al. 2017; Eggenberger et al. 2021), there is some mass limit (∼30M_⊙) below which massive stars evolve to the RSG phase and above which they evolve to the WR phase, possibly with a narrow band of masses where an RSG may become a WR. (For a somewhat different view, see Limongi & Chieffi 2018.) This upper mass for evolution to the RSG phase is set by the Eddington limit: stars that are more massive (and hence luminous) than some value cannot cross the Hertzsprung–Russell (H-R) diagram to cool, RSG-like temperatures, as opacities in their outer layers increase as the star evolves to cooler temperatures. At some point outwards radiation pressure approaches the pull of gravity, and the star belches out its cooler layers. This phase is often identified as the luminous blue variable (LBV) phase (see, e.g., Lamers & Fitzpatrick 1988). This process is revealed in the H-R diagram as the Humphreys–Davidson limit, the band of upper luminosity first characterized by Humphreys & Davidson (1979). This extends from the highest luminosities and effective temperatures diagonally to $\mathrm{log}L/{L}_{\odot }\sim 5.6$ and T_eff ∼ 15,000 K, and then horizontally (see Figure 22.1 in Lamers & Levesque 2017). In this classical, single-star picture, stars with masses above the RSG cutoff lose enough mass through stellar winds and (for the most massive stars) the LBV phase that the bare core of the star is exposed, and a WR star is born. In this picture, RSGs are the He-burning descendants of 9–30M_⊙ stars and have ages of 8–35 Myr (Ekström et al. 2012), while WRs come from higher-mass stars and have ages of 2–6 Myr (Georgy et al. 2012a). Even when binary evolution is included, WR stars have ages <10 Myr (Eldridge et al. 2008), still shorter than most RSGs.

Maeder et al. (1980) argued that the observed ratio of RSGs to WRs showed a very strong gradient with metallicity, increasing by two orders of magnitude within the Galaxy as one went from galactocentric distances of 7–9 kpc to 11–13 kpc, with a presumed decrease in metallicity. Furthermore, this trend seemed to be consistent with extension to the lower metallicities of the Large and Small Magellanic Clouds (LMC and SMC). Although the scaling of mass-loss rates with metallicity was unknown at the time, it was understood that these stellar winds were driven by radiation pressure in highly ionized metal lines; thus, mass-loss rates should be lower in stars of lower metallicity. This would mean that stellar winds would be more effective at producing WRs at higher metallicities. The Humphreys–Davidson limit was unrecognized at the time, and it was thought that RSGs were an earlier phase in the lives of most WRs, with stars with masses of 15–60M_⊙ passing through both phases. Thus, Maeder et al. (1980) proposed that the observed metallicity dependence in the RSG/WR ratio was due to the effect that different mass-loss rates would have on the relative lifetimes of the RSG and WR phases. However, even with the modern understanding that RSGs and WRs came from different mass ranges, the influence of metallicity on mass-loss rates still provided a reasonable explanation, as one might expect that the luminosity (mass) limit for RSGs would be higher in lower-metallicity environments, and hence the RSG/WR ratio would vary with metallicity.

Later, when stellar evolutionary models with scaling of mass-loss rates by metallicity were available, Massey (2003) found that (surprisingly) these early stellar evolution models really did not predict much of an expected change in the RSG/WR ratio at all. In Figure 1 we show how poorly the data and models compared at the time. Still, those models included neither rotation nor binarity, and the presumption was that as the models improved, so would the agreement, but that has not been the case. Eldridge et al. (2008) found that even by including binaries that the expected relationship in the RSG/WR ratio with metallicity remained relatively flat compared to the observed distribution (see their Figure 6). Similar disagreement can be seen in Figure 16 of Eldridge et al. (2017), which used an improved version of the binary models as well as a better estimate of the observed RSG/WR ratio provided by the first author of the present paper. Inclusion of rotation in the single-star models, and accounting for the effect of binarity in producing WRs, continued to show flatter distributions in the expected RSG/WR ratios with metallicity than the observations (Dorn-Wallenstein & Levesque 2020; Stanway et al. 2020).

Zoom In Zoom Out Reset image size

However, much has improved observationally in the decades since the Maeder et al. (1980) and Massey (2003) studies. New surveys for both RSGs and WRs have been carried out among the Local Group galaxies, with the result that we now have samples whose completeness are well understood (Section 2). Thus, we use this opportunity to present significantly improved values for the RSG/WR ratios as a function of metallicity (Section 3). We then go on to compare these ratios to two of the most used sets of evolutionary models, carefully describing our procedures so that others may follow in our footsteps with other evolutionary models, available either now or those of the future (Section 4). We summarize our results and discuss the implications in Section 5.

In order to compare the observed RSG/WR ratios with those predicted by the evolutionary models, we must make sure that the range of physical properties in our sample is well enough understood to correctly identify the corresponding phases in the tracks. The models give us the effective temperatures, luminosities, and surface compositions as a function of time, and in order to make valid comparisons, we must match these criteria to the physical properties of our observational sample. As we show below, this is relatively straightforward for the RSGs, as our samples have well-defined temperature and luminosity limits. However, for WR stars this is more complicated, as WRs are recognized observationally on the basis of emission lines formed in their stellar winds, and not simply on the basis of effective temperatures or luminosities. Fortunately, the surface chemical abundances, combined with temperature, provide the key for recognizing this phase in the evolutionary models.

Surveys for RSGs and WRs are now complete for the SMC, LMC, M31, and M33. In this section, we will describe the source lists and how we ascertain that we are identifying the RSG and WR phases in the models in the same way we have identified these stars observationally.

2.1. RSGs

The RSG populations of these four galaxies have all been recently identified using K_s versus J − K_s color–magnitude diagrams (CMDs). Early work by, e.g., Massey & Olsen (2003), Neugent et al. (2012b), and Drout et al. (2012) used either optical or near-infrared (NIR) CMDs to select RSG candidates in nearby galaxies. Radial velocity measurements were then used to remove foreground stars. However, a remaining concern was the contamination of the fainter stars in their sample by LMC asymptotic giant branch (AGB) stars. Contamination by AGBs could only be avoided by using a high-luminosity cutoff (typically $\mathrm{log}L/{L}_{\odot }\gt 4.9$), as first recommended by Brunish et al. (1986). Subsequently, the removal of foreground stars was greatly simplified by the availability of Gaia proper motions and parallaxes, as used by Aadland et al. (2018) to obtain a clean LMC RSG sample, and the AGB contamination problem was solved by Yang et al. (2019), who demonstrated that AGBs and RSGs in the SMC could be distinguished on the basis of their J − K_s colors, in accordance with the earlier theoretical and observational work of Cioni et al. (2006a, 2006b) and Boyer et al. (2011). Using these techniques, luminosity-limited samples of RSGs have now been identified over the whole of the LMC (Neugent et al. 2020a) and the two spirals M31 and M33 (Massey et al. 2021), and we adopt these source lists here.

Although Yang et al. (2019) identified a sample of RSGs in SMC, we have had to redo the selection here, for the following two reasons. The first concerns the lower limit used in J − K_s . Yang et al. (2019) used a lower limit cutoff, which paralleled the sloped upper-limit cutoff used to separate RSG and AGBs. As Neugent et al. (2020a) argue, there is no physical basis for a sloped lower limit, and using one has the potential to overlook RSGs with hot companions at the bright end. The second problem is that Yang et al. (2019) did not convert the photometry to physical properties, which we require for our analysis here. (The recent identification of RSGs in M31 and M33 by Ren et al. 2021 also suffered from these issues, and so we adopt the Massey et al. 2021 list for those galaxies). Thus, in the Appendix of this paper, we present our own source list of SMC RSGs.

Recognizing the RSG phase in the evolutionary tracks is relatively straightforward as long as we have a clear understanding of the luminosity and temperature criteria used in our observations. As emphasized in Massey et al. (2021), it makes little sense to talk about "the number of RSGs" in a galaxy without specifying the luminosity limit of the list. In addition, we must also be careful to match the temperature regime for what we are calling RSGs with the models. The J − K_s cutoffs that have been used for the source lists differ slightly from galaxy to galaxy, due to the change in the typical RSG temperature shifting to warmer temperatures at lower metallicities. Elias et al. (1985) first noted the shift toward earlier spectral types among the RSGs as one went from Galactic to LMC to SMC, and correctly attributed this to the changes in the Hayashi limit (Hayashi & Hoshi 1961) with metallicity. (See, e.g., Section 3.2 in Levesque 2017.)

In Table 1 we list the corresponding temperature ranges as a function of luminosity used in the four source lists. There does remain one uncertainty in relating this to the effective temperatures given by the evolutionary models. The conversion of dereddened J − K_s colors to effective temperatures comes about through the extension of the one-dimensional (1D) marcs atmosphere models (Gustafsson et al. 1975, 2008; Plez et al. 1992) to the low surface gravities appropriate for supergiants ($\mathrm{log}g[\mathrm{cgs}]\sim 0$), as described in Levesque et al. (2005). However, as noted by Levesque et al. (2005, 2006), and subsequent studies (see, e.g., Davies et al. 2013; Neugent et al. 2020a), there is a systematic difference between the effective temperatures derived from fitting the TiO bands using the MARCS models and fitting dereddened broadband colors, in the sense that the spectroscopic temperatures are cooler. (The size of this discrepancy using (V − K)₀ is 50–150 K according to Levesque et al. 2005, 2006). This problem has been attributed to the static 1D nature of the models (Levesque et al. 2005). Davies et al. (2013) argue that that the photometric (warmer) temperatures are more reliable as the continuum radiation is formed at deeper layers. However, the TiO spectroscopic temperatures are in better agreement with the evolutionary model temperatures (see, e.g., Figure 7 in Davies et al. 2013), although this is admittedly dependent upon the choice for α, the mixing-length parameter (Maeder & Meynet 1987; Chun et al. 2018). The TiO temperature determinations are also immune to uncertainties in the reddening corrections, which is an important consideration given the presence of circumstellar material around RSGs (Massey et al. 2005). A recent study by Taniguchi et al. (2021) reexamined the RSG temperature scale using a novel scaling of well-established red giant temperatures and the line-depth ratio method applied to Fe i lines, finding excellent agreement with the TiO method and the evolutionary models. (Note that their study is independent of the marcs models). In their Appendix, Neugent et al. (2020a) looked specifically at the differences between the temperatures derived from their TiO band fitting and their (J − K)₀ temperatures. The offset is about 200 K (see their Figure 7), again in the sense that the TiO band temperatures are smaller (cooler) than the ones derived from photometry. Thus, in our comparison below, we will use both the nominal (warmer) values in identifying the RSG phase in the evolutionary models and also a corrected version, found by subtracting 200 K from the photometric (J − K)₀ temperatures. Based on the fact that the latter is in better agreement with where stars fall on the evolutionary tracks, these would appear to be the more appropriate ones to use in identifying the RSG phase in the models.

Table 1. Temperature Selection of RSGs

Galaxy	Relations	Ref.
SMC:		1
	For $\mathrm{log}L/{L}_{\odot }\lt 4.32,{T}_{\mathrm{phot}}\lt 5613-279.1\mathrm{log}L/{L}_{\odot }$
	For $\mathrm{log}L/{L}_{\odot }\geqslant 4.32,{T}_{\mathrm{phot}}\lt 4407$
	For $\mathrm{log}L/{L}_{\odot }\lt 4.80,{T}_{\mathrm{phot}}\gt 5223-279.4\mathrm{log}L/{L}_{\odot }$
	For $\mathrm{log}L/{L}_{\odot }\geqslant 4.80,{T}_{\mathrm{phot}}\gt 3880$
LMC+M33-outer:		2,3
	For $\mathrm{log}L/{L}_{\odot }\lt 4.28,{T}_{\mathrm{phot}}\lt 5736-360.3\mathrm{log}L/{L}_{\odot }$
	For $\mathrm{log}L/{L}_{\odot }\geqslant 4.28,{T}_{\mathrm{phot}}\lt 4195$
	For $\mathrm{log}L/{L}_{\odot }\lt 4.78,{T}_{\mathrm{phot}}\gt 5326-360.0\mathrm{log}L/{L}_{\odot }$
	For $\mathrm{log}L/{L}_{\odot }\geqslant 4.78,{T}_{\mathrm{phot}}\gt 3604$
M33–inner+middle:		3
	For $\mathrm{log}L/{L}_{\odot }\lt 4.77,{T}_{\mathrm{phot}}\lt 5831-359.7\mathrm{log}L/{L}_{\odot }$
	For $\mathrm{log}L/{L}_{\odot }\geqslant 4.77,{T}_{\mathrm{phot}}\lt 4114$
	For $\mathrm{log}L/{L}_{\odot }\lt 4.78,{T}_{\mathrm{phot}}\gt 5324-359.7\mathrm{log}L/{L}_{\odot }$
	For $\mathrm{log}L/{L}_{\odot }\geqslant 4.78,{T}_{\mathrm{phot}}\gt 3604$
M31:		3
	For $\mathrm{log}L/{L}_{\odot }\lt 4.47,{T}_{\mathrm{phot}}\lt 5767-388.7\mathrm{log}L/{L}_{\odot }$
	For $\mathrm{log}L/{L}_{\odot }\geqslant 4.47,{T}_{\mathrm{phot}}\lt 4030$
	For $\mathrm{log}L/{L}_{\odot }\lt 4.73,{T}_{\mathrm{phot}}\gt 5249-388.7\mathrm{log}L/{L}_{\odot }$
	For $\mathrm{log}L/{L}_{\odot }\geqslant 4.73,{T}_{\mathrm{phot}}\gt 3412$
For all galaxies:		1,2
	Tspect = Tphot − 200

Note. T_phot refers to the effective temperature determined from the dereddened J − K_s photometry. T_spect refers to the effective temperature determined from spectroscopic analyses using the TiO band. See Section 2.1 of the text.

References. (1) This paper, (2) Neugent et al. (2020a), (3) Massey et al. (2021).

Download table as:  ASCIITypeset image

We do note that despite the uncertainty in the effective temperatures derived from the J − K_s photometry, this has a negligible effect on the derived luminosities. As described by Neugent et al. (2020a), the effect is on the order of 0.05 dex, comparable to the uncertainty in the values due to the uncertainties in the reddening corrections.

2.2. WRs

2.3. Spatial Coverage

In selecting our samples of RSGs and WRs, we must also make sure that the samples cover the same spatial extent. For the Magellanic Clouds, this was straightforward, as we had the luxury of designing the RSG studies to match the WR survey regions. Figures 2 and 3 show the RSG survey region (yellow circle) and the WR survey fields (green boxes).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figures 4 and 5 we show the survey regions for M31 and M33. The areas surveyed for the RSGs are outlined in yellow, while the surveyed regions for WRs are outlined in green. The WR survey regions are actually about 5% smaller than what are shown, in that there are gaps between the eight CCDs that made up the Mosaic detectors (see, e.g., Neugent & Massey 2011; Neugent et al. 2012a). Because we are also interested in investigating how the RSG/WR ratio changes with metallicity, and because we expect metallicity to be correlated with galactocentric distance ρ within the planes of these galaxies, it makes sense to restrict ourselves to the survey regions that are complete to a particular ρ. In M31 (Figure 4) the red oval denotes ρ = 0.75, where ρ has been normalized to an isophotal radius corresponding to μ_B = 25 mag arcsec², denoted R25. ⁷ We see that ρ = 0.75 falls just within the green Mosaic fields that denote the WR survey. For M33 (Figure 5) our survey for WRs extends all the way out to ρ = 1.0 thanks to the smaller angular size of the galaxy and the more favorable inclination. ⁸

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In studies of individual stars in other galaxies, crowding is always a concern, even within the Local Group. The presence of a close line-of-sight OB star of comparable luminosity to the WR star would lessen our detection limit. This issue was examined carefully for M33 by Massey & Johnson (1998) and Neugent & Massey (2011), who concluded that in M33 this was an issue only in the crowded cores of NGC 604 and NGC 595 and that these regions were surveyed for WRs using HST by Drissen et al. (1993). The same holds for the LMC, where the crowded R136 region of 30 Dor required HST to identify the individual O2–3/WN5–6 (Massey & Hunter 1998). The faintest WRs have M_V ∼ −3 (Crowther 2007), and thus any line-of-sight companion that would mask the detection of the WR would have to be very luminous, likely M_V < −5, i.e., a supergiant of some type. Such stars are extremely rare. As for the RSGs, Neugent (2021) showed that only a few percent of RSGs would be coincident with even B dwarfs, and such instance would not, in any event, mask the detection of the RSG in the NIR.

For M31 there is one additional concern. Massey et al. (2021) found that there were far fewer matches between their NIR-selected sample of RSGs and the optical Local Group Galaxy Survey (LGGS; Massey et al. 2006, 2016) than in M33. This could be due to extinction within the plane of M31, causing fewer optical detections on the far side of M31's disk. Massey et al. (2021) argue against this interpretation but it remains a concern. We therefore follow Neugent (2021) in considering just the subsample of fields where there are a high percentage (>50%) of matches. As we will see below, the results in terms of the RSG/WR ratio are essentially indistinguishable between the full sample and this subsample.

2.4. Star Formation Concerns

In this section, we present the counts of both the RSGs and WRs in these galaxies. We then compute the ratios and their uncertainties employing a Bayesian framework, following Dorn-Wallenstein & Levesque (2020). Finally, we discuss the likely initial metallicities of these stars in preparation for comparing the observed ratios to those predicted by the evolutionary models in Section 4.

3.1. The Number of RSGs

3.2. The Number of WRs

In Table 2 we list the number of WRs within our four galaxies based on the following surveys.

The source list of WRs for the SMC comes from Massey et al. (2003). There are 12 WRs known, and the recent deep survey (summarized in Neugent et al. 2018) failed to reveal any more. We show their distribution in comparison with the RSGs in Figure 6.

Zoom In Zoom Out Reset image size

For the LMC, the source list comes from Neugent et al. (2018), who list 154 stars. Of these, two (BAT99–93 and BAT99–105) were "downgraded" from O2–3If/WN6 stars to O2–3If stars, and we do not include them here. Another 10 are classified as O2–3If/WN5–6 and are mostly found in the R136 cluster at the heart of 30 Dor. As argued above, they are unlikely to be recognized as a WRs in the evolutionary models and so we exclude them here as well. However, we have retained the other two types of "slash" stars, the Ofpe/WN9s, and the WN3/O3s, as evolutionary models should recognize both types as WRs, as discussed above. We are left with 142 WRs. Removing stars within 30 Dor and Constellation III leaves 111 WRs. In Figure 7 we show the distribution of the 142 WRs in comparison to the 720 RSGs with $\mathrm{log}L/{L}_{\odot }\geqslant 4.5$ dex and indicate the exclusion regions around 30 Dor and Constellation III. Note the large number of RSGs in the Constellation III region with few WRs. These RSGs are likely the result of the 12–16 Myr burst found by Dolphin & Hunter (1998) mentioned earlier.

Zoom In Zoom Out Reset image size

For M31 our source list comes from Neugent et al. (2012a), who list 154 WRs. One additional WR star was subsequently found by Shara et al. (2016). (Our inspection of the survey images used by Neugent et al 2012a reveals that this star went undetected as it unfortunately fell in one of the narrow gaps between the CCDs in the Mosaic camera that they used.) Of these 155 WRs, 7 have ρ > 0.75 and are not included in our counts here. None of these are classified as O2–3I/WN5–6. We are therefore left with 148 WRs in our comparison sample. The distribution of RSGs and WRs are shown in Figure 8. Restricting ourselves to only the regions with low extinction (i.e., with good matches between the RSGs and the LGGS), we find 69 WRs.

Zoom In Zoom Out Reset image size

For M33, the counting is also a bit complicated. Neugent & Massey (2011) lists 206 WRs and six additional WR candidates. Neugent & Massey (2014) confirmed five of these as WRs; the sixth proved to be an Of-type star. Subsequently, Massey et al. (2016) confirmed another WR, J013344.05+304127.5 at ρ = 0.13. They also note that the star J013418.37 + 303837.0, classified as Ofpe/WN9 by Neugent & Massey (2011), is actually Var 2, one of the original Hubble-Sandage (Hubble & Sandage 1953) stars, now called LBVs. (This confusion came about due to the previous lack of good coordinates for the Hubble-Sandage variables.) They therefore discounted it as a WR, but here we retain its WR status as argued above. The first two authors of the present paper are involved in a deep search for additional WRs in M33, particularly WN3/O3s; so far, they have been able to confirm one new star, J013254.12+302313.7 (ρ = 0.70), as a WN3/O3. Thus, the total number of known WR stars in M33 is 213; all have ρ ≤ 1.0.

From this M33 total we must remove several likely hydrogen-burning WRs. Conti & Massey (1981) noted the similarity of several WRs in M33's brightest H ii regions (NGC 604, NGC 595, and NGC 592) to the overly bright WRs in R136 and NGC 3603 that we now classify as O2–3/WN5–6 stars. Currently available spectra do not have the necessary signal-to-noise ratio to accurately classify these M33 stars as O2–3/WN5–6, as the classification requires us to see absorption as well as emission. However, it is likely that several of them fall into this category. Inspection of the Neugent & Massey (2011) list allowed us to identify four stars (J013332.97+304136.1, J013311.85+303852.7, J013332.82+304146.0, and J013432.11+304705.8) whose brightness and spectral types suggest they may fall into this category. All four are located in M33's brightest three H ii regions as originally found by Conti & Massey (1981). ¹⁰ Thus the total number in our M33 comparison sample is 209. We show the distribution of RSGs and WRs in Figure 9 and demark the inner, middle, and outer zones.

Zoom In Zoom Out Reset image size

3.3. The Ratios

In Table 3 we list the RSG/WR ratios along with their associated errors. In many previous studies (see, e.g., Massey & Johnson 1998; Massey & Olsen 2003; Neugent et al. 2012a) we have used the ratio of the observed numbers directly and assumed that the errors were dominated by stochastic processes, i.e., that if we observed N objects today, at some other time, we might observe $N\pm \sqrt{N}$ such objects. Thus, the uncertainty σ_R in the ratio R = N_A /N_B is then computed using ${\sigma }_{R}=R\sqrt{\displaystyle \frac{{\sigma }_{{N}_{A}}^{2}}{{N}_{A}^{2}}+\tfrac{{\sigma }_{{N}_{B}}^{2}}{{N}_{B}^{2}}}$. With ${\sigma }_{A}=\sqrt{{N}_{A}}$ and ${\sigma }_{B}=\sqrt{{N}_{B}}$, this simplifies to ${\sigma }_{R}=R\sqrt{{N}_{A}^{-1}+{N}_{B}^{-1}}$. However, Dorn-Wallenstein & Levesque (2020) have recently argued one can do better than this by borrowing from the X-ray astronomy community and using a Bayesian framework. This avoids such issues as when N_B = 0 and the ratio is infinite, but finite values are also well within a 3σ error. Details are given in Section 3 of Dorn-Wallenstein & Levesque (2020), and we adopt their methodology here, and refer to the Bayesian-derived "true" ratio $\hat{R}$. They adopt the Markov Chain Monte Carlo package emcee (Foreman-Mackey et al. 2013) to sample the posterior probability distribution:

$$\begin{eqnarray}&&\begin{array}{l}p(\hat{R}| {N}_{A},{N}_{B})\propto {\displaystyle \int }_{{\hat{N}}_{B}}d{\hat{N}}_{B}{\hat{R}}^{\phi -1}{\hat{N}}_{B}^{2\phi -1}\\ \ \ \times \,\displaystyle \frac{{\hat{R}}^{{N}_{A}}{\hat{N}}_{B}^{({N}_{A}+{N}_{B})}{e}^{-{\hat{N}}_{B}(\hat{R}+1)}}{{N}_{A}!{N}_{B}!},\end{array}\end{eqnarray} \tag{ 1 }$$

where ${\hat{N}}_{B}$ is a nuisance parameter that they marginalize over corresponding to the "true" number of stars in the denominator, and ϕ is a parameter to ensure good coverage, which they set to 1/2. The $\hat{R}$ values, and their associated uncertainties, are given in Table 3 and will be what we use for comparing with the values derived from the evolutionary tracks.

Table 3. Number Ratios as a Function of Metallicity

Galaxy	log(O/H)+12	Z	${\hat{R}}_{\mathrm{RSGs}/\mathrm{WRs}}$
			L/L⊙ ≥ 4.2	L/L⊙ ≥ 4.5	L/L⊙ ≥ 4.8
			(9M⊙)	(12M⊙)	(15M⊙)
SMC	8.00	0.003	${41.0}_{-10.0}^{+14.0}$	${23.3}_{-5.7}^{+8.2}$	${8.4}_{-2.1}^{+3.0}$
M33–Outer a	8.29	0.006	${9.5}_{-1.1}^{+1.3}$	${4.4}_{-0.5}^{+0.6}$	${1.8}_{-0.3}^{+0.3}$
LMC(total)	8.38	0.007	${8.7}_{-0.7}^{+0.8}$	${5.0}_{-0.4}^{+0.5}$	${2.2}_{-0.2}^{+0.2}$
No 30Dor/III b	8.38	0.007	${9.1}_{-0.9}^{+1.0}$	${5.0}_{-0.5}^{+0.6}$	${2.3}_{-0.2}^{+0.3}$
M33–Middle c	8.41	0.008	${8.3}_{-1.0}^{+1.2}$	${4.2}_{-0.5}^{+0.6}$	${1.6}_{-0.2}^{+0.3}$
M33–Inner d	8.72	0.015	${5.9}_{-0.7}^{+0.8}$	${3.1}_{-0.4}^{+0.5}$	${1.5}_{-0.2}^{+0.2}$
M31 e	8.90	0.030	${12.5}_{-1.0}^{+1.1}$	${5.1}_{-0.4}^{+0.5}$	${2.4}_{-0.2}^{+0.2}$
Low Ext. 50%	8.90	0.030	${12.8}_{-1.5}^{+1.8}$	${5.5}_{-0.7}^{+0.8}$	${2.5}_{-0.3}^{+0.4}$

Notes.

^a M33–Outer: 0.5 ≤ ρ ≤ 1.0. ^b After removing stars within 10' of 2α₂₀₀₀ = 05:38:42.40, δ₂₀₀₀ = −69:06:03.4 (30 Dor), and within 45' of α₂₀₀₀ = 05:31:24.33, δ₂₀₀₀ = −66:54:45.0 (Constellation III). ^c M33–Middle: 0.25 ≤ ρ < 0.5. ^d M33–Inner: ρ < 0.25. ^e ρ ≤ 0.75.

Download table as:  ASCIITypeset image

We note that the distinction between the Bayesian-derived "true" value (based on the expected rational numbers) and the actual observed values (based on the integer values) differ only very slightly. The same is true for the uncertainties. The largest difference is of course for the situation with the smallest values. For instance, the Bayesian value for the RSG/WR ratio in the SMC for RSGs with $\mathrm{log}L/{L}_{\odot }\geqslant 4.8$ is $\hat{R}={8.4}_{+3.0}^{-2.1}$. This is essentially indistinguishable from the "classical" ratio and uncertainties calculated from the values in Table 2, 8.1 ± 2.5.

3.4. Metallicities

As discussed in the Introduction, the previous data showed a steeply falling RSG/WR ratio with increasing metallicity, while the (old) evolutionary models predicted a much flatter relation was expected. A causal examination of the values in Table 3 shows that the new ratios calculated here no longer span the two orders of magnitudes shown in Figure 1. Rather, they range over less than a single order of magnitude. Nor are they monotonic with metallicity. How do the new data now compare with modern evolutionary models?

For this comparison, we will use the current generation of the single-star Geneva evolutionary models as well as the Binary Population and Spectral Synthesis (BPASS) evolutionary models, as described below. We have selected these models because they have been shown to be very successful in other critical comparisons with massive star properties as described in the Introduction, and yet somehow failed to match the older RSG/WR ratio data. However, in the hopes that present and future workers will compare their models to our data in a similar fashion, we describe in some detail here how we identify the RSG and WR phases in the evolutionary tracks in such a way that they match what we do observationally.

The evolutionary models allow us to compute the lifetimes of RSGs and WRs as a function of initial mass. In order to turn these into an RSG/WR number ratio, we scale these lifetimes using the assumed slope of the IMF adopting Γ = −1.35 (see Section 2.4). We will begin by comparing with single-star models and then add in the complication of binarity.

Current versions of the Geneva models are so far available at four metallicities, Z = 0.014 (Ekström et al. 2012), 0.006 (Eggenberger et al. 2021), 0.002 (Georgy et al. 2013), and 0.0004 (Groh et al. 2019). The Z = 0.0004 models are not useful for our purposes, as they are significantly lower in metallicity than any of the galaxies we are considering here. At each metallicity and initial mass, there are models computed with both no initial rotation and with an initial rotation that is 40% of the breakup speed. A real ensemble of stars will likely have a range of rotation rates and so we have linearly interpolated between these values, i.e., a f_rot of 1 corresponds to the rotating tracks (initial rotation speed of 40% of the breakup speed), while an f_rot of 0 corresponds to the nonrotating tracks. Note that intermediate values of f_rot are approximate and were used to simply linearly interpolate between the results of the synthetic populations, and their values are not exact, as discussed in more detail in Section 2.4 of Dorn-Wallenstein et al. (2017).

We show the comparison between the predictions of the Geneva models with the RSG/WR ratios in Figure 10. The plots on the left show the results when we compare the results adopting the effective temperatures derived photometrically. The agreement is far better than in previous studies (compare with Figure 1), but the models predict too small a value at all luminosities. However, adopting the temperature limits corrected to the spectroscopic values, as shown on the right, is very good even at the lowest RSG luminosity we consider (4.2 dex), and excellent at the higher luminosities. We cannot predict what will happen at the higher metallicity corresponding to M31 (Z = 0.030) as modern Geneva models are not yet computed for this high a metallicity, but older Geneva models showed a flattening at higher metallicities, as shown in Figure 1. We note that the decrease in the RSG/WR ratio with an increase of metallicity is due both to an increase in the number of WRs and to a decrease in the number of RSGs. For the WRs, the increase is due both to the fact that the minimum mass for forming WR stars decreases with metallicity and to the fact that the WR lifetimes increase with metallicity. For the RSGs, the decrease is due to the fact that at higher luminosities mass loss results in the higher luminosity RSGs evolving back to the blue.

Zoom In Zoom Out Reset image size

Binary interaction is the other channel for driving the evolution of massive stars, by either donating or stripping off mass. Although the single-star evolutionary models have enjoyed remarkable success in passing the observational tests we (and others) have devised over the years, there is no escaping the fact that a substantial fraction of O-type stars are found in binary systems. The value is at least 35%–40% (Garmany et al. 1980; Sana et al. 2013) and may be 70% or even higher if long-period systems are included in the analysis (Aldoretta et al. 2015). The fraction of these systems that are close enough to interact during the star's lifetime is arguable, with estimates ranging up to 100% (Sana et al. 2012). A major complication is accounting for stellar mergers: as the more massive star expands during its evolution, it can engulf its neighbor, leaving behind no obvious sign (see, e.g., de Mink et al. 2014).

An essential tool for investigating the importance of binarity in massive star evolution is the BPASS evolutionary models (Eldridge et al. 2008; Eldridge & Stanway 2009; Eldridge et al. 2017; Stanway & Eldridge 2018). These models have been used for a variety of astrophysical phenomena; recent examples include comparing the expected and observed gravitational-wave detection rates (Eldridge et al. 2019), resolving the puzzle of an overly massive black hole (Eldridge et al. 2020), and identifying a pathway for creating hydrogen-poor superluminous supernovae (Stevance & Eldridge 2021). Their intent, however, is primarily in exploring the effects of including binarity on stellar populations (see, e.g., Dorn-Wallenstein & Levesque 2020; Stanway et al. 2020).

We show the comparison between the BPASS v2.2.1 models (Stanway & Eldridge 2018) ¹¹ and the observations in Figure 11. The binary fraction f_bin is again based on a linear approximation between the results of the BPASS single-star models and those containing a 100% binary fraction. The plots on the left show the comparison using the photometric temperature limits, while the plots on the right show the effect of applying the 200 K correction to bring these into accord with the spectrophotometric values. The data suggest a modest initial binary fraction (50%) for most metallicities. Taken at face value, the two extreme metallicity cases, SMC and M31, imply a much smaller initial fraction of binaries is needed. We discuss alternative interpretations below.

Zoom In Zoom Out Reset image size

We know that binarity and rotation can have similar effects on the predictions of evolutionary models, such as the relative number of WC- and WN-type WR stars (see, e.g., Dorn-Wallenstein & Levesque 2020). Here we see a quite different behavior. Including rotation shifts the RSG/WR ratios downwards but does not really alter the shape (Figure 10). On the other hand, adding in binarity not only lowers the curves but also changes the shape (Figure 11). Although binarity does affect the number of RSGs, the primary effect is in creating more WRs. Furthermore, the WRs created by binary evolution tend to have longer lifetimes than those created by single-star evolution because they have smaller masses. Even for higher-mass stars that would have become a WR star without binary evolution, binary interactions will cause the star to become a WR star sooner, so more of its evolution is during a WR phase. The relative importance of the binary channel becomes less important at higher metallicities, as the higher mass-loss rates that come with larger metallicities lowers the mass limit for becoming a WR. These lower-mass WRs have longer lifetimes during the WR phase, and so the net effect of binarity is to flatten the RSG/WR ratio as a function of metallicity.

We have determined the RSG/WR ratio for four Local Group galaxies, covering an order of magnitude in metallicity, from Z = 0.003 to Z = 0.030. We find that there is only a modest variation in this quantity with metallicity, contrary to the original suggestion by Maeder et al. (1980) ¹² and earlier estimates (e.g., Massey 2003). However, in the larger sense, Maeder et al. (1980) has certainly proven to be right, in that the RSG/WR ratios provide a very exacting test of modern stellar evolution theory.

In making the comparison with the Geneva and BPASS evolutionary models, we have been careful to match what we are using to define the RSG and WR phases in the models with what we are using for identifying these stars observationally. For the WRs, this means ignoring the handful of WRs that are likely still hydrogen-burning classified as O2–3If/WN5–6 stars and found in the largest H ii regions. For the RSGs, we have used both a luminosity limit in defining the number of RSGs and used the same temperature criteria in the models as we did observationally. When we do this, both the Geneva single-star models and the BPASS binary models do a good job of matching the observed ratios, as long as we adopt the temperature criteria corresponding to the spectroscopic effective temperatures rather than that derived from photometry. This is consistent with the recent analysis by Taniguchi et al. (2021) but runs contrary to the conclusion by Davies et al. (2013). When dealing with stars with extended atmospheres, such as RSGs, the definitions of radius and effective temperatures may differ between evolution interior models and atmospheric models. The fact that the spectroscopic temperatures agree better with the placements of stars in the HRD relative to the tracks is therefore encouraging.

Another factor to consider is that what we identify as RSGs and WR stars in the evolutionary models is, to some extent, uncertain. For RSGs, the model stars become large and cool, but it is unclear if their theoretical surface temperatures mean exactly the same thing as what we measure observationally as the temperatures of real RSGs. It is possible that we are missing some of the RSGs at lower metallicity because they remain hotter than the temperature limits imposed here. The correct boundary conditions for RSGs are also unknown, and the BPASS and Geneva models both use different methods to evaluate the surface of the RSG phase of evolution.

For WR stars, the luminosities of the evolutionary models mean the same thing as what we derive observationally from quantitative analysis of the spectra, but the meaning of their theoretical surface temperatures is more difficult to relate to what we derive from such studies. Here the fact the stars must be hydrogen free allows them to be easily identified. Uncertainty occurs when we begin to consider the lowest luminosity for a star to be a WR star. A low-mass and low-luminosity helium star will still eventually undergo core collapse but may not have a sufficiently high luminosity to drive optically thick winds and produce the characteristic emission WR lines, as discussed by Götberg et al. (2018). At solar metallicity this luminosity limit is $\mathrm{log}(L/{L}_{\odot })\sim 4.9$. Recently Shenar et al. (2020) suggested that this limit increases at lower metallicities due to the changing opacity of the stellar interiors, making it more difficult to drive stellar winds. This has not been taken into account in this study; however, Stanway et al. (2020) showed that including such a limit did decrease the number of stars observed as WR stars in a stellar population, which would improve agreement here between the BPASS models and the SMC numbers (see their Figure 4).

We note this is only important for binary populations because the single-star Geneva models only have the highest-mass stars becoming WR stars. With a binary population a much wider range of stellar masses become WR stars and without an accurate minimum WR star luminosity imposed, these models overpredict the number of WR stars at lower metallicities. If such a varying limit is correct, then there may be a number of products of binary interactions in low-metallicity environments that are difficult to identify.

How robust is our result to variations in the star formation rates? As emphasized earlier (Section 2.4), because of the difference in ages between RSGs and WRs, the expected RSG/WR ratio depends upon our assumptions about the constancy of the SFR. Let us briefly explore this, using as our example the BPASS models after adopting the spectroscopic temperatures and $\mathrm{log}L/{L}_{\odot }\gt 4.5$ for the RSGs. Figure 12 shows the effect of a steadily increasing or decreasing SFR over the past 40 Myr. For a 30% change in the SFR, the results are indistinguishable from the assumption of a constant SFR.

Zoom In Zoom Out Reset image size

Perhaps a more realistic approximation for changes in the SFH is shown in the upper panel of Figure 13. We simulated changing the SFHs using a Gaussian process, such that star formation is correlated on ∼5 Myr timescales and varies with an amplitude of 10%–40%. We used 200 random simulated SFHs; time is logarithmic to show the smooth variations on Myr timescales but are stochastic on longer timescales. When we applied this to the models, we obtained the comparison shown in the bottom panel of Figure 13. Although it broadens the expected ratios at a given metallicity, it does not change the results.

Zoom In Zoom Out Reset image size

Although the agreement with the models is very good, there are things we can learn from the discrepancies. First, in terms of the Geneva models, we find that the agreement is best when using an intermediate initial rotation rate. They adopted 40% of the breakup speed as their initial rotation based upon the observational work of Huang & Gies (2006) and Huang et al. (2010). But we actually expect there to be a distribution of initial rotation velocities as stars are formed. The agreement with the RSG/WR ratios suggests an intermediate value, but of course, these models do not include the effect of binarity. We will note that based on a more limited metallicity range than is currently available, Neugent et al. (2012a) also found that the relative number of WC and WN stars fell between the predictions of the rotating and nonrotating models (see their Figure 10) although this should be reexamined now that we have more complete data.

In terms of the BPASS models, we see really good agreement for most of the points at an intermediate initial binary fraction (about 60%–70%), except for the low- and high-metallicity examples (SMC and M31), which suggest a low binary frequency. Metallicity could plausibly affect the binary fraction in stellar populations through modifying the opacity of protostellar disks and hence the fragmentation of molecular clouds. This is challenging for hydrodynamical simulation to evaluate, and outcomes remain mixed in their predictions for any trend in binary fraction of massive stars with metallicity (see, e.g., Tokovinin & Moe 2020 for discussion). Or, it could be that (for the SMC), as noted above, our luminosity limit for when to include WR stars is too low.

Neugent (2021) found that the intrinsic binary fraction of RSGs is correlated with metallicity, with a higher binary fraction at larger metallicities; see also Dorn-Wallenstein & Levesque (2018). This is the reverse of the behavior seen in lower-mass solar-type stars, which are more likely to lie in binaries at low metallicities (Moe et al. 2019; Mazzola et al. 2020) and also show a dependence on α-enhancement, hinting at the complexity of disk fragmentation physics. This observed trend in RSGs could explain the SMC being consistent with a low binary fraction, but not M31. Note, too, that the observed binary fraction of WRs appears to be independent of metallicity (Bartzakos et al. 2001; Foellmi et al. 2003a, 2003b; Neugent & Massey 2014), although past mergers may affect those statistics.

Low apparent binary fractions at both ends of the metallicity range probed could hint that more than one aspect of the population is changing, with a different metallicity dependency. In this paper and in the context of the BPASS evolutionary models used, the term binary fraction is a simplification of a more complex parameter space and does not account for the variation in binary properties. In addition to an initial binary fraction, BPASS must also assume initial distributions in binary separation and mass ratio between the components. In BPASS v2.2 (used here), these are based on the empirical distributions of Moe & Di Stefano (2017), which have been calibrated on local stellar populations but have not been evaluated at high or low metallicities. A metallicity-dependent change in the parameter distributions favoring slightly wider binaries, or those with lower-mass ratios, would decrease the probability of interaction between binary components (rendering the stars effectively single) and thus have the same apparent effect on stellar number counts as a lower initial binary fraction. Further observations and modeling of binary populations in both high- and low-metallicity environments are required to determine whether such an explanation may plausibly explain the results for M31.

It should also be recalled in passing that while the BPASS models model binary interactions in detail, they do not include the evolutionary effects of higher-order multiple systems or of ongoing stellar rotation, which we know are also important, particularly at low metallicities.

Observationally, we can also make improvements. For instance, Shara et al. (2016) discovered a WR star in M31 that was not found as part of the Neugent et al. (2012a) study. The star is described as "heavily reddened," although no actual numbers are given. Is this star indicative of a missing population of WRs? We noted earlier that this WR was missing from the Neugent et al. (2012a) survey only because it fell within a gap between the mosaic CCD camera they used for their survey, but the question is still worth investigating. Neugent & Massey (2019) argues why a missing population of WRs is highly unlikely, but prompted by the current results, the first two authors plan a deeper survey in several spots in M31 to check out this possibility. (We note again that the RSG/WR ratio found using the entire sample and that found for the regions with the lowest extinction agree within the uncertainties.) Another improvement observationally would be to include additional galaxies, but here the options are rather limited. It would be practical to extend this work to include the galaxy NGC 6822, which has low metallicity (log(O/H)+12 = 8.25 or Z = 0.005; Pagel et al. 1980) and is known to contain four WRs (Massey & Johnson 1998). Yang et al. (2021) have identified RSGs in this galaxy using a unique method, but no studies have yet determined luminosities or temperatures for these objects; such a study involving several of the present authors is underway (T. Dimitrova et al. 2021, in preparation) However, other Local Group galaxies either have too few WRs to be statistically useful or else have other issues. For instance, IC 1613 has a very attractive low metallicity (log(O/H)+12 = 7.80 or Z = 0.0018, Bresolin et al. 2007), but contains just a lone WR, of WO type. Neugent & Massey (2019) cite this as an example of a WR that has almost certainly been formed by binary evolution, given the low-metallicity environment and the resulting expectations that stellar winds would be ineffective. ¹³ In contrast, the Local Group galaxy IC 10 contains several dozen WRs (Massey & Armandroff 1995; Massey & Holmes 2002; Tehrani et al. 2017) and has a metallicity similar to that of NGC 6822 (Garnett 1990). A preliminary accounting of IC 10's RSG population is given by Dimitrova et al. (2020). However, it is clearly undergoing a galaxy-wide starburst phase (e.g., Massey & Armandroff 1995; van den Bergh 2000). Thus, its RSG/WR ratio would not be a very useful test of models without a clearer understanding of the recent evolution of its SFR.

Using the RSG/WR value as a test of stellar evolution models is particularly challenging because of the difference in mass ranges leading to the RSG and WR stages. We have previously explored the use of the relative number of WC and WN stars both for the Geneva and BPASS models (see, e.g., Massey 2003; Neugent et al. 2012a; Eldridge et al. 2017; Stanway et al. 2020). However, one of the most useful observational tests has been the hardest to realize, namely the accurate measure of the relative number of blue and red supergiants, or, more specifically, the number of O-type stars to WRs (see, in particular, predictions in Dorn-Wallenstein & Levesque 2018, 2020; Stanway et al. 2020). The difficulties of a meaningful census have been discussed elsewhere (see, e.g., Massey 2010; Massey et al. 2017a), but observational progress continues to be made on that front.

In summary, the long-standing problem of the mismatch of the RSG/WR ratio and the evolutionary models has been resolved satisfactorily. The difference between this study and previous comparisons has been in the improved observational material rather than changes in the evolutionary models. We have utilized the complete surveys for WRs and RSGs that have become available over the years and also made use of comprehensive evolutionary models that include rotation (Geneva) and binarity (BPASS). We note that other modern evolutionary models are available, such as the MESA Isochrones and Stellar Models (Paxton et al. 2011, 2013, 2015; Choi et al. 2016; Dotter 2016), the FRANEC models (Chieffi & Limongi 2013; Limongi & Chieffi 2018), and the long-awaited Bonn evolutionary models (Szécsi et al. 2020). All of these include rotation, but not binarity. In general, we are unaware of detailed comparisons that have been made with these models with observations of massive stars, such have been carried out over the past decade for the Geneva and BPASS models. We would encourage their adherents to test these models against the RSG/WR ratios using the same methodology we have described here.

The lead author, P.M., would like to dedicate this paper to his former thesis advisor, Peter S. Conti, who recently passed away, and who would have enjoyed this study; it was also he who excitedly called the Maeder et al. (1980) paper to P.M.'s attention when it first appeared.

Lowell Observatory sits at the base of mountains sacred to tribes throughout the region. We honor their past, present, and future generations, who have lived here for millennia and will forever call this place home. Similarly, K.F.N., T.Z.D.-W., and E.M.L. would like to acknowledge that they work on the traditional land of the first people of Seattle, the Duwamish People past and present, and honor with gratitude the land itself and the Duwamish Tribe.

We are indebted to Georges Meynet for his thoughtful comments on an early draft of the manuscript. An anonymous referee also provided comments that helped us improve the paper. We also thank Maria Drout for help with the analysis of Gaia data in connection with earlier papers. This work was partially supported by the National Science Foundation (NSF) through AST-1612874. Although we do not include any new observations in the present paper, we wish to acknowledge that the two decades of work that have led to the present improvements were made possible thanks to the excellent observing facilities at Las Campanas Observatory (Magellanic Cloud WR surveys and spectroscopic followups), the Kitt Peak National Observatory (M31/M33 WR surveys), the MMT (M31/M33 WR spectroscopic followups), the United Kingdom Infrared Telescope and the Canada–France–Hawaii Telescope (M31/M33 RSG surveys). We are grateful to the generous and continuous support by the time allocation committees of the University of Arizona Observatories and the National Optical Astronomy Observatories. Similarly our knowledge of the RSG content of the Magellanic Clouds would not be possible without the 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 NSF.

Some of our illustrations were based on photographic data obtained using Oschin Schmidt Telescope on Palomar Mountain. The Palomar Observatory Sky Survey was funded by the National Geographic Society. The Oschin Schmidt Telescope is operated by the California Institute of Technology and Palomar Observatory. The plates were processed into the present compressed digital format with their permission. The Digitized Sky Survey was produced at the Space Telescope Science Institute under US Government grant NAG W-2166.

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.

Facilities: Swope (SITe No. 3 imaging CCD - , E2V imaging CCD) - , Magellan:Baade (MagE spectrograph) - , Mayall (Mosaic-1 wide-field camera) - , MMT (Hectospec multi-fiber spectrograph - , Blue Channel) - , UKIRT (WFCam NIR wide-field camera) - , CFHT (WIRCam NIR wide-field camera) - , CTIO:2MASS. -

In order to identify the RSGs in the SMC, we used the same technique Neugent et al. (2020a) recently applied to the LMC. We began by selecting objects from the VizieR version of the 2MASS point-source catalog (Cutri et al. 2003). We used a search radius of 3° centered on α₂₀₀₀ = 01: 08: 00.0 δ₂₀₀₀ = −73:10:00 in order to match the area searched for WRs (Neugent et al. 2018); see Figure 2. We used only the stars with photometric quality flags of "AAA" and artifact contamination flags of "000." We restricted the sample to stars with K_s ≤ 13 and J − K_s ≥ 0.5, much fainter and bluer than we expect to find RSGs at the distance of the SMC (59 kpc; van den Bergh 2000), as shown below. The resulting CMD is shown in Figure 14(a).

Zoom In Zoom Out Reset image size

As emphasized in earlier papers, one of the concerns is contamination by foreground stars. We therefore used Gaia proper motions and parallaxes to separate SMC stars from those in the foreground. The procedure is described in detail in Neugent et al. (2020a). As in that study, we retain stars with ambiguous results or incomplete Gaia data. Owing to the timing of our work, we relied primarily on the Gaia Data Release 2 (DR2, Gaia Collaboration et al. 2018), but as mentioned below, in some cases we also consulted Data Release 3 (DR3, Gaia Collaboration et al. 2021). Out of 15,005 stars in the photometry list, 9419 (62.8%) are considered to be certain members, 4903 (32.6%) are foreground, 415 (2.8%) have ambiguous results, and 268 (1.8%) lack Gaia complete Gaia data. We show the effects of removing the foreground stars in Figure 14(b). As in the LMC, the vast majority of the contaminants are at the warmer temperatures, outside the region we will use to define the RSGs.

Of greater concern is the separation of RSGs and AGBs. As discussed above, Yang et al. (2019) showed that RSGs and AGBs could be readily separated in a (J − K_s , K_s ) CMD, and we have followed that procedure here. Figure 15 shows the distribution of stars in such a CMD after foreground stars have been removed, where we have indicated the region we have adopted for the RSGs, as well as identifying the AGB sequences defined by Boyer et al. (2011). We list the photometric criteria adopted in Table 4.

Zoom In Zoom Out Reset image size

Table 4. Adopted and Derived Relations for the SMC

	Relation	Source
Adopted Distance:
	59 kpc (DM = 18.9 mag)	1
Reddening Relations:
	AK = 0.12AV = 0.686E(J − K)	2
	E(J − K) = AV /5.79	2
RSG Photometric Criteria:
	10.6 < Ks ≤ 13.5: Ks ≥ Ks0 and Ks ≤ Ks1	3
	Ks ≤ 10.6: J − Ks ≥ 0.812 and Ks ≤ Ks1	3
	Ks ≤ 9.0 and (J − Ks ) ≤ 1.8: J − Ks ≥ 0.812	3
	Ks0 = 24.00 − 16.50(J − Ks )	3,4
	Ks1 = 27.50 − 16.50(J − Ks )	3,4
Adopted Extinction:
	AV = 0.75	3
Conversion of 2MASS (J, Ks ) to Standard System (J, K):
	K = Ks + 0.044	5
	J − K = (J − Ks + 0.011)/0.972	5
Conversion to Physical Properties (Valid for 3500–4500 K):
	Teff = 5592.5 − 1656.0(J − K)0	3
	BCK = 5.495 − 0.73697 × Teff/1000	3
	K* = K − AK	⋯
	Mbol = K* + BCK − DM	1
	$\mathrm{log}L/{L}_{\odot }=({M}_{\mathrm{bol}}-4.75)/-2.5$	⋯

References. (1) van den Bergh (2000), (2) Schlegel et al. (1998), (3) This paper, (4) Cioni et al. (2006a), (5) Carpenter (2001).

Download table as:  ASCIITypeset image

Finally, we need to covert this photometry to the physical properties of effective temperatures and luminosities. As discussed above, we rely upon the marcs models (Gustafsson et al. 1975, 2008; Plez et al. 1992) computed with surface gravities appropriate for RSGs ($\mathrm{log}g\sim 0.0$) as described in Levesque et al. (2005). Following Neugent et al. (2020a, 2020b), we found a simple functional relationship between effective temperatures and the dereddened (J − K)₀ colors. Note that the model colors were computed on the standard Bessell & Brett (1988) system and that we must first transform our 2MASS colors to this system using the relationships described by Carpenter (2001), as given in Table 4. The previously analyzed RSGs have average A_V values of 0.75 (see Levesque et al. 2006 and discussion in Neugent et al. 2020a, 2020b), and we used this value to deredden the transformed colors, using the relationships give by Schlegel et al. (1998). The exception was a handful of stars that are too bright to be AGBs, but redder than most of the RSGs; as in previous papers we assume these are RSGs with extra circumstellar reddening (Massey et al. 2005) and correct for extinction by following the reddening line back to the mean relationship between J − K_s and K_s . As noted in earlier papers, the luminosity we derive is actually quite insensitive to the reddening, as the extinction in K is only 12% that in V. The relationships and transformations are all summarized in Table 4. The relationships between our color criteria, and the effective temperatures, have already been given in Table 1. The typical uncertainty in our values are 150 K in T_eff and 0.05 dex in $\mathrm{log}L/{L}_{\odot }$. Of course, as described above, we suspect that the actual "evolutionary" temperatures are 200 K cooler than these.

The final source list of SMC RSGs is given in Table 5. There are 1741 RSGs brighter than our K_s = 13.0 limit. This corresponds to a completeness limit of $\mathrm{log}L/{L}_{\odot }=3.7$, much less luminous than we are concerned with here.

Table 5. RSG Source List for the SMC

											Spectral Type
2MASS	αJ2000	δJ2000	Ks	σKs	J − Ks	σJ−Ks	Gaia a	Teff b , c	$\mathrm{log}L/{L}_{\odot }$ b , d	Type	Ref.	Other ID	Comments
J00450138–7350513	00 45 01.383	−73 50 51.34	11.170	0.019	0.918	0.032	0	4200	4.06	⋯	⋯	⋯	⋯
J00450313–7255156	00 45 03.135	−72 55 15.61	10.314	0.023	0.935	0.033	0	4200	4.39	G8 Ib	5	[GDN2015] SMC079	⋯
J00450456–7305276	00 45 04.566	−73 05 27.69	8.071	0.023	1.137	0.033	0	3850	5.19	M2 I	1	[M2002] SMC 005092	RV var Ref 9
J00450482–7340132	00 45 04.824	−73 40 13.28	12.209	0.027	0.791	0.036	0	4450	3.71	⋯	⋯	⋯	⋯
J00450746–7327417	00 45 07.464	−73 27 41.79	10.300	0.025	1.008	0.035	0	4050	4.36	⋯	⋯	⋯	⋯
J00450758–7358102	00 45 07.581	−73 58 10.21	10.868	0.025	0.923	0.035	0	4200	4.18	G7 Iab-Ib	5	[GDN2015] SMC080	⋯
J00450816–7325382	00 45 08.165	−73 25 38.22	12.131	0.029	0.849	0.038	0	4350	3.71	⋯	⋯	⋯	⋯
J00450831–7254146	00 45 08.312	−72 54 14.64	12.273	0.032	0.920	0.042	0	4200	3.62	⋯	⋯	⋯	⋯

Notes.

^a Membership flag based on Gaia information: 0 = probable member; 1 = uncertain member; 2 =i ncomplete Gaia data. ^b Computed assuming Av = 0.75 mag. ^c Typical uncertainties 150 K. ^d Typical uncertainties 0.05 dex.

References. (1) Levesque et al. (2006), (2) Massey et al. (2007), (3) Levesque et al. (2007), (4) Levesque et al. (2014), (5) Dorda et al. (2018), (6) This paper and references therein, (7) Neugent et al. (2020a), (8) Neugent et al. (2019), (9) Dorda & Patrick (2021), (10) Patrick et al. (2020), (11) Elias et al. (1980), (12) Beasor et al. (2018).

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

Download table as:  DataTypeset image

In Table 5 we have included spectral types for a variety of sources. Levesque et al. (2006) reclassified a number of the RSGs previously identified by Elias et al. (1985) and Massey & Olsen (2003) and fit their spectra to refine the effective temperature scale and determine luminosities. Subsequently, Massey et al. (2007) and Levesque et al. (2007) identified 12 highly unusual RSGs in the Magellanic Clouds, five of which were in the SMC. These stars underwent large swings in spectral types, changing from early K type to M4 or later on a timescale of a few years. In their cool state, these stars are much later in type than expected for their host galaxies. In their search for candidate Thorne–Żytkow objects, Levesque et al. (2014) examined these and other stars with similarly high J − K_s colors for abundance anomalies, leading to their finding of HV 2112 as the first such candidate. When we cross-referenced our list with the recent spectroscopy study of Dorda et al. (2018) we noticed two additional possible spectrum variables. The star 2MASS J01093824–7320024 ([M2002] SMC 069886) was classified as K5-M0 I by Levesque et al. (2006), but as a very late M-type star (M4 Iab) by Dorda et al. (2018). It was subsequently classified as M2 I by Neugent et al. (2020a). The star 2MASS J00485182–7322398 ([M2002] SMC 011939 may be a similarly case, as Levesque et al. (2006) called it a K2 I while Dorda et al. (2018) classified it as M2 Ia-Iab.

The main discrepancy we noted between researchers was that Dorda et al. (2018) classified a number of stars in our list as late G type. Four of these had previously been classified by Levesque et al. (2006), whose classifications range from K1 I to K2–3 I. Here we merely note this systematic issue; resolving this is well beyond the scope of the present paper. We will further note that one of the stars (2MASS J00382421–7410196) that was originally in our RSG list was classified by Dorda et al. (2018) as an early M dwarf (M1 V). Its SMC membership was ambiguous from DR2. Comparison with the improved astrometry from DR3, however, confirms it is likely foreground, and we have eliminated it both from our table and from our counts.

We also included spectral types by Neugent et al. (2019, 2020a) as part of their search for RSGs with binary companions, and noted if radial velocity variations were detected by Dorda & Patrick (2021). Spectral types of the recently studied stars in NGC 330 by Patrick et al. (2020) were also included. We note that cross-references to names and previous data are not intended to be complete; however, the 2MASS designations can readily be used in SIMBAD to find additional information.

We noticed that very few spectroscopically confirmed RSGs are missing from our list. 2MASS 01000932–7208441 ([M2002] SMC 48122) was classified by Levesque et al. (2006) as K1 I, but its 2MASS photometry is flagged as unreliable in J and K_s . Although the RSG+B binary 2MASS J00464984–7313525 ([M2002] SMC 7618) was originally identified as probable foreground from the Gaia DR2 astrometry, DR3 fixed the problem as the parallax went from 0.99 ± 0.25 to 0.13 ± 0.08, and we have retained the star in our source list. Finally, 9 of the 16 NGC 330 RSGs classified by Patrick et al. (2020) are missing, as crowding in the extreme center of the cluster compromised the 2MASS photometry. By contrast, there are 260 stars in Table 5. Thus, the requirement that we accept only the best 2MASS photometry in identifying RSGs in the SMC has had a very minimal effect on the statistics.