Red Supergiant Candidates for Multimessenger Monitoring of the Next Galactic Supernova

We compile a catalog of 578 highly probable and 62 likely red supergiants (RSGs) of the Milky Way, which represents the largest list of Galactic RSG candidates designed for continuous follow-up to date. We match distances measured by Gaia DR3, 2MASS photometry, and a 3D Galactic dust map to obtain luminous bright late-type stars. Determining the stars' bolometric luminosities and effective temperatures, we compare to Geneva stellar evolution tracks to determine likely RSG candidates, and quantify contamination using a catalog of Galactic AGB in the same luminosity-temperature space. We add details for common or interesting characteristics of RSG, such as multi-star system membership, variability, and classification as a runaway. As potential future core-collapse supernova (SN) progenitors, we study the ability of the catalog to inform the Supernova Early Warning System (SNEWS) coincidence network made to automate pointing, and show that for 3D position estimates made possible by neutrinos, the number of progenitor candidates can be significantly reduced, improving our ability to observe the progenitor pre-explosion and the early phases of the core-collapse supernova.


INTRODUCTION
Red Supergiants (RSGs) are one of the last evolutionary stages of massive stars, and have been directly linked to Type-IIP/L supernovae (SNe) through matching SN locations with archival pre-explosion images (Smartt et al. 2004(Smartt et al. , 2009;;Smartt 2015;Davies & Beasor 2020;Beasor & Davies 2017;Li et al. 2007;van Dyk et al. 2017;Van Dyk et al. 2012); a comprehensive summary of can be found in Van Dyk (2017); Smartt (2015).The connection between RSGs and core-collapse supernovae (CCSNe) links two stellar phases between a hugely disruptive explosion, which enables studies of stellar evolution as well as the mechanisms that cause CCSN.Furthermore, as progenitors of SNe, RSG properties impact the formation of stellar mass black holes and neutron stars, themselves sources of gravitational waves through possible eventual mergers.
The next Galactic CCSN will represent a once-in-a-generation opportunity to study important progenitor-SN connections in unprecedented detail.Such CCSNe are rare, but they allow us to study neutrinos (e.g.Scholberg 2012), gravitational waves (GW; e.g.Ott 2009;Kotake 2013), and nuclear gamma rays (e.g.Gehrels et al. 1987;Horiuchi & Beacom 2010).These observables allow us to ★ E-mail:healys@vt.eduprobe far more completely inside the stellar photosphere.The neutrinos are key since they are emitted from the proto-neutron star deep within the progenitor.They are also emitted before the SN starts and the neutrinos can be used to point to the location in the sky with an error circle of a few to ten degrees (Beacom & Vogel 1999;Bueno et al. 2003;Tomàs et al. 2003;Mukhopadhyay et al. 2020).These make neutrinos ideal as an early warning trigger for the SN and is the backbone of multi-messenger observation strategies for CC-SNe.In addition, the final stage of Si burning in massive stars emits neutrinos with ∼ 5 × 10 50 erg (Arnett et al. 1989), which can be detected for progenitors within a few kpcs (Odrzywolek et al. 2004;Kharusi et al. 2021;Asakura et al. 2016), acting as another layer of early trigger.Various aspects, including predictions, implementations, and detectability of the multi-messenger signals of Galactic and nearby CCSNe have been explored over the years (Adams et al. 2013;Nakamura et al. 2016;Kharusi et al. 2021).
To most effectively capture this opportunity, awareness of the properties and spatial distributions of RSGs is crucial.With such knowledge, the rapid identification of the SN progenitor would become more realistic, e.g., mirroring the strategies in the gravitational wave community with lists of galaxies for follow-up searches (White et al. 2011; and more recently Dálya et al. 2018).The effective use of a pre-compiled target list would allow, among others, monitoring of the earliest light curves of the progenitor/SN, which are crucial for accurate reconstruction of the CCSN evolution and determining progenitor properties (e.g., Tominaga et al. 2011).
However, our knowledge of RSGs remains incomplete.In particular, the spatial and luminosity distributions are not well known.Spectra from massive stars are limited, and even more so are those with assigned spectral types.While the temperature scale of Galactic RSGs has been a topic of study for decades (Johnson 1964(Johnson , 1966;;Flower 1975Flower , 1977;;Lee 1970), and has been updated in a series of recent papers- Levesque et al. (2005) and Levesque et al. (2006) compared the strength of the TiO bands in the optical to that of MARCS's model atmospheres-the debate about the precise temperatures of Galactic RSGs is still ongoing (Taniguchi et al. 2021;Taniguchi & Winered Team 2021).It is also known that the most evolved RSGs have substantially higher levels of circumstellar extinction, but the dust sizes and distributions are not well modeled.Since dust estimations directly affect luminosity calculations, this can lead to errors in a star's bolometric luminosity, mass, and stellar radius (Beasor & Davies 2016).
While the extensive information that can be obtained from Galactic CCSNe cannot be matched by extragalactic ones, it is worth noting that RSG samples in the Magellanic clouds, M31, and M33 are estimated to be near-complete.The number of known RSGs in the LMC and SMC was initially a few tens (Feast et al. 1980;Catchpole & Feast 1981;Wood et al. 1983;Pierce et al. 2000) and increased to a couple hundred (Massey 2002;Massey & Olsen 2003;Yang & Jiang 2011, 2012;Neugent et al. 2012;González-Fernández et al. 2015) but the more recent work of Yang et al. (2019Yang et al. ( , 2021) ) identified an approximately 90 percent complete sample with 1405 (SMC) and 2974 (LMC).Similarly, work done in M31 and M33 initially selected 437 (M31) and 776 (M33) candidates (Massey et al. 2006(Massey et al. , 2007(Massey et al. , 2009;;Drout et al. 2012) with 255 of those confirmed to be RSGs by Massey & Evans (2016).However, Ren et al. (2019) found a lack of known faint stars motivating the identification of a near-complete sample containing 5498 (M31) and 3055 (M33) RSGs (Ren et al. 2021).
Mapping the Milky Way's RSGs has been attempted in the past.For example, those compiled with a greater focus on pre-explosion studies, including Messineo & Brown 2019 with 889 late-type bright stars of which ∼382 are RSG candidates, and Messineo (2023) with 203 bright late-type stars of which 20 are candidate RSGs.Also, lists exist that are focused more on post-explosion studies, including Nakamura et al. 2016's 212 stars andMukhopadhyay et al. 2020's 31 stars.We build on the methods of Messineo & Brown 2019 but design our methodology to minimize losing a true RSG even at the expense of keeping more non-RSG contaminates.To further augment our list, we include a second method for search.We also compile or estimate additional stellar characteristics, e.g., variability, radius, multi-star system, cluster membership, runaway status, presence of a significant magnetic field, mass-loss rate, and classification as a runaway star.We present a systematic list of Galactic binaries.Finally, we look into our final samples' angular separation compared to the angular resolution anticipated from future CCSN neutrino detections in preparation for use in multi-messenger astronomy.Further discussion of Galactic RSG surveys in the literature and comparisons to our work is in Section 5.
A comprehensive list of evolved massive stars in the Milky Way galaxy will have numerous broader benefits beyond multi-messenger astronomy.For example, it can be used to study the effects of massive star evolution, binary influence, as well as metallicity.RSGs trace stars with masses from about 9 to 40 M ⊙ , i.e., ages of 4 to 30 Myr (Ekström et al. 2012;Chieffi & Limongi 2013).The binary frequency of RSGs can be used to test interaction rates when compared to unevolved massive stars.At the point of writing this paper, approximately a dozen Galactic RSG binaries are confirmed (Neugent et al. 2019).The main factor in making these improvements is having a relatively complete Galactic RSGs population that also spans various metallicities.
This paper is organized as follows.In Section 2, we compile our list of bright and cool stars and measurements of their distance, photometry, and dust extinction.We determine our most confident sample using luminosity, stellar evolution tracks, and note stars with typical RSG characteristics in Section 3. In Section 4, we discuss specific objects of interest, the sample's spatial distribution, the range of estimated radii, and mass loss.The implications for the list and coordinated use with multi-messenger astronomy (MMA) for Galactic CCSNe and massive star astronomy are discussed in Section 5 6.We conclude in Section 7.

METHODS
We aim to determine whether a star is a RSG by using its effective temperature and bolometric luminosity.We first collect a sample of luminous late-type stars, which we define by spectral type K or M and luminosity class I. Based on spectral types, T eff is estimated with the temperature scale for RSGs determined in Levesque et al. (2005).From there, we collect photometry measurements, model dust extinction, and estimate the intrinsic bolometric magnitude and T eff , which are crucial to determining whether the star is more likely an RSG or contamination from a non-RSG.

Data Selection
We consider two data selection methods, each designed to create a starting list of stars with RSG-like properties.The first "compilationbase" method focuses on selecting stars with predetermined spectral types matching the expectations of RSGs, i.e., a collection of K-M types stars with luminosity class I.However, spectral classifications of Galactic massive stars are relatively small in size.To remove this limitation, we use a second "Gaia-based" method, which focuses on obtaining stars whose position on the HR (or pseudo HR) diagram matches the expectations of RSGs.To this end, we utilize Gaia photometry and color to place a cut on Gaia DR3 stars to obtain stars that compare to K-M stars with luminosity class I.We describe in detail the two methods below and summarize total counts in Table 1.
In an effort to be as complete as possible, we also included catalogs that were compilations of others (Humphreys 1978;Jura & Kleinmann 1990).Skiff (2014) is the largest such catalog, containing 1,058,791 entries, and not all entries can be found in Simbad 1 .Entries in the catalog have been converted to the MK type where possible.
We restrict ourselves to K and M type stars with at least one classification of or including luminosity class I.We do not restrict the inclusion of spectra based on location in the sky.Special care was made to ensure that if there were repeats in original source classifications between these catalogs or with the earlier list, only one entry was included.A portion of the list can be found in Table 2.
Not all observations were unique and a single star could have a dozen classifications, so we implemented a method for adopting spectral types.For catalogs with well-studied objects, we adopted spectral types based on the classification assigned in those catalogs.We gave preference to Levesque et al. (2005) as they used both moderate-resolution optical spectrophotometry and the MARCS stellar atmosphere models fits to determine the spectral type.Thereafter, we prioritize the catalogs with well-studied classifications in Dorda et al. (2016a), Dorda et al. (2018), Jura & Kleinmann (1990), Elias et al. (1985), andHumphreys (1978).
If none of the previous works provide spectral types, we adopt a mean spectral type weighted by uncertainties.This keeps in mind that classifications derived from fewer than the 8 TiO bands or from bands that have less characteristic behavior for massive stars and those who have observational limitations, like Keenan & McNeil (1989), generally have larger uncertainties between ± 1-2 subtypes (usually ± 0.5 subtypes for classification done by visual inspection of the whole 8 TiO band spectrum).We include White & Wing (1978), who, along with stating previous spectral classification, finds the mean spectral type from the 8-color indices, which use a photoelectric system defined by eight narrow bands between 0.7 and 1.1 m to determine two-dimensional spectral classification.We do not adopt 2-D spectral types but use them to inform mean values for entries with many classifications.This is accomplished by converting K0-M10 to a numerical scale, calculating the weighted average, and then rounding to the nearest half-integer, except for special cases with 4 or more entries where a single spectral type appears a majority of the time.
Initial celestial positions were taken directly from the source paper, 1 These stars are still retained in our analysis.but after combining all catalogs, some variance in RA and Dec was observed, partially due to the different epochs used.In order to improve the positions, we matched the stars listed in Table 2 to Simbad based on the alias of the star and took the corresponding position with epoch J2000.0, also shown in Table 2.In summary, our compilation-based Method provided 2,051 stars with either K or M spectral type and at least one classification of luminosity class I.

Method 2: Gaia-based sample
In an effort to expand our list further, we utilized the vast stellar database of Gaia (DR2/DR3) and a subset of our compilation-based sample to select RSG candidates.We use Gaia's G passband (roughly 330-1050 nm) magnitude and BP-RP (330-680 nm; 630-1050 nm) colour band for rough estimates of luminosity and effective temperature, respectively.While the absolute G magnitude does not provide total bolometric luminosity, we use it to indicate where a star would be if we had included all magnitudes from other bands; the Gaia G band works well for this purpose as it spans near-ultraviolet to near-infrared.Larger differences between passbands correlate with cooler stellar temperatures.The star's temperature can be inferred as the larger the difference between passbands, the cooler the star.These are shown in Fig. 1.Here, the color indicates the log of stellar density.
Table 2. Summary of late-type bright stars compiled from both method whose distance uncertainty meet the requirements detailed in §2.2.For stars from the compilation-based method, we list all spectral types recorded along with the mean adopted spectral type.Spectral types from (Houk 1978) for stars in the Gaia-based sample are included with the adopted spectral type set equal to the spectral type.Alias, Gaia DR3 IDs, RA, Dec, parallax, and distance are included along with the T eff determined using the adopted spectral type.The full table is publicly available at https://github.com/SNEWS2/candidate_list.Davies et al. (2007), 25 Davies et al. (2008), 26 Keenan & McNeil (1989), 27 Alexander et al. (2009), 28 Bidelman (1957b), 29 Bidelman (1957a), 30 Halliday (1955), 31 Medhi et al. (2007), 32 Ginestet et al. (1999), 33 Ginestet et al. (1997), 34 Fawley & Cohen (1974), 35 Garmany & Stencel (1992), 36 Humphreys (1970a), 37 Humphreys et al. (1972), 38 White & Wing (1978), 39 Marco & Negueruela (2013), 40 Boulon (1963), 41 Westerlund (1987), 42 Marco et al. (2014), 43 Currie et al. (2010), 44 González-Fernández & Negueruela (2012), 45 Humphreys (1970b), 46 Figer et al. (2006)  4 r_med_geo: geometric and photometric distance estimate from Gaia DR3's parallaxes (Bailer-Jones et al. 2021) We used a subset from our compilation-based sample to identify the region in the colour-magnitude diagram (CMD) corresponding to late-type bright stars.For this purpose, we used stars whose RSG status is more confident, either low uncertainty or how well-analyzed it is.This includes four Galactic clusters: RSGC1, RSGC2 (Stephenson 2), RSGC3, and NGC 7419, which cover the range of expected luminosity as well as include masses ≃ 9M ⊙ to ≲ 25M ⊙ .These clusters are some of the most well-studied in terms of RSG population; most stars included in these regions have follow-up spectra, allowing us to restrict our set further to those that are or have a follow-up observation.Ultimately we use 156 stars compiled from Levesque et al. (2005), Humphreys (1978), Figer et al. (2006)  The reference RSGs were matched to Gaia DR3 (DR2) objects with both G magnitude and estimated extinctions, which reduced the size to 30 (52) stars.The DR2 matches are shown by red stars in Fig. 1 and reveal the region where RSGs populate.The strong bimodality of the RSGs on the CMD presented in Fig. 1 likely results from the nature of stellar metallicity distribution and its effects on stellar evolution.Studies of the metallicity of the horizontal branches of globular clusters (using Gaia DR2 data; Gaia Collaboration et al. 2018b andZinn 1985) and halo stars (using Gaia DR2 data; Gaia Collaboration et al. 2018b and TGAS data with RAVE and APOGEE;Bonaca et al. 2017) show the distribution is double-peaked.It is also known that the evolutionary tracks of stars are highly dependent on metallicity (El Eid et al. 2009;Langer 2012;Maeder 2009).The double-peaked metallicity distribution of stars before the RSG phase and the impact metallicity has on stellar evolution suggests the distribution of stars in the CMD is similarly bimodal as seen in Fig. 1.

Alias
Including both regions, we derive a cut to select Gaia stars in the CMD that are likely to be RSGs.While different cuts can be found for different samples, we choose to use the cut derived from DR2 for both Gaia datasets as few of our reference RSGs have estimated extinction in DR3.This limitation is primarily driven by the nature of bright stars and the high uncertainties they would likely have when methods for determining extinction values are based on training sets that exclude massive stars, as is the case for Gaia DR2 and DR3 (Delchambre et al. 2023;Bailer-Jones et al. 2013;Andrae et al. 2022).As our cut determined from Gaia DR2 is derived from absolute G mag and BP-RP colours, it is applicable to any dataset containing those measurements.It includes all in the DR3-matched reference set.Our data cut is: and shown by a thick solid line in Fig. 1.Stars above the cut are bright and cool enough to be candidates RSGs. Figure 1 shows how all our reference RSGs are above the cut.All stars are retained regardless of uncertainty; we discuss the accuracy of distance in §2.2.The cut makes a homogeneous sample in terms of extinction and magnitude.
Still, those values are not used for later analysis, as our goal is to determine bright late-type stars homogeneously.The datasets were merged, preserving only unique values.
Our cut in the CMD gives a degree of confidence in the luminosities of selected stars.Thus, we only need to take stars with appropriate spectral types.However, our CMD provides only rough estimates for spectral types.Thus, we obtained spectral types from the Henry Draper (HD) catalogs.Over multiple publications, the HD catalogs (1918-1924), Henry Draper Extension (HDE;1925-1936), andHenry Draper Extension Charts (HDEC;1934) determined spectral classification and rough positions for 272,150 stars.These spectral classifications are based on the Harvard system, which has slightly different notations for spectral types than the MK system.For consistency, we use the Michigan Catalog of HD (MCHD) , which contains ∼161,000 HD stars that were reclassified to the MK system using the spectra from HD, HDE, and HDEC (Houk 1978).
Matching was done for all stars to available M and K spectral type objects from MCHD providing an additional, though not entirely unique from the compilation-based method, set of 3,499 candidates.
The overlap with the compilation-based method is <100 stars, reinforcing that our Gaia-based method provides us with valuable new candidate RSGs.

Distance related cuts with Gaia
We utilized the Gaia DR3 catalog, which contains over 1.8 billion sources (Gaia Collaboration et al. 2016Collaboration et al. , 2022) ) to provide reliable parallaxes for both the compilation-based and Gaia-based samples, including uncertainties for each star.After converting coordinates from J2000 to J2016, Gaia matches were found within a 2" radius of each late-type star, producing objects for ∼ 96% of the compilationbased method and, by design, 100% of the Gaia-based method.
The Renormalised Unit Weight Error (RUWE), described in the technical note (Lindegren, L. et al. 2018), is recommended by the Gaia team as a proxy for a good astrometric solution to a star.The unit weight error (UWE) dependence on the magnitude and color of the source means that for consistent estimation, it needs to be renormalized depending on color and magnitude.RUWE is less certain for objects without color information.Lindegren, L. et al. (2018) recommends using RUWE ≤ 1.4 as a 'good' solution.However, since our list is entirely made of bright objects, they will inherently have more noise, so we choose a cut of RUWE ≤ 2.7 (see also Messineo & Brown 2019).
The Bayesian distance estimator R_med_geo is the recommended package for Gaia stellar distances (Luri et al. 2018;Bailer-Jones et al. 2021).Originally proposed in Bailer-Jones ( 2015) and further studied in Astraatmadja & Bailer-Jones ( 2016), estimating distances with Gaia becomes an inference problem that is best solved using a simple exponential decreasing space prior.With this method, non-positive parallaxes can still be used, and a bias correction is not required.We use observed parallaxes when estimating the galactocentric coordinates as other position information, such as proper motions, are considered, as suggested in Gaia Collaboration et al. (2018a).
Furthermore, distance errors are crucial for us to separate RSGs from impostors.The Gaia collaboration defines two errors: external which includes both random and systematic errors, and internal which is the formal errors reported in Gaia DR3 but only quantify the consistency of measurements 2 .The former is the total error, while the latter is not.We calculated the external parallax uncertainties as recommended by applying Eq. ( 2) to the internal uncertainties found in DR3.The external error of the parallax,   , is defined as where depending on the G mag, the values of  and   are either  ≲ 13 :  = 1.08,   = 0.021 or  ≳ 13 :  = 1.08,   = 0.043.With these calculations, we look for relative error, /  , and RUWE such that we retain only objects whose errors will not cause significant uncertainties in later derived quantities.Previous cuts from Messineo & Brown (2019) restricted relative error > 4 to keep the difference between 1/( −  0 ) 3 and r_med_geo to less than 5% throughout the entire sample.In an effort to maximize final numbers while balancing quality, comparisons were made between different restrictions, shown in Fig. 2. Allowing objects whose relative parallax error > 2 produces a sample whose maximum difference between distances derived from parallax and r_med_geo is > 5%, but on average, the difference is still ∼ 0.7% with a median value of ∼0.2%, yet almost doubles our sample size in the most favorable regions while preventing the exponentially decreasing space prior to influence r_med_geo to the degree that it would be dependent on parallax and the Galactic model used by Bailer-Jones et al. ( 2018).4,863 unique values are retained after applying the distance cut.

Photometry
We utilize near-infrared photometry from the Two Micron All Sky Survey (2MASS).2MASS finished collecting 25.4 tbytes of data in February 2001, including raw images of 99.98% of our celestial sphere.The bands measured were the near-infrared (NIR) J (1.25 m), H (1.65 m), and K s (2.16 m) of which errors were calculated and produced for the data set which includes images from almost the entire sky.All images were taken through either Mount Hopkins, Arizona's and Cerro Tololo, Chile's 1.3m diameter telescopes both with 7.8s time accumulated for each sector (Skrutskie et al. 2006).
Using ESA pre-crossmatched best neighbour and good neighbourhood list, available 2MASS Point Source Catalogue (Skrutskie et al. 2006) matches were found for > 99% of stars that had a Gaia DR3 ID, providing a near-complete sample of JHK s measurements and their associated uncertainties.The K s magnitudes range from -2 to 16.3 mag with a median value of 4.3 mag.Where available, we also include 2MASS values for the blue, visual, or red magnitude of the associated optical source.

Dust
In Messineo et al. (2005), dust extinction of RSGs was calculated by using the extinction power law with an index of 1.9 and the intrinsic colors expected for each spectral type.However, RSGs in general show excess circumstellar extinction to those of OB stars in the same regions (Massey et al. 2005), causing the extinction ratios to differ from those of the main sequence or other stars.This increase in extinction for brighter massive stars is also seen in Neugent et al. (2020), which takes direct measurements of   for a sample of RSGs in M31.
While a linear relationship between intrinsic color, ( − ) 0 , and effective temperature works well for RSGs, metallicity in the models causes discrepancies as well as unlikely large extinction values.We also have to take into account that a large portion of our sample's  eff is based on a mean spectral type.Deriving dust based on  eff and intrinsic magnitude means not only would we need to propagate the uncertainties of  eff and observed photometry, but we would also have to deal with the dependence between the two when calculating the uncertainties of intrinsic magnitudes.This method would keep us from treating our uncertainties as uncorrelated, complicating our error propagation.For these reasons and those described in the next section, we use Galactic dust maps to estimate reddening while keeping the extinction ratios from Messineo et al. (2005).

3D Dust Map
In recent years, the development of interstellar dust extinction maps has improved to 3 dimensions spanning many kiloparsecs.We use the combination mwdust4 created by Bovy et al. (2016) and later updated to replace Green et al. (2015) with Green et al. (2019).While there are others based on modeling stellar photometry (e.g., Marshall et al. 2006;Sale et al. 2014;Green et al. 2019), which would be ideal for our use with 2MASS photometry, none provide results for the entire sky that we need.
The combination of Bovy et al. (2016) starts with Marshall et al. (2006) based on 2MASS passbands and augments with Drimmel et al. (2003) based on dust distribution fit to COBE DIRBE data.Bovy et al. (2016) unified the projection of Marshall et al. (2006), Green et al. (2015), and Drimmel et al. (2003) to HEALPix; note however the resolutions are variable.Collectively, Bovy et al. (2016) provides full sky coverage and preserves small-scale structures in the dust extinction.
As our sample focuses on RSGs, we modified some constants used to transform between bands which are used in the maps so that they fit the extinction transformation equation found in Messineo et al. (2005), specifically where R V = 3.1.Massive stars, however, can have additional attenuation due to local extinction.For example, studies have found values of 3.6 (Lee 1970;McCall 2004)  In the absence of excess circumstellar medium (CSM) or other peculiar reddening law, 3.1 still gives good agreement with near-UV for Galactic RSGs (Massey et al. 2005).Additional dust is a major source of systematic uncertainty, but we can obtain some estimates of their impacts.In the sample of Galactic RSGs from Levesque et al. 2005, approximately 15% of those with estimations for V band extinction from spectrophotometry differ significantly-here quantified as impacting our classification of a star as a RSG or not (see §3.3)-from estimations based on R V = 3.1.While their use of moderate-resolution SEDs significantly reduces the chance of discounting high-extinction RSGs, dust attenuation in the NIR is less than that of the optical or lower wavelengths, minimizing our sensitivity to extinction estimations and improving our ability to retain RSG with excess reddening.For these reasons, we conservatively use R V = 3.1, giving us a final sample of stars with either approximately correct or underestimated luminosity.While this makes our sample biased against RSGs with extra extinction, any higher value of R V would only increase the luminosity of stars and hence increase our sample size.
As Marshall et al. (2006) uses A K s = 0.089A V , we make the appropriate correction to match the above.The Green et al. (2019) map is intended to provide reddening in a similar unit as E(B−V) SFD (Schlegel et al. 1998) which can be covered to E(B−V) by multiplying 0.884, the constant found in Schlegel et al. (1998) and then moving through the equations above to get A Ks .While Drimmel et al. (2003) is normalized to E(B−V) SFD , mwdust appropriately converts to A V .
Since mwdust does not provide uncertainties, we estimate them by comparing different maps.We used dustmaps, a python package that can be used to query several commonly used maps of interstellar dust, including 2D maps such as Schlegel et al. (1998), Planck Collaboration et al. ( 2014), Lenz et al. (2017), and 3D maps such as Marshall et al. (2006) and Green et al. (2015).If Marshall et al. (2006) was the map with the largest impact on final mwdust extinction value, we pulled the uncertainties from dustmaps as they were included in the dustmaps module.For Green et al. (2019), values were not directly provided by dustmaps, so we followed the method in Green et al. (2015) and took the uncertainty to be half the difference between the 84th and 16th percentiles.While Drimmel et al. (2003) does not provide uncertainties, they point to López-Corredoira et al. (2002) who estimates the uncertainties in the K s band extinction to be less than 0.015, so all objects where Drimmel et al. (2003) largely determines extinction value have the error set to 0.015 for simplicity.Even though this reduces the confidence of calculated uncertainty for relevant objects, it is < 25% of the final sample.These uncertainties also required the matching conversion applied to the extinction above.
The values derived from 3D Galactic dust extinction maps are noticeably different from those calculated using intrinsic colors from T eff .We summarize the details in §4.7.

RESULTS
We determine our most probable RSGs by using the parameter space of log(L/L ⊙ ) vs T eff .To this end, we first determine the luminosity and T eff , then compare them to simulated evolutionary tracks for massive stars, and finally determine the most likely regions by considering the contamination from AGB stars.

Luminosity
Using the extinction method described in §2.4,we estimate the absolute magnitude in the K s band based on the de-reddened 2MASS K s magnitude and the distance moduli derived from the Bailer-Jones et al. ( 2018) geometric distance estimates.Bolometric corrections determined by Levesque et al. (2005) were matched to stars based on their spectral type and then applied to the absolute K s magnitude, yielding the bolometric magnitude.From the bolometric magnitude, we determined bolometric luminosities relative to the sun, log(L/L ⊙ ).By design, all uncertainties can be treated as independent, so we propagate errors of bolometric luminosities using those of T eff , 2MASS photometry, dust extinction, and geometric distance estimates.

Comparison to Evolutionary Tracks
With  eff and luminosity, we can compare the distribution of stars in the H-R diagram against evolutionary tracks predicted for massive stars.In Fig. 3, we show our sample against tracks of Ekström et al. (2012) for rotating stars computed at solar metallicity (z = 0.014) over the mass range 7-25 M ⊙ .For solar metallicity, anything above 25 M ⊙ is not cool enough to truly be called RSGs, and anything below 7 M ⊙ would be too small.While these evolutionary tracks do not account for binary evolution, the tracks are still appropriate as binary evolution does not greatly affect the location of RSG on the HR-diagram; Levesque (2017) succinctly showed this in their Fig.8.3 which shows a comparison between the single-star Geneva and BPASS binary evolutionary tracks (Eldridge & Stanway 2016;Eldridge et al. 2017) (Fig. 7 of Levesque (2018) also supports this).
The tracks for stars of 7 and 9 M ⊙ are particularly uncertain as their evolutions are not well understood, and this can be seen as these lines cover wider ranges in luminosity and are different shapes from those of higher mass.If these lines are used as cuts, it is possible that we would exclude true RSGs from our list.To take this uncertainty into account, we follow Messineo & Brown (2019) and use a diagonal cut base on the 7 M ⊙ track but shifted down and right to the ascending track.This line, plotted in thick purple in Fig. 3, gives us a limit for how faint or red a star could be and still be a candidate RSG: where log(T eff ) varies from 3.54 to 3.6 (or M4 to K1, based on Levesque et al. 2005).As RSGs candidates should be on or between the stellar tracks, with uncertainties taken into account, we can define approximate regions for RSGs.In the following, we determine this in more detail by also considering contamination.34.7% (0.0%) 5.4% E (Green) 44.8% (38.8%) ∼0%

Identifying regions
The main contamination for RSGs is AGB stars.AGBs have lower initial masses and spend more time on the main sequence; however, in their extreme cases, those that are oxygen-rich or carbon-rich increase in luminosity at a warmer temperature than the average AGB and will cross into the areas that we deem regions for RSG.To quantify this, we compare the distribution of Galactic AGBs and RSGs in the HR diagram.We sourced from the literature Galactic AGBs surveys that include values for both T eff and Luminosity (or M bol ) (McDonald et al. 2011;Blum et al. 2003;Groenewegen 2022;Milam et al. 2009) to quantify how much each region is contaminated.We include a list of extreme AGBs (Groenewegen 2022), including both C-rich and Orich, as these would most likely be incorrectly selected as RSGs.The Extreme AGBs list includes a small fraction of AGBs supplemented from the Magellanic Cloud.Altogether, we compile 377 AGBs to compare to the reference RSG sample compiled in §2.1.2.The results are shown in Fig. 4. Based on these, we define five regions.The percentages of AGBs and RSGs in these regions are given in Table 3.The five regions are: • Region A (Blue): defined by log( L L ⊙ ) ≥ 4.76.It has little contamination: only 1.7% of the reference AGBs are in region A, while 55.4% of the reference RSG population is here.Compared to stellar tracks, region A contains high mass RSGs with ≥ 15M ⊙ .
• Region B (Purple): defined by log( L L ⊙ ) between 3.92 and 4.76 and whose spectral type is earlier than M4, T eff ≥ 3535 K.It is still RSG-rich (39.2% of the reference RSGs) but has more contamination than region A (13.9% of reference AGBs).Region B contains intermediate mass RSGs with masses between 15 and 9 M ⊙ , which fall along well-defined stellar tracks.
• Region C (Red): defined by log( L L ⊙ ) between 3.36 and 3.92 and is bluer than Eq. ( 5).It contains intermediate-mass AGBs (4.9% of AGBs), and as discussed below, lower mass 9M ⊙ RSGs should be in this region but biased against observation (0% of RSGs).
• Region D (Orange): defined by log( L L ⊙ ) between 3.92 and 4.76 and anything T eff ≤ 3535 K. Compared to region B it contains a larger fraction of AGB contamination (34.7% of reference AGBs) and a smaller fraction of RSGs (5.4% of reference RSGs).
• Region E (Green): defined by log( L L ⊙ ) between 3.36 and 3.92 and redder than Eq. ( 5) or anything with log( L L ⊙ ) < 3.36.It contains a large fraction of the reference AGBs (44.8%) and none of the reference RSGs.Region E stars are either too red or too dim to be RSGs of any mass.
The reference RSGs along with our reference AGBs illustrate an issue in massive star studies, which is that we lack stellar tracks certain enough and photometric characteristics that could distinguish between 7-9M ⊙ RSGs and contamination from lower mass stars and extreme AGBs (Siess 2006;Poelarends et al. 2008).The best options beyond modeling and photometry also fail in certain cases; observation of spectral lines sensitive to surface gravity fail for cool emission spectra, high-resolution follow-up is not realistic logistically for large data sets, and use of maser emissions requires circumstellar masers to distinguish between supergiant and giant stars (Levesque 2017).Even in this list of well-studied RSGs, almost none have masses below 9 M ⊙ even though these should be more common than higher-mass RSGs.This results from the prevailing selection criteria, similar to that in Massey & Olsen (2003), which applies a conservative lower limit and filters out lower mass RSGs.Ultimately, we similarly remove Region C even though it should contain low-mass RSGs.

The Catalog
We show in Fig. 5 our five regions and our sample combining the compilation-based and Gaia-based methods.We keep stars in Regions A and B for our final RSG list.This results in 134 and 436 stars, respectively.However, we can also include stars that satisfy the region cuts after consideration of errors.We include distance, photometry, dust extinction, and effective temperature as errors.These stars are listed in their original region identifier, meaning C, D, or E (see Appendix A for a complete column listing of the catalog).Still, they are included as "likely" RSGs in the final catalog due to uncertainties.This adds another 62 stars, primarily from Region C.

Characteristics Flags
We add additional flags that give extra insight into how likely an object is an RSG, going beyond the required luminosity and temperature cuts.The statistics of these flags are shown in Table 4.

Clusters
Due to the dense molecular cloud that massive stars are born in, RSGs often form as part of a cluster or OB associations, so RSGs searches are often restricted to a specific cluster.We thus search the literature to see if each star is a member or candidate member of any Milky Way cluster.We include the cluster name, any radial velocity derived from cluster dynamics, and any cluster details, like age, approximate distance, and mass, into the catalog if available.The main clusters found include Stephenson 2 (RSGC2), PER OB1, and CAR OB1.

Multi Star Systems
Like OB associations and clusters, the dynamics of the molecular clouds are appropriate for forming binary systems.Considering the significant binary fraction for main-sequence massive stars, once a massive star reaches mass transfer at the Hayashi limit, it is likely that this will cause an increase in the binary fractions of massive stars.To what extent it increases depends on the theoretical treatment of RSGs.It is almost certain that the limitation of observations and the rarity of RSGs are the causes of so few known binaries.
We thus include a flag for stars that have evidence of interactions with a companion.The catalog contains a column to flag binary and multi-star systems and a column detailing if the spectra show a binary (Spectral Binary), if Gaia DR3 has identified it as a multiple object system (non_single_star), or if the literature defines it as such.There is also a corresponding reference column that will give the relevant reference.

Variability
RSGs are expected to be radial pulsators with luminosity-dependent pulsational periods (Stothers 1969(Stothers , 1972;;Heger et al. 1997;Guo & Li 2002) leading to photometric variation.Mass-loss rates have been measured up to 10 −5 M ⊙  −1 (Beasor & Davies 2016) and are another driver for variability.Spectroscopic variability is affected by the dynamics of their outer layers, like extended atmosphere and optical depth effects near the Hayashi limit or convection, pulsation, and mass loss.While the dominant driver responsible for this variation is unknown, it is an important characteristic of RSGs.We, therefore, included a flag for any object which is designated as a variable star, either in the literature, AAVSO, ASAS-SN (Shappee et al. 2014;Jayasinghe et al. 2020), or through Gaia DR2/DR3's photometric variability flag.The latter uses statistics, time-series, and data mining analysis on Gaia G, G BP , G RP photometry as well as parallaxes and positions (Holl et al. 2018) and using statistical and machine learning methods built from a global revision of major published variable star catalogs (Eyer et al. 2023)).

RSGs which are not in Gaia DR3
Even though Gaia has a large magnitude range G ∼ 3 to 21 mag, for any star with G mag < 5, the brightness could lead to saturation and subsequent systematic errors in parallax.This happens since Gaia astrometric measurements use the stars' magnitudes in part to self-calibrate, so this becomes a problem as bright stars are scarce (Lindegren, L. et al. 2018).Also, the extreme widths of massive stars can limit the accuracy of Gaia's parallax measurement.For stars with magnitudes bright and radii wide enough to have large enough errors to be removed from our pipeline but are generally accepted as RSGs in literature, we reintroduce them into our sample.There are 12 stars in this category.It includes  Del,  Cep, S Per, VX Sgr, and NML Cyg who are in Gaia DR3 but whose measurements are saturated and accuracy is reduced below our threshold, and alf Ori (aka Betelgeuse), alf Her, alf Sco, eps Peg,  Cep, and lam Vel who are not in Gaia DR3 due to some combination of their distance, brightness, and radii.We also include VY CMa, whose brightness and distance should not lead to over-saturation but is excluded as its radius, one of the widest RSGs in our galaxy, hinders accurate distance measurements.
As the above list is not measured with Gaia DR3, we do not use our pipeline to determine their bolometric magnitude or luminosity.Instead, we used values determined throughout the literature.These papers are denoted in the catalog, and a full list can be found in Appendix B).While we tried to fill as many columns as possible for each non-Gaia object, this was not possible for some columns, especially those from Gaia.A list is provided in a single reference column for objects whose values come from multiple sources.

MY Cep
The use of stellar tracks and reference RSGs/AGBs to determine our regions for possible candidates resulted in MY Cep being removed from our list, not because it does not have a high enough luminosity, but because its T eff is too cool.However, MY Cep is generally acknowledged as a RSG (Levesque et al. 2005;Mauron & Josselin 2011;Humphreys et al. 2020) and it makes sense to keep it when considering the uncertainty in T eff .Blaauw (1956a,b) discovered a sample of OB stars with significantly higher space velocities than their surroundings, denoted runaway stars (RWs).RWs were first theorized to originate from binary systems whose primary undergoes SN and shoots the secondary out, creating their peculiar velocities.However, this has been disproven as a primary mechanism (Gies & Bolton 1986).Other explanations hypothesize their existence from dynamical ejection (Leonard & Duncan 1990) 2011) which interestingly is both the only known K-type and whose proposed velocity evolution favors the nondominant mechanism of ejection due to SN (Comerón & Figueras 2020)).All of these are included in the final sample.

Presence of Magnetic field
While the extent to which magnetic fields can contribute to total mass loss in RSGs is unknown, perturbation of fields with strength greater than several Gauss could have effects through the production of Alfven waves (Hartmann & Avrett 1984;Charbonneau & MacGregor 1995;Vidotto & Jatenco-Pereira 2006).There have been several Galactic RSGs with measured magnetic fields, including ∼ 1 G in Betelgeuse (Aurière et al. 2010;Tessore et al. 2017).Observations of surrounding circular polarization from Zeeman-splitting suggest a surface magnetic field as strong as 4 G for VX Sgr (Vlemmings et al. 2005), ∼ 1 G for 32 Cyg, and ∼ 2 G for  Vel (Grunhut et al. 2010).Tessore et al. (2017) detected ∼ 1 G from CE Tau and  Cep, and several G in  1 Her similar to the magnetic fields of AGBs.
Measurements of VY Cma show disagreement with estimations from H 2 O maser observation ranging from 90 to 180 mG (Vlemmings et al. 2002) and estimations based on linear polarization suggesting a lower limit of ∼ 10 G (Shinnaga et al. 2017).Discussion detailing possible reasons for this disagreement are presented in Shinnaga et al. (2017).Non-detections of post-RSG stars have also been made, such as yellow SG (YSG)  Cas (Tessore et al. 2017), which suggests either a very weak magnetic field or its dissipation, giving insight into the stellar evolution of RSGs that evolve back to warmer temperatures.

Spatial Distribution
The current collection of Galactic RSGs is still highly incomplete, with little known about the total spatial distribution (see, e.g., Davies et al. 2009;Messineo et al. 2016).Even with our catalog being the largest collection, it only extends out 12.9 kpc when the radius of the Milky Way is ∼ 30 kpc. Figure 6 shows the spatial distribution of our catalog.In the right panel, we compare with the spatial distribution of other catalogs in the literature: Nakamura et al. (2016) uses less stringent luminosity cuts but still reaches approximately the same distance, and Mukhopadhyay et al. (2020) where the inner Galactic Center's RSGs and Blue SGs are represented but only goes out to < 1 kpc.Note that we see a more prominent dearth of stars around the solar system, driven by the saturation of Gaia photometry (see §4.1.1).
The driving factor for the lack of stars beyond 12.9 kpc comes down to the availability of spectra.Even in Skiff (2014), one of the largest collections of Galactic stellar spectra, only ∼ 0.3% meets the most basic spectral characteristics for RSGs.As the determination of spectral type from spectra for RSGs involves either labor, limitations to sample size, or reduction of confidence, RSGs searches have been limited by pointing in the direction of OB associations or clusters, which generally have a higher population of RSGs.However, for inner Galactic supergiants, only ≈ 2% are associated with stellar clusters (Messineo et al. 2017).The spectra of RSGs are not only observationally limited but also lack well-defined spectral standards.This is especially true of K-type stars, which are often broken into late and early K-type stars rather than subtypes.
Even after locating candidates, intrinsic characteristics that help determine whether it is a RSG (e.g., pulsation properties and chemical abundances) are not easy to obtain.As discussed in §3.3, reliable confirmation is also difficult because the colors of RSGs are not unique and match those of giant late-type stars, specifically from low masses to super-AGBs of 9-10 M ⊙ .These issues all lead to the number of Galactic late-type stars of class I being less than ∼1000, and, when not considering our new catalog, ∼400 RSGs known throughout various surveys.Major catalogs like Humphreys (1978) lists 92, Elias et al. (1985) lists 90, Levesque et al. (2005) analyzed the spectra of 62, Jura & Kleinmann (1990) lists ∼135, even though more than ∼5000 RSGs are predicted by Gehrz (1989).

Radii range
The expanse of RSG radii makes measurements challenging as it limits the observational effectiveness of interferometric determinations and requires sufficiently accurate distances.For the several dozen stars with well-measured radii, comparisons to estimates from L = 4R 2 T 4 eff show good agreement (Van Belle et al. 2009;Wittkowski et al. 2017Wittkowski et al. , 2012;;Arroyo-Torres et al. 2013, 2015).Solving the equation for the radius of Betelgeuse, ∼ 855 R ⊙ , compares well to interferometric measurements (Dolan et al. 2016) ∼ 890 R ⊙ , especially as there is variation depending on observed wavelength, treatment of asymmetries, among other things (Townes et al. 2009).Based on this simple but effective relationship, we provide radii for all stars in our catalog, resulting in a range of 10 1.38 to 10 3.28 R ⊙ .

Mass loss rates
For massive stars, more than half of their mass loss happens after the main sequence, with RSGs losing between 10 −7 to 10 −3 M ⊙ yr −1 (Mauron & Josselin 2011;van Loon et al. 2005).Observations of CSM around type II SNe suggest that mass loss during the final stages of the RSG phase is accelerated.Förster et al. (2018) surveyed 26 SN II within hours of their discovery and obtained optical light curves from High Cadence Transit Survey (HiTS) (Martínez-Palomera et al. 2018) that when compared to detailed models suggest density profiles consistent with > 10 −4 M ⊙ yr −1 .Confined dense CSM was confirmed to surround SN2013fs that was estimated to be ejected during the final ∼1 yr prior to explosion at a rate ∼ 10 −3 M ⊙ yr −1 (Yaron et al. 2017).However, with the dominant mechanism for significant mass loss unknown and explanations incomplete, derivations of general relations are restricted to observations and limited modeling.While there are some discrepancies between different mass loss rate relations, we estimate the mass loss of our RSG candidates based on van Loon et al. (2005), which gives log 10 ( )= −5.5 based on observations of a small sample of M-type stars and the assumption of a dust-driven wind model, most applicable to dustenshrouded RSGs and oxygen-rich AGB stars.We find a range of 10 −7.76 to 10 −3.87 M ⊙ yr −1 in our sample.

Dust Extinction Method Comparison
In Section 2.4, we discuss our reasoning for departing from previous methods, which use observationally determined intrinsic colors for different spectral types of massive stars and working backward to determine their dust extinction.However, when a comparison is made with the previous method, our 3D dust map, and other dust maps, we see differences significant enough for us to outline here.We compare six different methods of estimating dust: Schlegel et al.  (2023).We also explore how the calculations of bolometric luminosity would shift depending on the method used for extinction.Figure 7 shows comparisons of dust corrections, spanning 3D dust maps, mwdust, and the use of intrinsic colours of massive stars (ICMS).The close alignment to the 1-to-1 line, shown dashed in pink, for Marshall and Green results from mwdust predominately using those two maps for their healpix extinction map.We see that mwdust produces a smaller spread across the board and generally has smaller differences from the average values as well.In Fig. 8 we show how they impact our determined stellar luminosities.Except for the earliest 2D maps of Schlafly & Finkbeiner (2011) and Planck Collaboration et al. (2014), the luminosities determined by our dust method remain very similar to those from more modern dust maps.
We also look into how our final estimates for bolometric luminosity compare to those in the literature; their values are more confident as they generally focus on a small sample of RSG candidates and directly determine their characteristics.Using our reference RSGs from §2.1.2,we compare the changes that would result from using DR2 geometric distance (r_est Bailer-Jones et al. 2018) instead of DR3's.The mean, median, and mode values for the difference between the reference RSG values and those estimated using our method and DR3 are all ≥ 15% smaller than those from DR2.The general trend of DR3 based values being closer to their reference values reinforces our use of distances primarily from r_geo_med but supplemented with r_est.However, both methods systematically underestimate luminosity, with a few outliers generally caused by areas of high dust.Stars which have the required luminosity and T eff are very likely to be bright or brighter than our limits as systematics taken into account would shift them higher, validating the inclusion of stars whose estimated values are not high enough but whose uncertainty places them in one of the two retained regions.

COMPARISON WITH PREVIOUS WORKS
Collecting the Milky Way's RSGs has been a necessary challenge, particularly for two areas of astrophysics: stellar astronomy and studies of massive stars pre-explosion (Messineo & Brown 2019;Messineo 2023) and astroparticle focused more on neutrinos and the explosion itself (Nakamura et al. 2016;Mukhopadhyay et al. 2020).The motivations and approaches behind all efforts vary considerably.Messineo & Brown (2019) focused on finding candidate RSGs in Gaia DR2.They compiled Galactic spectral catalogs whose objects are in Gaia DR2 to estimate the effective temperature (T eff ) and determine reliable distance measurements, and match them with infrared and optical measurements from various sources.Their combination of multiple bands with effective temperature allowed them to estimate stellar bolometric magnitude in two ways: by using photometric measurements and by integrating under the SED.This resulted in 889 late-type stars and 43 highly-probable RSGs.
The release of Gaia DR3 and associated data products presented new opportunities to categorize Galactic RSGs.A new list of latetype bright stars, including RSG and bright AGBs was developed using Gaia DR3 Apsis (Delchambre et al. 2023;Bailer-Jones et al. 2013) and Gaia DR3 GSP-Phot and GSP-Spec parameters of known K-and M-type stars of Class I luminosity by Messineo (2023).In addition to the previous work, there are 203 new entries of late-type bright stars with 15 S-type, 1 S/C, 9 C-rich, and 20 confirmed new RSGs with 6 having bolometric magnitudes brighter than the AGB limit.2020) heavily focused on the potential of liquid scintillators neutrino detectors to localize presupernova neutrino signals, a task which may be achieved at close distances.Thus, their list is restricted to ≲ 1 kpc, for a total of 31 containing red and blue massive stars.No analysis is done beyond compilation to determine the accuracy of CCSN progenitor status.As the purpose and method for listing possible CCSN progenitors candidates were relatively consistent, there was an overlap of the most nearby candidates.
By comparison, our work resulted in 578 Milky Way candidate RSGs.To obtain this enhancement, we begin by following the method of Messineo & Brown (2019) but modify it in a few key ways to tailor our analysis for use with multi-messenger astronomy.For example, we do not restrict our sample to those in the Galactic plane, and we place an emphasis on having more RSGs at the expense of having slightly more AGB stars contaminating our list, similar to Nakamura et al. (2016) and Mukhopadhyay et al. (2020) but to a lesser degree.Furthermore, we work with the updated Gaia DR3, though we do not yet include Gaia DR3 spectra or Apsis parameters as Messineo (2023) does.Our method for determining a single spectral type for objects with multiple designations considers the increase in data and follows a more strict prescription than that of Messineo & Brown (2019).We also implement a new method that adds 33 unique highly probable and likely RSGs.This method uses a sample of known Kand M-type stars included in Gaia DR3 to determine a Gaia G band photometry cut-off.It is combined with external spectra, working similarly to Messineo (2023)'s use of Gaia Apsis parameters and spectra.
Even though we obtain a larger number of RSG candidates-with potentially more contaminants-it is worth pointing our we are far from complete.Given the Milky Way's core-collapse supernova rate of a few per century (Tammann et al. 1994;Rozwadowska et al. 2021) and the duration of the RSG phase for the most common 8 ⊙ star being ∼ 0.5 Myrs (Eldridge et al. 2008), the number of RSGs at any given time in the Milky Way is approximately ∼ 5 × 10 5 /100 ∼ 5000 (Gehrz 1989).

