© The Authors 2025

1. Introduction

Massive stars (M ≳ 8 M_⊙) represent the minority of stars in the local Universe (Salpeter 1955; Chabrier 2005). Despite their low relative number, their influence on the interstellar medium and the chemical evolution of their host galaxies is enormous (e.g., Maeder & Meynet 2000; Langer 2012; Nomoto et al. 2013; Smith 2014). Furthermore, they are among the most luminous stars in the Universe, making them useful tools as distance indicators (Kudritzki et al. 1999, 2003), and the main stellar contributors to the brightness of high-redshift galaxies (e.g., Dijkstra & Wyithe 2007). Despite their influence and importance for many branches of astrophysics, the physical properties and processes that govern their evolution are not yet fully understood.

As for less-massive objects, massive stars spend most of their lifetimes in the main sequence (MS), a phase in which stars burn their core hydrogen while slowly becoming cooler and larger in size. Upon termination of the MS, stars undergo a rapid reconfiguration of their cores, which translates into an important reduction of the star’s surface temperature and a significant increase in their sizes. In the classical Hertzsprung-Russell (HR) diagram, the lack of massive stars with intermediate effective temperatures (between those characterizing the stars on the MS and the red supergiant phase, respectively) caused by this quick transition is known as the Hertzsprung gap.

Traditionally, blue supergiants (BSGs; B-type stars with luminosity classes I and II) were considered to be post-MS objects following MS O-type stars. However, the discovery of an overdensity of BSGs by Fitzpatrick & Garmany (1990) within the Hertzsprung gap raised questions about the true extension of the main sequence of massive stars and the true nature of BSGs. One possibility that would solve the still persistent uncertainty is considering that a significant fraction of BSGs are MS objects. This would imply a reconsideration of the MS extension in the evolutionary models. In this regard, the length of the MS depends not only on the stellar mass and the initial rotation rates of the stars (Meynet & Maeder 2005; Heger et al. 2005; Georgy et al. 2013) but also on different proposed internal mixing processes that can provide additional hydrogen to the star’s core (Maeder 1981, 2009; Schootemeijer et al. 2019; Martinet et al. 2021), thus extending the MS. Furthermore, massive stars are affected by strong stellar winds, which can also affect their lifetimes through the effect of mass loss (Kudritzki & Puls 2000; Vink et al. 2010; Smith 2014). Another possibility is that additional evolutionary channels contribute to the observed overdensity of BSGs. These may include stars evolving blueward, such as those that undergo a blue loop after the red supergiant (RSg) phase (Stothers & Chin 1975; Ekström et al. 2012), but also stars that have evolved through binary interactions (Sana et al. 2012; de Mink et al. 2013; Menon et al. 2024). In this entangled context, observational constraints are the key to improving the situation.

Since the work by Blaha & Humphreys (1989) and Fitzpatrick & Garmany (1990) using photometric data of stars in the Galaxy and the Large Magellanic Cloud, only the study conducted by Castro et al. (2014) using a statistically significant sample of massive stars with derived spectroscopic parameters from the literature has attempted to better constrain the width of the MS. Although their sample comprised ≈400 Galactic OB stars and their results evidenced the extension of the MS, the presence of observational biases and the lack of homogeneity in their analysis prevented them from quantitatively providing the location of the terminal-age MS (TAMS) for masses greater than 20 M_⊙. Consequently, high-quality statistically significant samples of BSGs homogeneously analyzed and for which additional spectroscopic information is available (including spin rates, binary status, and surface abundances) are key to improving the current situation, which is the main goal of this work. Following Martinet et al. (2021), we define the TAMS by the position in the MS band where the evolutionary tracks reach their lowest effective temperature before momentarily returning to a hotter region (the so-called hook) at core hydrogen exhaustion (see also Castro et al. 2014).

In this paper, we benefit from the spectroscopic observations gathered within the IACOB project (Simón-Díaz et al. 2020) and from the astrometric distances delivered by the Gaia mission (Gaia Collaboration 2016, 2023; Babusiaux et al. 2023) to investigate three empirical aspects that can shed light on the extension of the MS and the location of the TAMS in the 12−40 M_⊙ range of initial masses. These are the distribution of stars in the HR diagram in terms of density and relative number, the spin-rate properties, and the binary status.

This paper is organized as follows. Section 2 presents the spectroscopic sample of stars and data used to carry out this work. Section 3 describes the methodology used to evaluate the completeness of our sample. In Sect. 4 we present the HR diagram of the stars, their empirical properties evaluated in this work, our proposed TAMS compared to the most-used evolutionary models, and a comparison with the results of Castro et al. (2014). In Sect. 5 we discuss the location of the TAMS in the HR diagram and the rotational properties of the sample, as both are connected with further implications for massive star evolution. Section 6 presents the main conclusions and the follow-up work.

2. Sample and data

2.1. Spectroscopic sample

The spectroscopic sample of this work comprises 876 O3- to B6-type stars that cover the ranges of 49−14 kK in effective temperature (T_eff), 4.2−1.8 dex in surface gravity (log g), and 6.0−4.3 dex in luminosity (log(L/L_⊙)). Of them, 678 were selected from the sample of late O- and early B-type stars presented in de Burgos et al. (2023) as stars putatively descended from stellar objects born as O-type stars (in the context of single-star evolution). This sample was later enlarged with additional observations of stars of similar types to increase its completeness (see Sect. 3). The remaining 198 objects were selected from Holgado et al. (2018, 2020, 2022). In total, ≈45% of the sample corresponds to stars classified as B-type supergiants and bright giants and ≈35% to O-type stars. The remaining ≈20% corresponds to B-type stars classified as giants to dwarfs.

All stars in the sample have magnitudes of B_mag ≲ 11, and their optical spectra were collected from the IACOB spectroscopic database and the ESO public archive. In particular, they were acquired using either FIES (Telting et al. 2014, 2.5 m NOT), HERMES (Raskin et al. 2011, 1.2 m Mercator), or FEROS (Kaufer et al. 1997, 2.2 m MPG/ESO) high resolution (R = 25 000 − 85 000) optical (3800−7000 Å) echelle spectrographs. The average signal-to-noise ratio is ≈130 (measured at 4500 Å; for more details, see de Burgos et al. 2023).

Our sample includes likely single stars and single-line spectroscopic binaries (SB1) in which only one set of spectral lines is detected and where the Doppler shifts indicate the presence of a companion. The sample excludes hypergiants, classical Be-type stars or stars showing strong evidence of being surrounded by a disk, and double- or higher-order spectroscopic binaries (SB2+) where the visual inspection of the spectra shows two or more separated components associated with two or more companions. As also indicated in de Burgos et al. (2024a), the exclusion of these objects is due to the inherent restrictions imposed by the models employed in the quantitative spectroscopic analysis (see Sect. 2.3).

2.2. Astrometric distances

We adopted the distances quoted in de Burgos et al. (2023) or obtained them using the same approach. In general, most of them were taken directly from the geometric distances quoted by Bailer-Jones et al. (2021) or, for the few brightest objects, from the parallaxes (ϖ) provided by van Leeuwen (2007) using the Hipparcos data (Perryman et al. 1997). The former also includes a direction-dependent prior and a 3D extinction map of the Galaxy that are used in the determination of the distances. On average, our adopted distances have an uncertainty of ≈150 pc.

2.3. Spectroscopic parameters

We applied the same methodology as in de Burgos et al. (2024a) to derive the spectroscopic parameters for the O9–B6 stars not analyzed there that are included in the present sample. Briefly, the line-broadening analysis was performed using the iacob-broad tool (see Simón-Díaz & Herrero 2014), which allows one to obtain estimates of the projected rotational velocity (v sin i). For the remaining spectroscopic parameters used here, namely T_eff and log g, the analysis was carried out using a grid of unclumped model atmospheres (that is, not considering wind inhomogeneities) computed with the non-local thermodynamic equilibrium model atmosphere and the Fast Analysis of STellar atmospheres with WINDs line synthesis code (FASTWIND v10.4.7, Santolaya-Rey et al. 1997; Puls et al. 2005, 2020; Rivero González et al. 2011). Although we refer to Sect. 3 in de Burgos et al. (2024a) for additional details on the properties of the grid, the main advantage of our analysis comes from the use of a statistical emulator of synthetic spectra, which is combined with a Markov chain Monte Carlo (MCMC) method to obtain the probability distribution function for the parameters mentioned above. For the 198 O-type stars from Holgado et al. (2018, 2020, 2022), v sin i was also obtained using the iacob-broad, while T_eff and log g were obtained using the iacob-gbat/FASTWIND tool (Simón-Díaz et al. 2011; Sabín-Sanjulián et al. 2014). We refer to those works for more details.

2.4. Fundamental parameters

