© G. Holgado et al. 2022

1 Introduction

Stellar rotation affects the global physical characteristics and the evolution of massive stars in many different ways (Maeder & Meynet 2000; Langer 2012). As a massive star evolves, it is expected that its surface rotational velocity will be modified due to different processes. Generally speaking, in single stars the equatorial velocity always decreases monotonically, and the net braking rate of the stellar surface mostly depends on the combined effect of angular momentum losses because of stellar winds and the efficiency of internal transport of angular momentum from the core to the envelope. In binary systems – especially in cases where the two components are close enough to interact at some point in their lives – the situation becomes a bit more complex. There are several additional effects that do not operate in single stars (such as tides and mass transfer events) that are capable of both decelerating or accelerating each of the individual components (de Mink et al. 2013, and references therein).

In order to advance our understanding of the global rotational properties of massive stars – as well as the impact of this important parameter on their internal structure and evolution –, numerous studies have obtained and analyzed spectroscopic data sets of difference size and quality (Conti & Ebbets 1977; Penny 1996; Howarth et al. 1997; Huang & Gies 2006a,b; Dufton et al. 2006; Daflon et al. 2007; Hunter et al. 2008; Penny & Gies 2009; Fraser et al. 2010; Bragança et al. 2012; Simón-Díaz & Herrero 2014; Ramírez-Agudelo et al. 2013, 2015). Some of them have also compared the gathered empirical information with the general predictions of evolutionary models and population synthesis computations that include rotation (Simón-Díaz et al. 2010; Sundqvist et al. 2013; Granada & Haemmerlé 2014; Martins et al. 2017; Markova et al. 2014, 2018; Keszthelyi et al. 2017, 2020). Others have used the obtained distributions of projected rotational velocities (v sin i) to establish empirical constraints to our theories of massive star formation (Wolff et al. 2006; Mokiem et al. 2006; Huang et al. 2010; Rosen et al. 2012).

Following (Holgado et al. 2018, 2020), this paper is a continuation of a series devoted to the study of a sample of more than 400 Galactic O-type stars, for which the IACOB and OWN surveys (last described in Simón-Díaz et al. 2020b; Barbá et al. 2017) have been gathering high-quality, multi-epoch spectroscopic observations over the last few decades. In this paper, we focus on presenting a re-evaluation of our empirical knowledge about the global spin-rate properties of the Galactic main sequence stars with masses in the range ~15–80 M_⊙.

To this aim, we used the outcome of a homogeneous (and objective) quantitative spectroscopic analysis, performed with a battery of semi-automated tools based on state-of-the-art techniques (Simón-Díaz et al. 2011; Simón-Díaz & Herrero 2014), and the stellar atmosphere code fastwind (Santolaya-Rey et al. 1997; Puls et al. 2005), as well as information about the spectro-scopic binary status of the complete sample of stars obtained from the available multi-epoch data set (Holgado 2019). In addition, we benefited from the theoretical predictions provided by two different single star evolutionary codes (Brott et al. 2011; Ekstrõm et al. 2012), and the results of the binary population synthesis simulation presented by de Mink et al. (2013).

This paper is intended to broaden the scope of the works by Simón-Díaz & Herrero (2014) and Markova et al. (2014), also improving and superseding the information presented in some important reference studies such as Conti & Ebbets (1977), Penny (1996), and Howarth et al. (1997). In addition, it can be considered as an extension of the investigation of the rotational velocities of a sample of presumably single and spectroscopic binary O-type stars in the 30 Doradus region of the Large Magellanic Cloud (LMC), which was performed by Ramírez-Agudelo et al. (2013, 2015) in the framework of the VLT-FLAMES Tarantula Survey (VFTS, Evans et al. 2011).

The structure of this paper is as follows. The sample, observations, and spectroscopic analysis tools used are briefly described in Sect. 2. Results, in Sect. 3, constitute the core of the paper, where we discuss the global properties of the sample and the interdependence between v sin i and evolution. Section 4 evaluates the empirical results, taking into account the trends that evolutionary models produce, with Sect. 5 focusing on the initial rotation distribution in O-type stars. Concluding remarks and future prospects are included in Sect. 6.

2 Sample definition and methodology

The final sample of targets considered for this work was drawn from an initial sample of 415 Galactic O-type stars, for which the IACOB and OWN projects have presently available ~3900 high-resolution (R = 25000-85 000), multi-epoch spectroscopic observations. To avoid too much repetition, we refer the reader to Holgado et al. (2020) for a detailed description of how this sample was built. There, we also describe its main characteristics in terms of spectral type and luminosity class coverage, as well as its completeness when compared with the list of O-type stars quoted in version 4.1 of the Galactic O-star catalog (GOSC, Maíz Apellániz et al. 2013).

From this initial sample, we excluded all those stars clearly identified as double-line spectroscopic binaries (SB2s, 113 objects) or Wolf-Rayet (4), as well as all the targets presenting features in their spectra that are typically associated with Oe (6), or magnetic (7) stars, for which a successful fitting with our fastwind set of models was not attainable¹. The main reason supporting this decision was this inherent limitation of our analysis strategy (see below) to provide reliable stellar parameters for these stars.

We also benefited from the multi-epoch character of the IACOB and OWN surveys to flag the stars in the remaining sample that could be clearly identified as single-line spectroscopic binaries (SB1s). To this aim, following the strategy presented in Holgado et al. (2018), we obtained radial velocity (v_rad) estimates for all available spectra. Then, after computing the associated dispersion – σ(v_rad) – and the peak-to-peak amplitude – ∆v_rad – of all measurements, we identified as clear SB1s the stars that fulfilled the criteria described in Appendix A (see also Holgado 2019).

As a result, we ended up with a final sample of 285 stars, for which we were able to obtain their projected rotational velocities (v sin i) and the spectroscopic parameters. Among them, 55 stars were flagged as SB1s. The other 230 stars are referred to as likely single (LS), although some of them might still be undetected spectroscopic binaries and/or SB1 stars with low amplitude v_rad variations. In this regard, we remark that, by the time of finalizing this paper, we had only one epoch available for 54 out of the 230 stars in the LS subsample.

Details about the tools we used and the methodology we followed to perform the quantitative spectroscopic analysis of our working sample of 285 O-type stars can be found in Holgado et al. (2018, 2020) and references therein. In brief, our analysis strategy, applied to the best signal-to-noise ratio (S/R) spectrum per star, can be summarized in two main steps.

First, we obtained estimates for the two line-broadening parameters (v sin i and v_mac) using the iacob-broad tool (Simón-Díaz & Herrero 2014) and following the guidelines presented in Holgado et al. (2018). In particular, we based our v sin i determination on the O iii λ5592 line whenever possible (see, however, notes in Sect. 3), and we assumed a radial-tangential definition of the macroturbulent broadening profile (as suggested by Simón-Díaz & Herrero 2014). In addition, we benefited from the possibility offered by iacob-broad to compare the v sin i values resulting from the Fourier Transform (FT) and goodness-of-fit (GOF) analysis, as an assessment of the reliability of the results regarding this quantity.

Next, we performed a HHe analysis² using the semi-automated tool iacob-gbat (Simón-Díaz et al. 2011; Holgado et al. 2018), which has incorporated a grid of synthetic spectra computed with the fastwind (Santolaya-Rey et al. 1997; Puls et al. 2005; Rivero González et al. 2012) stellar atmosphere code. As a result, we obtained estimates and associated uncertainties for the stellar effective temperature (T_eff), surface gravity (log g), helium surface abundance (Y_He), microturbulence (ξ_t), and wind-strength Q-parameter. The complete results of this analysis were presented in Holgado (2019), and a subset of the derived parameters can be found in Holgado et al. (2020, see also Table F.1).

In this paper, we mainly make use of three of the parameters determined spectroscopically: v sin i, T_eff, and log g³. Indeed, instead of log g, we use the quantity $L:=T_{\text{eff}}^4/g$, which allows us to locate the stars in the so-called spectroscopic Hertzsprung–Rusell diagram (sHRD, Langer & Kudritzki 2014).

Similarly to the case of Holgado et al. (2020), we preferred to use this diagram instead of the traditional Hertzsprung–Rusell diagram (HRD), given that 45 (out of the 285) stars in the sample present a value of renormalized unit weight error (RUWE) > 1.4 in the Gaia EDR3 release (Lindegren et al. 2018; Gaia Collaboration 2021), and hence the reliability of available parallaxes for these stars still needs to be put under quarantine. As discussed in Langer & Kudritzki (2014), the sHRD can be considered as an analogous version of the HRD, but with the benefit of allowing the comparison of observations and evolutionary models regardless of distance and extinction constraints, since it is purely based on parameters determined spectroscopically.

Comparison with previous studies: Sample size and available empirical information

As indicated in Sect. 1, one of the initial drivers of the study presented in this paper was to provide an updated and improved empirical overview of the spin-rate properties of Galactic O-type stars. Table 1 allows us to evaluate the main improvements that we have achieved with respect to several commonly referenced studies found in the literature. We cover four different aspects: (a) the size of the investigated sample of stars; (b) whether the study uses v sin i estimates decontaminated from the effect of the so-called macroturbulent broadening; (c) the availability of T_eff and log g estimates; and (d) whether the study makes any specific study comprising the spectroscopic binaries that are present in the sample under study.

As can be concluded from examination of Table 1, the empirical information provided in (and discussed throughout) this paper implies a clear improvement – in at least two of the aspects mentioned above – with respect to all quoted studies about rotational velocities in Galactic O-type stars. More specifically, our working sample of stars is not only ~20–40%o larger⁴ than those considered by Conti & Ebbets (1977), Penny (1996), or Howarth et al. (1997), but we have also been able to: (a) supersede their derived v sin i with new (more robust) estimates of this quantity; (b) better identify spectroscopic binaries; and (c) routiney incorporate information about T_eff and log g in all stars in the sample that are not identified as SB2s.

Regarding the more recent works by Simón-Díaz & Herrero (2014) and Markova et al. (2014), our paper can be seen as a natural extension of these two works. Here, we investigate a considerably larger sample of stars, with a better coverage in terms of spectral type and luminosity class⁵, and with freshly incorporated information about the spectroscopic binary status for a large fraction of them.

3 Results

3.1 General overview

Table F.1 summarizes the compiled empirical information of interest for our working sample of 285 Galactic O-type stars⁶. We repeat the same data provided in Holgado et al. (2020) regarding spectral classification (as quoted in GOSC), spectroscopic parameters – T_eff, log g, log (L/L_⊙) –, and whether the star was identified as line profile variable (LPV) or SB1, using the set of available multi-epoch spectra. In addition, we complement this information with the new listed derived line-broadening parameters v sin i and the v_mac used for this work⁷, also indicating the diagnostic line used in the IACOB-BROAD analysis. Similarly to the table presented in Holgado et al. (2020), stars in Table F.1 are grouped by luminosity class (LC) and ordered by spectral type (SpT).

Figures 1 to 3 provide a first general overview of the location of the sample of investigated stars in the sHRD and the associated v sin i distribution. We refer the reader to Holgado et al. (2020) for a thorough discussion about the general distribution of the sample in the sHRD, including possible explanations for the scarcity of stars close to the zero-age main sequence (ZAMS) in the mass range ~30–80 M_⊙ (marked in gray in Fig. 1). Also, while we do not consider it necessary to provide additional information about the reached level of accuracy and reliability of the measured effective temperatures and surface gravities (since this is already properly discussed in Holgado et al. 2018), we indicate below some information of interest regarding our v sin i measurements. In addition, a discussion of some potential caveats and limitations that could still affect the reliability of the estimates of this quantity in the low- v sin i regime can be found in Sect. 3.2.3.

As indicated in Sect. 2, for the sake of homogeneity, we tried to mainly base our v sin i determinations on the O III λ5592 line. This was possible for almost 90% of the stars in the sample (see top panel in Fig. 2). There were, however, a few late and early O-type stars for which the O III line was weak and we had to rely on alternative lines. In most ot those cases, we were able to still use a metal line (either Si iii 24552, Niv 26380, or N v λ4603/20⁸), allowing us to reach 95% of the stars. But, for some very fast rotators (37 stars, 5% of the total sample) all metal lines appeared too diluted to provide reliable results, and we had to rely on Hel lines. Although the use of a He i line may present caveats in the case of slow rotators (especially in analyses based on medium and low-resolution spectra, or when an important nebular contamination affects the He i lines, see, e.g., Ramírez-Agudelo et al. 2013), we note that almost all the stars for which we needed to use a He i line have a v sin i value above 200 km s⁻¹. As a result, in these stars the effect of Stark broadening in the determination of v sin i can be considered as negligible. Also, none of the analyzed He i lines presented any nebular contamination.