MULTI-MESSENGER ASTRONOMY
Multi-messenger astronomy is based on coordinating observations across multiple messengers-from electromagnetic radiation to gravitational waves and neutrinos-and on the interpretation of the joint results to optimize the science attained.With coordinated efforts, both detections and non-detections lead to useful information on the physics behind different astrophysical phenomena.A key to successful coordination for transient events such as CCSNe is in aligning the observatories in time.To this end, advanced warning of an impending CCSN is highly useful.
Neutrinos will be prolifically produced during stellar CC, and since they arrive minutes to days before the rise of the electromagnetic emission, provide the desired early alert (see, e.g., Kharusi et al. 2021).Also, neutrino production rapidly rises during the last stages of nuclear burning, which can provide an advanced warning of the CC itself.One early established neutrino alert system is the SN Early Warning System (Antonioli et al. 2004).SNEWS provides an early warning for a Galactic CCSN by tracking coincidences between neutrino experiments.In addition to providing coincident alerts, the successor to SNEWS, SNEWS 2.0, will aim at providing as much  useful information as possible to capture the multi-messenger counterparts of the CC (Kharusi et al. 2021).In particular, the distance and position of the source will be key parameters that can be estimated with the neutrino data.The snewpdag5 software will be in charge of these calculations and used to release automated information in case of a CCSN alert.
SNEWS also coordinates with citizen astronomers, for instance, with the American Association of Variable Star Observers (AAVSO) to monitor potential CCSNe candidates long before they collapse, allowing more detailed observations on the development of stars as they evolve towards CC.The alert follow-up is being facilitated with the Recommender Engine For Intelligent Transient Tracking (Sravan et al. 2020).This is particularly useful as peculiar behavior prior to CC has been documented in several SNe with spurious observations.For example, the progenitor of SN 2009ip underwent a series of mini explosions prior to SN (Smith et al. 2022).It could also allow us to check for evidence of pre-SN mass loss, more accurately relate explosion energy with progenitors' initial mass, and define the relationship of light curve plateau brightness, temperature evolution of SN emission, and progenitor radius.