Reliable distances fulfilling σ_ϖ/ϖ ≲ 0.15 in the full sample (> 97% of the total; see Bailer-Jones et al. 2021) were used to derive stellar radii and luminosities for the 876 stars. To do this, we first downloaded the optical B and V magnitudes available from Mermilliod (2006) and the infrared J, H, and Ks magnitudes from Skrutskie et al. (2006, 2MASS). Then, a synthetic spectral energy distribution (SED) corresponding to the solution obtained in the spectroscopic analysis was calculated in FASTWIND. Following an MCMC approach in which E(B − V) and R_V were varied, the observed and synthetic SEDs of the star were compared, allowing us to obtain the total V-band extinction (A_V) with an average uncertainty of 0.07 mag. Combining A_V and the apparent V magnitudes with the distances mentioned above, we computed the absolute V magnitudes. Then, following Eq. (1) in Herrero et al. (1992), we obtained the stellar radii. In a last step, we used the Stefan-Boltzman law ($L=4\pi R^2\sigma T_{\mathrm{eff}}^4$) to obtain the luminosities, for which we achieved an average propagated uncertainty of 0.1 dex.

2.5. Radial velocity measurements

We obtained radial velocity (RV) estimates and peak-to-peak amplitude values (RV_pp) in 605 stars for which multi-epoch spectra in the IACOB spectroscopic database are available. For most O-type stars, we adopted the values from Holgado et al. (2018). The methodology followed there consisted of a cross-correlation between the observed spectra and the best-fitting synthetic FASTWIND spectra for a given set of diagnostic lines. For most O9- and all B-type stars in the sample, our estimates were obtained following the methodology described in de Burgos et al. (2020). In this case, individual lists of several diagnostic lines are selected and optimized for each sub-spectral type. Then, a combination of line-fitting and sigma-clipping routines allowed us to obtain an average estimate for each spectrum. Alternatively, and only for stars with v sin i ≳ 150 km s⁻¹, the estimates were derived using a cross-correlation technique similar to Holgado et al. (2018). Peak-to-peak amplitudes were obtained simply as the difference between the largest difference of all the individual values.

3. Completeness analysis

In de Burgos et al. (2023), a preliminary analysis of the completeness of the available spectroscopic sample was presented, including a first analysis of the number of observed and missing stars to be included in follow-up observing campaigns (see Table 1 of that work). The main goal there was to corroborate that the sample had a statistically significant number of stars up to a certain distance from the solar neighborhood. In de Burgos et al. (2024a), the study focused on presenting the methodology and estimates of some of the spectroscopic parameters of a sample of ≈500 stars, but did not include any analysis of the completeness of that sample. This work includes ≈200 additional stars analyzed using the same methodology plus ≈200 stars from Holgado et al. (2018, 2020, 2022), and we evaluated our completeness against the same reference sample: the Alma Luminous Star catalog (Pantaleoni González et al. 2021, ALS III; Pantaleoni González et al., in prep.).

Figure 1 shows our analyzed sample (green and blue circles) along with missing objects included in the ALS III catalog (gray circles) on a Gaia color-magnitude diagram (CMD). The top dashed black line starting at the 20 M_⊙ evolutionary track and extending several magnitudes of extinction separates our stars at the approximate age at which O-type stars evolve into B-type stars. As shown in de Burgos et al. (2023), all stars above this line should, in principle, be O-type stars or their descendants when evolving as predicted by single-star evolutionary models, independently of extinction.

For our completeness evaluation, we selected only those stars in the ALS III with B6 spectral types or earlier and with a B_mag < 11, as quoted in the SIMBAD astronomical database (Weis & Bomans 2020). This was made to limit the catalog to stars with the same spectral types and magnitude range as in our sample (see also below). We also removed hypergiant, Be, and SB2+ stars identified as such in previous works, as they are not in our sample (see Sect. 2.1). Additionally, we removed missing stars identified as classical Be-stars in the BeSS (Neiner et al. 2011) and BeSOS (Arcos et al. 2021) databases, plus some other stars identified as such in SIMBAD. Last, we removed stars with a distance larger than 4000 pc in Bailer-Jones et al. (2021) or the ALS III catalog and those with a poor distance determination, which mostly corresponds to very bright stars with a large parallax uncertainty (σ_ϖ/ϖ ≳ 0.15).

Table 1 summarizes the completeness of our sample relative to the ALS III catalog. It includes information on the number of stars and the observed percentages relative to the total for four different distances up to 4 kpc. In addition, stars are separated by their B_mag into three groups. The percentage of analyzed stars with B_mag < 9 consistently reaches ≈100% throughout all distances considered. For stars with B_mag < 10, the percentage decreases from ≈90% to ≈70% at 4 kpc (below 1 kpc, the total number is not relevant due to the lack of faint stars). Here, we must also take into account that, despite the high level of completeness indicated by the ALS III authors for stars up to 5 kpc and B_mag ≲ 16, the further the distance, the higher the possibility of highly reddened stars along the line of sight not included in the catalog. In this regard de Burgos et al. (2023) showed that up to 2 kpc, sources within the original sample of O9–B9 type stars could reach up to 8−9 mag of extinction and would still have B_mag < 16, reassuring a good level of completeness in the ALS III up to that distance.

Fig. 1. Gaia CMD including the stars in the sample and other missing stars included in the ALS III catalog (Pantaleoni González et al., in prep.). All the stars have Bmag < 11 and distances below 4 kpc. The two diagonal dashed black lines mark the reddening line of a 20 M⊙ star at the zero-age MS (ZAMS, bottom) and at 30 kK (top). Both are extended up to a Δ(AV) = 2. Stars shown in green include those for which both the use of the maximum and minimum distances result in the star above the top reddening line. Those for which minimum distance locate them below the same reddening line are indicated in blue. Stars missing from the ALS III catalog are shown in gray. Four evolutionary tracks corresponding to 9, 15, 20, and 40 M⊙ are included for reference. They were downloaded from the MESA Isochrones & Stellar Tracks online tool (MIST; see Dotter 2016; Choi et al. 2016; Paxton et al. 2013, 2015, for references purposes) for solar metallicity, no initial rotation, and AV = 0.0. The average uncertainty in the y-axis is indicated in the legend, while it is negligible in the x-axis.

Based on these results, we concluded that the completeness level achieved relative to the ALS III catalog is particularly good for stars with B_mag < 9, and it is sufficiently high up to 3 kpc (> 60%) to study evolutionary aspects of massive stars, such as the extension of the MS. However, we considered a more conservative distance of 2.5 kpc to evaluate this and other aspects. In particular, this distance provides a confident 70% completeness for all stars brighter than B_mag = 11. We also note that the percentages of stars analyzed are approximately similar for all the spectral types and for both hemispheres, minimizing the possibility of biases in that regard (see Appendix A).

4. Results

4.1. The Hertzsprung-Russell diagram: Drop in density

One of the most robust empirical evidence for the location of the TAMS is the expected drop in the density of stars that populate the HR diagram as a function of T_eff. As explained in Sect. 1, this is a consequence of the rapid evolution of massive stars once they reach the end of the MS, the hydrogen is exhausted in their cores, and core nuclear burning stops. In this section, we investigate the occurrence of this feature.

Figure 2 includes three different HR diagrams showing our sample of O3–B6 spectral type stars for which we account with reliable fundamental parameters. Each one considers stars up to a different distance. The central panel shows stars located up to 2500 pc, that is, to our considered best distance in terms of completeness, and includes a total of 670 stars. The right panel extends the distance to 4000 pc and displays 852 of the 876 stars of the complete sample. The left panel is limited to stars up to 1500 pc and includes 221 stars.

Fig. 2. Hertzsprung-Russell diagrams showing the stars in our sample. From left to right, each panel limits the sample to stars within 1500 (221 stars), 2500 (670 stars), and 4000 pc (852 stars) of distance, respectively. The stars with v sin i > 100 km s−1 are indicated with cyan circles, whereas the rest are shown in orange. The average uncertainties are indicated with error bars in the lower-left corner of the panels. Each panel displays four luminosity ranges included for easier interpretation of the other figures. Our proposed location for the TAMS is marked with a purple line in all panels (see Fig. 3 and Sect. 4.1). The evolutionary tracks are taken from the MIST web-tool for solar metallicity and no initial rotation. The dashed green line on the right panel corresponds to a second-order polynomial fit of the points in the tracks where vcrit reaches 220 km s−1. A contour density mesh based on a Gaussian kernel is included to represent the density of points. Fig. 3. Number of stars with respect to their effective temperature. The histograms in the top four rows of subpanels separate stars within the indicated range of luminosities using different colors. Cyan bins indicate fast-rotating stars with v sin i > 100 km s−1. Each histogram displays the total number of stars and of fast-rotating objects (FR). Each column separate stars by their distances, as indicated in the top panel. The adopted position of the TAMS for each luminosity range is included with a vertical dashed line (see Sect. 4.1). The bottom row of subpanels displays the corresponding cumulative distribution function following the color code of the histograms.