The bottom panel of Fig. 2 shows the good agreement found between the two analysis strategies provided by iacob-broad, with the comparison of v sin i(FT) and v sin i(GOF) estimates. Differences are smaller than 10 km s⁻¹ in ~85% of the sample and, except for a few stars, are not larger than 15 km s⁻¹ for the rest. This result can be used as an assessment of the reliability of the derived v sin i values, since they were obtained using two complementary analysis strategies (see, however, Sect. 3.2.3).

In view of this result, we decided to use v sin i (GOF) estimates and to establish 5–15 km s⁻¹ as the typical (formal) uncertainty associated with the measurement of this quantity from our high-quality spectroscopic data set. This is also the reason why we use a bin size of 20 km s⁻¹ in the v sin i histograms presented hereafter throughout the paper. See Appendix E for a review on the significance of this selection.

Fig. 1 Spectroscopic HR diagram indicating the location of the 285 Galactic O-type stars investigated in this work. Different colors and symbols are used to highlight luminosity classes. Individual uncertainties are included as error bars. Evolutionary tracks and position of the ZAMS from the non-rotating, solar-metallicity models by Ekström et al. (2012), plus the log g = 3.7 dex line of constant gravity, are also included for reference purposes. The gray shadowed area highlights a region of the diagram which is void of stars (see main text and Holgado et al. 2020, for further explanations).

Fig. 2 Rotational properties of the sample. Top: Histogram of υ sin i values resulting from the iacob-broad analysis of the 285 O-type stars in our working sample, highlighting the diagnostic line considered for the analysis. Bottom: Comparison of the v sin i estimates resulting from the FT and GOF analysis strategies. The diagonal dotted lines indicate differences between both quantities of ±10 and ±20 km s−1.

Fig. 3 Histogram (top) and cumulative distribution (bottom) of υ sin i estimates for the 285 O-type stars investigated in this work, separating results for the various luminosity class subgroups (same color code as in Fig. 1). Ranges in v sin i labeled with letters “A” to “F” are used in Fig. 5 and explained in Sect. 3.3.

3.2 The empirical v sin i distribution of Galactic O-type stars

3.2.1 Properties of the global v sin i distribution

Figure 3 depicts the information about the v sin i distribution of our sample of 285 Galactic O-type stars in the form of a histogram (top) and a cumulative distribution function (bottom). In addition to the global distribution, we also provide the resulting distributions separated by luminosity class groups. This helps us to establish a first link between our results and those presented in previous investigations about the rotational distribution of O-type stars.

Inspection of these two panels allows us to point out two main features already highlighted in the seminal work by Conti & Ebbets (1977), and later on reproduced by other similar studies, such as those by Penny (1996) and Howarth et al. (1997) for the case of Galactic stars, or Penny & Gies (2009) and Ramírez – Agudelo et al. (2013, 2015) for the case of stars in the Magellanic Clouds. On the one hand, the global distribution shows a clear bimodal character, with a slow velocity peak comprising most of the stars (in our case, ~75% of the total sample, if we consider that this part of the distribution ends in υ sin i ~ 150 km s⁻¹), and a tail of fast rotators, reaching υ sin i values up to ~450 km s⁻¹. On the other hand, the ranges in v sin i covered by the various luminosity class groups are different, with the two most extreme cases found between the dwarfs (LC V) and the supergiants (LC I). While the former are distributed over the whole range, the latter are mostly concentrated within the low-velocity peak feature of the distribution. This feature is well-documented also for other spectral classes (e.g., Royer et al. 2002; Dufton et al. 2013, for A- and B-type stars, respectively).

Several studies have proposed that for less massive stars, the origin of this bimodality may lie in their ability to maintain, or not, the presence of a disk of material around themselves, producing an interaction between the stars and their accretion disks. This interaction is able to remove angular momentum and prevent the spin up, and is normally simplified as the concept of “disk locking” (see Bastian et al. 2020, and references therein). In the cases of stars born with masses in the range covered by the O-type stars, the bimodal character of the υ sin i distribution (mainly present in the dwarf sample, and initially speculated to reflect the initial spin distribution, see Fig. 1) could not be convincingly explained in the early works by Conti & Ebbets (1977). However, more recent studies have proposed that this feature is a natural consequence of the effect of binary interaction during massive star evolution (de Mink et al. 2013, see also Sect. 4.2). Therefore, taking into account the large percentage of O-type stars that are expected to evolve as part of a binary system (Sana et al. 2012), the incorporation of information about the spectroscopic binary status of individual targets in any global empirical study of the rotational properties of massive stars has become of crucial interest.

3.2.2 Likely single stars versus single-line spectroscopic binaries

To the best of our knowledge, there is only one comprehensive study that performs a detailed comparison of the rotational properties of a statistically meaningful sample of single and spectroscopic binary stars in the O-star domain. We refer to the work by Ramírez-Agudelo et al. (2015) in the 30 Doradus region of the LMC. No similar study has been carried out yet for Galactic O-type stars.

Thanks to the multi-epoch character of the IACOB and OWN surveys, we can take the first steps in this direction. As mentioned in Sect. 2, the compiled spectroscopic data set allowed us to identify about 170 spectroscopic binaries (including 113 SB2 systems and 55 clearly detected SB1 stars) within our initial sample of 415 Galactic O-type stars. The rest of the stars in the sample (230, excluding 17 peculiar stars) have been flagged as likely single. While the SB2 systems will not be further considered for this work (see Sect. 2), throughout this paper we also investigate similarities and differences that can be found between ourLS and SB1 subsamples.

In view of this, Fig. 4 depicts the v sin i distributions resulting from these two subsets of stars separately. As illustrated by the lower panel of Fig. 4, the distribution of likely single stars is very similar to the combined one, meaning that they clearly dominate the global distribution (as expected from the relative percentage of stars in each subsample). In contrast, the comparison of the LS and SB1 samples highlights clear differences. In fact, a Kuiper’s test (Kuiper 1960) indicates that there is a very low probability (below 1%) that the two observed distributions are randomly drawn from the same parent population.

More specifically, the bimodal character of the global v sin i distribution is not so clearly defined in the case of the SB1 sample. Indeed, one could assume that there is only one component that, compared to the main (low v sin i) component of the LS star distribution, is wider (by a factor of ~2) and it is centered at a somewhat larger v sin i (120 vs. 65 km s⁻¹). In addition, there is a lack of SB1 stars with projected rotational velocities larger than 300 km s⁻¹.

Interestingly, similar results were obtained by Ramírez-Agudelo et al. (2015) when they compared the observed distributions of the presumed single stars and the primary components of those spectroscopic binaries comprising the O-type star population of the 30 Doradus region of the LMC. Following de Mink et al. (2013), these authors propose that the observed differences are most likely produced by a combination of effects taking place along the various stages of early evolution of high-mass binaries (such as tides, mass and/or angular momentum transfer, and magnetic braking; see also additional notes in Sect. 4.2). For a more detailed comparison of our resulting v sin i distribution with the one obtained by Ramírez-Agudelo et al. (2013), we refer the reader to Appendix C.

Fig. 4 Comparison of the υ sin i distributions of the likely single (LS) and clearly detected SB1 subsamples, presented as histograms (top and middle panels, respectively) and cumulative distributions (bottom panel). Stars in the LS sample, for which we only have one spectrum, are highlighted in gray. Ranges in υ sin i labeled with letters “A” to “F” are used in Fig. 5 and explained in Sect. 3.3.

Fig. 5 Location in the sHRD of our working sample of stars as a function of the various v sin i ranges indicated in Figs. 3 and 4 (see also Sect. 3.3). The LS and SB 1 subsamples are highlighted in black and red, respectively, and their relative percentage in each panel is included for reference. Further information about other features presented in this figure can be found in the caption of Fig. 1.

3.2.3 The low- v sin i regime

In previous sections, we have indicated that our revised v sin i estimates can be considered as more robust and reliable than those provided by previous works (e.g., Conti & Ebbets 1977; Howarth et al. 1997; Penny 1996). This is mainly a result of the combined use of a high-quality spectroscopic data set (in terms of spectral resolution and S/R) and a line-broadening analysis technique, allowing us to disentangle the effect of the rotational and the so-called macroturbulent broadening components. However, following Sundqvist et al. (2013), Simón-Díaz & Herrero (2014), and Markova et al. (2014), we are aware that our analysis methodology might still present some limitations in the low-v sin i regime.

This potential failing of the considered line-broadening analysis technique clearly stands out from inspection of the v sin i distribution presented in Fig. 3: there is a notorious deficiency of O stars with spectroscopically inferred v sin i below ~40–50 km s⁻¹ (particularly in the case of stars with LCs I, II and III). Taking into account that our instrumental set-up imposes a lower limit in the measured v sin i of 5–10 km s⁻¹, plus the statistical effect of the inclination angle (i), one would expect a larger number of stars in the low- v sin i region of the distribution.

As extensively discussed in Simón-Díaz & Herrero (2014) and references therein, this long-standing problem – which was already present (to a larger extent) in the works by Conti & Ebbets (1977); Howarth et al. (1997), or Penny (1996) – has been only partially alleviated by accounting for the effect of macroturbulence on the determination of projected rotational velocities. Other sources of line-broadening seem to still be hindering the suitability of available methods to reach actual v sin i measurements below ~40–50 km s⁻¹ in some specific situations.

While there is not yet a definitive answer as to what is the main culprit of this methodological limitation, some empirical hints seem to point toward so-called microturbulent broadening. As pointed out by Gray (1973), classical atmosphere microturbulence also produces zeroes in the FT at (high) frequencies associated with low values of v sin i that may be falsely identified as the zeroes associated with the rotational broadening below a certain v sin i limit. In this context, Simón-Díaz & Herrero (2014) showed that microturbulence values in the range 5–25 km s⁻¹ can lead to Fourier space minima at frequencies corresponding to v sin i ~ 10–40 km s⁻¹. The good correlation between the Fourier and GOF methods shown in Fig. 2 indicates that the latter is probably also affected by this unaccounted broadening.

Some additional notes about this line of research can be found in Appendix D, where we also summarize the main results obtained in the context of this paper, which could help to shed light on this problem. Hereafter, throughout the paper, we take this warning into account for any interpretation of results. In particular, we consider (as a first order approach) that any v sin i estimate between ~40–50 km s⁻¹, could actually correspond to an upper limit of this quantity, especially in the upper part of the sHRD, or in the case of stars for which a high value of v_mac has been obtained (following the results presented in Fig. 4 of Simón-Díaz et al. 2017).

3.3 The v sin i distribution of Galactic O-type stars in the sHRD

In most previous studies about the rotational properties of O-type stars, spectral types and luminosity classes were used as a proxy of mass and evolutionary status, since direct information about the stellar parameters of the investigated samples was not available. Having access to T_eff and log g estimates allows us to take one step further in the description and interpretation of results (see, e.g., Markova et al. 2014).

To serve as a first example, Fig. 5 shows (in six separated panels) how stars covering different ranges of the global v sin i distribution populate the sHRD. In each panel, we differentiate those stars identified as SB 1. The considered ranges (indicated in Figs. 3 and 4 with letters “A” to “F”), were specifically selected to separate several regions of interest within the two main components of the global (bimodal) v sin i distribution. In particular, regions “A”, “B”, and “C” (depicted in the top panels of Fig. 5) cover the low-velocity peak of the distribution, ranging up to 150 km s⁻¹ and including ~75 % of the stars in the working sample. Regions “D”, “E”, and “F” (bottom panels) correspond to the tail of fast rotators. Although the assumed boundaries for the various regions within each of the two main components could be initially seen as somewhat arbitrary, they were actually chosen to illustrate several features of interest that emerge when considering how the stars in each v sin i range are distributed in the sHRD.