Directional information from neutrinos
There are two main methods to obtain angular pointing information with the neutrino data: a direct direction reconstruction using anisotropic interactions and triangulation.The first exploits an interaction between the neutrino and a target with some intrinsic directionality with a detector that can exploit such information.For example, the electron scatters preferentially along the neutrino direction in neutrino-electron elastic scattering.Cherenkov detectors can track the direction through the light that is released when the electron is kicked.With this method, Super-Kamiokande can give the direction within 3-8 degrees for a SN at 10 kpc (Mukhopadhyay et al. 2020).In relation to our catalog, neutrino pointing can reduce the possible potential stellar candidates on average within the error circles of 3-8 degrees down to 12-41 stars, respectively.Other neutrino interactions either do not retain very well directional information or have smaller cross sections.
The second method is triangulation (or, more precisely, multilat-eration), which uses the relative arrival times of the neutrino burst through multiple detectors in different geographical locations to determine from where those neutrinos originate (Vogel 1999;Beacom & Vogel 1999;Mühlbeier et al. 2013).The time of flight through the Earth for a neutrino is on the order of milliseconds.The success of the triangulation approach, therefore, depends on detectors registering samples of neutrino events large enough to achieve sufficient burst timing precision, or to have a very low background so as to precisely tell which is the first signal event.It also benefits from wide geographical spread among the contributing detectors.Estimates for the precision of triangulation under various assumptions and as a function of distance and location have been studied in various forms (e.g.2020)) where the most favorable conditions give a 1 area of a few percent of the sky.Once a trigger is received, SNEWS's network runs the information through a pointing algorithm, producing a trigger packet of a skymap of direction probabilities and distance estimations (See §6.2).Direction information from anisotropic interactions, calculated by individual experiments, can be incorporated into a final probability skymap, which can then be used to prioritize follow-up observation, especially if one is able to focus on a finite list of high-probability candidates, maximizing the prospects of obtaining both pre-SN observations and capturing the early rise of the SN with better sensitivity and cadence.
In Fig. 9, we show simulations of the triangulation error box estimated with the snewpdag package and overlaid with our RSG list.We considered four stars ([W60] C10, X-46, CD-34 11794, and HD 303250) in the catalog covering different positions on the sky, distances, and stellar clustering.We fix the CC time epoch to illustrate both sky location dependence and the distance dependence.Simulations using SNEWPY (Baxter et al. (2022)) of the CCSN models in Burrows & Vartanyan (2021), with the corresponding stellar masses, have been used for this study.We apply the method described in Coleiro et al. (2020) to inverse beta decay (   →  + ) samples at four detectors (Super-Kamiokande, JUNO, IceCube, and SNO+) to estimate the time delays.This method to estimate the time delay uncertainties by combining the observed neutrino lightcurves was chosen as it does not rely on any model template.The time uncertainties obtained range from 1ms, for distances up to ∼1kpc, to ∼10 ms for 12-13 kpc.We see that in all cases, the angular resolution significantly reduces the number of possible candidates, and further reduction can be achieved with the use of a distance estimate.