Regardless of the distance, the distribution of the stars in the three panels shows an apparent lower density of objects below 23−21 kK. To better characterize this feature, Fig. 3 shows the number of stars against their effective temperature on a grid of histograms that separate stars in different distances and luminosity ranges as indicated in the panels (the latter ranges are also included in Fig. 2). Looking at the histograms for stars within 2500 pc (middle column), those panels evidence a consistent drop in the relative number of stars below 22 kK, being less steep only in the top panel of stars in the highest luminosity range. As is discussed in Sect. 5 (see also Sect. 4.3), we interpret the drop in density as a consequence of stars leaving the MS and therefore outlining the location of the TAMS.

We used the associated cumulative distribution function (CDF) of the central bottom panel of Fig. 3 to establish a uniform criterion for the position of the drop in T_eff. In particular, we took into account the increase in the CDFs at ≈0.15 with increasing T_eff, which is particularly steep for stars with log(L/L_⊙) < 5 dex, where also the number of objects is larger. We decided to adopt this value and the associated T_eff in each luminosity range to define the drop. These four values of T_eff are marked with a dashed vertical line in the central histograms and are also copied into the other histograms with stars up to 1500 and 4000 pc. We paired each T_eff with the average log(L/L_⊙) of each range to analytically describe the TAMS using a linear fit, which is shown in the three HR diagrams of Fig. 2 with a purple line. The fit has the form: log(L/L_⊙) = 0.47 T_eff −5.42 [dex], where T_eff is given in kK.

Figure 3 shows that for stars with log(L/L_⊙) < 5.0 dex the adopted position of the drop is consistent for stars within 1500 pc and 4000 pc. Regarding the stars with log(L/L_⊙) > 5.0 dex up to 1500 pc, we argue that despite the higher degree of completeness, the statistical number of stars is not sufficient to show the drop. For stars with log(L/L_⊙) > 5.35 dex within 4000 pc, it can be seen that the drop is somewhat shifted toward cooler temperatures. This could be the result of missing a non-negligible number of stars that are intrinsically fainter on the hot side of the drop, as mid-B-type stars are intrinsically brighter than their early-B and O-type companions (in other words, caused by the Malmquist bias; Malmquist 1922).

Regarding the distribution of the stars excluded from our sample (see Sect. 2.1), we expect Be stars to only populate the lower panels of Fig. 3, since classic Be stars mostly comprise early-B MS objects with luminosity classes III to V, and to do it evenly with T_eff (see, e.g., Arcos et al. 2018). On the other hand, SB2+ systems represent an important fraction in the O-star domain. For example, Holgado et al. (2022) showed that ≈30% of their sample of O-type stars are SB2 systems (see also Barbá et al. 2017; Sana 2017), while in de Burgos et al. (2023, see Fig. 6) we found that the percentage among early- to mid-B-type supergiants drops below ≈10%. The lack of SB2+ systems could for instance explain the apparently lower number of objects above 30 kK shown in Fig. 2 (see also Sect. 5.2). In addition, the lower number of objects around 32 kK with log(L/L_⊙) > 5.0 dex (more clearly shown in Fig. 3) could be explained considering that O9-type stars, which populate that parameter space, gather the largest number of detected SB2 systems. However, a proper diagnostic of this feature would also require a better understanding on how binary interaction products populate the HR diagram. We also note that the left-skewed distributions shown in Fig. 3, particularly for stars with log(L/L_⊙) < 5.0 dex, partially result from the effect of stars reaching different maximum T_eff values along the ZAMS.

4.2. Spin-rates of the sample

Another observational characteristic proposed to delineate the location of the TAMS is the drop of the upper envelope of the v sin i distribution as a function of T_eff (see Vink et al. 2010; Brott et al. 2011). Figure 4 shows our distribution of v sin i values against T_eff separating stars in the same luminosity ranges as in Fig. 3. To qualitatively evaluate the aforementioned feature, we included all the stars regardless of their distance. The figure indicates the location of the ZAMS (gray areas), the proposed empirical TAMS (dashed colored lines), and outlines the critical rotational velocity (v_crit; solid green lines), that is, the equatorial velocity at which the centrifugal force balances gravity, as a function of T_eff (see below). We also include the median v sin i values for stars within bins of 4 kK (black-edge colored dots), and the associated 25% and 75% (solid black lines) and 10% and 90% (dashed gray lines) percentiles. Although a detailed discussion of the connection with the TAMS is presented in Sect. 5.1.1, here we highlight the three most interesting features of Fig. 4 (see also Appendix C).

Fig. 4. Projected rotational velocities of the stars in our sample against their effective temperature. In this case, the sample is not limited to any distance. The panels separate stars into different luminosity ranges as in Fig. 3, following the same color code. Solid black error bars and central circles indicate the 75% and 25% percentiles and the median values of the stars in bins of 4 kK. The dashed gray error bars correspond to the 90% and 10% percentiles. The dotted line at 100 km s−1 separates fast-rotating stars. The location of the TAMS is also included with a vertical dotted line. The solid green lines correspond to the critical velocity calculated from the models in Fig. 2. Stars showing double sub-peak emission in Hα are marked with additional lime-colored crosses (see Appendix B). The gray area corresponds to the temperature range covered by the theoretical ZAMS in each luminosity range.

A first feature is the separation of the stars into two main components (see also Fig. 10 of de Burgos et al. 2024a): a low-v sin i densely populated group of stars (colored differently) and a less populated group of fast-rotating objects (cyan circles; indicated also in Fig. 2). Following de Burgos et al. (2024a), we adopt the value of v sin i = 100 km s⁻¹ that best separates both components. This separation is well known to be present in the O-star domain and extends toward the early-B supergiant domain (see, e.g., Holgado et al. 2022; de Burgos et al. 2023, and references therein). Using as reference those stars within 2500 pc, we find a decreasing fraction of fast-rotating objects from ≈50% for stars with luminosities below 4.65 dex, to ≈35% for stars between 5.0 and 4.65 dex, and to ≈30% for stars with luminosities higher than 5.0 dex. On average, 38% of the stars located on the hot side of the TAMS belong to the fast-rotating component. A further discussion of both components can be found in Sects. 5.1.2 and 5.3.

Another feature is the absence of fast-rotating objects below a certain T_eff, located very close to the TAMS, which is especially evident in the middle panels. This can also be seen in Fig. 3, where fast-rotating stars are indicated in cyan bins in all histograms. In addition, the T_eff that outlines the lack of these objects seems to shift with luminosity by a few thousand kelvins, with fast-rotating stars extending beyond the TAMS (in evolutionary terms) toward lower luminosities and vice versa. As shown in Sect. 5.3, this behavior is similar to what is predicted in some evolutionary models of stars born with half of their critical velocity.

A third characteristic is the consistent decrease of the maximum v sin i reached by stars with decreasing T_eff within each luminosity range. In fact, the upper envelopes closely follow v_crit, which is shown with a solid green line (shown also in the right panel of Fig. 2) and has been calculated from the same models used in Fig. 2 employing Eq. (2) in Langer (1998). All these features are discussed in more detail in Sect. 5.1.

4.3. Comparison with evolutionary tracks

Three publicly available evolutionary calculations widely used in the literature are those of Brott et al. (2011), Ekström et al. (2012), and Choi et al. (2016). Each uses different assumptions, prescriptions, or approximations to control internal physical processes including, among others, angular momentum transport from the core to the surface, rotational mixing, or core overshooting. Figure 5 displays the evolutionary tracks of these models without initial rotation (top panel) and with considerable initial rotation¹ (v_ini/v_crit = 0.4; bottom panel). To compare our proposed empirical TAMS (dashed-dotted purple line) with the equivalent TAMS of evolutionary models, we performed a polynomial fitting to those T_eff − log(L/L_⊙) points on the coolest side of the hook toward the end of the MS (see Sect. 1 and Martinet et al. 2021). Since our empirical TAMS mixes stars with different rotations, the top panel includes for reference an additional TAMS where we removed all fast-rotating stars (dashed-dotted red line; hereafter SR-TAMS). This results in an almost vertical line at T_eff = 22.65 kK (σ = 0.11 kK). Complementary, the bottom panel includes an alternative empirical TAMS when we limit out sample to only fast-rotating objects (dashed-dotted yellow line; hereafter FR-TAMS). In this case, the second-order polynomial fit follows the form: log(L/L_⊙) = −0.001 ${T_{\mathrm{eff}}}^2$ + 0.154 T_eff +1.71 [dex], where again T_eff is in kK.

Fig. 5. Hertzsprung-Russell diagrams showing different evolutionary tracks for solar metallicity together with different predictions for the location of the TAMS as indicated in each legend. The top panel includes tracks with no initial rotation, and in the bottom panel the tracks correspond to vini/vcrit ≈ 0.4. Intervals of the same age difference are marked with different symbols. The gray density mesh in the background is based on a Gaussian kernel and includes our sample stars up to 2500 pc. Evolutionary tracks are limited up to the RSg phase (if applicable).

