Astrophotometric search for massive stars in the Milky way. Confronting Random Forest predictions with available spectroscopy

TL;DR

This study develops a Gaia-based machine learning method to identify bright massive star candidates in the Milky Way, achieving high completeness and purity, and validates it through spectroscopic analysis, potentially increasing known bright massive stars by 50%.

Contribution

Introduces a Gaia DR3 and literature data-driven Balanced Random Forest approach for identifying bright massive star candidates with validated spectroscopic confirmation.

Findings

Achieves ~80% completeness in candidate selection.

Purity ranges from 51% to 85% depending on probability threshold.

Spectroscopically confirms 300 stars as B3 or earlier, including 107 new stars.

Abstract

Massive stars play a significant role in different branches of astronomy, from shaping the processes of star and planet formation to influencing the evolution and chemical enrichment of the distant universe. Despite their high astrophysical significance, these objects are rare and difficult to detect. With Gaia's advent, we now possess extensive kinematic and photometric data for a significant portion of the Galaxy that can unveil, among others, new populations of massive star candidates. In order to produce bonafide bright (G magnitude $<$ 12) massive star candidate lists (threshold set to spectral type B2 or earlier but slight changes in this threshold also explored) in the Milky Way subject to be followed up by future massive spectroscopic surveys, we have developed a Gaia DR3 plus literature data based methodology. We trained a Balanced Random Forest (BRF) with the spectral types from the compilation by Skiff et al. (2014) as labels. Our approach yields a completeness of $\sim80\%$ and a purity ranging from $0.51 \pm 0.02$ for probabilities between 0.6 and 0.7, up to $0.85 \pm 0.05$ for the 0.9-1.0 range. To externally validate our methodology, we searched for and analyzed archival spectra of moderate to high probability (p $>$ 0.6) candidates that are not contained in our catalog of labels. Our independent spectral validation confirms the expected performance of the BRF, spectroscopically classifying 300 stars as B3 or earlier (due to observational constraints imposed in the B0-3 range), including 107 new stars. Based on the most conservative yields of our methodology, our candidate list could increase the number of bright massive stars by $\sim$50\%. As a byproduct, we developed an automatic methodology for spectral typing optimized for LAMOST spectra, based on line detection and characterization guiding a decision path.

Tables

Table 1. Table 1: Number of classifications excluded

Spectral type (labels)	Number
Binaries	2456
em	6220
OB	16363
OBe	721
Non MK and uncertainty	16345

Table 2. Table 2: hyperparameters tested for the Random Forest and the best value

Features	Tested parameters	Best parameter
n_estimators	[500, 800, 1200]	1200
max_features	[sqrt, log2]	log2
max_depth	[10, 20, 30, None]	30
sampling_strategy	[majority, not majority, all]	majority

Table 3. Table 3: Performance metrics for the test sets

	Precision	Recall	F1 Score
Later than B2	0.968 $\pm$ 0.002	0.878 $\pm$ 0.005	0.921 $\pm$ 0.003
B2 or earlier	0.48 $\pm$ 0.01	0.80 $\pm$ 0.01	0.599 $\pm$ 0.009

Table 4. Table 4: Categories of possible peculiarities:

Category	Peculiarity
Uncertainty	:, ?
Miscellaneous	comp, k, p, pec, !, var, v, w, wl, wk
Line-broadening	(n), n, nn, nnn, [n], s, sh, shell
Emission lines	e, em, (e), [e], eq, neb, ((f)), (f), f, ((f*)), (f*), f*, c, f?p, (n)fp, nfp, z, q
Carbon/Nitrogen enhanced	N, C

Table 5. Table 5: Minimum and maximum values for different features.

Feature	Min value	Max value
$J$	5.531	14.69
$H$	5.502	14.69
$K_s$	5.462	14.70
$G$	5.665	14.97
$G_{BP}$	5.641	15.19
$G_{RP}$	5.611	14.73
$J$-$H$	-0.125	0.433
$J$-$K_s$	-0.154	0.661
$J$-$G$	-2.426	0.3933
$J$-$G_{BP}$	-3.246	0.4882
$J$-$G_{RP}$	-1.434	0.2395
$H$-$K_s$	-0.077	0.261
$H$-$G$	-2.853	0.4913
$H$-$G_{BP}$	-3.668	0.5923
$H$-$G_{RP}$	-1.856	0.3399
$K_s$-$G$	-3.086	0.5268
$K_s$-$G_{BP}$	-3.892	0.6251
$K_s$-$G_{RP}$	-2.087	0.3714
$G$-$G_{BP}$	-0.8725	0.2467
$G$-$G_{RP}$	-0.1592	1.046
$G_{BP}$-$G_{RP}$	-0.2591	1.829
${\sigma }_{\overline{\omega }}$	0.01088	0.4444
$\overline{\omega }/{\sigma }_{\overline{\omega }}$	-2.518	161.8
$\overline{\omega }$	-0.1438	5.356
${\mu }_{\alpha }$	-20.03	12.75
${\mu }_{\delta }$	-17.82	11.18
${\sigma }_{{\mu }_{\delta }}$	0.01095	0.4344
${\sigma }_{{\mu }_{\alpha }}$	0.00973	0.473
$ASM$	0.01693	0.7232
$AEN$	0.000	3.907
$RUWE$	0.6365	20.98
$AGOFAL$	-9.306	265.8
$IPDGOFHA$	0.00265	0.1691
${\sigma }_{AEN}$	0.000	1.564e+04

Equations

$$\text{Recall}=\frac{TP}{TP+FN}$$

$$\text{Precision}=\frac{TP}{TP+FP}$$

$$\text{F1 Score}=2\times \frac{\text{Precision}\times \text{Recall}}{\text{Precision}+\text{Recall}}$$

$$TNR=\frac{TN}{TN+FP}$$

$$Dom=recall-TNR$$

$$G\text{-}mean=\sqrt{\text{Recall}\times TNR}$$

$$IB{A}_{\alpha }(M)=(1+\alpha \cdot Dom)\cdot M$$

Full text

11institutetext: 1Departamento de Astronomía, Universidad de La Serena, Avenida Juan Cisternas 1200, La Serena, Chile

2European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching, Germany

3Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, the Netherlands

4Instituto de Física y Astronomía, Universidad de Valparaíso, Avenida Gran Bretaña 1111, Valparaíso, Chile

5Gemini Observatory, NSF’s NOIRLab, Av. J. Cisternas 1500 N, 1720236 La Serena, Chile

11email: nicolas.monsalves@userena.cl

Astrophotometric search for massive stars in the Milky way. Confronting Random Forest predictions with available

spectroscopy