Distance estimate with neutrinos
It was shown in Segerlund et al. (2021) that the   event rate in the early stages of the neutrino burst is related to the CCSN distance.Two methods proposed have been implemented into snewpdag: the first compares the observed number of events ( 50 ) to the expected value weighted over the initial mass function (the "IMF method"), and the second (the "f Δ method") relies on a distance independent parameter, f Δ = N(50)/N(100-150), and its linear relation with  50 (Horiuchi et al. 2017).Here, N(50)/N(100-150) is the ratio of the expected number of events in the two respective time windows.We use both approaches to evaluate the distance uncertainty for our selected four candidate stars of the catalog, using a detector of the size of Super-Kamiokande or JUNO.We consider the same CCSN models as in §6.1 and include the systematic uncertainty due to the model dependency of the methods.A further systematic uncertainty arises from the time uncertainty of §6.1, which affects how the time windows (and therefore their event counts) are reckoned.This systematic effect has been studied by shifting the neutrino light curve in simulations with respect to the true bounce time according to the expected time uncertainty.The results, shown in Fig. 10, show the complementarity between the two methods: the f Δ method, dominated by statistical uncertainty, works better at smaller distances, while the IMF method outperforms at larger distances, albeit dominated by systematic uncertainties.Using this information on top of the triangulation skymap, the number of RSGs is further reduced in the best of cases to a few stars.For areas of high stellar density, the reduction is to less than 33% of the candidates falling within the 90% confidence level sky region, as seen in Fig. 9.