The most evident cases correspond to panels A and F, which comprise the two extreme regions of the v sin i distribution. On the one hand, the upper v sin i boundary of region A has been specifically selected to show the lack of stars with measured v sin i below ~35 km s⁻¹ above the 32 M_⊙ track. We note that the value of 35 km s⁻¹ is used here only for that particular reason, as stars with low v sin i are described throughout the paper as those with v sin i< 40–50 km s⁻¹. On the other hand, we find it interesting to remark how the (19) stars with v sin i > 300 km s⁻¹ (panel F) are mostly concentrated in the lower-mass region of the O-star domain (except for two targets⁹, all of them are located below the 40 M_⊙ evolutionary track). Furthermore, no SB 1 is detected among this latter sample, even though we have a minimum of five epochs for most of them (see Table F.2).

A more careful inspection of the various panels in Fig. 5 allows us to highlight another group of stars which deserves further attention. We refer here to the bunch of (seven) stars closer to the upper part of the ZAMS (i.e., those with T_eff > 45 000 K and log(L/L_⊙) ~ 4.1 dex, see also Table F.3), which are almost exclusively found in panel B (except for the only SB 1 star in this group, located in panel C). This result was already pointed out by Markova et al. (2014), who found suggestive evidence that stars with initial masses ≳50 M_⊙ rotate with velocities not exceeding 26% of their critical rotational velocities when they are close to the ZAMS. We note, however, that this sample of seven stars is still not statistically large enough to provide a robust confirmation of the lack of fast-rotating stars in this region of the sHRD, especially taking into account that we are actually measuring v sin i and not the surface equatorial velocity.

Also, it can be noticed that the upper right region of the sHRD (corresponding to the more evolved O-Supergiants) is void of stars in panel E, and only a few stars in this region are found in panel D. This result is coherent with Fig. 3, where it is shown that the relative percentage of O stars with luminosity classes I and II drops drastically above v sin i~ 150 km s⁻¹, and becomes zero above ~250 km s⁻¹ (see also the specific location of the sample of LC I and II stars in the sHRD in Fig. 1).

Last, but not least, Fig. 5 offers us an interesting overview of the relative percentage (and distribution) of stars labeled as SB 1 in each of the considered v sin i ranges. In particular, the relative percentage of SB1 stars reaches its maximum (~30–35%) in panels C and D (v sin i = 80–200 km s⁻¹), and decreases significantly below and above these limits. This is just a consequence of the different shape of the v sin i distributions of the LS and SB1 samples (see Fig. 4).

The unexpected distribution of stars in panel A obviously reminds us that sometimes our presently available methods are still not suitable to provide reliable v sin i estimates in the low-v sin i regime (see notes in Sect. 3.2.3 and Appendix D). Other than that, the results summarized above are expected to provide important empirical constraints for any attempt to investigate, from a theoretical point of view, the impact of high-mass star formation and (single and binary) evolution on the spin-rate properties of main sequence O-type stars (see further notes in Sect. 4).

Fig. 6 v sin i–Teff diagram including our working sample of stars. The LS and SB1 subsamples are highlighted in black and red, respectively. The v sin i = 35, 150, and 300 km s−1 values are highlighted for reference purposes (see Sect. 3.4).

3.4 v sin i dependence with T_eff and log (L/L_⊙)

Figures 6 and 7 offer a complementary view of the information presented in Fig. 5 and described in the previous section. In particular, the type of diagram presented in Fig. 6 (either using spectral types or directly using T_eff in the x-axis¹⁰) has traditionally been used as a first approach to investigate the behavior of the (projected) rotational velocity of stars in the high-mass domain as a function of their evolution (see, e.g., Conti & Ebbets 1977; Fraser et al. 2010; Vink et al. 2010; Markova et al. 2014; Simón-Díaz & Herrero 2014; Keszthelyi et al. 2017).

Thanks to the increased statistics, combined with the direct spectroscopic determination of the effective temperatures, our sample offers a more complete and robust overview (compared to previous attempts) of how Galactic O-type stars are distributed in this diagram. For example, in the case of Vink et al. (2010), which used the long sample of O stars investigated by Howarth et al. (1997, see also Table 1), effective temperatures were obtained by converting spectral types, using SpT–T_eff calibrations. In addition, the considered v sin i estimates were still affected from so-called macroturbulent broadening. More recently, Markova et al. (2014) used a similar approach to the one utilized here to get estimates of v sin i and T_eff; however, their sample only included 39 O-type stars (see Table 1) and, in their own words, was importantly affected by some selection effects.

The v sin i-log L diagram (Fig. 7) has not been so commonly used in previous works, but we also include it here because of its utility for the correct interpretation of some specific features shaping the distribution of stars in the v sin i - T_eff diagram. In addition, it also provides interesting clues to aid in the endeavor to identify the origin of the still remaining methodological limitations hindering a reliable determination of v sin i in the low-velocity regime.

Indeed, Fig. 7 clearly shows that there exists a strong correlation between the lower limit in v sin i imposed by our analysis methodology and the quantity log L. As discussed in more detail in Appendix D (see also Sect. 3.2.3), this result indirectly supports the hypothesis that in low-v sin i stars, we might be misidentifying zeroes in the FT that are produced by the line-broadening effect of microturbulence, as if they were associated with rotational broadening.

Once we are aware of this correlation, we can better interpret the curious shape defining the lower boundary of the distribution of points in the v sin i-T_eff diagram (Fig. 6). As seen in Fig. 1, the range in Teff covered by stars with a given value of log £ becomes smaller as we move down in the sHRD. If we combine this characteristic distribution of our sample of O-type stars in the sHRD with the abovementioned correlation, we can understand why the lower v sin i limit in the v sin i-Teff diagram reaches a minimum at Teff ~ 35 000 K, and tends to increase toward both higher and lower effective temperatures. This is also the reason why stars in Panel A of Fig. 5 only populate the region below log(L/L_⊙)~3.8dex.

Beyond that, the most noticeable feature that stands out from inspection of Fig. 6 is the concentration of stars with T_eff ≳ 45 000 K around the v sin i range between ~50 and ~75 km s⁻¹ (except for one SB 1 star which has a somewhat larger v sin i). Interestingly, when combined with the information presented in Fig. 5 (see panels B and C), it becomes clear that this appendix in the v sin i-Teff distribution corresponds to the bunch of stars located closer to the ZAMS around the 85 M_⊙ track. A similar distribution appears in Berlanas et al. (2020) for a sample of Cygnus OB2 OB stars. This lack of fast rotators at the hot end of the distribution – assuming that it is not an inclination effect – together with the clear decrease in the upper v sin i boundary of the distribution of points in the v sin i–log L diagram (Fig. 7), could be interpreted as empirical evidence of the increasingly important braking effect produced by stellar winds, as more luminous stars are considered. However, as we will show in the next sections, this is not necessarily the correct interpretation, especially if we take into account, following de Mink et al. (2013), that a high fraction of the stars in the sample could actually be post-interacting binaries (even if they are currently detected as likely single stars). The work by Berlanas et al. (2020) also proposes a similar direction of investigation, with alternative evolutionary channels proposed to explain these objects.

Fig. 7 v sin i–log(L/L⊙) diagram including our working sample of stars. The LS and SB1 subsamples are highlighted in black and red, respetively. The v sin i = 35, 150, and 300 km s−1 values are highlighted for reference purposes (see Sect. 3.4).

3.5 v sin i dependence on mass and evolution

In order to better investigate how the distribution of projected rotational velocities in the O-star domain depends on mass and evolution, we decided to go beyond the traditional use of the Teff-v sin i diagram (Fig. 6) and split the global sample of studied stars into three subsamples according to their location in the sHRD. By doing this, we could still keep a good level of statistical significance in each of these three subgroups, with a minimum of 75 objects in the less populated group. There have been similar studies for very large samples of lower-mass stars (e.g., Zorec & Royer 2012, for A-type stars).

We first used the 32 M_⊙ non-rotating evolutionary track computed by Ekström et al. (2012) to roughly separate the sample in two mass ranges. We chose this boundary because stars with higher masses are expected to develop winds strong enough to produce significant losses of angular momentum, even during the early phases of evolution along the main sequence. Then, we further divided the high-mass sample in two, using the log g = 3.7 dex line of constant gravity, which allowed us to roughly separate the expected younger dwarf stars from the more evolved giants and supergiants (see Fig. 1).

The leftmost panels of Fig. 8 show the location of these three groups of stars in the sHRD, where we again highlight the stars identified as SB Is in red. These panels are complemented with the corresponding empirical v sin i distributions, separating once more the likely single and SB1 subsamples (middle and rightmost panels).

In each of these latter panels, we indicate three v sin i values with vertical lines. These correspond to the rough boundary (50 km s⁻¹) below which our empirical v sin i measurements can be affected by some type of methodological limitations (see Sect. 3.2.3), the approximated limit we established in the global v sin i distribution between the low-velocity component and the tail of fast rotators (150 km s⁻¹, see Sect. 3.2.1), and the v sin i value above which no SB1 stars are detected (300 km s⁻¹, see Sect. 3.2.2).

Hereafter, we denote the two high-mass subgroups as “a” and “b”, respectively, while group “c” corresponds to the sample of O stars below the 32 M_⊙ evolutionary track. In a single star evolution context, stars in group “a” can be regarded as progenitors of the stars in group “b”. Equivalently, the v sin i distribution of group “a” is expected to become the v sin i distribution of the group “b” by the effect of stellar evolution. Again in this context, group “c” could be considered as a separated group in terms of evolution. However, as illustrated by Ahumada & Lapasset (2007), Sana et al. (2012), Schneider et al. (2014), and Wang et al. (2020), one must keep in mind that this naive interpretation of the location of stars in the sHRD in evolutionary terms becomes a bit more complex if we take into account that a high percentage of stars in our working sample might have suffered from some type of interaction with a companion, even if they are not presently detected as a spectroscopic binary. Therefore, in this section, we restrict ourselves – as a first step – to highlighting the main characteristics of the various v sin i distributions presented in Fig. 8 from a purely empirical point of view. We refer the reader to Sect. 4 for a further discussion and interpretation of results taking into account our present theoretical knowledge about single and binary star evolution of high-mass stars.

Roughly speaking, it becomes clear from inspection of the various panels in Fig. 8 that the individual v sin i distributions of the LS and SB 1 samples comprising each of the three subgroups retain similar properties to the corresponding ones resulting from the global sample (see Sect. 3.2.2). In particular, we refer to the bimodal character of the v sin i distributions associated with the LS samples, and the more uniform distributions found in the cases of the SB1 samples.

To complement the information presented in Fig. 8 and help to better identify additional similarities and differences between the various v sin i distributions, we summarize some statistical quantities of interest, extracted from each of the distributions, in Table 2. In all cases, we provide separated information for two v sin i ranges, using 150 km s⁻¹ as boundary.

The first remarkable result is the striking similarity¹¹ of the v sin i distribution associated with the LS samples comprising groups “a” and “b”. In both cases, there is a clear separation between the low- and fast-usim components at ~150 km s⁻¹. In addition, the statistical properties of both components in the two groups (in terms of relative percentage of stars, associated mean υ sin i value, and standard deviation, see Table 2) are basically the same. Also, the corresponding v sin i distributions for the SB 1 samples are remarkably similar (although with a much lower statistical significance), especially for the covered range in v sin i, extending from ~50 to ~200 km s⁻¹ (except for one star in group “a”, which is very close to the boundary with group “c”).

This empirical result has far-reaching implications that will improve our understanding of the importance that angular momentum losses because of stellar winds, as well as several proposed mechanisms that could transport angular momentum from the stellar core to the envelope, have on the evolution of surface rotational velocities in stars with masses above 30 M_⊙ (see Sect. 4.1). In this regard, if we assume that there is an evolutionary connection between both groups, the abovementioned result would imply that no important braking is detected in any of the two main components of the v sin i distribution as stars evolve from region “a” to region “b”.

Indeed, the mean v sin i value associated with the low-velocity component of the LS distributions is slightly higher in the case of group “b”, something that is not expected from and evolutionary point of view. However, this result is most likely a consequence of the methodological limitation associated with our v sin i measurements. As commented in Sect. 3.2.3, this limitation may affect the relative number of stars populating the first two bins of the distributions. The lack of stars produces a shift in the associated mean v sin i value of the low-velocity component to somewhat larger values, especially in the case of stars in group “b”.