Starting with the top panel, the comparison with the models by Ekström et al. (2012) and Choi et al. (2016) shows that in the 15−30 M_⊙ range, both predict a hotter effective temperature for the TAMS. However, both agree well below the 15 M_⊙. In contrast, the TAMS predicted by Brott et al. (2011) is closer to ours up to the 25 M_⊙ track, but significantly extends the width of the MS above it. We note that this extended MS, while also present in the other evolutionary calculations, only becomes noticeable above 40 M_⊙. Unfortunately, the number of stars in our sample is not sufficient to draw any decent comparison. We also note that for models computed without initial rotation, the difference between our TAMS and SR-TAMS is not substantial (see below).

In the bottom panel of Fig. 5, our TAMS only coincides with those of Ekström et al. (2012) and Choi et al. (2016) for masses ≲15 M_⊙. Interestingly, our FR-TAMS overlaps well with both sets of models up to ≈25 M_⊙, and for models by Ekström et al. (2012), the agreement continues up to 40 M_⊙. However, models by Ekström et al. (2012) for this initial rotation also predict an important spin-down of the stars near the MS. Therefore, an adequate comparison requires further considerations regarding the evolution of the spin-rate properties of the stars. The models of Brott et al. (2011) extend the width of the MS toward cooler temperatures with respect to our TAMS and FR-TAMS in the entire range of covered masses.

We remark here that our definition of the TAMS has been adopted as the temperature that includes 15% of the coolest observed objects, but this depends, for example, on the number of objects located on the hot side of the TAMS. Therefore, the above comparison is affected by this definition. Despite this, the general comparison with the different models indicates that there is no definitive set of models that can reproduce our empirical TAMS in the full range of masses (see Sect. 5.3, where we discuss the reasons for the observed discrepancies).

4.4. Distribution of SB1 objects in the HR diagram

As explored in McEvoy et al. (2015) using a sample of BSGs in the 30 Doradus region, the lack of SB1 systems in the HR diagram can also provide additional constraints on the location of the TAMS. Moreover, the lack of binaries on the cool side of the TAMS can provide clues about the nature of those objects located beyond the MS.

We followed the guidelines presented in Simón-Díaz et al. (2024; see also Holgado et al. 2018) to separate the observed peak-to-peak amplitudes of the radial velocity caused by the intrinsic variability of the stars from true SB1 systems. In particular, the work by Simón-Díaz et al. (2024) establishes RV_pp thresholds with the effective temperature and the spectroscopic luminosity² for this separation (shown in Appendix D), in addition to several visual guidelines and a minimum number of multi-epoch data to avoid erroneous detections.

We evaluated the presence of SB1 systems in our sample limited to a distance of 2500 pc. Figure 6 shows the location in the HR diagram of the 494 systems with three or more available spectra, where we identified 176 SB1 systems (green circles). It can be seen that they are spread mainly within the MS band of the HR diagram, with only a few stars beyond the position of the TAMS. A further discussion of this difference is included in Sect. 5.2.

Fig. 6. Similar HR diagram as in Fig. 2 but limiting the sample to stars within 2500 pc with three or more available multi-epoch spectra. Green circles highlight those stars identified as SB1.

4.5. Comparison with Castro et al. (2014)

This result represents an important improvement over the previous reference work by Castro et al. (2014), who used a compilation of spectroscopic observations and results from the literature to show the distribution of 439 massive OB stars, for the first time, in a spectroscopic HR (sHR) diagram, substituting the luminosity of the traditional HR diagram by the spectroscopic luminosity. Comparing our results with those of Castro et al. (2014), we first point out that despite the fact that part of their sample populates the same parameter space as ours (see de Burgos et al. 2024a), their sample of stars was not created homogeneously and included targeted groups of stars of interest for specific studies. Our sample also benefits from being analyzed in a homogeneous self-consistent manner, rather than being composed of heterogeneous results available in the literature.

Additionally, observational biases due to the different intrinsic brightnesses of the stars are likely present in their work, whereas our volume-limited sample with high degree of completeness solves this issue. Furthermore, while the sHR diagram is used in Castro et al. (2014), here we make use of the HR diagram, combined with additional information on the spin-rate properties and binary fraction of the stars (Sects. 4.2 and 4.4). Moreover, we consider double the number of stars in Castro et al. (2014) up to 2500 pc within the same T_eff range. We also note that in Castro et al. (2014), the surface gravities were not corrected by centrifugal forces (see Herrero et al. 1992; Repolust et al. 2004), shifting fast-rotating stars toward apparently higher luminosities in the sHR diagram.

Despite the differences, Castro et al. (2014) found a drop in the density of stars below 25 M_⊙ at T_eff ≈ 20 kK, which was also interpreted as the TAMS. However, for stars with 25−40 M_⊙ the drop was located at T_eff ≈ 10 kK. As shown in Figs. 2 and 3, our results indicate that the drop occurs throughout the full range of luminosities and masses considered here, and it is located approximately at the same T_eff (24−22 kK). The overdensity of stars shown in Castro et al. (2014) beyond 20 kK, which is not observed in any of the panels of Fig. 2, is likely the result of a selection bias introduced by considering any star with available parameters from the literature rather than a volume-limited sample. In particular, since mid-B-type stars are intrinsically brighter in the optical compared to early-B- and O-type objects, the aforementioned overdensity is likely caused by mid-B-type stars located at larger distances.

5. Discussion

5.1. Spin-rate properties of the sample

Rotation in massive stars plays an important role in their entire evolution (Meynet & Maeder 2000; Heger & Langer 2000; Maeder & Meynet 2005; Yoon & Langer 2005; Langer 2012). It can enhance internal mixing processes that provide the stellar core with additional hydrogen, allowing stars to extend their time in the MS. Rotation can also trigger increased mass-loss rates, which can affect their evolution and final fates. Therefore, understanding the observed spin-rate properties of massive stars is necessary to better constrain their evolution.

As indicated in Sect. 4.2, two main observable features are known to characterize the spin rates of luminous O- and B-type stars. One is the bimodal distribution of the projected rotational velocity (e.g., Conti & Ebbets 1977; Howarth et al. 1997; Dufton et al. 2013; Ramírez-Agudelo et al. 2013; Holgado et al. 2022; de Burgos et al. 2024a), which includes a main component of slowly rotating stars and a secondary more extended tail of fast-rotating objects. The latter have been proposed to be the products of binary interactions (e.g., Ramírez-Agudelo et al. 2013, 2015; de Mink et al. 2014; Holgado et al. 2022; Britavskiy et al. 2023; de Burgos et al. 2024a). The second feature is the lack of fast-rotating stars below a certain temperature (see Howarth et al. 1997; Ryans et al. 2002; Vink et al. 2010; Simón-Díaz & Herrero 2014; McEvoy et al. 2015; de Burgos et al. 2024a). As indicated by different studies, this characteristic could be interpreted as the empirical location of the TAMS (see Vink et al. 2010; Brott et al. 2011; Martinet et al. 2021; de Burgos et al. 2024a). We explore this possibility further.

5.1.1. The lack of fast-rotating stars beyond the TAMS

Previous representations of the v sin i distribution against T_eff have shown the presence of a clear drop from 300−200 km s⁻¹ to 70−30 km s⁻¹ of the most rapidly rotating objects around 22−20 kK (see Vink et al. 2010; McEvoy et al. 2015; de Burgos et al. 2024a). The distributions shown in Fig. 4 for stars with log(L/L_⊙) < 5.35 dex also evidence this feature.

On the one hand, massive stars leaving the MS should undergo an expansion of their envelopes and, in turn, also decrease their observed surface rotational velocities. However, this expansion is not sufficient to drastically reduce spin rates within the investigated T_eff range where the drop occurs (which is of ≈200 km s⁻¹; see, e.g., the second panel of Fig. 4).

On the other hand, stars with expanding envelopes and high spin rates near the TAMS can get too close to their critical rotational velocities (see, e.g., Meynet & Maeder 2006; Meynet et al. 2007). It is predicted that before reaching that point, stars expel enough material into the medium (e.g., Kurfürst et al. 2014), triggering a spin-down of the star’s rotation (e.g., Georgy et al. 2011; de Mink et al. 2013). As shown in Fig. 4 (see also the right panel of Fig. 2), the critical rotation that we extract from the models approximates the upper envelope of the v sin i distribution, particularly in the two central panels. We note that, except for the top panel, we could, in principle, expect to find some fast-rotating objects that have not yet reached their v_crit.

A third possibility is related to mass loss. In this respect, Vink et al. (2010) suggested that the bi-stability nature of the winds (see Pauldrach & Puls 1990; Vink et al. 1999, 2000) can lead to an important increase in the mass-loss rates at ≈22 kK, and to a subsequent enhancement of the angular momentum loss (known as the bi-stability braking). However, in de Burgos et al. (2024b) (see also Rubio-Díez et al. 2022; Krtička et al. 2024; Bernini-Peron et al. 2024) it is shown that such an increase in mass-loss rates (by a factor of ten to twenty) is not observed³. Therefore, we argue that the lack of fast-rotating stars is not due to the latter proposed mechanism.