Implications for MMA
Obtaining detailed observations of a CCSN, from pre-explosion and neutrino burst through to shock breakout and the SN, would provide a treasure trove of data.Our RSG list would help realize this in the next Galactic CCSN by providing likely targets in the event of a neutrino alert in the future.
To date, pre-explosion observations of CCSN have been few and far between, but they have revealed interesting, unexpected behaviors.For example, the progenitor of SN2009ip underwent a series of mini explosions prior to SN (Maza et al. 2009;Miller et al. 2009;Smith et al. 2010) and SN1987a had a blue supergiant progenitor (Sonneborn et al. 1987) contrary to the standard expectation for a RSG progenitor.Pre-explosion images have also revealed the so-called "red supergiant problem" wherein Type IIP SN progenitors appear to have a maximum mass of ∼ 17 ⊙ , which is smaller than the most massive RSGs seen locally (Smartt et al. 2009).Various possible solutions have been suggested.For example, the most massive RSGs may not be exploding as Type IIP SNe, necessitating a revision of our understanding of how RSGs end their lives, with deep implications for the formation rate of black holes and the influence of SN feedback in galaxy evolution (e.g.Horiuchi et al. 2011;Crain et al. 2015;Kochanek 2015).Alternatively, the discrepancy could derive from our lack of understanding of stellar physics such as mass-loss (Ekström et al. 2012;Georgy 2012), or larger than assumed systematic effects, for example, in bolometric corrections (Davies et al. 2013) or circumstellar extinction (Walmswell & Eldridge 2012;Beasor & Davies 2016), causing an underestimation of the progenitors mass.Having detailed pre-explosion monitoring data, in both cadence and multiple bands, would help in all the above debates.
As discussed above, the detection of GWs, multi-wavelength photons, and neutrinos provides a unique opportunity to probe the mechanisms of CC; however, the current setup of telescopes and detectors have limitations.Most Galactic SNe should be observable in the optical and virtually all should be observable in the near-IR (Adams et al. 2013).There is even a one-in-three chance for the SN to be visible to the naked eye.Neutrino detections easily cover the entire Milky Way galaxy.The detection by Super-Kamiokande with Gadolinium of a CCSN within 8.5 kpc from the Sun will have a pointing accuracy of ∼ 3 degrees or better.Note that the region spanning a radius of 8.5 kpc includes 99% of our full catalog.The detectability of GWs depends strongly on the uncertain core rotation.Nakamura et al. 2016 has shown that for an initially non-rotating core, GW detection is possible for sources up to the Galactic center if there is a coordinated observation with GW detectors and neutrino telescopes.
Statistically, the most likely distance to the next CCSN is close to the Galactic Center (Adams et al. 2013).For example, the modeling of the progenitor and Galactic dust as a double-exponential spatial distribution indicates that within 5.8 kpc of the Galactic center corresponds to 62 percent of the Galactic CCSNe rate (Adams et al. 2013;Nakamura et al. 2016).While estimates from Galactic distributions of SN remnants (SNR) (Ahlers et al. 2009;Mirizzi et al. 2006) show some variation, a radius of ∼6 kpc consistently contains ≥ 50%.In our sample, this results in 41 candidates.On the other hand, CCSN closer to the Sun will have better multi-messenger prospects, as we summarized in the previous paragraph.Here, we have 594 stars within 8.5 kpc from the sun (excluding the stars added based on uncertainty), i.e., within the distance regime where neutrinos can provide ∼ 3 degrees of pointing accuracy and detection prospects for GWs and infrared photons are the most optimistic.Thus, the most promising objects are those within both regions; all 41 stars of the former set are also in the latter.Assuming that only M-type stars are approaching explosion following from the stellar evolutionary paths used in Section 3.2, 36 highly promising stars are retained within the overlapping spatial regions.
Further studies are guaranteed.Within the sample of 36 promising candidates, greater than > 70% have observations in more than 3 electromagnetic bands.Nevertheless, effort must be made to expand the available high-quality images to more wavelengths and maintain up-to-date follow-up spectra.