Again (see also Sect. 3.3), inspection of the v sin i distributions of stars in group “c” (and particularly that of the LS sample) allows us to confirm that the abovementioned methodological limitations do not significantly affect our v sin i estimates in the low-luminosity O-type stars. Contrarily to the case of the higher-mass samples, in which there is a very low percentage of stars with v sin i <40–50 km s⁻¹, the sum of the first two bins of the v sini distribution of the LS sample comprising group “c” reaches almost 30%. This feature, together with: (a) the existence of a clearly separated group of LS stars populating the v sin i range between ~250 and 450 km s⁻¹ in the lower-mass sample, which is not present at higher masses; and (b) the different v sin i limit above which no SB1 stars are detected in groups “c” (~300 km s⁻¹) and “a” (~200 km s⁻¹), seems to indicate that there is a dependence on mass in the efficiency of the various physical processes, which may affect the rotational properties of high-mass stars in the O-star domain along their main sequence evolution.

Fig. 8 Rotational properties of our working sample of Galactic O-type stars separated in three groups (from top to bottom), comprising stars covering different ranges in mass and evolutionary status. Left panels: location of each subgroup in the sHRD; see also Fig. 1 for additional information about other features presented in each sHRD. Middle and right panels: Corresponding v sin i distributions for the LS and SB1 samples, respectively, within each group. In each panel, we indicate the v sin i values of 50, 150, and 300 km s−1 with vertical dashed lines, for reference purposes. Similarly to Fig. 4, the stars in the LS sample, for which we only count on one spectrum, are indicated in gray in the corresponding v sin i distributions.

4 Comparison with well-established theoretical predictions

In this section, we provide further insights on the interpretation of the distribution of our working sample of O-type stars in the various diagrams presented in Sect. 3. To this aim, as a guideline we use some of the theoretical predictions resulting from well-established studies that investigate the importance that single and binary star evolution can have on the rotational properties of massive stars (e.g., Brott et al. 2011; Ekstrõm et al. 2012; de Mink et al. 2013; Wang et al. 2020). Being aware of the existence of other more recent model computations obtained with other single and binary stellar evolution codes, such as PARSEC (Bressan et al. 2012), FRANEC (Degl’Innocenti et al. 2008), BASTI (Hidalgo et al. 2018), or BPASS (Eldridge et al. 2008), we decided to mainly focus on the abovementioned computations because they are well established within the community and have been the object of previous comparisons.

As a disclaimer, we remark that we do not intend to tight-fit the empirical distributions using the available sets of single star evolutionary models and binary population synthesis simulations. Rather, we use the predicted trends to better understand the empirical properties of the sample under study and, at the same time, to evaluate possible observational constraints that our empirical results could provide to presently available (as well as future) theoretical developments.

4.1 Single star evolutionary models

Figure 9 shows the behavior of the two sets of single star evolutionary models including the effect of rotation that have been commonly used by the massive star community in the last decade. Namely, those of Brott et al. (2011) and Ekstrõm et al. (2012), hereafter quoted as the Brott+11 and Ekstõm+12 models. In particular, for the case of the latter models, we chose those computed with an initial rotational velocity (v_ini) being 40% of the critical velocity (v_crit) at the ZAMS. For consistency, among the models available in the Brott+11 database¹², we selected those assuming a v_ini ~ 220 km s⁻¹.

The top panels in Fig. 9 present the evolutionary tracks of both set of models in the sHRD. For reference purposes – and to better identify the three regions in the sHRD discussed in Sect. 3.5 (see also Fig. 8) – we also depict the log g = 3.7 dex line of constant gravity, and the ZAMS and the 20, 32, and 85 M_⊙ evolutionary tracks corresponding to the Ekström+12 models without rotation. The bottom panels in the same figure show the predicted evolution of equatorial velocity (v_eq, without inclination effect), using T_eff as a proxy of time (see also the corresponding empirical diagram in Fig. 6).

In all panels, tracks are colored depending on the value of v_eq, and the dashed red line marks the point on each track where v_eq has dropped by 20% of its v_ini. We also highlight, as a gray shadowed area, the region in the sHRD close to the ZAMS, which is void of observed O stars (Holgado et al. 2020). Knowing the specific location of this region in the v_eq–T_eff diagram is of critical interest when we compare model predictions and the distribution of observational points. If this particular region of the v_eq–T_eff diagram is void of stars, it could mean that either stars with these properties do not exist, or they have escaped from our optical spectroscopic surveys.

As expected, Fig. 9 shows that the models for single massive stars predict a braking effect that is stronger for the most luminous stars because of stronger winds and enhanced mass loss (Brott et al. 2011; Ekström et al. 2012). It also shows that the Ekström+12 models slow down faster than the Brott+11 models, the effect being more clearly seen in the lower panels. This is mainly a consequence of the specific treatment of an internal magnetic field in each of the two evolutionary model computations, which in turn causes the consequent reinforcing and/or limitation of angular momentum transport within the star (see Martins & Palacios 2013; Keszthelyi et al. 2017, 2020, and references therein).

The model predictions presented in Fig. 9 are also complemented with the corresponding distribution of empirical points in each of the considered diagrams. From first inspection of the lower panels of this figure, it becomes already clear that – although the v_eq–T_eff (or v sin i–T_eff) diagram has been used elsewhere to contrast model predictions about the spin evolution of massive stars with observational data – the comparison of empirical and modeled values using exclusively this diagram is complex. In particular, we have to take into account that stars with a given T_eff and different luminosities (and thus masses) can appear mixed in this diagram, and are not necessarily placed near their correct corresponding evolutionary track.

Neither of the two sets of evolutionary models with rotation is able to reproduce the lack of stars with v sini values between ~75 km s⁻¹ and the considered v_ini at the hot side of the diagrams, an effect that is more evident in the Ekström+12 models. Moreover, Brott+11 models with v_ini ~200 km s⁻¹ cannot explain the existence of low-velocity stars across the whole O-star domain, as they do not slow down sufficiently before reaching temperatures below 30 000 K.

These results seems to indicate that single star evolutionary models with an initial equatorial velocity of ⪅150 km s⁻¹ might be more convenient to explain the spin-rate properties of a high percentage (75%, as stated in Sect. 3.2.1) of the stars in our sample. However, this consideration will still leave us with an unsatisfactory explanation for the stars with rotational velocities above the selected v_ini.

In Fig. 10, we show the distribution of empirical points in the v_eq–log (L/L_⊙) diagram, along with the predicted behavior of the evolutionary tracks for the same set of evolutionary models considered for Fig. 9. For a better comparison, we present a separated diagram for each observational subsample on the three areas of the sHRD defined in Fig. 8 (from top to bottom), as well as for the Ekström+12 (left) and Brott+11 (right) models. In each panel, we highlight the section of the evolutionary tracks that corresponds to the represented area with a thicker colored line. This is the part of the evolutionary track that has to be compared with the subsample of observational points depicted in each panel. The other parts of the tracks, represented as thin dashed lines, are only included for reference; they do not cross the corresponding area in the sHRD (see Fig. 9).

As in previous figures, we also mark the sections of the evolutionary tracks that correspond to the area near the ZAMS void of stars with a thick gray line. Lastly, we note that in this figure, the temporal evolution goes from right to left, and we recall that the quantity represented in the x-axis is different for the model (v_eq) and the observations (v sin i); thus, the projection effect must be taken into account for a correct interpretation of the information presented in this figure.

At first glance, the empirical distributions do not coincide with the expected theoretical ranges in most cases. In panels a and b – namely stars with M > 32 M_⊙ at gravities higher and lower than log g= 3.7 dex, respectively – most of the observational points appear at rotational velocities much lower than those predicted by the Brott+11 models (right column). Although there is a projection effect involved, the high number of points clustering around v sin i ~ 65 km s⁻¹ (see also Fig. 8 and Table 2), allows us to discard this effect as a possible explanation of the discrepancy. This result reinforces our previous statement that – if interpreted in terms of Brott+11 model predictions – an initial velocity of ~150 km s⁻¹ (or even lower) seems to be more appropriate to explain the rotational properties of more than 70% of the O stars in the high-mass region.

The situation is somewhat different – but still not very encouraging – for the case of Ekström+12 models (left panels). As previously commented (see also bottom left panel in Fig. 9), the specific physical ingredients considered by Ekström et al. (2012) in their evolutionary model computations lead to an important braking of the stellar surface soon after the star leaves the ZAMS, especially in the case of stars with masses above ~30 M_⊙. In addition, as illustrated by the top left panel in Fig. 9, the Ekström+12 evolutionary tracks depart from their non-rotating counterparts, this effect becoming more important for increasing masses. Both effects together make it possible to explain the existence of stars in region “a” with log (L/L_⊙) in the range ~4.0–4.2 dex, with both low and intermediate values of v sin i (even if models with v_ini/v_crit = 0.4 are considered). Also, the 32 and 40 M_⊙ evolutionary tracks reach region “b” with equatorial velocities that are compatible with the low-v sin i component of the empirical distribution of points. However, there is still a large amount of stars in regions “a” and “b” with a location in the v_eq–log (L/L_⊙) diagram that is incompatible with the considered Ekström+12 models.

The top and middle panels of Fig. 10 also allow us to contrast which of the two model predictions regarding the spin evolution of stars in the mass range between ~30 and ~80 M_⊙ is in better agreement with the observations. As already commented in Sect. 3.5 (and illustrated by Fig. 8), the relative percentages of stars in regions “a” and “b” with intermediate values of v sini (50-100 km s⁻¹) are quite similar. This is in contrast with the strong braking of the stellar surface predicted by the Ekström+12 models, but in relative fair agreement with the more moderated spin down characterizing the Brott+11 models.

In this context, we must remember that there is empirical evidence indicating that our line-broadening analysis strategy might be failing for some stars for which we obtained v sin i estimates between ~40–50 km s⁻¹ (see Sects. 3.2.3 and 3.4); however, as illustrated by Fig. 7, this methodological limitation should affect the stars in the low-v sin i boundary of the distribution in both regions “a” and “b” in a similar way, and thus we expect the abovementioned statement to still be valid once we are able to overcome this methodological issue.

Furthermore, the strong braking effect present in Ekström+12 models makes it challenging to find a straightforward explanation – even if we assume a higher value of v_ini/v_crit – for the existence of evolved stars (i.e., those in region “b”) with log (L/L_⊙) ≳ 4.0 and v sin i values in the range 50–300 km s⁻¹. Indeed, this conclusion is not affected by the limitations of our line-broadening analysis strategy since we are referring to v sin i estimations well above the problematic range.

This latter result could be considered as a strong empirical argument to support the various decisions taken by Brott et al. (2011) regarding the treatment of the internal angular momentum transport in the evolutionary model computations. In particular, the activation of the effect of internal magnetic fields (via the Taylor-Spruit dynamo) on this transport helps to compensate the braking of the stellar envelope produced by the presence of strong stellar winds in the more massive stars.

However, we must emphasize that this is not necessarily the last word on this matter, since the mass loss recipes assumed in the model computations also have important consequences for the predicted evolution of surface angular momentum in massive stars (Keszthelyi et al. 2017). For example, these authors showed that, decreasing the usually considered mass loss rates by 30% (from Vink et al. 2000), they could make a model without an internal magnetic field result in a similar evolution of v_eq as a model with a magnetic field with the original values of mass-loss rate. Interestingly, Björklund et al. (2022) has recently provided a new prescription for steady-state, radiation-driven mass loss from hot, massive stars, also indicating that for the case of O stars, the predicted mass loss rates are lower by about a factor of three than the rates typically used in previous stellar evolution calculations. It would therefore be interesting to perform a similar comparison to the one presented in Fig. 10, but assuming a 30% of reduction in the mass loss rates incorporated in the Ekström+12 models (and a smaller value of v_ini/v_crit, see Sect. 5).

This problem of reconciling observations and models was already pointed out by different previous studies. Markova et al. (2018), using 53 Galactic O-type stars and the same sets of models, found that none were able to reproduce the helium and nitrogen surface enrichment appropriately when considering only the expected behavior of the v sin i evolution included in single star rotating evolutionary models. Also, Martins et al. (2017) found difficulties in explaining the chemical properties of 15 late O giants at low metallicity using the single star evolution models.