N. Monsalves 112[\scalerel*

—](https://orcid.org/0000-0002-4129-8195)2[\scalerel*

—](https://orcid.org/0000-0002-4129-8195)   

A. Bayo 2[\scalerel*

—](https://orcid.org/0000-0001-7868-7031)2[\scalerel*

—](https://orcid.org/0000-0001-7868-7031)   

M. Jaque Arancibia 1[\scalerel*

—](https://orcid.org/0000-0002-8086-5746)1[\scalerel*

—](https://orcid.org/0000-0002-8086-5746)   

J. Bodensteiner 332[\scalerel*

—](https://orcid.org/0000-0002-9552-7010)2[\scalerel*

—](https://orcid.org/0000-0002-9552-7010)   

A. G. Caneppa 4[\scalerel*

—](https://orcid.org/0009-0003-3496-2611)4[\scalerel*

—](https://orcid.org/0009-0003-3496-2611)   

P. Sánchez-Sáez 2[\scalerel*

—](https://orcid.org/0000-0003-0820-4692)2[\scalerel*

—](https://orcid.org/0000-0003-0820-4692)   

R. Angeloni 5[\scalerel*

—](https://orcid.org/0000-0001-7978-7077)5[\scalerel*

—](https://orcid.org/0000-0001-7978-7077)

Massive stars play a significant role in different branches of astronomy, from shaping the processes of star and planet formation to influencing the evolution and chemical enrichment of the distant universe. Despite their high astrophysical significance, these objects are rare and difficult to detect. With Gaia’s advent, we now possess extensive kinematic and photometric data for a significant portion of the Galaxy that can unveil, among others, new populations of massive star candidates. In order to produce bonafide bright (G magnitude ¡ 12) massive star candidate lists (threshold set to spectral type B2 or earlier but slight changes in this threshold also explored) in the Milky Way subject to be followed up by future massive spectroscopic surveys, we have developed a Gaia DR3 plus literature data based methodology. We trained a Balanced Random Forest (BRF) with the spectral types from the compilation by Skiff (2014) as labels. Our approach yields a completeness of $\sim 80\%$ and a purity ranging from $0.51\pm 0.02$ for probabilities between 0.6 and 0.7, up to $0.85\pm 0.05$ for the 0.9–1.0 range. To externally validate our methodology, we searched for and analyzed archival spectra of moderate to high probability (p ¿ 0.6) candidates that are not contained in our catalog of labels. Our independent spectral validation confirms the expected performance of the BRF, spectroscopically classifying 300 stars as B3 or earlier (due to observational constraints imposed in the B0-3 range), including 107 new stars. Based on the most conservative yields of our methodology, our candidate list could increase the number of bright massive stars by $\sim$50%. As a byproduct, we developed an automatic methodology for spectral typing optimized for LAMOST spectra, based on line detection and characterization guiding a decision path.

\nolinenumbers

1 Introduction

It is commonly agreed that the galactic stellar mass distribution peaks at $\sim 0.5M_{\odot}$, and that its shape has a certain universality (Bastian et al., 2010). However, the slopes towards the higher and lower mass ends are still under debate. Despite their lower numbers in our Galaxy, the effect of massive stars ($M>8M_{\odot}$) on their surroundings is critical. These stars play a fundamental role in astrophysics, from shaping the processes of star and planet formation (Bally & Zinnecker, 2005) to providing us with key knowledge of the distant universe (Eldridge & Stanway, 2022). In terms of direct outcomes, massive stars can produce a wide range of phenomena, such as Type II supernovae (Smartt, 2009), neutron stars and black Holes (Heger et al., 2003), gamma-ray bursts (Izzard et al., 2004; Nakar, 2007), and due to their high multiplicity fraction (Sana et al., 2012; Moe & Di Stefano, 2017; Barbá et al., 2017), they are sources of gravitational waves (Schneider et al., 2001).

Observational classification of stars regarding their mass cannot be directly achieved (as the mass is not an observable). However, classification systems like the MK (Morgan & Keenan, 1973) aim to map observables (spectral types) to masses. This system is anchored with a set of standard stars, and classification relies on comparing spectral morphology with these standards. The system focuses on the blue part of the optical spectrum, which historically were the regions of highest sensitivity in photograph plates (Pickering, 1890; Walborn, 1972; Gray & Corbally, 2009; Sota et al., 2011, 2014). Regarding massive stars, the 8 $M_{\odot}$ boundary corresponds to spectral type B2 for dwarfs (V), B5 for giants (III), and B9 for supergiants (Ia/Ib, Morgan 1951; Gray & Corbally 2009; Ma´ız Apellániz et al. 2024).

Focusing mostly on the Milky Way, extensive efforts have been made to identify and spectroscopically classify massive stars. These include, among others, targeted studies (Sota et al., 2011; Simón-D´ıaz et al., 2011a, b; Ma´ız Apellániz et al., 2013; Sota et al., 2014; Simón-D´ıaz et al., 2015; Li, 2021), specialized compilations such as Pantaleoni González et al. (2021), and broader spectral type compilations (Skiff, 2014). Combined, these works yield a total of $\sim$7500 stars spectroscopically classified as B2 or earlier.

This low figure has been subject of a heated debate for several decades. Kroupa & Weidner (2003) reported an apparent lack of massive stars in the field when compared to their relative number with respect to late-type stars in clusters. This finding is supported by observational studies such as Rybizki & Just (2015), however, due to the low number and non-homogeneous distribution of high-mass stars, the slope of the mass function above 10 $M_{\odot}$ is an extrapolation. One possible explanation for the apparent deficit of field massive stars is provided by the Integrated Galaxy initial mass function theory (Kroupa et al., 2024), but see Cerviño et al. (2013b), and Cerviño et al. (2013a) for statistical counterarguments. It is therefore necessary to obtain more complete and robust censuses of high-mass stars.

Astronomy is undergoing a transformative phase with an unprecedented influx of data from large-scale surveys. The Gaia mission with its Data Release 3 (Gaia DR3, Gaia Collaboration et al. 2023) has delivered six-parameter astrometric data for nearly 900 million stars. This information will soon be complemented by upcoming spectroscopic surveys, for instance, the 4-metre Multi-Object Spectroscopic Telescope (4MOST) (de Jong et al., 2016; de Jong et al., 2019), WHT Enhanced Area Velocity Explorer (WEAVE) (Jin et al., 2024), Sloan Digital Sky Survey-V (SDSS-V) (Kollmeier et al., 2019). Additionally, VLT/MOONS (Gonzalez et al., 2022) will be equipped with 1000 fibers, significantly more than any existing near-infrared spectrograph. Therefore, to fully and efficiently exploit these extensive capabilities, it is essential to inform the surveys with the most robust candidate lists.

Within this context, Zari et al. (2021) presented a comprehensive catalog of hot, luminous stellar candidates with estimated spectral types from O to A in the Milky Way. An immediate goal of the work was to provide robust candidates for SDSS-V follow-up (Kounkel et al., 2023). The candidate catalog consisted of $\sim$990000 sources with $G<16$ mag selected with Gaia Early Data Release 3 (Gaia eDR3; Gaia Collaboration et al. 2021) photometry and astrometry and Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) photometry. The sample was designed to be more complete than pure with a reported purity of 20% for stars with temperatures consistent with spectral type B7V.

In this work, we present a supervised machine learning approach designed to improve the yield of massive star candidates with spectral types B2 or earlier in the Milky Way. The method builds on and updated version (fully Gaia DR3) the selection criteria proposed by Zari et al. (2021). In order to validate our approach, we mined spectroscopic archival data for our candidates. As a byproduct, we derived a decision path for the classification of low-resolution spectra of massive star candidates. The path relies on the equivalent width and the presence or absence of specific spectral lines.

This paper is organized as follows: in Section  2, we describe the parental catalog along with its characteristics, and we describe the sources of spectroscopic labels. In Section  3, we explain the methodology employed to identify moderate to high probability candidates. In Section 4 we present the BRF results and the moderate to high probability candidates of spectral type B2 or earlier. In Section 5, we discuss the validation of our methodology, including external spectroscopic diagnostics and the comparison with databases and catalogs specialized to massive stars. Finally, in Section 6, we present the conclusions and summarize the main results of this work.

2 The data

In the following we present the different sets of data gathered to develop our classification flow, namely, the source of attributes and lables.

2.1 OBA GDR3: Parental catalog

The Zari et al. (2021) catalog consisted of sources in Gaia Data Release 2 (Gaia DR2, Gaia Collaboration et al. 2018 ) crossmatched with 2MASS (see Section 2.1 of Zari et al. 2021). Posterior to the matching, the authors updated the Gaia photometry and astrometry to those of the Gaia Early Data Release 3 (Gaia eDR3, Gaia Collaboration et al. 2021). Their selection criteria filtered stars brighter than $G<16$ mag, and then applied simple photometric cuts to select stars of approximate spectral type OBA.

To define our parental catalog, we took advantage of the availability of Gaia DR3 and the existing Gaia DR3–2MASS crossmatch111ESAC’s external crossmatch catalog service. We applied their selection function to this updated release. From now on, we will refer to this catalog as OBA GDR3, which contains 1016743 entries, $\sim$3% more than OBA GR2. This difference maps naturally to the difference in “seed” catalog as Gaia DR3, as a whole, contains $\sim 7\%$ more sources than Gaia DR2.

Figure 1 presents a color-color diagram of the entire $G<16$ mag population in Gaia DR3, together with the OBA GDR3 subsample corresponding to a slightly updated version of Fig. 2 in Zari et al. (2021).

2.2 OBA GDR3: Features description

We used photometric, astrometric, and color information as features in our analysis. Regarding photometry, the Gaia bands ($G$, $G_{BP}$, and $G_{RP}$) covered the optical domain and the 2MASS ones ($J$, $H$, and $K_{s}$), the near-infrared. In addition, we considered 15 color indexes resulting from the combination (differences) of pairs of photometric bands. The astrometric parameters considered were parallax ($\overline{\omega}$), proper motion in right ascension ($\mu_{\alpha}$) and declination ($\mu_{\delta}$) with their respective errors ($\sigma_{\overline{\omega}}$, $\sigma_{\mu_{\alpha}}$, $\sigma_{\mu_{\delta}}$).

Finally, we incorporated six quality indicators that measure the goodness of fit of the astrometric solution, such as the astrometric_sigma5d_max (ASM), astrometric_excess_noise (AEN), Renormalised Unit Weight Error (RUWE), astrometric_gof_al (AGOFAL), ipd_gof_harmonic_amplitude (IPDGOFHA), and astrometric_excess_noise_sig (AENS). Extreme values of these features have been interpreted in the literature as sings of multiplicity (Belokurov et al., 2020; Dodd et al., 2024). For more detailed information on these parameters, refer to the Gaia Data Release 3 documentation.222GDR3 Documentation.

2.3 OBA GDR3: Spectral type compilation

The aim of this work is to identify reliable early type candidates based on the features presented in the previous section. We tackle this problem as a supervised classification one with a taxonomy based on spectral types. Obviously, not all objects in the OBA DR3 catalog have spectral type determinations, but a subset do and those will compose the “universal truth” that we want to propagate with our methodology. The spectral types will be referred to as labels from now on, and we used two sources for them: the spectroscopic catalog from Skiff (2014), hereafter referred to as Skiff and wider and very commonly used database, Simbad (Wenger et al., 2000). We note that Simbad is a manually curated survey that may include mismatches and fall behind the current literature; however, it is the most complete stellar database available at present.

To assign labels to the entries of the OBA GDR3 catalog, we performed two independent crossmatches: to Simbad and Skiff. In both cases, we keep every entry from OBA GDR3 and employed a 2″radius for coordinate matching.

Regarding multiple counterparts, over 95% of targets in Simbad had a one to one correspondence in the defined radius, and hence, in case of multiple-match, we opted to use as “true” counterpart, the closest target in terms of angular distance. The situation was drastically different in the matching with Skiff, where over 70% of the targets have indeed multiple counterparts. This was expected, as Skiff, by construction, does not prevent multiple entries for the same astrophysical object, in Sec. 2.3.2 we present our approach to deal with duplications.

2.3.1 Numerical encoding

The labels provided in Skiff, and Simbad are alphanumerical quantities that we proceeded to encode in numerical tuples to preserve the ordered relation of the labels (an O stars is closer in properties to a B star than to an A and so forth). These encodings are used in Section 3 to train the model and revised in Section 5.3 to assess possible limitations.

For the first element of the tuple, the letter, O to M, was mapped to integer values in increments of 10, ranging from 0 to 60. If a subtype was reported in either catalog, its integer value (0–9) was added. Otherwise, in cases where a subtype was not present, an integer was sampled from a uniform distribution of possible values and then added. For instance, a star classified only as type B was initially encoded as 10. We then sampled a random integer from 0 to 9; if the result was 1, the final encoded label was set to 11, corresponding to B1. We noticed that this may introduce noise in the labels; however, this uncertainty with the subspectral types only affects 1.6% of the labels.

In addition, whenever uncertain values or ranges (e.g., B9/A0, F1-F3) were reported, the individual encoding for the extremes were computed and used as the inclusive limits for a distribution to be uniformly sampled. It must be noted that this sampling artificially “produces” spectral types for which no spectra of standard stars may exist (Ma´ız Apellániz et al., 2024), however, the classifier (presented in Sec. 3) will be trained on a lightened binary taxonomy “earlier or later than B2”, and hence this (ease of computation motivated) artifact can be neglected. Finally, when the uncertainty in the spectral classification spanned more than one main spectral type (e.g., OB), the label was saved as a known label but excluded from the training, while targets with reported composite spectral types (such as kA2hA9mF2333This corresponds to A2 for the Ca K-line subtype (k), A9 for the hydrogen-line type (h), and F2 for the metallic-line type (m) ) were discarded.

Although Table 4 provides the full encoding of the spectral types here described, when training the Random Forest classifier we simplified our labels to two: 0 for encoded spectral type ¿ 12 (‘negative” class), and 1 for encoded spectral type $\leq$ 12 (‘positive” class).

Since the different elements used as attributes among themselves (photometry and kinematics) and labels are not coetaneous, for objects exhibiting extreme variability this could pose some source of noise in the modeling. However, we expect the number of such cases to be negligible. For example, the fraction of extreme variability cases reported in the VISTA Variables in the Via Lactea (VVV) survey (Minniti et al., 2010) is below $10^{-5}\%$ (Lucas et al., 2024).

2.3.2 Addressing repeated entries

While the (2″radius) crossmatch of OBA GDR3 data with Skiff may yield multiple entries for the same astronomical sources, either due to different sources or to multiple spectral classifications of the same source, the crossmatch with SIMBAD does not. Thus, in order to faithfully compare both crossmatched datasets against each other in the following subsection, we first performed an internal match within the Skiff-correlated set with respect to its astrometric coordinates, as to identify and group multiple records that describe the same object. The match threshold criterion was a maximum radius of 2″of angular distance. For each group, we computed the mean and standard deviation of their spectral type encodings. Specifically, grouped entries with standard deviation values greater than $1.0$, $3.0$, and $5.0$ (i.e., $1$, $3$, and $5$ units of spectral subtype) were discarded, resulting in three filtered sets. This standard deviation will hereafter be referred to as $\sigma_{\text{literature}}$.

2.3.3 Comparing the two sources of labels

We compared Simbad and Skiff in detail to assess the differences between these two databases. Fig. 2 shows the number of spectral types recorded in Skiff and Simbad, considering the previously discussed $\sigma_{\text{literature}}$. We collected spectral type information for approximately 55000 stars, which represented $5\%$ of the OBA GDR3 sample. In all cases, Simbad contains more stars with spectral type information than Skiff. For O stars in Skiff, the number increases by $30\%$ when the $\sigma_{\text{literature}}$ threshold rises from 1 to 3 and by $5\%$ when it increases from 3 to 5, with a similar trend observed for other spectral types. Based on this, we adopt a $3-\sigma_{\text{literature}}$ threshold for grouping Skiff data, achieving a balance between robustness and sample size.

The complete description of the encoding for the luminosity class can be found in Appendix A. On Fig. 3, however, we do incorporate this information in the comparison between Simbad and Skiff to offer a complete view. As can be seen in the figure, both databases exhibit a similar distribution, with dwarfs being the most abundant class. However, there is a discrepancy in the number of subgiants (class IV), due to uncertainties in the luminosity classifications in Skiff (e.g., IV/V or V/IV), which are encoded as intermediate between these classes (see Section 2.3.1). When combined with the binning process used in the histogram, these uncertainties result in a higher contribution for subgiants.

Figure 4 compares the labels identified for the same object in the two databases for those cases where the luminosity class is in agreement. We included various dispersion metrics for comparison: the Mean Absolute Error (MAE), Root Mean Square Error (RMSE), Coefficient of Determination ($R^{2}$), Pearson Correlation Coefficient ($r$), and Spearman’s Rank Correlation Coefficient ($\rho$), available in Pedregosa et al. (2011). The MAE and RMSE were 0.36 and 0.91, respectively, with differences remaining below one spectral subtype unit. The other metrics were very close to 1, confirming a strong correlation between the two databases and indicating that they were highly comparable.

We identified a few outliers above the bulk of the data and investigated them using both the SIMBAD and Skiff databases. In most cases, we could not retrieve the reference listed in SIMBAD. In a few others, we found misreferences in Skiff or more updated spectral classifications in one of the two databases. We did not exclude them from our analysis, since all these stars are later than B2 and represent only a very small fraction of the labeled data.

Given this strong correlation, we ultimately decided to use the Skiff database for our analysis. This choice is motivated because approximately $\sim$1000 objects in Simbad lack references for their spectral types; $\sim$7000 objects potentially suffer from suboptimal astrometry, with a positional precision of $8$–$12\arcsec$ compared to the $\sim 1\arcsec$ coordinate precision guarantee in Skiff. Finally, Skiff ensures spectroscopic origins for its compilation, whereas Simbad does not. Although it was possible to combine classifications from both databases to increase the sample size, this approach could have introduced additional dispersion or variance due to the inclusion of objects classified through different techniques.

2.3.4 Simplification of encoded labels and adopted threshold

We simplified the encoded labels into more global and representative categories. Since the fundamental truth sample is meant to provide ’standard examples’ of each class, we filtered out phase-dependent spectral types (e.g., OBe), stellar spectra with reported uncertainties in the classification, and non-MK classifications (e.g., Wolf-Rayet stars, white dwarfs, planetary nebulae, and carbon stars). Additionally, ’em’ designations were removed, as they lack information about the primary spectral type. Lastly, we excluded the $\sim$2000 systems reported in the catalog as multiples, due to the limited number of cases (only 2%) in which spectral types for both components are provided. This decision could be seen as at odds with the fact that multiplicity is known to scale with primary mass. However, the uncertainty introduced by the lack of individual spectral types (and hence across the period - mass-ratio parameter space) for these systems was deemed detrimental to the training process.

The ”boundary” between massive and not massive stars is not a simple cut in terms of spectral type for B stars, as the luminosity class has to be considered. However, more than half of the available encoded labels lack this information. Promoting purity over completeness, we set B2 as the threshold, as this spectral type corresponds to a massive star across all luminosity classes (Ma´ız Apellániz et al., 2024).

3 Methodology: Balanced Random Forest

In this section, we used the simplification of the encoded labels introduced in section 2.3 and the features detailed in section 2.2 to address the classification of massive stars as a supervised learning problem. We applied a Balanced Random Forest (BRF), suitable for handling imbalanced training datasets (Chen et al., 2004). We trained the BalancedRandomForestClassifier method from the imbalanced-learn555BalancedRandomForestClassifier using the simplification of the encoded labels from Skiff within the OBA GDR3 sample and applied the classification to the unlabeled data.

3.1 Preprocessing Features

We removed sources with NaN values in any feature or duplicate $J$, $H$, and $K_{s}$ photometry. The percentage of duplicate sources using the tmass_oid identifier was $0.013\%$. Then, we defined the minimum and maximum values for each feature by calculating the 1st and 99th percentiles of the encoded labeled data. We only selected the pool of new candidates that fell within these boundaries. This ensured that all unlabeled candidates were within the specified feature restrictions, as shown in Table 5.

We aimed to quantify how representative the OBA encoded labeled data are with respect to the pool of new candidates for each feature. In Fig. 5, we present a similarity test using the Wasserstein distance666https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wasserstein_distance.html between these samples. In Fig. 6, we show the cumulative number of sources per bin of G magnitude.

The similarity test was normalized between the maximum and minimum distance for each bin of magnitude. The most distinct feature across all magnitude bins was AENS, followed by $\sigma_{\overline{\omega}}$, the second most important feature for the model (see Section 4). This feature reached a minimum at $G=11$ mag and began to increase for both bright and faint stars. From Fig 6, we observed that for bright stars ($G=10$ mag), the number of OBA encoded labeled sources at $G=10$ mag is limited, which may explain the larger distance in this magnitude bin. For fainter stars, starting around $G\approx 12$, there was no noticeable increase in the cumulative number of OBA encoded labels. This suggested a bias in the labeled sample toward brighter stars, making extrapolation to fainter magnitudes uncertain. To be more conservative, we selected stars brighter than $G_{\text{mag}}\leq 12$ as a compromise between completeness and purity. The final sample consisted of 99388 stars up to $G=12$ mag, of which 35465 had spectral types and 63923 were unlabeled.

3.2 Metrics

We evaluated the model’s performance using Recall, Precision, F1-score, Geometric Mean, and Index Balanced Accuracy (IBA). In the text, we refer to completeness and purity as recall and precision, respectively. The last two metrics, designed for imbalanced datasets, were taken from the imbalanced-learn package (Lemaître et al., 2017) 777https://imbalanced-learn.org/stable/references/metrics.html. The IBA incorporates any metric ($M$) and allows emphasis on specific classes through a weighting parameter ($\alpha$). These metrics helped assess model construction, final classification, and comparison with other models. The metrics are presented as follows:

[TABLE]

Where TP are the number of true positives, FP the false positives, FN the false negatives and TN true negatives. See subsections 3.3 and 4.1 for the use and results of these metrics.

3.3 Building the model

The random forest (RF) is a tree-based algorithm that consists of an ensemble of decision trees (Breiman, 2001). Each tree is constructed using a random subsample with replacement (a bootstrap sample) from the training set, and the final prediction is determined by a majority vote among all the trees.

The structure of a decision tree consists of a root node, internal nodes, and terminal (leaf) nodes. Data is classified based on decision rules inferred from features at different levels, from the root node (level 0) to the terminal nodes (maximum depth). A terminal node is reached when all training examples within it belong to the same class. These components are defined by hyperparameters: n_estimators indicates the number of decision trees, max_features denotes the subset of features used to determine the best split at each node, and max_depth sets the limit for recursive division.

Traditional RFs show high performance in different problems, however they underperform on highly imbalanced datasets since bootstrap samples often underrepresent the minority class. Instead, we used the Balanced Random Forest (BRF), a variation of RFs that modifies the sample distribution to ensure equal class representation in each tree. We selected the BRF from the Imbalanced-learn module, which offers various under-sampling strategies for bootstrapping: majority, not minority, not majority, all, and auto. We tested independent strategies for binary classification: ’majority’, ’not majority’, and ’all’.

To optimize BRF construction, we used Grid Search CV with stratified k-fold cross-validation. This method systematically explores hyperparameter combinations from a predefined grid. It divides the data into ’k’ folds, using ’k-1’ folds for training and the remaining fold for testing. We used a 5-fold stratified cross-validation and selected the IBA combined with $Gmean^{2}$ as the scoring metric, setting $\alpha=0.1$, as recommended in García et al. (2012). Table 2 presents the tested and optimal hyperparameters.

4 Results: Candidates from BRF

In the following, we use the built model and the preprocessed features to train a BRF and extend the classification to the unlabeled data. We also present the performance metrics of the algorithm and the moderate to high probability massive star candidates.

4.1 Performance of the BRF

We defined the BRF using the best combination of hyperparameters and evaluated the performance of the model roughly following Sánchez-Sáez et al. (2021). We created 20 random permutations of the training and test sets using StratifiedShuffleSplit from scikit-learn, allocating $80\%$ of the data for training and $20\%$ for testing in each permutation. This corresponded to 28372 sources in the training set and 7093 in the test set per split. Across the full dataset, 4393 sources were classified as B2 or earlier, and 31072 as later than B2. Subsequently, we trained 20 BRFs, each on a separate splitting iteration using the training sample, and evaluated their performance with the respective test sets. We used a probability threshold of 0.6, above which stars are classified as spectral type B2 or earlier. We did not include oversampling techniques, as the ratio between the majority and minority classes is approximately seven, significantly lower than the $\sim$42 reported in the original SMOTE paper (Chawla et al., 2002).

Figure 7 summarizes the results obtained for the 20 random permutations for each test set. TP, TN, FP, and FN values are averaged over these permutations, with the error defined as their standard deviation. We obtained a TN fraction of $\sim 88\%$ and an FP fraction of $\sim 12\%$. The TP fraction was $\sim 80\%$, while the FN fraction was $\sim 20\%$. We measured a standard deviation of $0.47\%$ for stars with spectral type later than B2, lower than that of stars with spectral type B2 or earlier.

Figure 8 shows the variations of the relevance of each feature in the deciding power of the 20 models. When the randomness of the different models is taken into account, the relevance of most features are indistinguishable, with two exceptions: the parallax and the parallax-over-error.

The H–G color was the least informative for the model’s performance, although it traces the slope of the stellar spectrum. We explored this effect using the spectral energy distributions (SEDs) of 100 nearby sources ($<$1kpc) with high predicted probabilities of being of spectral type B2 or earlier. These SEDs were fitted using VOSA (Bayo et al., 2008), allowing extinction to vary freely within the catalog provided limits. The results showed that H–G is highly degenerate with extinction, which explains its low relevance with respect to effective temperature.

We trained the final BRF classifier using a single data split, motivated by the small variation observed in the confusion matrices across different $k$-folds. All features were used to train the model, which was then applied to classify the OBA candidates in GDR3. We used $80\%$ of the data for training and $20\%$ for testing, applying the best hyperparameters from Table 2. Fig. 9 shows the kernel density estimates (KDE) for our test sets by metric. The left panel presents the BRF probability, while the right panel displays its distribution across spectral types. The TN peak was near 0.15, transitioning smoothly from lower probabilities before distorting as it approached 0.5, where the model had the highest classification confusion. For stars with probabilities above 0.6, TP and FP clearly overlapped. In the right panel, the metrics were associated with the spectral type, with the TN peak corresponding to A0 stars, while the FP peak was centered around B5 stars.

The peak in FP around spectral type B5 might be explained by the subtleties in spectral classification. The main spectroscopic difference between a B5 and a B2 star is the disappearance of silicon lines (very weak absorption features), accompanied by a strengthening of the Mg II 4481 line. Given that the model is informed only with photometry (“even worse”, broad band photometry), the model has not been exposed to data that traces such subtle differences and, therefore struggles to distinguish between these spectral types.

Table 3 presents the performance metrics for the test sets. We observe high recall for stars with spectral types later than B2 and B2 or earlier. While precision remains high for stars later than B2, it is lower for B2 or earlier stars, indicating a high fraction of false positives. Fig. 10 shows how precision improves at higher probability thresholds. In addition, as noted in Section 2.3, the classification of massive stars based on spectral type is a working approximation, typically set at B2 or earlier. However, some authors, such as Zinnecker & Yorke (2007) (see Table 1), have considered B3 or early stars as massive. Therefore, we also presented the results using this alternative threshold in Appendix B. Finally, this trade-off between high completeness and lower purity can be addressed, as the final classification of massive stars requires spectroscopic confirmation.

4.2 Moderate to high probability massive stars candidates

From the sample of 99 388 stars with $G\leq 12$ mag, we classified the 63 923 unlabeled candidates using the BRF. We show the galactic distribution of our candidates in Fig. 11. The colors indicate the probability of classification as a star with spectral type B2 or earlier. We observed that the highest probabilities are more concentrated near the Galactic plane, even without RA and Dec coordinates as input for the algorithm.

Figure 12 shows the probability distribution for Em, OB, and OBe stars, which are the most abundant simplified encoded labels that were not included in the training or test sets. The distributions for OB and OBe stars are skewed toward higher probabilities, reflecting their predominance as early-type stars. In contrast, the Em distribution shows two peaks at low and high probabilities, suggesting the presence of both early- and late-type stars with emission. These trends suggest that the simplification of the encoded labels captured the global phenomenon without an evident systematic pattern for phase-variable spectral types.

Other spectral classes not considered in the training were excluded due to the feature space cuts applied during sample construction. As a result, we do not have a statistically significant number of objects from these classes to evaluate model performance reliably. In addition, assessing the classifier’s performance on Be stars is dificult because the Gaia and 2MASS photometry, the labels, and the spectroscopic data are non simultaneous. Since emission lines are time dependent features, this prevents us from drawing reliable conclusions. However, we expect Be stars to show distinct patterns in the Gaia feature space (Radley et al., 2025).

From the $\sim$ 100000 stars analyzed, we identified 9089 moderate to high probability massive star candidates (probability ¿ 0.6) that were not included in our encoded labeled dataset. According to Table 3, this threshold corresponds to an expected precision of $\sim 0.5$ for identifying massive stars. Based on the FP and TP rates shown in Fig. 9, we expect the FP to be centered around B5 stars. This sample is more than twice the number of massive stars used for the development of the model, significantly expanding the set of potential massive stars and representing a valuable target for spectroscopic follow-up.

5 Discussion

In the following we discuss the validation of our methodology incorporating external spectroscopic diagnostics, and the characteristics of our sample with respect to the literature.

5.1 Multiplicity considerations

The process of simplification of the numerical encoding of labels discussed in Section 2.3.4, implied the removal of multiple systems. However, as previously mentioned, multiplicity is a very common feature of massive stars (Offner et al., 2023). In order to assess the performance of our (singled source trained) algorithm in detecting multiple systems with a massive primary, we analysed the predictions of the model for a sample of 29 confirmed binaries from the Skiff catalogue. The number of systems was limited, as we only considered multiple systems that were not used in the previous analysis.

Figure 13 shows the spectral type of the primary stars against the BRF probability, with the spectral type of the secondary indicated by color. Systems with mass ratio ($q\sim 1$) have similar spectral types and are therefore expected to exhibit similar colors. We observe that most systems with a massive primary and q$\sim$ 1 are identified by our methodology as massive indeed. On the other hand, the recovery efficiency is not $100\%$ as two systems with a B primary and a close to 1 mass ratio lie just below threshold probability. With this small test we can assess that, although we introduced a bias removing multiples from the training set, this does not seem to have a very strong effect (baring the small numbers of the experiment) in correctly predicting multiples with a massive primary.

5.2 Homogeneous validation dataset: LAMOST spectroscopy

We searched for spectral information in the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST, Zhao et al. 2012) to validate the classifier performance. This survey can obtain 4 000 spectra in one exposure for a limiting magnitude of 19 in r band and a spectral resolving power $R\sim 1800$ around the g band. It has a Field of View (FoV) larger than 20 square degrees and a fiber scale of 3.3″on the focal plane. The last public data release (DR10) 888http://www.lamost.org/dr10/ contained low-resolution spectra for 11100139 stars, covering approximately $-10^{\circ}<\text{DEC}<89^{\circ}$, making it a valuable database for large-scale spectroscopic studies.

5.2.1 Automatic spectral classification

We found 1066 spectra in LAMOST within a 5″radius for the moderate to high probability massive star candidates, considering all matches. To enhance the precision of our selection, we only considered counterparts located at the intersection of the initial 5″radius and a 5″wide band along the star’s proper motion direction. We discarded 106 spectra that showed anomalous offsets inconsistent with the expected proper motion direction. The final sample contained 880 spectra, from which we identified 585 unique sources. For each source with multiple observations, we retained the spectrum with the highest signal-to-noise ratio (SNR).

To derive spectral types for this sample, we developed an automatic methodology for classifying the low-resolution LAMOST spectra, based on a ground-truth sample of 50 massive standard stars of luminosity class V with spectral types O, B, and A collected from Didelon (1982); Sota et al. (2011); Ma´ız Apellániz et al. (2016); Martins (2018). The ground-truth spectra were obtained using the High-Efficiency and high-Resolution Mercator Echelle Spectrograph (HERMES; Raskin et al. 2011) and the Fiber-fed Extended Range Optical Spectrograph (FEROS; Kaufer et al. 1999), with resolutions of 85000 and 48000, respectively.

We degraded the high-resolution spectra of standard dwarfs to match LAMOST’s observational capabilities using the specutils package (Earl et al., 2024)999https://specutils.readthedocs.io/en/stable. We applied convolution based smoothing to the spectrum using a Gaussian kernel, reducing the resolution to R$\sim$1300, which is the expected resolution for LAMOST at He i $\lambda 4471$ (see Figure 2 of Zhu et al. (2006)), an important line for distinguishing O and B stars. We resampled the flux to the LAMOST dispersion grid while preserving the integrated flux and added normal noise to compensate for the increase in SNR caused by the resolution reduction. We used a Wasserstein distance to compare the SNR distributions before and after adding noise. The distance decreased from 0.008 to 0.001, making the SNR of both distributions more similar. In Fig. 14, we show the spectrum of the massive standard star HD 242926, classified as O7Vz (Sota et al., 2011), as observed with HERMES, along with its degraded version at LAMOST resolution. For comparison, we also include the LAMOST spectrum of the O7V((f))z star HD 277878 (Roman-Lopes & Roman-Lopes, 2019).

We selected 18 spectral lines used in Kyritsis et al. (2022) and tested which lines were detected in the standard dwarfs with low resolution. For the detection of the lines, we combined the methodology presented in Bayo et al. (2011) and Kyritsis et al. (2022). The first approach detected the continuum of the line in a completely automatic way, while the second one used static spectral windows, taking into account the surroundings of each line.

We present an example of the detection of the He i $\lambda 4471$ line in Fig. 15, the same procedure was applied to each line. We defined two fixed spectral regions (R1 and R2), in terms of central wavelength and width, to model the continuum, following Kyritsis et al. (2022), and a third one to trace the line (hereafter, the ”line region”). An initial estimate of the continuum was derived from a linear fit to the R1 and R2 regions. This continuum was then refined by applying a one-sigma clipping. Once the final local continuum was identified, we assessed the significance of a detection by comparing the flux within the line region to the dispersion of the continuum. We subtracted the continuum and considered a detection whenever the peak of the line (defined as the 95 percentile of the flux distribution) was above 5 sigma.

We calculated the equivalent width (EW) of the line if it was detected. We considered that the line follows a Gaussian profile and calculated the associated $\sigma$. We performed three types of line fitting to characterize different scenarios: isolated line, blended line, and core emission line using a single Gaussian fit, a double Gaussian fit, and double inverted Gaussian fit, respectively. We selected the best fit using the Bayesian Information Criterion (BIC). This criterion is useful for model selection, giving priority to models with fewer parameters.

The line region was defined as $\pm 4\sigma$ around the center of the fitting model. We found that this threshold provided the best separation between the He i $\lambda$4471 and Mg ii $\lambda$4481 lines in early-type stars when these lines are not blended. For consistency, the same criterion was applied to all other lines. Finally, we calculated the EW over the measured flux (not the fitted line), within the line region.

We applied the automatic methodology to characterize line detection in low-resolution spectra of dwarf massive standard stars. For a given spectral line, its detection can be influenced by the SNR of the spectra. To quantify the impact of SNR on line detection, we selected SNR values sampled from the distribution of moderate- to high-probability candidates and assessed the corresponding detections in the low-resolution spectra of massive dwarf standard stars.

We performed 20 iterations, each using a different SNR value. In each experiment, we calculated the detection rate for all lines across the sample. The resulting matrix, shown in Fig. 16, reports the mean and standard deviation of the detection rates for each spectral type and line. This results are highly compatible with the results from Kyritsis et al. (2022). The He i lines are clearly identified in B-type stars, while the He ii lines are systematically observed in O-type stars. The Si lines remain challenging to detect in both works.

From this matrix, we selected the following lines for classification: He i $\lambda$4471, Mg ii $\lambda$4481, He i $\lambda$4144, He ii $\lambda$4541, He ii $\lambda$4686, Fe ii $\lambda$4233, He i + He ii $\lambda$4026. These lines are consistent with the expected spectral types, and based on them, we created an empirical decision path to classify the LAMOST spectra. It is important to note that we are only able to spectroscopically confirm stars of type B3 or earlier, due to the difficulty of detecting silicon lines at this resolution using an automatic approach, as illustrated in Fig. 16 and in agreement with Kyritsis et al. (2022).

Figure 23 summarizes the empirical decision path using automated line detection to classify O and early B stars. The first selection focused on He ii $\lambda 4686$ and He ii $\lambda 4541$, as these lines were consistently present in O stars. To minimize false positives due to the SNR in B-type stars, we applied a strict detection criterion for both lines. We then used the detection and EW ratio of He i $\lambda 4471$ and Mg ii $\lambda 4481$, which is useful to classify helium-normal stars later than B3 (Gray & Corbally, 2009), with a ratio near 1 for B6-B7 stars (Ram´ırez-Preciado et al., 2020). The Mg ii $\lambda 4481$ line was consistently detected from B3 to A1, while He i $\lambda 4471$ became undetectable from B9 onward and was almost always absent in A1 stars. We combined the detection of He i + He ii $\lambda 4026$, which is always present in O and B stars, with the detection of He i $\lambda 4144$, primarily observed up to B3. For stars later than A1, we used the Fe ii $\lambda 4233$ line, characteristic of A1 stars, and kept the designation ”later” since we did not include spectra of standard stars of later types.

5.2.2 Estimated spectral types for the validation data-set

We applied the automatic methodology and empirical decision path to classify the LAMOST spectra of moderate to high probability massive star candidates. Fig. 17 presents the classification percentages, while Fig. 18 illustrates examples of spectra classified with an SNR in the g band detection in SDSS of approximately 250. We classified $57.2\%$ of the sample as massive stars. Approximately $19.6\%$ were misclassified as mid B-type stars (B4–B7). The results could improve by considering the luminosity class; however, developing an automated methodology for luminosity classification is beyond the scope of this work. We refer to Negueruela et al. (2024) for the standard classification of B stars. Additionally, $16\%$ of the sample corresponds to stars later than type A, representing the lower limit of the model’s misclassification.

Figure 19 illustrates the number of stars classified into different spectral type groups, separated by the probability assigned by the BRF. At the top of each bar, we present the detection rate of massive stars, calculated as the percentage of the O-B0 and B0-B3 groups relative to all other spectral types within each probability interval. The detection rates align with the expected precision of $0.564\pm 0.009$ for detecting B3 or earlier spectral types, as shown in Table 3. We observe that the detection rate increases with higher probabilities, while the total number of classified stars decreases. This reduction in moderate to high probability candidates is due to the selection function of the LAMOST survey, which does not cover the Galactic center, where a higher concentration of massive star candidates is expected.

5.3 Comparison with state-of-the-art massive star data bases and catalogs

We selected dedicated studies focused on the spectroscopic classification of massive stars to evaluate the performance of our BRF algorithm. Specifically, we examined the IACOB spectroscopic database (Simón-D´ıaz et al., 2011a), GOSSS (Ma´ız Apellániz et al., 2013), Chini et al. (2012), Li (2021), and the ALS compilation (Pantaleoni González et al., 2021). Within a 2″search radius, we identified 124, 36, 17, 29, and 831 stars, respectively, that were not included in our training or test sets.

Figure 20 presents the probability distribution for different studies. The GOSSS and Li (2021) samples are more concentrated at high probabilities, as they consist exclusively of O-type stars. In the GOSSS sample, we classified all stars with a probability greater than 0.5, whereas in the Li (2021) sample, we misclassified two O9 stars. The comparison with the Chini et al. (2012), showed the poorest classifier performance, with most misclassifications resulting from intermediate B-type stars present in this sample. Among the five misclassified stars, three were B4 V stars, one was a B2 V star, and one was a B1 III star. In the IACOB and ALS samples, we classified approximately 10% and 20% of the stars, respectively, with probabilities below 0.5. However, not all stars in these studies had luminosity classifications, and some belonged to later B spectral types. In general, the probabilities were high, with a missed detection rate consistent with the algorithm’s expected metrics.

6 Summary and Conclusions

In this study, we present a moderate to high-probability catalog of massive star candidates with estimated spectral types B2 or earlier in the Milky Way. We defined our parental catalog using stars with spectral type O to A with astrometry and photometry from Gaia DR3, combined with photometry from 2MASS. We defined a subsample spectroscopically confirmed from the literature to serve as a source of labels. Since the labeled sample is biased toward brighter stars, we restricted our sample to $G<12$ mag. In addition, we discarded possible outliers by discarding objects with any feature outside of the 1-99 percentile range.

A BRF classifier was trained on the encoded labeled data, simplifying the taxonomy to “B2 or earlier” or “later than B2” spectral types. The trained model was used to predict spectral types in this binary fashion on the unlabeled sample. The performance of our classifier shows high completeness of $0.8$ and a purity ranging from $0.51\pm 0.02$ for probabilities between 0.6 and 0.7, up to $0.85\pm 0.05$ for the 0.9–1.0 range.

We mined spectroscopic archives for data on our moderate to high-probability candidates (estimated “B2 or earlier spectral type” with $P>0.6$). The most homogeneous subsample obtained consisted of 534 objects with available spectra from the LAMOST survey. We developed an automatic methodology to detect the presence and characteristics of spectral lines on a sample of dwarf massive standard stars. Based on the extracted line(s) information, we designed a decision path capable of differentiating stars with spectral types O–B0, B0–B3, B4–B7, B8–A1, and later than A1. Applying this scheme to the 534 validation set, we obtained a confirmation rate of $\sim$54%. We further explored the dependence of this confirmation rate on the probability cut form the model, reaching 83% for $P>0.9$. From the validation set (534 source), to the best of our knowledge 107 are completely new discoveries.

As a result of our work, we obtained a catalog of $\sim$9000 moderate to high probability $(P>0.6)$ “B2 or earlier” bright ($G<12$ mag) stars. From those, 5158 are not present in known massive stars databases (or general databases as Simbad). The remaining, 3931 were not present in our homogeneous set of encoded labels but were reported in other works in the literature. The number of new candidates presented in this work (5158) is therefore comparable to the currently known number of stars spectroscopically confirmed with spectral type B2 or earlier, 7500 and 5793 for G ¡ 12 mag respectively (Simón-D´ıaz et al., 2011a; Chini et al., 2012; Ma´ız Apellániz et al., 2013; Skiff, 2014; Pantaleoni González et al., 2021; Li, 2021). Assuming our metrics hold for the full sample, we expect to increase the known population of massive stars by approximately 50%.

Acknowledgements.

We thank the anonymous referee for the helpful comments that improved the quality of this work. We acknowledge the use of TOPCAT (Taylor, 2005) and the following Python libraries used in this research: pandas (pandas development team, 2020), matplotlib (Hunter, 2007), numpy (Harris et al., 2020), ray (Moritz et al., 2017), scikit-learn (Pedregosa et al., 2011), imblearn (Lemaître et al., 2017), astropy (Astropy Collaboration et al., 2013, 2018, 2022), scipy (Virtanen et al., 2020), and specutils (Earl et al., 2024). We acknowledge support from the National Agency for Research and Development (ANID) through the Scholarship Program DOCTORADO BECAS CHILE/2021 - 21211323. A.B acknowledges support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2094 – 390783311. Based on observations made with the Mercator Telescope, operated on the island of La Palma by the Flemish Community, at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. Based on observations obtained with the HERMES spectrograph, which is supported by the Research Foundation - Flanders (FWO), Belgium, the Research Council of KU Leuven, Belgium, the Fonds National de la Recherche Scientifique (F.R.S.-FNRS), Belgium, the Royal Observatory of Belgium, the Observatoire de Genève, Switzerland and the Thüringer Landessternwarte Tautenburg, Germany. Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences. 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.

Appendix A Complementary spectroscopic codification

For spectral classifications that contain additional information beyond the spectral type and subtype, such as luminosity class and spectral peculiarities, we provide complementary data in the catalogue. This information is included only for targets that were not discarded during the encoding process described in Section 2.3.1. Each luminosity class, from hypergiants to subdwarfs, was assigned an integer value from 0 to 5 (as hypergiants and supergiants were grouped with the value of 0). For intermediate classes (e.g. II-III, II/III) the mean value, rounded to one decimal was adopted (for instance 2.5 was adopted for III/IV). In addition to the tuple, reported spectral peculiarity information was propagated by means of five, not mutually exclusive, flags (uncertain, miscellaneous, line-broadening, emission, and nitrogen/carbon enriched, see Table 4)).

Appendix B Effects on the spectral type threshold in the BRF

The limit of massive stars can be different depending on the work and on the luminosity class, as we mention in Section 2.3.4. During the body of the work “B2” was used as a threshold. However, in this appendix we show the impact that a slight change in threshold to “B3” would have in our results. To this effect, we trained a new BRF following the same procedure described before but using B3 as the boundary between the two classes. The confusion matrix of this new model is shown in Fig. 21. In addition, we also performed the same flow described in the text for external validation with LAMOST spectroscopy. The final results of this external validation are shown in Fig. 22.

Even though the difference in performance of the algorithm is not significant when comparing the two confusion matrices (Figs 7 and 21), the external validation procedure seems to indicate an underestimation of the contamination in the 0.6-0.7 probability bin for the case where B3 was used as a threshold. Therefore, thresholding with ”B2”, provides more robust results.

Appendix C