SUMMARY AND FUTURE WORK
Using two complementary methods, one involving a thorough literature survey and the other using Gaia data, we compiled a catalog of 578 Milky Way RSGs.This is the largest catalog of its kind in the literature.To differentiate between RSGs and non-RSG contaminants, the most problematically extreme AGBs, we determined the bolometric luminosity and T eff .We compared them to stellar evolutionary tracks and Galactic AGB surveys.In this context, dust extinction is a major systematic uncertainty, which we estimated using mwdust, a 3D dust map, to estimate dust extinction for each star to reduce and prevent correlations in uncertainties.The accuracy of this method was tested by comparing both the estimated extinctions and final luminosities to several alternative methods.Consideration of uncertainties adds another 62 stars that enter the RSG parameter space when individual errors are considered, resulting in a total of 640 candidates.
Along with the information collected to determine their luminosities, the catalog contains details of each star's known characteristics, including evidence of binary interactions, magnetic fields, variabil-ity, and cluster membership.Using the RSG catalog compiled, we explored the RSG's spatial distribution, stellar radii range, and massloss rates.Finally, we explored the impacts of the RSG catalog for multi-messenger observations of the next Galactic CCSN.The intense neutrino burst from the CC of a massive star provides a natural alert, but our RSG catalog can help suggest locations for observations of a CCSN.Typically, a dozen or so of the RSGs reside within the typical pointing uncertainty given by the neutrino alert.
Future work includes both improvement and expansion of the current catalog.As extinction is likely our greatest source of error, it seems clear that a better understanding of the complexity of Galactic ISM, including the completion of Galactic 3D dust maps and the excess of CSM around RSGs, such as a more robustly defined R V and associated extinction ratios, would not only improve luminosity estimation but increase the likelihood of including dusty RSGs.Achieving the collection of more candidates in Region C is also important for future work.Issues of AGB contamination could be further mitigated with improvements in infrared spectra and flux values, like those that could be provided by the Roman Space Telescope scheduled to launch in 2027.The subsequent infrared CMDs can be combined with knowledge of an object's variability to reduce the ambiguity between AGBs and RSGs.Inspecting the new infrared CMDs can also provide an avenue for finding previously unknown candidates for further study.Expansion of available wavelength and SED fits will allow us to build a more robust profile for each of our candidate RSGs.
While we wait to detect neutrinos, the plan is to improve the characterization of candidates in our list so that post-explosion analysis will require as few assumptions and be as complete as possible.Along with our coordinated efforts to photometric variability monitoring via AAVSO, which is done primarily in the optical band, it would be important to expand the number of pre-explosion images across a greater range of wavelengths for each candidate.Mid-infrared photometry, in particular, will be useful to determine the dusty CSM around the RSGs, for which detailed knowledge of the environment will improve future analyses.This can be accomplished with the Roman infrared space telescope and its 2.4-meter mirror and Wide Field Instrument.Also important would be starting a campaign of follow-up spectroscopic observations to determine the variability, confirm our estimated luminosities, and collect spectra of the stellar progenitors as close to the explosion as possible.These efforts will also aid in determining fundamental properties like surface gravity and chemical composition to enhance our ability to trace stars from the zero-age main sequence to core collapse.Some next steps are more immediate; for example, while a RSG is statistically the most likely next CC progenitor, other types of massive stars must also cause CCSNe, such as stripped-envelope SNe.This motivates the inclusion of a broader set of stars to complete further the target list (e.g., Rosslowe & Crowther 2015).Continuation or initiation of monitoring more broadly CC candidates will provide more insight into each.
The next Galactic CCSN will be a once-in-a-generation opportunity to collect exquisite multi-messenger observations.With a more complete pre-assembled target list, the completeness of such observations will be improved and would likely shed new insights and potential surprises into stellar and CCSN physics.