Regarding stars in region “c”, the main problem arises from the fact that neither of the two sets of models is able to explain the broad range in v sin i that characterizes the empirical distribution. In this case, angular momentum losses due to stellar winds are much less relevant, and thus the only way to be able to reconcile single star evolutionary model predictions with the empirical distribution of points would be to consider that stars in this mass range are born covering a very broad range of initial spin rates (from almost zero and up to ~60% of the critical velocity at the ZAMS, see Fig. 11). Indeed, to some extent, a similar problem is also present in the upper mass region, especially if we assume that angular momentum losses because of stellar winds do not result in an important braking effect of the stellar envelope (and thus the observed values of v sin i).

In view of the discussion above, one would be tempted to explain the various features characterizing the observational distributions presented in Figs. 6, 7, and 10 by means of a combination of Brott+11 models with different vini. However, this would result in a very ad hoc distribution of initial rotational velocities. For example, an anticorrelation between the maximum initial rotational velocity and the initial mass would be needed, in order to reproduce the lack of massive stars with high rotational velocities close to the ZAMS and the large number of lower-mass stars with high rotational velocities. In addition, as already commented, in the 20–40 M_⊙ range, we would need a much broader range of initial spin rates than in the case of stars with higher masses (see also further notes in Sect. 5). Last but not least, there are some features in the distributions that could be only explained by assuming that some stars are spinning up along their evolution, something that is not possible in the case of single star evolution.

Fig. 9 Predicted behavior of the equatorial rotational velocity in the single star evolutionary models, with rotation computed by Ekström et al. (2012, left panels) and Brott et al. (2011, right panels). See Sect. 4.1 for an explanation of the selected values of initial rotational velocity in each set of models. Cross symbols represent time steps in each evolutionary tracks, colored accordingly to the corresponding veq. The dashed red line marks the point in each track with a 20% reduction of veq from vini, and the yellow line indicates the log g = 3.7 dex line of constant gravity. Observations are represented by gray points. The shaded gray area highlights the area in the sHRD close to the ZAMS, which is void of observed stars. Top panels: Considered evolutionary tracks in the sHRD, where we also depict the ZAMS and the 20, 32, and 85 M⊙ evolutionary tracks corresponding to Ekström+12 models without rotation to better identify the three regions in the sHRD (marked as “a”, “b”, and “c”, respectively) discussed in Sects. 3.5 and 4.1 (see also Figs. 8 and 10). Bottom panels: Considered evolutionary tracks in the veq − Teff diagram (or v sin i–Teff for the case of the observations).

Fig. 10 Predicted behavior of the evolutionary tracks computed by Ekström et al. (2012, left panels) and Brott et al. (2011, right panels) in the Oeq-log (L/L⊙) diagram. From top to bottom, the sections of the evolutionary tracks, which are located within the three regions (“a”, “b”, and “c”, respectively) indicated in Fig. 8 and described in Sect. 3.5, are highlighted with blue and green thicker lines. In addition, the gray thick part of the evolutionary tracks indicates that the corresponding section of the track lies inside the region of the sHRD close to the ZAMS where no O stars are detected. The stars in the observed sample that populate the three considered regions of the sHRD are also included in each panel, separately, for comparison (this time using v sin i values instead of veq). As in other figures, we indicate three v sin i values of interest (50, 150, and 300 km s−1) with vertical lines. Size represents log q value.

4.2 Exploring the effect of binary interactions

The present formation scenario of massive stars states that ~90% of them are born as binary¹³ systems (Sana et al. 2012), As a result, it is natural to think that the features in the various diagrams presented in previous sections that cannot be easily explained by means of single star evolutionary models may be produced by binarity effects. To evaluate this, we follow the ideas presented by Ramírez-Agudelo et al. (2015) in their study of the rotational properties of O-star binary systems in the 30 Doradus regions of the LMC, which, ultimately, are based on theoretical results obtained by de Mink et al. (2013).

The evolution of the individual spin rates in a binary system can be affected by several additional effects not present in the case of isolated stars, some of them giving rise to a more or less important acceleration of one of the two components (or even both under certain circumstances). For example, in close binaries (P ≲ 10 days), tidal interaction may synchronize rotation and orbit¹⁴ (either increasing or decreasing the equatorial velocities). But more importantly in the context of our study, mass (and angular momentum) transfer occurring after Roche-lobe overflow can spin up the initially less massive star to a very high rotational velocity, which can even reach critical velocity in some cases (see Petrovic et al. 2005; de Mink et al. 2013, and references therein). In this situation, the secondary star in the binary system does not only experience an increase in rotational speed, but it also modifies its luminosity and apparent age. Eventually, this mass-transfer event could modify the appearance of the binary system, leading to a fast-rotating, rejuvenated (apparently) single object in many cases (Ahumada & Lapasset 2007; Schneider et al. 2014, 2015, 2016; Wang et al. 2020). All of this simulated behavior should be considered with caution, since it relies on some physics that are not fully constrained. There are still different criteria when considering the amount of mass and angular momentum that is lost in collisions and mass transfer events (Wang et al. 2022), as well as in the generation of a magnetic field and its braking properties in mergers events (Schneider et al. 2019).

The binary population synthesis simulations performed by de Mink et al. (2013) showed that binary interaction offers a natural explanation for the bimodal characteristic of the global rotational distribution commonly found in the O-star domain (Sect. 3.2.1). By assuming a flat initial rotational distribution between 0 and 200 km s⁻¹ – and a continuous star formation – they found that, after several Myrs, the effects of evolution, tides, mass transfer, and mergers modify this initial distribution to another one, including a dominant low-velocity component and a tail of fast rotators. In this context, the simulations by de Mink et al. (2013) predict that 19% and 11% of the stars will develop rotational speeds above 200 and 300 km s⁻¹, respectively. These numbers can be compared with the 17% and 7% resulting from our empirical distribution of projected rotational velocities. In this context, we find it interesting to indicate that the total percentage of detected spectroscopic binaries in our sample of 415 stars is ~41% (55 SBls and 113 SB2s). Indeed, it could reach ~50% if we roughly consider half of the LS stars with a single spectrum as hidden binaries. These numbers must be taken into account when comparing with the simulation from de Mink et al. (2013), who assumed an initial binary fraction of 70%. That would imply that about 20% of the present single stars should have been born as binaries. Also, more in-depth future studies connecting the complete evolutionary path followed by massive stars should take into account the increasingly available empirical information about detected spectroscopic binaries in more evolved stages, such as B supergiants (Simón-Díaz et al. 2020a, and in prep.), red supergiants (e.g., Neugent 2021; Dorda & Patrick 2021), or WR stars (e.g., Dsilva et al. 2022), which routinely obtain much lower binary fractions.

Aside from projection effects, we can expect some differences between our empirical results and de Mink’s simulations due to processes favoring a faster spin down than modeled, such as the formation of magnetic fields during mass accretion (as mentioned by de Mink et al. 2013). Evidence supporting such processes (not necessarily a magnetic field) is found in the observation of super-synchronization but clearly subcritical rotation in OB+WR systems (Shara et al. 2020). Moreover, we note that the de Mink et al. (2013) simulation includes SB2 stars that have not been considered in our observational sample, and that the difference in metallicity – LMC for de Mink’s study – could lead to a higher number of high-spin objects in their simulations. Considering all these uncertainties, the agreement seems to be surprisingly good. The scenario proposed by de Mink et al. (2013) provides a more natural way to explain the distribution of stars in the v sin i–T_eff and v sin i–log(L/L_⊙) diagrams (Figs. 9 and 10, respectively) than the assumption of a very ad hoc initial spin-rate distribution (see Sect. 4.1). Despite considering that the most straightforward solution for bimodality seems to be the effect of binary interaction and its spin-up effect on stars, it is not possible to totally rule out the existence of stars that are born with an extremely high rotation naturally, although as assumed in Sect. 4.1, it may seem a more artificial solution.

The second empirical feature that we evaluate is the clear dominance of apparently single stars in the tail of fast rotators (especially above 250–300 km s⁻¹ in the lower-mass domain, and 150–200 km s⁻¹, for higher masses, see Fig. 8). This result agrees with the hypothesis that the tail of fast rotators is mostly populated by post-interaction products of binaries (de Mink et al. 2013) that, in addition – as indicated by de Mink et al. (2011, 2014) – will often appear to be single because the companion tends to be a low-mass, low-luminosity star in a wide orbit. Alternatively, they became single stars after a merger or disruption of the binary system during the supernova explosion of the primary.

On a related matter, we consider, as a third characteristic to be interpreted in the context of binary evolution, the apparent upper limit in v sin i of ~300 km s⁻¹ found in the sample of clearly detected SB1 stars (see Fig. 4). This phenomenon is also present, although not as abrupt, in the study by Ramírez-Agudelo et al. (2015), in spite of the different metallicity.

In principle, our methodology for distinguishing between LS and SB1 stars (see Appendix A) could be dismissed as a plausible factor, making this results to be spurious. Being aware that the detection of small changes in v_rad due to orbital motions becomes more complex when the rotational broadening becomes larger, we see no obvious reason indicating why the value of v sin i = 300 km s⁻¹ could represent a limitation. Indeed, it is also interesting to note that, although we identify some SB1 stars among the lower-mass sample (see Fig. 8) with v sin i in the range ~200–300 km s⁻¹, there is only one clearly detected SB1 star with v sin i in this range in regions “a” and “b”.

We could also claim that this result is just a consequence of the limited number of epochs available for the sample of fast rotators. However, as indicated in Table F.2, we count on a minimum of five epochs for most of them (except for two stars)¹⁵.

Following de Mink et al. (2014), our Sb1 star sample can be assumed to be dominated by pre-interaction objects – in other words, previous to mass exchange – so the only phenomenon that may be altering the rotational velocity beyond that expected by single star evolutionary models is tidal interactions. If we look at the results provided by de Mink et al. (2013), these tidal interactions have several visible effects on their simulated distribution, such as a tail of accelerated stars that reaches just 300 km s⁻¹. Hence, the effect of tidal interactions could explain why the maximum speed reached by our sample of clearly detected SB1 stars, if mostly pre-interaction, is precisely that value.

The effect of tidal interactions and orbital-rotation synchronization may also be the origin of the difference between the distributions of single and SB1 stars (Fig. 4), the latter being characterized by a wider low-v sin i peak and a somewhat higher mean value of v sin i when compared to the former (see Sect. 3.2.2). Furthermore, de Mink et al. (2013) also indicate that tidal effects prevent stars in a binary system from slowing down excessively, something which could also explain a lower frequency of stars with v sin i < 40–50 km s⁻¹ found among the SB1 sample (compared to the LS sample).

All these three characteristics are also found in the work by Ramírez-Agudelo et al. (2013, 2015), although with small differences (e.g., they have higher v sin i velocities, also for binary components) and without considering stellar masses and effective temperatures. These authors delved into the empirical existence of the bimodal distribution of v sin i and the differences between apparent distributions of single and binary components (mainly primaries). They also used the works by de Mink et al. (2013, 2014) for the interpretation of their results, finding good agreement between observations and simulations. Our study reproduces their findings. Interestingly, however, in comparison with Ramírez-Agudelo’s study, ours represents not only a different metallicity but also a different environmental origin, as the LMC sample comes from a single cluster, 30 Dor. This reinforces that the origin of the bimodality and the lower limit in v sini among SB1 as compared to apparent single stars is, at least in the first order, independent of metallicity and environment.

The fourth and last characteristic we want to investigate is the specific location in the sHRD of those stars in the sample with v sin i > 300 km s⁻¹ (Sect. 3.5 and Fig. 8). As already discussed, the origin of such a high velocity is probably the binary interaction. Mass accretion or merger during the evolution of the binary systems may produce a rejuvenation and spin up of the initial secondary star or the final product (if no additional mechanisms such as the formation of magnetic fields play a role, de Mink et al. 2013, 2014). They may even appear as blue stragglers in a cluster when compared to the rest of the population (Ahumada & Lapasset 2007; Schneider et al. 2014, 2015, 2016; Wang et al. 2020). The works of Ahumada & Lapasset (2007) and Schneider et al. (2016) find that the number of blue stragglers is larger, and their rejuvenation more intense, for smaller masses and more evolved clusters. This, together with a faster evolution at higher masses, could explain the concentration of very fast rotators that we see below 32 M_⊙.