Given these considerations, the lack of fast-rotating objects could be the result of the statistically and relatively lower number of stars expected beyond the TAMS, that is, the chances of observing these objects is lower. This could be further supported by the following. Considering that for stars with log(L/L_⊙) = 5.35−4.65 (where we account for more stars), the histograms of Fig. 3 show the largest difference in the relative number of stars when crossing the TAMS, which we locate at an average of 22.5 kK. For stars up to 2500 pc, there are five times fewer stars between those located on the hot and the cool sides of the TAMS within 5 kK. Considering that the proportion of fast-rotating objects located within each side with respect to the total is 1:3, this exercise would result in only six fast-rotating objects below 22.5 kK. We account for two objects that meet these conditions, which we consider sufficient to support the aforementioned possibility. We note that we do not account for the fact that some of the fast-rotating stars may not be seen as such because of a projection effect. Furthermore, if fast-rotating stars are the product of binary interaction, it is expected to occur mainly when at least one of the components of the binary system is still on the MS (see de Mink et al. 2013, 2014). Therefore, it would be difficult to observe these objects beyond the TAMS. Additionally, fast-rotating objects evolving near or beyond the TAMS will eventually reach their v_crit, thus becoming slowly rotating stars. Nonetheless, one cannot rule out the possibility of unconsidered effects or forces that can trigger moderate braking toward the end of the MS. This can include, for example, the effect of magnetic fields (see Keszthelyi et al. 2020, 2022).

In conclusion, we argue that the lack of fast-rotating stars below a given temperature is unlikely to be caused by enhanced mass-loss rates and may have its origin in two possible compatible explanations. One, that beyond the TAMS, the statistical number is not sufficient to observe them; and two, that if fast-rotating objects result from binary interaction within the MS, then it would be difficult to observe them beyond the TAMS. Either way, we argue that the drop of fast-rotating objects can also be used to outline the end of the MS.

5.1.2. The evolution of vsini along the MS

Figure 4 provides the largest v sin i distribution of Galactic OB stars to date, which is combined with the largest T_eff and L coverage of stars homogeneously analyzed. This allows us to confidently extend our discussion of how the stars within each of the two v sin i components may evolve along the MS.

Interestingly, the median v sin i values of the two top panels in Fig. 4 show similar trends toward lower effective temperatures. The values are approximately constant until 30 kK and then slowly decrease toward the cool end. This trend is less obvious in the third panel, but may still be present, and it is certainly not the case of the bottom panel, which we discuss separately in Appendix E. We also note that the first median value may not be so representative, as it overlaps with the position of the ZAMS, where the number of stars is significantly lower than for the other mean values. Excluding that one, our results strongly suggest a slow-braking process occurring along the MS, which has already been empirically supported by the works of Holgado et al. (2022), de Burgos et al. (2024a). Furthermore, the number of stars on the hot side of the TAMS with v sin i ≲ 100 km s⁻¹ (i.e., the low-v sin i, component of the spin-rate distribution) represents 65% of the total number. This large fraction may indicate that massive stars are born with low or mild initial rotations (not more than 0.2 of the critical velocity), as indicated also in the aforementioned works. However, we note that since we do not account for projection effects, that percentage may represent a lower limit.

In Sect. 4.2 we showed the moderate decrease of the highest v sin i values with decreasing T_eff until the TAMS is reached (see Fig. 4; see also Sect. 5.1.1). In general, these values follow the limit imposed by the critical velocity, indicating that fast-rotating objects are ubiquitous within the MS. In this regard, the 90%-percentiles associated with the median values shown in Fig. 4 also provide a better statistical representation of this fact. We also note here that the percentile corresponding to the first median value is not representative, whereas the rest of the values monotonically decrease toward lower temperatures. We noticed that the 90%-percentiles associated with the top panel (stars with log(L/L_⊙) = 5.7−5.35 dex) are more separated from the v_crit than in the other panels. A possible explanation for this could be related to the importance of mass loss with increasing luminosity and its connection with angular momentum loss (Maeder 2009). In this regard, de Burgos et al. (2024b) showed that the difference in mass-loss rates between BSGs with log(L/L_⊙) < 5.0 dex and those with log(L/L_⊙) > 5.3 dex can be more than a factor of ten. If this is the case, stars with higher luminosities would lose angular momentum at a higher rates than their less luminous counterparts.

5.2. The impact of binaries

An important aspect to our understanding of massive stars is the impact that binary interaction has on the interpretation of the observed properties. This is particularly relevant as a result of the large proportion of objects that are predicted to interact with a companion during the MS (see Sana et al. 2012; de Mink et al. 2014).

Regarding the detection of SB1 binary systems, it is important to remark that there is no definitive method to distinguish RV variations produced by stellar companions from those caused by intrinsic stellar variability (see Bowman 2020, for a review). Moreover, this variability can have different physical origins, including pulsations that originate in the stellar interiors or strong stellar winds that originate in their atmospheres (see, e.g., Simón-Díaz et al. 2010, 2024; Aerts et al. 2017, 2018; Burssens et al. 2020). Therefore, the percentage of detected SB1 systems is highly dependent on the adopted threshold used to separate the effect of intrinsic variability from the variations linked to orbital motion. In this work, we followed the approach of Simón-Díaz et al. (2024), which provides a study of the characteristic amplitude of RV variations produced by stellar pulsations in the O-stars and BSGs domains, allowing us to more confidently indicate the likely presence of a stellar companion.

In Sect. 4.4 we present the distribution of detected SB1 systems in an HR diagram, limited to a distance of 2500 pc. Figure 6 shows the proposed TAMS and the lack of SB1 systems beyond it, on the cool side. By evaluating the fraction of SB1 systems on each side of the TAMS, we obtain 39% within T_eff = 30−22.5 kK and 15% within T_eff = 22.5−15 kK, highlighting an important decrease when crossing the TAMS. Interestingly, the average percentage is similar to the 23 ± 6% found by Dunstall et al. (2015) in the 30 Doradus region of the Large Magellanic Cloud. Furthermore, our results approach those found by McEvoy et al. (2015) in the Tarantula Nebula, where SB1 systems were not detected below a given temperature. In this work, we show that SB1 systems are also present beyond the TAMS, but their number is very low. We conclude that the lower percentage of SB1 systems beyond the TAMS can be used to infer its location with respect to T_eff. We also argue that the few SB1 systems beyond the TAMS likely correspond to stars heading to the RSg phase, since this phase practically eliminates all detectable binaries beyond the TAMS (see de Mink et al. 2014).

One of the predicted types of binary interaction involves that one of the stars in the system reaches a more evolved state in which its larger size leads to a mass transfer into the companion. The donor can then become a compact object (i.e., a black hole or neutron star) or a stripped helium star. This interaction may result in an SB1 system, a merge, or in the primary star being kicked by supernova explosion of the secondary. The latter two cases do not produce an SB1 system. In either case, this type of interaction is less likely to occur outside the MS, as it would require both stars to have similar masses to reach the end of the MS at the same time (see Marchant & Bodensteiner 2024). Therefore, one would expect the number of SB1 systems to decrease significantly beyond MS, not only in relative terms (caused by the fast post-MS evolution), but also as a result of this less probable scenario. This situation is supported by our findings. We refer to Simón-Díaz et al. (in prep.) for a detailed evaluation of the evolution of the fraction of SB1 and SB2+ systems across the MS.

Our results also indicate that the percentage of SB1 systems within T_eff = 31−22.5 kK decreases in the upper half of the HR diagram. In particular, from 40% for stars with log(L/L_⊙) < 5.0 to 31% for stars above that luminosity. This case could be partially explained by the higher detection threshold of RV (up to 25 km s⁻¹; see Appendix D) toward our highest luminosities, which would eliminate some true SB1 systems.

Similarly as for SB1 systems, the presence of a drop in SB2+ systems on the cool side of the TAMS could provide additional evidence of its location. Although SB2+ systems were intentionally excluded from our sample, and therefore, we do not count with reliable quantitative spectroscopic analysis of each individual component, we can still use to some extent the available information on their spectral types as a proxy of the corresponding effective temperatures.

The upper and middle histograms of Fig. 7 show the distribution of SB2+ and SB1 systems with respect to their spectral type⁴. In this case, SB2+ systems are also limited to a distance of 2500 pc. SB2+ systems are limited to B-type supergiants (a primary component of luminosity classes I and II) and O9-type stars (of any luminosity class), or with a FW3414(Hβ) < 7 Å. As shown in de Burgos et al. (2023), the quantity FW3414(Hβ) is defined as the difference between the width of the Hβ line measured at three-quarters and one-quarter of the depth of the line and can be used as a proxy of ℒ to select BSGs independently of their spectral classification and v sin i. The selected SB2+ systems include those identified in de Burgos et al. (2023), plus some others found later from the acquisition of new spectra. The SB1 systems are limited to stars in our sample with spectral types O9–B6.