Figure 1 .
Figure 1.Log of the density distribution of Gaia DR2 stars in the colourmagnitude diagram of absolute Gaia G passband magnitudes (near-ultraviolet to near-infrared) and BP-RP colour.Overlaid are 52 stars in our reference RSG sample (see selection criteria in §2.1.2),plotted as red stars.Our selection cut is indicated by the solid thick dark red line (see Eq. (1)).Stars above the cut in the green-blue shade are retained as they have temperatures and luminosity similar to those in our reference RSG sample.

2Figure 2 .
Figure 2. r_med_geo distance derived by Bailer-Jones et al. (2021) using Milky Way model, against parallactic distances from direct inversion of the parallaxes.Those with relative error > -1 and RUWE < 2.7 are in green, and relative error > 2 and RUWE < 2.7 are in dark blue.The cyan objects have relative error > 4 and RUWE < 2.7.The pink dashed line shows the line of zero difference.

Figure 3 .
Figure 3. Luminosity vs. T eff of stars in both methods, colored by their intrinsic K s magnitude.Stellar tracks from models at solar metallicity and including rotation from Ekström et al. (2012) are shown for different ZAMS masses, as labeled.The thick purple diagonal line recreates the ascending track from the 7 M ⊙ progenitor, shifted to account for uncertainties in ≲ 7 − 9 M ⊙ mass stellar tracks.