Fig. 11 Properties of the supposed youngest stars in the sample. Left: sHRD highlighting a subsample of 118 LS O-type stars with log g > 3.7 dex. The symbol size depends on the value of the ratio v sin i/vcrit, as indicated. Right: υ sin i/vcrit distributions (as normalized histograms) for the indicated subsample, separating those stars located above (upper panel) and below (lower panel) the 32 M⊙ track. The vertical dashed lines mark several values of interest for the v sin i/vcrit ratio (see text for explanation). The vertical dotted line at v sin i/vcrit = 0.06 marks the approximate limit below which the v sin i determination can be affected by methodological limitations in those stars located above the 32 M0 track).

5 A final note on the initial spin-rate distribution

The ratio of initial equatorial velocity (v_ini) against the critical velocity (v_crit) is an influential input parameter for any evolutionary model of massive stars (see review by Maeder & Meynet 2000, and references therein). As a consequence, having access to reliable empirical knowledge about the initial spin-rate distribution of main sequence O- and B-type stars is of utmost importance for understanding the importance of rotation on the evolution and final fate of these important stellar objects. In addition, it provides a key constraint for our theories of massive star formation, and is an essential ingredient in population synthesis models.

There are several empirical studies in the literature that focus on the study of the main characteristics of the v_ini/v_crit distribution associated with stars covering different mass ranges in the high-mass domain and located in different environments (e.g., Wolff et al. 2006; Martayan et al. 2006; Huang et al. 2010). Of particular interest for comparison with our study are the works performed by Wolff et al. (2006) and Markova et al. (2014). In their investigation of the spin-rate properties of a sample of newly formed stars spanning a range in mass between 0.2 and 50 M_⊙, Wolff et al. (2006) included 44 O-type stars. These authors, who define v_obs as the projected rotational velocity corrected statistically from the projection effect, find that the average ratio v_obs/v_crit for stars between¹⁶ 25 and 50 M_⊙ is 0.2, and 0.13 for stars between 8 and 25 M_⊙. We should note, however, that in their study, Wolff et al. (2006) (a) neglected the effect of radiation pressure on the estimation of v_crit (that may be important for the first group), (b) did not consider possible binarity effects, (c) did not make any attempt to disentangle the effect of macroturbulent broadening, and (d) as commented above, did a statistical correction to account for the effect of the inclination angle.

More recently, Markova et al. (2014) extended the investigated mass range (now also disentangling the effect of macroturbulence on the line-broadening analysis) by studying 31 apparently single O stars, with 18 of them having initial masses between 50 and 80 M⊙. They find that above 50 M⊙ the ratio v sin i/v_crit does not exceed 0.26 (and is mostly below 0.17), whereas they find a limit of 0.50 between 22 and 50 M⊙, with most stars rotating at less than 30% of their critical velocity. Indeed, only two stars in their sample exceed 200 km s⁻¹.

The sample analyzed in this paper allows us to review and considerably extend the works performed by Wolff et al. (2006) and Markova et al. (2014) thanks to, among other things, the more complete mass range coverage, and the identification of SB1 stars.

To construct an estimated v sini/vcrit distribution¹⁷ from our data, we need to select the v sin i of the stars that have evolved as little as possible since birth. To this aim, we selected stars in our sample having log g > 3.7 dex, using gravity again as a proxy for evolution. In Sect. 3.5 we showed that the v sin i distribution tends to maintain shape during evolution (see also top and middle panels in Fig. 8). This further reinforces our consideration of selecting stars up to the log g = 3.7 dex line of constant gravity to approach an initial rotation distribution while trying to keep statistics as large as possible in the resulting distribution. Lastly, we discard SB1 stars, as these have most likely suffered variations in their vini since the beginning of their evolution because of the influence of their companions (see Sect. 4.2).

Although the v sin i values of the subsample described above should approximate in some way the initial velocity distribution of Galactic O stars, we must keep in mind that – taking into account the ideas presented in Sect. 4.2 – its similarity with the real distribution will largely depend on the net effect of the binaries on the sample, and the absence of a strong and fast braking effect from the ZAMS to the position of these stars.

The left panel of Fig. 11 presents the location of the above-mentioned subsample in the sHRD, whereas the right panel shows the associated v sin i/v_crit histograms, further separating the investigated subsample in two (considering the stars located above and below the evolutionary track of 32 M⊙).