Fig. 7. Histograms of the number of spectroscopic binaries against spectral type or log10(Teff), limited to a distance of 2500 pc. The top panel include SB2+ systems for which the spectral type is taken from the primary star (see also Sect. 5.2 for further details). The middle panel is similar to the top one but including the SB1 systems detected from the RVpp variability. The bottom panel include the same stars as in the middle one, but we used the log(Teff) instead of the spectral type. In the top two panels, the position of the TAMS is highlighted with a patch of purple lines, derived from the range of spectral types of stars within 1 kK of the average Teff of the TAMS, which is shown in the bottom panel as a dashed purple line.

Interestingly, SB1 and SB2+ systems follow a similar overall distribution with respect to their spectral types: a high number of systems from O9 to B1 spectral types, followed by a steep drop at the B2 spectral type. In both cases, this drop roughly coincides with the empirical location of the TAMS (patched areas with diagonal purple lines). To provide a rough idea of the comparison between T_eff and spectral types, the bottom histogram of Fig. 7 shows this equivalency, where the drop of SB1 systems very well matches the position of the TAMS at log₁₀(T_eff) = 4.35 dex (≈22.5 kK; dashed purple line).

These results also provide evidence that the number of SB2+ systems in the BSG domain before reaching the TAMS is lower than for O9 stars. In this regard, given their relatively low number, the lack of SB2+ systems in the sample should not produce important biases with respect to the distribution of stars in Figs. 2, 3. However, the difference between O9- and B0-type stars could partially explain the apparent lower density of objects around 32 kK (see Sect. 4.1). We also notice the larger number of B3–B4 SB2+ systems compared to SB1 and we refer to Simón-Díaz et al. (in prep.) for further discussion.

From this exercise, we conclude that the drop in SB2+ systems coincides with that of the SB1 systems and could also be used as an indicator of stars reaching the TAMS. In fact, we also expect the number of SB2+ systems to decrease beyond the MS as a result of (at least) one of the stars in the system either turning into a compact object or interacting with the companion.

5.3. Implications for stellar evolution models

The drop in density in the HR diagram and the lack of fast-rotating stars and SB1 systems beyond that drop all support the location of a new empirical TAMS at the position described in Sect. 4.1. In Sect. 4.3 we show that this position implies an extension or contraction of the width of the MS depending on the set of evolutionary models considered. Moreover, we showed that the location of the TAMS depends on the selection of stars based on their v sin i. In this section, we further discuss the comparison with widely used publicly available stellar evolutionary models.

We argue in Sect. 5.1.2 that one possibility to explain the observed v sin i distribution and the larger fraction of stars with v sin i ≲ 100 km s⁻¹ is that massive stars are born with low to mild initial rotation (v_ini ≲ 150 km s⁻¹) and continue evolving while slowly decreasing their spin rates. This scenario is more compatible with the predictions of Brott et al. (2011) and Choi et al. (2016), where the effect of surface braking is much lower than in Ekström et al. (2012). Although the models of Ekström et al. (2012) and Choi et al. (2016) do not include an intermediate initial rotation (v_ini/v_crit ≲ 0.2) we continue our discussion assuming that the available models without initial rotation (top panel of Fig. 5) are more representative to compare the positions of the theoretical and empirical TAMS. In this regard, two important processes that significantly impact the width of the MS across all the masses considered here are the above-mentioned initial rotation and the amount of overshooting. Other mechanisms such as semi-convection, inflation, mass-loss rates, rotationally induced chemical mixing, and magnetic fields also influence the MS, some of which are restricted to a specific range of masses (e.g., inflation for masses ≳30 M_⊙; see Ishii et al. 1999; Köhler et al. 2015), or their influence also depends on additional assumptions (e.g., rotational mixing which relies on an efficient angular momentum transport; see Martinet et al. 2021). As shown by Keszthelyi et al. (2022), the effect of surface fossil magnetic fields can largely influence the evolution of stars and certainly deserves further consideration. However, the models of Brott et al. (2011), Ekström et al. (2012), and Choi et al. (2016) do not account for this effect; therefore, we decided not to include it in our discussion.

This leaves the main difference for the position of the TAMS in Fig. 5 set by the overshooting parameter (α_ov). In the models of Brott et al. (2011), this parameter is fixed to α_ov = 0.335, whereas in Ekström et al. (2012), it takes the value of α_ov = 0.10. For Choi et al. (2016), they used the exponential overshoot formalism, which for models without initial rotation is equivalent to α_ov ≈ 0.20⁵. In general, in all available models, the coefficients that regulate this parameter have been set to match either the empirical location of the drop in the relative number of stars (e.g., Ekström et al. 2012), or the drop in rotational velocities (e.g., Brott et al. 2011). The latter would be more risky under the aforementioned assumption of a low initial rotation, plus the additional uncertainty in the case that fast-rotating objects have their origin in binary interaction (see Sect. 5.2). Additionally, the range of available masses used to calibrate the coefficients in those models where overshooting is not mass-dependent is an additional source of discrepancy. This explains, for example, the better agreement between the models by Brott et al. (2011) and Ekström et al. (2012) for intermediate masses below ≈15 M_⊙.

Recent work by Martinet et al. (2021) used the observational information presented in Castro et al. (2014) to investigate the sizes of the convective cores and constrain the amount of overshooting. This and other works have also shown that the size of the cores is mass-dependent (e.g., Langer & Kudritzki 2014; Claret & Torres 2019; Tkachenko et al. 2020, 2024). Using our volume-limited sample and its high degree of completeness, we can contribute to a better definition of the width of the MS by adding new constraints to fine-tune the overshooting parameter and its dependence on the initial mass for evolutionary masses from 12 M_⊙ and up to 40 M_⊙. In this regard, we show in Sect. 4.3 that our empirical TAMS lies between the theoretical TAMS of Ekström et al. (2012) and Choi et al. (2016), and the one of Brott et al. (2011). Although more in-depth comparisons are certainly required, our results suggest higher overshooting values in the first two cases and a lower value above 25 M_⊙ in the latter case. Moreover, the lack of an overdensity of objects with 25−40 M_⊙ below ≈20 kK as shown in Castro et al. (2014) (see Sect. 4.5) also adds an important constraint for future models. Although some authors attributed the observed extension to the effect of inflation (e.g., Choi et al. 2016; Sanyal et al. 2017), this is no longer needed within that range of masses. Unfortunately, the number of stars in our sample is not sufficient to draw conclusions above 40 M_⊙.

In Sect. 1 we mentioned the possibility that stars undergoing a blueward evolution may also be present in our sample. If this is the case and some of these stars overlap with those used here to define the position of the TAMS, it could lead to its incorrect definition. We evaluated this possibility using the models by Ekström et al. (2012) and Choi et al. (2016)⁶. The models of Choi et al. (2016) only predict stars evolving blueward for v_ini/v_crit = 0.4 and initial masses below 12 M_⊙ and above 40 M_⊙, whereas in the models of Ekström et al. (2012), we find their presence for masses below 12 M_⊙ (which corresponds to stars undergoing blue loops) and above 32 M_⊙ or 20 M_⊙ for models without rotation or with v_ini/v_crit = 0.4, respectively. Assuming the slow initial rotation scheme, the models of Ekström et al. (2012) for masses ≳32 M_⊙ predict stars to quickly move blueward and then spend a non-negligible time (≈0.2−0.3 Myr) in the 30−20 kK T_eff range before reaching their final fates. If this is the case, then some of the stars in this T_eff range could contaminate the distributions shown in Fig. 3 for stars with log(L/L_⊙) ≳ 5.3 dex. To further evaluate the nature of these objects (whether they are evolving toward the RSg phase or coming from it), additional clues from the surface chemical composition and the pulsational properties are required (see Bowman et al. 2019).

As also mentioned in Sect. 1, the binary channel can also contribute to the observed distribution of single stars and SB1 systems used to delineate the TAMS (see, e.g., Menon et al. 2024). Unfortunately, our understanding of how different types of binary interaction spread over the HR diagram is still in an early stage. Nevertheless, we may expect that for those cases where the merging of two stars leads to a reconfiguration of the outcome core, the evolution follows the same path as for a more massive single stars, thus not altering the drop in relative number.

The initial stellar rotation also influences the extension of the MS predicted by evolutionary models (see, e.g., Meynet & Maeder 2000; Brott et al. 2011; Martinet et al. 2021). Although, as pointed out in Sect. 4.3, we cannot directly compare our empirical FR-TAMS with the theoretical ones, our large sample compared to previous studies (e.g., Vink et al. 2010; Castro et al. 2014; McEvoy et al. 2015) allows us to withdraw one conclusion. If fast-rotating stars retain their spin rates (regardless of whether they got it as a result of binary interaction, or they were born as such), then rotationally induced effects seem to shorten the width of the MS for higher masses and extend it for lower ones. This behavior could be similar to what predicted by Ekström et al. (2012) for v_ini/v_crit = 0.4. However, in that case, the models by Ekström et al. (2012) also predict an important spin-down on the initial part of the MS, which is not consistent with our observations.