Figure 4 .
Figure 4. Left: 377 Galactic AGB, including a set of extreme AGBs from the literature (McDonaldet al. 2011;Blum et al. 2003;Groenewegen 2022;Milam et al. 2009); see §3.3 for further details.In each of the regions A-E, we label the percentage of the AGB sample residing in that region.Right: our 156 reference RSGs(Figer et al. 2006;Negueruela et al. 2012;Clark et al. 2009;Marco & Negueruela 2013;Levesque et al. 2005;Humphreys 1978), compiled under the guidelines set in §2.1.2.The percentages indicate the fraction of the reference RSG residing in each of the regions A-E.

Figure 5 .
Figure 5. Luminosities vs. T eff values of the compilation-based method (arrows pointing right) and the Gaia-based method (arrows pointing left), for stars with relative error > 4 and RUWE < 2.7.Stars within Region A (blue) and Region B (purple) are highly probable RSGs.Region B contains some contamination from AGBs.Region C (Red) contains intermediate-mass AGBs and late M and K type stars with solar masses around 9M ⊙ .As discussed in §3.3, lower mass RSGs should be in this region.Region D (Orange) and Region E (Green) stars are either too red or too dim to generally be RSGs of any mass, so they are unlikely regions to find RSGs.The number of stars for each region, not accounting for repeats between methods, is shown in the legend in parentheses.

Figure 6 .
Figure 6.Left panel: Our final RSGs (teal stars) in Galactocentric coordinates compared to a population of inner Milk Way young stars (shaded based on stellar density from a 2D histogram using cells of 150x150 pc for y vs. x and 150x115 pc for z vs. x) on the upper main sequence selected from Gaia DR3 based on methods of Gaia Collaboration et al. (2023b) whose density distribution reconstructs the spiral arms of the Milky Way, specifically, Perseus arm, Orion spur, and Sagittarius arm.The Sun's location (8, 0, 0) is marked in white.Right panel: Our RSGs catalog (large black stars) compared to previous studies: Nakamura et al. (2016) (small light blue stars), Mukhopadhyay et al. (2020) (small green stars), and reference RSGs (pink).The Sun location (8, 0, 0) is marked in yellow.The number of stars is shown in the legend in parentheses.

Figure 7 .
Figure 7. Comparisons of 3D dust maps, mwdust, and using ICMS to other methods of determining Galactic extinction.The top row shows ICMS along the y-axis and the bottom row shows mwdust extinction converted to the K s band along the y-axis, with each being compared to 2D maps (Schlafly & Finkbeiner 2011; Planck Collaboration et al. 2014; Delchambre et al. 2023) and 3D maps (Sale et al. 2014; Green et al. 2019; Marshall et al. 2006) along the x-axes.Points are colored by the size of the difference in the y-axis estimated extinction and the average extinction value across all methods.

Figure 8 .
Figure 8. Same as Fig. 7 but showing how different dust correction methods impact our determined luminosity.
(a) Calculations based on CD-34 11794 (263 • , -34 • , and 4.1 kpc) as the progenitor.15 stars lie within the 90% credible region with 2 of those matching the combined distance estimation and angular resolution.(b) Calculations based on HD 303250 (161 • , -58 • , and 2.4 kpc) which lies within cluster CAR OB1 as the progenitor.85 stars lie within the 90% credible region, with 20 of those matching the combined distance estimation and angular resolution.(c) Calculations based on X-46 (113 • , -22 • and 7.6 kpc) as the progenitor.105 stars lie within the 90% credible region, with 2 matching the combined distance estimation and angular resolution.(d) Calculations based on [W60] C10 (83• , -67 • , 12.9 kpc) as the progenitor.187 stars lie within the 90% credible region, with 3 of those matching the combined distance estimation and angular resolution.

Figure 9 .
Figure 9. Confidence level localisation skymaps derived from triangulation, using the simulation of four CCSN events corresponding to the four stars shown in Fig. 10.Overlaid is our RSG list (light gray stars), with those inside the 90% angular error as dark purple open circles.RSGs, which also match the distance estimated from neutrinos, are shown by dark purple stars.The four panels show different sources at various RA, Dec, and distances.

Figure 10 .
Figure10.On the right, we see the absolute distance uncertainties in kpc for our four selected stars, and for each of the two distance estimation methods.On the left, the relative distance uncertainty as a function of the distance.Star markers show the method (statistics plus systematics) error for each of the stars.The color bands correspond to the uncertainties when adding the systematics related to the timing uncertainty.

Table 1 .
Progression of sample size moving through steps building up to calculations of bolometric luminosity.

Table 3 .
Percentages of the reference RSGs and AGBs in the five regions labeled A∼E.For each region, the percent of AGBs known to be extreme, either carbon or oxygen, is also stated in parentheses.Note that percentages in region C are likely severely underestimated as there is a bias against confirming RSGs or AGBs (see §3.3).

Table 4 .
Breakdown of characteristics for highly probable RSGs included in the final catalog separated into Region A (138 stars) and Region B (460 stars).Close RSGs are those whose proximity, extended radius, and/or brightness lead to poor or no measurements in Gaia and were manually added to our catalog as described in §4.1.1.