Both histograms show a peak near v sin i/v_crit = 0.1, with the lower-mass stars extending to lower values (probably as consequence of the limitations in obtaining low v sin i values for stars of high luminosity, see Sects. 3.2.3, 3.4 and Appendix. D). But the most obvious feature in the two histograms is the existence of low- and high-velocity components clearly separated by a gap (between ~0.175 and 0.25 for stars with M > 32 M⊙, and between 0.325 and 0.45 for stars below this mass.

Following the discussion in Sect. 4.2, we associate these gaps to the effect of binary evolution. On the one hand, the region with v sin i/v_crit < 0.175 would be dominated in both cases by single stars or wide binaries that have suffered little evolution, and thus are representative of the initial spin-rate distribution. On the other hand, following de Mink et al. (2013, 2014), the fast rotating components are most likely populated by post-interacting binaries. Therefore, these high-velocity components, which indeed comprise a small percentage of the investigated sample, should not be associated with the initial spin-rate conditions of stars at birth.

The difference in the upper boundary of the low-velocity components (and thus the location of the gap) in the high- and low-mass distributions may have different interpretations. One possibility would be to attribute it to the influence of the stellar winds. Surface angular momentum losses because of stellar winds are expected to play a more important role in the higher-mass domain, thus producing some additional braking when compared to lower-luminosity stars. However, we recall that, in Sect. 4.1, we showed that stellar winds do not seem to be very effective at producing an important braking of the stellar envelope as the more luminous O-type stars evolve. It could also be due to a small mass dependence of the covered range of the vini/vcrit ratio, with decreasing values for increasing initial stellar mass. This would be the opposite effect to the one presented by Wolff et al. (2006), but it agrees with the results by Huang et al. (2010) for 220 high-gravity B stars. They find a decreasing range of v_eq/v_crit from 0.8 for M = 2 M⊙ stars to 0.2 on the frontier with O stars. Lastly, the low-mass sample could still have a small contamination of SB1 stars between v sin i/v_crit ~ 0.2-0.3, as a significant number of stars in this range have less than two epochs (see Fig. 8).

Whatever the reason for the differences between the two v sin i/vini distributions, and even considering projection effects, we conclude that O stars in their unaltered early stages have velocities between v sin i/v_crit – 0.1-0.175. If we rule out a fast and extreme braking effect just next to the ZAMS, not present in any of the two sets of evolutionary models considered here, O-type stars are typically born with rotational velocities near v_ini/v_crit ~ 0.1 and not significantly exceeding v sin i/v_crit = 0.2 (in agreement with the results presented by Markova et al. 2014). Therefore, evolutionary models for single stars should not use higher values (except perhaps for some peculiar individual cases). Stars in binary systems will be spun up in the course of their evolution through binary interaction, populating the upper v_ini/v_crit values up to ratios ≤0.7.

6 Summary and conclusions

This article focuses on the results for the rotational analysis (v sin i and v_mac) of the 285 likely single and SB1 stars of the total 415 galactic O stars in the IACOB and OWN sample, excluding SB2 systems and peculiar stars. This represents a step forward from previous studies in terms of the number of stars (+43%), the separation between extra-broadening components, and the separation between different multiplicity statuses. We made use of spectroscopic parameters obtained for the sample in previous articles (T_eff, log g) and the multiplicity statuses to present a combined vision in the study of rotational properties. The sample was divided into subsamples to investigate rotation on each of them, and their relationship in terms of evolution. The youngest subsample was studied in particular depth with the intention of exploring the vini of the O-type stars. In each step we compared the empirical results with predictions by evolutionary models and population syntheses of O stars.

The main conclusions we reached with this work are:

• The bimodal appearance of the general v sin i distribution persist with respect to previous works, as well as the absence of SB1 stars with extremely high rotational speeds. We support binary interaction, as depicted by the de Mink et al. (2013) simulations, as the plausible origin of the bimodality feature, and the difficulty to detect binarity in the spun-up binary products, explained in de Mink et al. (2014), for the absence of very fast SB1 stars. All of this is in accordance with the Ramírez-Agudelo et al. (2013) results in the LMC.

• The extended width of the v sin i distribution of SB1 stars compared to the single distribution and its abrupt cutoff at 300 km s⁻¹ agrees with tidal effects, since we could consider them stars without probable mass exchange interaction episodes (de Mink et al. 2013).

• We found an unexpected scarcity of very low-rotation stars (<~40–50 km s⁻¹) that appear to be concentrated on the lower part of the sHRD. We claim microturbulence as the plausible origin, whose effects have not been completely eliminated, and which depends strongly on the location of the star in the sHRD (see Simón-Díaz & Herrero 2014).

• Stars with a very high rotation (>300 km s⁻¹) are grouped in the sHRD mostly below the 32 M⊙ evolutionary track. Their grouping in the sHRD may have its origin in the faster evolution in the higher-mass interaction products and/or mergers witnessed in the simulations of Wang et al. (2020), and the increase in blue stragglers for smaller masses (Ahumada & Lapasset 2007; Schneider et al. 2014, 2015).

• We do not detect a strong braking effect on the evolution of massive stars, even above the 32 M⊙ evolutionary track. This weak braking is in accordance with the Brott+11 models (Brott et al. 2011), instead of the strong and fast braking obtained in the Ekström+12 models (Ekström et al. 2012).

• Using the v_ini considered in most current single evolutionary models (v_ini/v_crit ~ 0.4), the presence of stars with v sin i > 300 km s⁻¹ cannot be explained. Binary interaction effects, as depicted in de Mink et al. (2013), could be a plausible explanation for the tail of the distribution and the slow rotation peak.

• We find that the empirical rotational properties of the “young” subsample of O-type stars is also close to a value v_ini/v_crit ≤ 0.2, in agreement with the trend for B stars presented in the literature (Huang et al. 2010). Considering that the empirical distribution does not changes dramatically with evolution, this could be representative of the initial v_ini/v_crit.

In short, the empirical results lead us to believe that the majority of massive stars are only able to lose between 10 and 20% of their rotation rate from birth. We consider that, in order to generate the observed distribution, the bulk of massive stars would need to be born with rotation rates closer to ~ 150 km s⁻¹. This corresponds to an initial rotation of ~20% of critical velocity. These results match with the rotational study of O-type stars in the 30 Dor region of the LMC carried out by Ramírez-Agudelo et al. (2015), with the interest of having a different metallicity, and our case being a study that comprises different clusters and field stars within the magnitude-limited sample of IACOB. For a much more productive comparison with modeled values, new sets of models and/or population synthesis with an initial velocity v_ini/v_crit = 0.2, and including binary interaction, should be generated.

Collectively, our data represent a solid empirical basis to investigate the reaction and dependence of evolutionary models to changes in rotational properties, as well as other ingredients, such as mass loss or the inclusion of magnetic fields, as investigated in Keszthelyi et al. (2017) and Keszthelyi et al. (2020). In the future, the general results we found with this work could be confirmed with an even more representative observed Galactic sample using data from the upcoming WEAVE survey (Dalton et al. 2012).

Acknowledgements

In memoriam of Dr. Rodolfo H. Barba, whose dedication to the astrophysics was an inspiration for all. We are extremely grateful to the referee of this paper, Zsolt Keszthelyi, for his valuable comments and for providing a constructive and helpful report to improve this work. G.H., S.S.-D. and A.H. acknowledge support from the Spanish Ministry of Science and Innovation (MICINN) through the Spanish State Research Agency through grants PGC-2018-0913741-B-C22 and the Severo Ochoa Programe 2020-2023 (CEX2019-000920-S). This work has also received financial support from the Canarian Agency for Economy, Knowledge, and Employment and the European Regional Development Fund (ERDF/EU), under grant with reference ProID2020010016. G.H. also acknowledge support from the Spanish Government Ministerio de Ciencia through grant PGC2018-95049-B-C22, and from the Science and Operations Department of the European Space Agency - Contract Number 4000126507-126507/19/ES/CM. Based on observations made with the Nordic Optical Telescope, operated by NOTSA, and the Mercator Telescope, operated by the Flemish Community, both at the Observatorio del Roque de los Muchachos (La Palma, Spain) of the Instituto de Astrofísica de Canarias. Based on observations at the European Southern Observatory in programs 073.D-0609(A), 077.B-0348(A), 079.D-0564(A), 079.D-0564(C), 081.D2008(A), 081.D-2008(B), 083.D-0589(A), 083.D-0589(B), 086.D-0997(A), 086.D-0997(B), 087.D-0946(A), 089.D-0975(A) This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This paper made use of the IAC Supercomputing facility HTCondor (http://research.cs.wisc.edu/htcondor/), partly financed by the Ministry of Economy and Competitiveness with FEDER funds, code IACA13-3E-2493. G.H. wants to thank N. Langer and the Bonn group for their helpful comments that improved this manuscript significantly, and M. Cervino for the precise corrections.

Appendix A Identification of SB1 stars in our working sample of O-type stars

In this appendix, we briefly outline the strategy we followed to identify the stars in our working sample that can be tagged as SB1 to a high level of confidence¹⁸. A more detailed description about this process can be found in Holgado (2019), where we extend the criteria already presented in Holgado et al. (2018) to the case of stars with broad line-profiles (i.e., stars with v sin i > 150 km s⁻¹).

As a first step, we performed a visual inspection of a set of diagnostic lines which allowed us to optimize the identification of SB2 systems and to flag (qualitatively) the stars for which we had more than one epoch as presumably constant (C), line-profile variable (LPV), and clear or most likely SB1 systems. We refer the reader to Holgado et al. (2018) for some illustrative examples of the various types of detected variability.

As a second step, once we had access to the outcome of the iacob-broad and iacob-gbat analyses of those stars not identified as SB2s (see Sect. 2), we performed a (more) quantitative reassessment of the abovementioned classification. In brief, the fastwind synthetic spectrum of the iacob-gbat best fitting model for each star (convolved with the corresponding v sin i, v_mac, and R) was used to improve the radial velocity determination of each associated individual (multi-epoch) spectrum by means of a simple (iterative) cross-correlation technique.

To this aim, we considered the subsequent initial set of lines included in the grid of fastwind models presently incorporated to iacob-gbat: He I λ4026, He I λ4387, He I λ4471, He I λ4713, He I λ4922, He I λ5015, He I λ5047, He I λ5875, He ii λ4200, He ii λ4541, He ii λ4686, and He ii λ5411. Each line was analyzed separately, and the associated spectral window selected automatically using the corresponding synthetic line. The lines that were in emission or absent in the observed spectrum were automatically eliminated from the very beginning on a star-by-star basis. Then we followed an iterative process in which the lines providing v_rad measurements that deviated more than 2σ from the mean were removed. Final values and associated uncertainties were then obtained from the mean and standard deviation of those individual v_rad measurements resulting from the remaining set of lines.

As a last step, we computed the dispersion – σ(v_rad) – and the peak-to-peak amplitude – ∆v_rad – associated with each set of radial velocity measurements for a given star. Then, as described below, we used either of these two quantities to reassess the flags assigned to the targets initially identified as C, LPV, and SB1. In particular, we followed two different criteria depending on the projected rotational velocity of the star:

• if v sin i ≤ 180 km s⁻¹, we flagged those cases with ∆v_rad in the range 5–20 km s⁻¹ and σ(v_rad) in the range 2.5–10 km s⁻¹ as LPV, while those cases with lower (higher) dispersion were flagged as C (SB 1), respectively;

• if v sin i > 180 km s⁻¹, we flagged those cases with a ratio ∆v_rad v sin i in the range 0.02–0.1 and σ(v_rad)/v sin i in the range 0.01–0.05 as LPV, while those cases with lower (higher) ratios were flagged as C (SB1), respectively.

If either of the two conditions was met, the object was reclassified. The indicator with the greatest variability always prevailed, that is, SB1, then LPV, and finally C.

This separation between stars with low and high projected rotational velocities accounts for the fact that subtle changes in the shape of the line profile due to, for example, pulsations, could led to a larger dispersion in the measured (multi-epoch) radial velocities in stars with broad profiles. This situation could be erroneously interpreted as the detection of signatures associated with spectroscopic binarity, if not correctly taken into account.

Figures A.1 and A.2 present a summary of the outcome of this second assessment of spectroscopic variability in the targets initially flagged as OLPV/SB1. The limiting values of ∆v_rad and σ(v_rad) – or their respective ratios with respect to v sin i – quoted above are indicated as vertical dotted lines.

Fig. A.1 ∆vrad (left) and σvrad (right) distributions for the O-type stars in our initial sample not identified as SB2s and having a v sin i≤ 180 km s−1. We indicate with different colors how the stars were finally classified as C/LPV/SB1 in the stars with more than one spectrum. Horizontal dashed lines indicate the thresholds considered to separate the stars in the three categories. For most stars, both indicators coincide. In case of divergence, the higher variability flag prevails (SB1, then LPV, and then C).

Fig. A.2 Similar to Fig. A.1, but for stars with a projected rotational velocity larger than 180 km s−1. We note that, in this case, we consider the ratio of ∆vrad (left) and σvrad (right) with respect to v sin i.

Appendix B Peculiar stars in the sample

As commented in Sect. 2, there are 13 peculiar stars in our sample (six Oe and seven magnetic, see Table F.4) that were not further considered in our study because of the inherent limitations of our analysis strategy to provide reliable stellar parameters for these stars. However, we could obtain υ sin i estimates for all of them. Figure B.1 shows the v sin i distribution for this small sample of stars, where the global v sin i distribution of the remaining LS and SB1 stars is also depicted in the background for reference.

Fig. B.1 Normalized v sin i distribution of the 13 peculiar stars in the sample (six Oe and sevenφlus three) magnetic stars, see text. The global (also normalized) v sin i distribution of the full working sample of 285 LS and SB1 stars is depicted in the background for reference.

Magnetic stars are expected to present a distinct v sin i distribution when compared with the global sample (ud-Doula & Owocki 2002; Grunhut et al. 2017). Even though the coupling between magnetic fields and winds reduces mass loss, the presence of a magnetic field induces an extra braking effect an order of magnitude more efficient than wind-momentum loss (Ud-Doula et al. 2009; Petit et al. 2013).

For the Oe, as explained in Negueruela et al. (2004), they represent the natural continuation of Be stars toward higher masses. There is a strong correlation between fast rotation and the Be phenomenon, manifested clearly when examining the Ω/Ω_crit factor, the surface angular velocity as a fraction of the critical breakup velocity.

The number of magnetic and Oe stars within the sample is low, but Fig. B.1 shows that most Oe stars have v sin i values above the main peak of the distribution (> 100 km s⁻¹). Only one Oe star has a v sin i<50 km s⁻¹, and the origin could be in the inclination effect. For the magnetic stars, the majority present values within the main peak of the v sin i distribution, and only two have an unexpected high value of v sin i (v sin i > 100 km s⁻¹). These are the companion stars of the black hole Cyg X-1 and the post-interacting binary ζ Pup (Howarth & van Leeuwen 2019). We note that the latter, ζ Pup, has a peculiar spectral type (“fp”), but the detection of its magnetic field has been considered spurious (Grunhut et al. 2017).

In addition to these seven magnetic stars with problems in the spectroscopic fitting, the sample includes three other stars recognized as having a detectable magnetic field: HD 57682, HD 37742, and HD 54879 (see e.g., Keszthelyi 2019). These three stars are included in the general sample as they do not exhibit any caveat in the automated fastwind model fit, and they do not represent a very significant sample within the total.

Appendix C Comparison with Ramírez-Agudelo et al. (2013)

Figure C.1 presents the histogram and the associated cumulative histogram with the v sin i distribution obtained by (Ramírez-Agudelo et al. 2013) for a sample of 216 likely single stars O-type in the 30 Doradus region, of the LMC. This is a representative distribution of v sin i values for stars in a different metallicity environment (Z_LMC ~ 0.5Z_⊙). Their v sin i values were derived using the same methodology as our work, separating the macroturbulent broadening contamination. However, the lower resolving power (R~8000) and more limited wavelength coverage of their spectra imposed some caveats in the line-broadening characterization of their sample of stars. For example, their resolution implied that below 40 km s⁻¹ any v sin i determination was uncertain. In addition, their reduced wavelength range force them to rely, in almost all cases, on He i and He ii lines (which are not optimal lines to determine v sin i in O-type stars, Holgado 2019).

Fig. C.1 v sin i histogram (up) and cumulative histogram (down) for Ramírez-Agudelo et al. (2013) and our sample. It includes the results for our sample of Galactic likely single O-type stars (black).

The main expected difference with respect to the Galactic case is that, at lower metallicities, the O-type stars are characterized by developing weaker winds (Vink et al. 2000; Mokiem et al. 2006). This reduces the amount of angular momentum that the star loses through the wind, which would finally produce a v sin i distribution with a higher average speed and a possible greater extreme rotation than in the Galactic case.

While the latter effect can be actually identified when comparing both distributions (the fast-rotators tail of the LMC distribution reach up¹⁹ to 600 km s⁻¹), no clear conclusions can be extracted in the low-velocity peak (below ~40 km s⁻¹, or even ~100 km s⁻¹ in some cases). This is because of the limitations that affect the v sin i measurements of the LMC distribution.

Appendix D Understanding whether microturbulence is the culprit of the scarcity of low-v sin i stars in the O-star domain

We have indicated on several occasions throughout this paper that there is increasing empirical evidence that our line-broadening analysis methodology might be failing under certain limiting circumstances in the low-v sin i regime (see also Sundqvist et al. 2013; Simón-Díaz & Herrero 2014; Markova et al. 2014). For example, in the histogram presented in Fig. 3, there is a notorious scarcity of O stars with spectroscopically inferred v sin i below ~40–50 km s⁻¹ (particularly in the case of stars with LCs, I, II and III). Also, panel A of Fig. 5 highlights that the only 21 stars in our working sample of 285 O-type stars with v sin i<35 km s⁻¹ are concentrated below the 32 M_⊙ tracks in the sHRD. Lastly, Fig. 7 clearly shows that there exists a strong correlation between the lower limit in v sin i imposed by our analysis methodology and the quantity log L.

Fig. D.1 Location in the sHRD of a subsample of 174 LS and SB 1 stars with v sin i < 100 km s−1, separating the targets with a vmac lower or higher than 50 km s−1.

None of these results are expected if we assume a random distribution of inclination angles i and we consider that more than 75% of the investigated stars have a v sin i < 150 km s⁻¹. Under these circumstances, one would roughly expect to find at least ~15% of the stars in the sample (independently of their location in the sHRD) with a v sin i lower than 35 km s⁻¹.

Following ideas presented in Gray (1973), one of the possible culprits that might be at the origin of this problem is the effect of microturbulence. This affects the shape of the line-profiles and their associated FTs (see notes in Sect. 3.2.3 and Simón-Díaz & Herrero 2014).

While a thorough investigation of this subject is planned for a separate paper, we provide a few preliminary results of interest below, which can be obtained using the already available information (including rough estimations of the microturbulent velocity resulting from the iacob-gbat analysis, not quoted in Table F.1, and available at Holgado 2019).

The first result is summarized in Fig. D.1. With this figure we want to highlight that there seems to exist a correlation between the quantity v_mac provided by iacob-broad and the detectability of stars with v sin i below 35 km s⁻¹ (see panel A in Fig. 5). To better illustrate this statement, we only include in the figure the stars with a v sin i estimate below 100 km s⁻¹ (following Simón-Díaz & Herrero 2014, above this value, the v_mac estimates are more uncertain and normally limited to upper limits). We separate the sample in two groups, considering v_mac = 50 km s⁻¹ as the limiting boundary. Interestingly, the region where stars for which a v_mac 50 km s⁻¹ has been measured roughly match the region where no stars with v sin i < 35 km s⁻¹ are found.

This result could be interpreted as empirical evidence indicating that we are still not adequately disentangling the contribution of the so-called macroturbulent broadening from the pure effect of rotation. Alternatively, it could indicate that there are other sources of non-rotational line-broadening that are still affecting our v sin i estimations, microturbulence being one of the possible sources. Indeed, taking into account that, by means of the way we are using iacob-broad²⁰, the quantity v_mac resulting from the GOF analysis encloses any source of non-rotational extra-broadening, a large value of v_mac might also indicate that the line-profile is affected by a relatively high microturbulence. If this is the case, the larger the microturbulence, the lower the frequency (or, equivalently, the higher the “equivalent” v sin i) at which the associated zero will be found in the FT.

Fig. D.2 ξt–log (L/L⊙) diagram including only those stars in our working sample with v sin i < 80 km s−1 and ξt available.

As commented in Holgado et al. (2018), a proper characterization of microturbulence in the O-star domain is still lacking due to the more limited number of metal lines available to determine accurate values of this parameter (e.g., using the curve of growth method), compared to the case of B-type stars. Despite this limitation, several studies have found suggestive evidence of the existence of a correlation between microturbulence – obtained from the analysis of the Hel lines – and luminosity class in these type of stars (e.g., Massey et al. 2013; Markova et al. 2014; Holgado et al. 2018).

Keeping in mind that the accuracy reached in the estimation of the microturbulent velocity (ξ_t) is not so good when using He lines as when metal lines are considered, we decided to use the outcome of iacob-gbat to investigate whether we find any interesting correlation between ξ_t and log (L/L_⊙). Figure D.2 summarizes the result of this investigation, where we only consider those stars in our working sample having a v sin i < 80 km s⁻¹. This sort of correlation seems to be present, although with a significant scatter.

Fig. D.3 Distribution in the v sin i–ξt diagram of the same set of stars considered in Fig. D.2.

Awaiting a more accurate determination of ξ_t in the stars depicted in Fig. D.2, we performed the following formal exercise to determine the location in the FT (in terms of “equivalent” v sin i) of the zero associated with a given value of ξ_t covering the range shown in Fig. D.2. To this aim, we calculated a grid of fastwind models including H, He, and O as explicit elements. The surface gravity was fixed to log g = 3.8 dex, and we covered a range in Teff from 25 to 50 kK. For each of the models, we computed a synthetic spectrum assuming two values of ξ_t in the formal solution (5 and 20 km s⁻¹), and each spectrum was convolved to a resolving power R=50000 and a v sin i = 5 km s⁻¹. We then used iacob-broad to perform the analysis of all the OIIλ5592Å synthetic lines resulting from the fastwind computations. Those cases in which the estimated v sin i(FT) was larger than 5 km s⁻¹ indicated that the first zero of the FT was actually associated with microturbulence and not with rotation (see also the outcome of a similar exercise in Simón-Díaz & Herrero 2014).

Interestingly, we found that, while in the case of the O iii synthetic lines with ξ_t = 5 km s⁻¹, the values of v sin i(FT) that were recovered were close to 10 km s⁻¹, much larger values (up to 30 km s⁻¹) were obtained when considering the synthetic profiles computed with a ξ_t = 20 km s⁻¹. In addition, we found that the resulting v sin i(GOF) estimates provided by iacob-broad (assuming a δ-function as intrinsic profile) were very similar to those corresponding to the erroneously estimated v sin i(FT) values. Therefore, this simple academic exercise served to illustrate that classical microturbulence could explain both the observed tendency in Fig. 7, as well as the good agreement found between v sin i(FT) and v sin i (GOF) in Fig. 2 in the low- v sin i range, even for cases in which v sin i is most likely affected by methodological limitations.

To reinforce this statement, in Fig. D.3 we present again the distribution of LS and SB 1 stars with v sin i < 80 km s⁻¹, but this time in the v sin i–ξ_t diagram. Roughly speaking, there seems to be an increasing tendency in the lower v sin i, which is achieved for increasing values of ξ_t.

In summary, although there have been advances in disentangling the pure effect of rotation from other sources of line-broadening to get more reliable estimates in υ sin i, we are still not getting the correct rotational measurements for very low values of v sin i. The inadequately disentangled effect, which is most likely produced by microturbulence, seems to be hampering our ability to reach actual v sin i estimates below ~40–50 km s⁻¹. We thus concur with the precautions raised by Sundqvist et al. (2013) regarding the interpretation of any statistical distribution of observed rotation rates of massive-star populations. However, this effect most likely does not affect our measurements above ~50–60 km s⁻¹ or the measurement of v sin i in stars with low microturbulent velocities.

Appendix E Impact of bin size

The bin size used throughout this study for the v sin i histograms is 20 km s⁻¹, following the reasoning presented in Sect. 3.1. Figure E.1 shows the effect of selecting different bin sizes, and how distributions could be misinterpreted as a result.

Fig. E.1 v sin i distribution for the whole sample using different bin sizes.

Appendix F Tables

$$L=T_{\text{eff}}^4/g_{\text{corr}},g_{\text{corr}}=g_{\text{sp}}+g_{\text{cent}},g_{\text{cent}}\approx \frac{{\left(V_{\text{rot}}\sin i\right)}^2}{R_{*}}$$

Columns are: (1, 2) Name and spectral type from GOSC, (3, 4) T_eff, and log L (from Holgado et al. 2020). From this work are v sin i and vmac quoted in Columns 5 and 6. (7) is the number of epochs available, and (8) whether the star was clearly identified as LPV/SB1 from the available multi-epoch observations.

$$L=T_{\text{eff}}^4/g_{\text{corr}},g_{\text{corr}}=g_{\text{sp}}+g_{\text{cent}},g_{\text{cent}}\approx \frac{{\left(V_{\text{rot}}\sin i\right)}^2}{R_{*}}$$

References

All Tables

All Figures

Fig. 1 Spectroscopic HR diagram indicating the location of the 285 Galactic O-type stars investigated in this work. Different colors and symbols are used to highlight luminosity classes. Individual uncertainties are included as error bars. Evolutionary tracks and position of the ZAMS from the non-rotating, solar-metallicity models by Ekström et al. (2012), plus the log g = 3.7 dex line of constant gravity, are also included for reference purposes. The gray shadowed area highlights a region of the diagram which is void of stars (see main text and Holgado et al. 2020, for further explanations).

In the text

Fig. 2 Rotational properties of the sample. Top: Histogram of υ sin i values resulting from the iacob-broad analysis of the 285 O-type stars in our working sample, highlighting the diagnostic line considered for the analysis. Bottom: Comparison of the v sin i estimates resulting from the FT and GOF analysis strategies. The diagonal dotted lines indicate differences between both quantities of ±10 and ±20 km s−1.

In the text

	Fig. 3 Histogram (top) and cumulative distribution (bottom) of υ sin i estimates for the 285 O-type stars investigated in this work, separating results for the various luminosity class subgroups (same color code as in Fig. 1). Ranges in v sin i labeled with letters “A” to “F” are used in Fig. 5 and explained in Sect. 3.3.
In the text

	Fig. 4 Comparison of the υ sin i distributions of the likely single (LS) and clearly detected SB1 subsamples, presented as histograms (top and middle panels, respectively) and cumulative distributions (bottom panel). Stars in the LS sample, for which we only have one spectrum, are highlighted in gray. Ranges in υ sin i labeled with letters “A” to “F” are used in Fig. 5 and explained in Sect. 3.3.
In the text

Fig. 5 Location in the sHRD of our working sample of stars as a function of the various v sin i ranges indicated in Figs. 3 and 4 (see also Sect. 3.3). The LS and SB 1 subsamples are highlighted in black and red, respectively, and their relative percentage in each panel is included for reference. Further information about other features presented in this figure can be found in the caption of Fig. 1.

In the text

	Fig. 6 v sin i–Teff diagram including our working sample of stars. The LS and SB1 subsamples are highlighted in black and red, respectively. The v sin i = 35, 150, and 300 km s−1 values are highlighted for reference purposes (see Sect. 3.4).
In the text

	Fig. 7 v sin i–log(L/L⊙) diagram including our working sample of stars. The LS and SB1 subsamples are highlighted in black and red, respetively. The v sin i = 35, 150, and 300 km s−1 values are highlighted for reference purposes (see Sect. 3.4).
In the text

Fig. 8 Rotational properties of our working sample of Galactic O-type stars separated in three groups (from top to bottom), comprising stars covering different ranges in mass and evolutionary status. Left panels: location of each subgroup in the sHRD; see also Fig. 1 for additional information about other features presented in each sHRD. Middle and right panels: Corresponding v sin i distributions for the LS and SB1 samples, respectively, within each group. In each panel, we indicate the v sin i values of 50, 150, and 300 km s−1 with vertical dashed lines, for reference purposes. Similarly to Fig. 4, the stars in the LS sample, for which we only count on one spectrum, are indicated in gray in the corresponding v sin i distributions.

In the text

Fig. 9 Predicted behavior of the equatorial rotational velocity in the single star evolutionary models, with rotation computed by Ekström et al. (2012, left panels) and Brott et al. (2011, right panels). See Sect. 4.1 for an explanation of the selected values of initial rotational velocity in each set of models. Cross symbols represent time steps in each evolutionary tracks, colored accordingly to the corresponding veq. The dashed red line marks the point in each track with a 20% reduction of veq from vini, and the yellow line indicates the log g = 3.7 dex line of constant gravity. Observations are represented by gray points. The shaded gray area highlights the area in the sHRD close to the ZAMS, which is void of observed stars. Top panels: Considered evolutionary tracks in the sHRD, where we also depict the ZAMS and the 20, 32, and 85 M⊙ evolutionary tracks corresponding to Ekström+12 models without rotation to better identify the three regions in the sHRD (marked as “a”, “b”, and “c”, respectively) discussed in Sects. 3.5 and 4.1 (see also Figs. 8 and 10). Bottom panels: Considered evolutionary tracks in the veq − Teff diagram (or v sin i–Teff for the case of the observations).

In the text

Fig. 10 Predicted behavior of the evolutionary tracks computed by Ekström et al. (2012, left panels) and Brott et al. (2011, right panels) in the Oeq-log (L/L⊙) diagram. From top to bottom, the sections of the evolutionary tracks, which are located within the three regions (“a”, “b”, and “c”, respectively) indicated in Fig. 8 and described in Sect. 3.5, are highlighted with blue and green thicker lines. In addition, the gray thick part of the evolutionary tracks indicates that the corresponding section of the track lies inside the region of the sHRD close to the ZAMS where no O stars are detected. The stars in the observed sample that populate the three considered regions of the sHRD are also included in each panel, separately, for comparison (this time using v sin i values instead of veq). As in other figures, we indicate three v sin i values of interest (50, 150, and 300 km s−1) with vertical lines. Size represents log q value.

In the text

Fig. 11 Properties of the supposed youngest stars in the sample. Left: sHRD highlighting a subsample of 118 LS O-type stars with log g > 3.7 dex. The symbol size depends on the value of the ratio v sin i/vcrit, as indicated. Right: υ sin i/vcrit distributions (as normalized histograms) for the indicated subsample, separating those stars located above (upper panel) and below (lower panel) the 32 M⊙ track. The vertical dashed lines mark several values of interest for the v sin i/vcrit ratio (see text for explanation). The vertical dotted line at v sin i/vcrit = 0.06 marks the approximate limit below which the v sin i determination can be affected by methodological limitations in those stars located above the 32 M0 track).

In the text

Fig. A.1 ∆vrad (left) and σvrad (right) distributions for the O-type stars in our initial sample not identified as SB2s and having a v sin i≤ 180 km s−1. We indicate with different colors how the stars were finally classified as C/LPV/SB1 in the stars with more than one spectrum. Horizontal dashed lines indicate the thresholds considered to separate the stars in the three categories. For most stars, both indicators coincide. In case of divergence, the higher variability flag prevails (SB1, then LPV, and then C).

In the text

	Fig. A.2 Similar to Fig. A.1, but for stars with a projected rotational velocity larger than 180 km s−1. We note that, in this case, we consider the ratio of ∆vrad (left) and σvrad (right) with respect to v sin i.
In the text

	Fig. B.1 Normalized v sin i distribution of the 13 peculiar stars in the sample (six Oe and sevenφlus three) magnetic stars, see text. The global (also normalized) v sin i distribution of the full working sample of 285 LS and SB1 stars is depicted in the background for reference.
In the text

	Fig. C.1 v sin i histogram (up) and cumulative histogram (down) for Ramírez-Agudelo et al. (2013) and our sample. It includes the results for our sample of Galactic likely single O-type stars (black).
In the text

	Fig. D.1 Location in the sHRD of a subsample of 174 LS and SB 1 stars with v sin i < 100 km s−1, separating the targets with a vmac lower or higher than 50 km s−1.
In the text

	Fig. D.2 ξt–log (L/L⊙) diagram including only those stars in our working sample with v sin i < 80 km s−1 and ξt available.
In the text

	Fig. D.3 Distribution in the v sin i–ξt diagram of the same set of stars considered in Fig. D.2.
In the text

	Fig. E.1 v sin i distribution for the whole sample using different bin sizes.
In the text