The comparison with the different models therefore suggests that a revision of the overshooting is required in all the cases to match our empirical TAMS, which is a conclusion also supported by asteroseismic studies of massive stars (see Burssens et al. 2023). In the case of Ekström et al. (2012) and Choi et al. (2016), the overshooting should be adjusted to further extend the width of the MS to cooler temperatures for masses between 15 to 30 M_⊙, whereas models by Brott et al. (2011) appear to require further adjustments.

6. Conclusions and future work

The question of whether BSGs are MS or post-MS objects arose over three decades ago as a result of the observed overdensity of BSGs situated beyond the theoretical location of the TAMS. To date, we know that a significant fraction of BSGs are likely hydrogen-burning stars and that a significant fraction of their observed population may be binary products. However, we do not have a clear definition of the location of the end of the MS, a situation that gets worse considering that, apart from single MS stars, additional evolutionary channels can also populate the BSGs domain, making a very difficult task to disentangle the different groups. A key to improving this situation requires high-quality statistically significant volume-limited samples of BSGs analyzed homogeneously.

In this work, we considered 670 spectroscopically observed Galactic O-stars and BSGs. Using the ALS III catalog as reference, we found that our sample comprises ≈60% of all massive stars within 2500 pc and B_mag < 11 that are born as O-type stars and reach the spectral type B6. This results in the most complete sample of Galactic massive stars spectroscopically analyzed to date. Combining the results of the quantitative spectroscopic analysis with Gaia distances, we were able to derive their absolute luminosities. This allowed us to represent the sample in the HR diagram and to analyze its properties in the context of stellar evolution.

Our results evidenced a drop in the density of stars below 25−22 kK and for luminosities between 5.7 and 4.3 dex (corresponding to initial masses between 12 and 40 M_⊙). We interpret this drop as the location of the TAMS. Moreover, this drop is still present when the volume of our sample is reduced to 1500 pc or increased to 4000 pc. This is an important improvement over previous studies that attempted to empirically locate the TAMS. We particularly improve the situation for initial masses above ≈25 M_⊙, where the location of the TAMS was suggested to be at T_eff ≈ 10 kK based on the lower number of objects in the sHR diagram.

We studied the projected equatorial velocities of the stars in our sample with respect to T_eff and log(L/L_⊙). We found an important lack of fast-rotating stars (v sin i > 100 km s⁻¹) on the cool side of the TAMS, almost outlining its position. The lack of these objects has been suggested to be the result of stars leaving the MS, or due to an important braking effect triggered by an increase in mass-loss rates at ≈22 kK. However, the latest findings on mass-loss rates in BSGs do not support the latter scenario. We concluded that the lack of fast-rotating objects can be a consequence of the lower relative number of stars expected beyond the TAMS due to rapid post-MS evolution, combined with the fact that those stars possibly reach their critical velocity near the TAMS, triggering a spin-down of their rotation. Furthermore, this scenario is compatible with the possibility that fast-rotating objects are products of binary interaction occurring during their MS phase. We also find that the T_eff at which fast-rotating stars are no longer present varies with log(L/L_⊙), with cooler fast-rotating stars toward lower log(L/L_⊙), even beyond the TAMS. We argue that this could be a consequence of rotationally induced effects that change the width of the MS.

Based on the distribution of v sin i values with T_eff, we argue that the evolution and the larger fraction of the slowly rotating component may be an indicator that massive stars enter the MS with low to mild initial rotations (v_ini/v_crit ≲ 0.2). This would be further supported under the assumption that the fast-rotating component mainly emerges from binary interaction.

Using multi-epoch spectra, the associated peak-to-peak RV variations, and an adequate diagnostic of the intrinsic variability of O-stars and BSGs, we were able to find 176 SB1 systems among 494 stars located within 2500 pc. We find an average of 39% of SB1 systems within the MS, which is significantly reduced to 15% beyond the TAMS, further delineating its position in the HR diagram. We found that the behavior of SB2+ systems with T_eff is qualitatively similar to that of SB1 systems.

We compared the location of the theoretical TAMS in the reference sets of evolutionary models of Brott et al. (2011), Ekström et al. (2012), and Choi et al. (2016) with its empirical location obtained in this work. Our results suggest a revision of the overshooting parameter for non-rotating models, with higher overshooting values in the first two sets of models and a lower value in the latter set. We show that the behavior of the FR-TAMS is closer to that displayed by the theoretical TAMS in the models of Ekström et al. (2012).

The results presented in this work provide strong empirical constraints for improving stellar evolution models. Particularly with respect to the width of the MS in the massive star domain. Although our results also agree with the possibility that fast-rotating stars may have binary origin, further investigation is required to confirm this hypothesis. Furthermore, stars evolving blueward may also populate the upper part of the HR diagram near the TAMS. Determining the surface abundances and analyzing the pulsation properties of the stars in our sample will help to find new clues regarding their evolutionary nature and will be addressed in a future work. Moreover, the fraction of spectroscopic binaries within our sample demands further investigation and will likely add important constraints on the frequency of binary interaction in massive star evolution.

Acknowledgments

The authors thank the anonymous reviewer for a constructive and positive report that has led to the improvement of the manuscript. The authors also thank Dominic Bowman for providing valuable comments on this manuscript. AdB, SS-D, and GH acknowledge support from the Spanish Ministry of Science and Innovation and Universities (MICIU) through the Spanish State Research Agency (AEI) through grants PID2021-122397NB-C21, PID2022-136640NB-C22, 10.13039/501100011033, and the Severo Ochoa Program 2020-2023 (CEX2019-000920-S). Regarding the observing facilities, this research is based on observations made with the Mercator Telescope, operated by the Flemish Community at the Observatorio del Roque de los Muchachos (La Palma, Spain), of the Instituto de Astrofísica de Canarias. In particular, 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. This work is also based on observations with the Nordic Optical Telescope, owned in collaboration by the University of Turku and Aarhus University, and operated jointly by Aarhus University, the University of Turku and the University of Oslo, representing Denmark, Finland and Norway, the University of Iceland and Stockholm University, at the Observatorio del Roque de los Muchachos, of the Instituto de Astrofísica de Canarias. Additionally this work is based on observations obtained with the FEROS spectrograph attached to the 2.2 m MPG/ESO telescope at the La Silla observatory (Chile). This work has used 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.

References

Appendix A: Further notes on the completeness of the sample

In Sect. 3 we evaluated the completeness level of the stars analyzed with respect to the ALS III catalog. Figure A.1 extends this by showing the observed and missing stars within 2500 pc, separated by their spectral type. We do this to further evaluate the potential existence of biases in the distribution of the stars shown in this work that may lead to misinterpretations of the results.

For most of the spectral types, the percentage of observed stars is higher than the missing ones. Only the B5 spectral type has more missing objects. The B0-type gathers more missing stars than any other group. The associated mean T_eff of the observed objects that corresponds to this spectral type is 28 kK. Interestingly, this corresponds to a temperature close to the TAMS, on the hot side. Having these stars would probably strengthen the position of the TAMS with respect to what already shown.

The dashed pattern included in the bins of Fig. A.1 indicates stars for which we have available spectra, but we were unable to provide reliable estimates of the parameters. This fact has not been taken into account in the completeness analysis of Sect. 3 as we only aim to evaluate the level of completeness of the analyzed stars regardless of whether we have observed them or not. However, we can see that for many of the spectral types, they represent an important fraction of the missing objects. These cases usually correspond to peculiar objects or with problematic features in the spectra (e.g., strong nebular emission) or even stars simply not yet analyzed. Based on these results, we conclude that the lack of observed or analyzed objects should not affect the interpretations made in this work.

Fig. A.1. Histogram of the number of stars within 2500 pc located above the top reddening line of Fig. 1 in bins of spectral type. As in that figure, we use green and gray colors to separate analyzed stars in our sample from missing stars listed in the ALS III catalog (see Sect. 3 for further details). The dashed pattern within the missing sources corresponds to observed stars that have not been analyzed.

Appendix B: Evidence of disks near critical velocity

Oe and Be stars that exhibit disks are not present in our sample because of the impossibility of analyzing these objects with conventional 1-D atmospheric models (see Sect. 2). Despite this, in de Burgos et al. (2024a) some objects exhibiting a small double sub-peak structure in the Hα line were analyzed, suggesting the presence of some material around them. Looking at Fig. 15 of de Burgos et al. (2024a), we observed that the location of these objects in the sHR diagram shifts with T_eff and log(ℒ/ℒ_⊙) in a similar way as does the drop of fast-rotating stars in the HR diagram⁷. To further explore this connection, we marked with lime crosses in Fig. 4 those stars flagged in de Burgos et al. (2024a) for having such sub-peak structure. Interestingly, we can see that most of them are located close to their v_crit. As explained in de Mink et al. (2013), a consequence of stars reaching near-critical rotation is the formation of a near-Keplerian disk, which helps stars to decrease their angular momentum. We suggest that our results indicate clear evidence of this phenomenon.

Appendix C: The minimum observed spin rates increasing with luminosity

The associated 10%-percentiles of the top three panels of Fig. 4, representative of the low-v sin i component (Sect. 4.2), clearly show an almost flat trend, with a slow decay of ≈20 km s⁻¹ from the hot- to the cool-end of temperatures in the top two panels. Although the connection between the minimum measured v sin i values with the spectroscopic luminosity (Holgado et al. 2022; de Burgos et al. 2024a) or the evolutionary mass (Markova et al. 2014) has already been shown, our results display for the first time this characteristic with respect to the luminosity. For completeness, the average values of the 10%-percentiles drop by ≈10 km s⁻¹ from the top to the third panel. However, as pointed out in Sundqvist et al. (2013) (see also Simón-Díaz & Herrero 2014), we may be overestimating our measured values in the low-velocity regime as a consequence of additional turbulent motions present in the star’s atmospheres (namely micro- and macroturbulence).

Appendix D: Peak-to-peak detection thresholds for SB1 systems

Figure D.1 shows the same stars with three or more available multi-epoch spectra and SB1 systems shown in the HR diagram of Fig. 6, but this time in an sHR diagram. The figure includes two subpanels showing the RV_pp – T_eff and log(ℒ/ℒ_⊙) – RV_pp detection thresholds (gray areas) used in Simón-Díaz et al. (2024) to separate intrinsic stellar variability from SB1 systems (see Sect. 5.2). As explained there, the latter corresponds to those cases in which the RV_pp value is above any of the two threshold limits.

Fig. D.1. spectroscopic Hertzsprung-Russell diagram of the stars in the sample within 2500 pc with three or more available multi-epoch spectra. Green circles correspond to systems identified as SB1. The bottom and right subpanels show RVpp against Teff and log(ℒ/ℒ⊙), respectively. The gray areas highlight the threshold limits for detection of SB1 systems (see Appendix D). The evolutionary tracks correspond to nonrotating Geneva models. The average Teff of our proposed TAMS is shown with a dashed purple line.

Appendix E: Spin rate properties of stars with log(L/L_⊙) = 4.65 – 4.3

The bottom panel of Fig. 4 has some characteristics that are different from the other panels. It is the only panel in which the median values are above the threshold limit of 100 km s⁻¹, resulting from the smaller fraction of stars in the low-v sin i component (see the corresponding panels of Fig. 3). Although this characteristic is not new (see, e.g., Hunter et al. 2008; Holgado et al. 2022), the higher number of fast-rotating stars allows us to better see the “diagonal" decrease of v sin i in the 25 – 17 kK range. This panel also shows that most of the stars beyond the TAMS are fast-rotating objects. However, we note the lack of slowly rotating stars in the 21 – 15 kK range, for which we do not find a specific reason, but it could be related to the selection of stars with spectroscopic luminosities with log(ℒ/ℒ_⊙) > 3.5 dex in de Burgos et al. (2024a). The inclusion of additional stars in that T_eff range would shift the location of the TAMS toward cooler temperatures, resulting in a lesser number of fast-rotating objects beyond the TAMS.

All Tables

All Figures

Fig. 1. Gaia CMD including the stars in the sample and other missing stars included in the ALS III catalog (Pantaleoni González et al., in prep.). All the stars have Bmag < 11 and distances below 4 kpc. The two diagonal dashed black lines mark the reddening line of a 20 M⊙ star at the zero-age MS (ZAMS, bottom) and at 30 kK (top). Both are extended up to a Δ(AV) = 2. Stars shown in green include those for which both the use of the maximum and minimum distances result in the star above the top reddening line. Those for which minimum distance locate them below the same reddening line are indicated in blue. Stars missing from the ALS III catalog are shown in gray. Four evolutionary tracks corresponding to 9, 15, 20, and 40 M⊙ are included for reference. They were downloaded from the MESA Isochrones & Stellar Tracks online tool (MIST; see Dotter 2016; Choi et al. 2016; Paxton et al. 2013, 2015, for references purposes) for solar metallicity, no initial rotation, and AV = 0.0. The average uncertainty in the y-axis is indicated in the legend, while it is negligible in the x-axis.

In the text

Fig. 2. Hertzsprung-Russell diagrams showing the stars in our sample. From left to right, each panel limits the sample to stars within 1500 (221 stars), 2500 (670 stars), and 4000 pc (852 stars) of distance, respectively. The stars with v sin i > 100 km s−1 are indicated with cyan circles, whereas the rest are shown in orange. The average uncertainties are indicated with error bars in the lower-left corner of the panels. Each panel displays four luminosity ranges included for easier interpretation of the other figures. Our proposed location for the TAMS is marked with a purple line in all panels (see Fig. 3 and Sect. 4.1). The evolutionary tracks are taken from the MIST web-tool for solar metallicity and no initial rotation. The dashed green line on the right panel corresponds to a second-order polynomial fit of the points in the tracks where vcrit reaches 220 km s−1. A contour density mesh based on a Gaussian kernel is included to represent the density of points. Fig. 3. Number of stars with respect to their effective temperature. The histograms in the top four rows of subpanels separate stars within the indicated range of luminosities using different colors. Cyan bins indicate fast-rotating stars with v sin i > 100 km s−1. Each histogram displays the total number of stars and of fast-rotating objects (FR). Each column separate stars by their distances, as indicated in the top panel. The adopted position of the TAMS for each luminosity range is included with a vertical dashed line (see Sect. 4.1). The bottom row of subpanels displays the corresponding cumulative distribution function following the color code of the histograms.

In the text

Fig. 3. Number of stars with respect to their effective temperature. The histograms in the top four rows of subpanels separate stars within the indicated range of luminosities using different colors. Cyan bins indicate fast-rotating stars with v sin i > 100 km s−1. Each histogram displays the total number of stars and of fast-rotating objects (FR). Each column separate stars by their distances, as indicated in the top panel. The adopted position of the TAMS for each luminosity range is included with a vertical dashed line (see Sect. 4.1). The bottom row of subpanels displays the corresponding cumulative distribution function following the color code of the histograms.

In the text

Fig. 4. Projected rotational velocities of the stars in our sample against their effective temperature. In this case, the sample is not limited to any distance. The panels separate stars into different luminosity ranges as in Fig. 3, following the same color code. Solid black error bars and central circles indicate the 75% and 25% percentiles and the median values of the stars in bins of 4 kK. The dashed gray error bars correspond to the 90% and 10% percentiles. The dotted line at 100 km s−1 separates fast-rotating stars. The location of the TAMS is also included with a vertical dotted line. The solid green lines correspond to the critical velocity calculated from the models in Fig. 2. Stars showing double sub-peak emission in Hα are marked with additional lime-colored crosses (see Appendix B). The gray area corresponds to the temperature range covered by the theoretical ZAMS in each luminosity range.

In the text

Fig. 5. Hertzsprung-Russell diagrams showing different evolutionary tracks for solar metallicity together with different predictions for the location of the TAMS as indicated in each legend. The top panel includes tracks with no initial rotation, and in the bottom panel the tracks correspond to vini/vcrit ≈ 0.4. Intervals of the same age difference are marked with different symbols. The gray density mesh in the background is based on a Gaussian kernel and includes our sample stars up to 2500 pc. Evolutionary tracks are limited up to the RSg phase (if applicable).

In the text

	Fig. 6. Similar HR diagram as in Fig. 2 but limiting the sample to stars within 2500 pc with three or more available multi-epoch spectra. Green circles highlight those stars identified as SB1.
In the text

Fig. 7. Histograms of the number of spectroscopic binaries against spectral type or log10(Teff), limited to a distance of 2500 pc. The top panel include SB2+ systems for which the spectral type is taken from the primary star (see also Sect. 5.2 for further details). The middle panel is similar to the top one but including the SB1 systems detected from the RVpp variability. The bottom panel include the same stars as in the middle one, but we used the log(Teff) instead of the spectral type. In the top two panels, the position of the TAMS is highlighted with a patch of purple lines, derived from the range of spectral types of stars within 1 kK of the average Teff of the TAMS, which is shown in the bottom panel as a dashed purple line.

In the text

Fig. A.1. Histogram of the number of stars within 2500 pc located above the top reddening line of Fig. 1 in bins of spectral type. As in that figure, we use green and gray colors to separate analyzed stars in our sample from missing stars listed in the ALS III catalog (see Sect. 3 for further details). The dashed pattern within the missing sources corresponds to observed stars that have not been analyzed.

In the text

Fig. D.1. spectroscopic Hertzsprung-Russell diagram of the stars in the sample within 2500 pc with three or more available multi-epoch spectra. Green circles correspond to systems identified as SB1. The bottom and right subpanels show RVpp against Teff and log(ℒ/ℒ⊙), respectively. The gray areas highlight the threshold limits for detection of SB1 systems (see Appendix D). The evolutionary tracks correspond to nonrotating Geneva models. The average Teff of our proposed TAMS is shown with a dashed purple line.

In the text