© The Authors 2023

1. Introduction

One decade ago, de Mink et al. (2013, 2014) performed a detailed theoretical evaluation of the impact that binary interaction could have on the spin-rate properties of massive O-type stars (i.e. main sequence stars with masses in the range of ∼20–80 M_⊙). This study was partly motivated by the necessity to provide an explanation for an empirical result already highlighted by Conti & Ebbets (1977), Wolff et al. (1982) and subsequently confirmed by some other authors (see further references in Simón-Díaz & Herrero 2014; Ramírez-Agudelo et al. 2013; Holgado et al. 2022); namely: the spin-rate distribution of any investigated (large) sample of O-type stars (even in different metallicity environments) is characterised by a main component, including stars spinning with projected rotational velocities (v sin i) below ∼150–200 km s⁻¹, and a tail of fast rotators reaching values of v sin i up to 400–600 km s⁻¹. Such a tail of fast rotators normally comprises ∼20–25% of stars in the considered samples.

By performing a specific simulation of a massive binary-star population typical for our Galaxy and assuming continuous star formation, de Mink et al. (2014) found that binary interaction during main sequence evolution could easily explain the existence of the empirically detected tail of fast-rotating O-type stars. In brief, mass (and angular momentum) transfer from the initially more massive star (donor) to the lower mass companion (gainer) is an efficient mechanism to spin up the latter which, under certain circumstances, may become the more massive component of the binary system after Roche lobe overflow (see also e.g. Packet 1981; Pols et al. 1991).

Indeed, by assuming the empirical distributions of mass ratios, orbital periods, and the corresponding binary fraction obtained by Sana et al. (2012) as input for their simulations, de Mink et al. (2013) ended up with a fraction of fast rotators (i.e. assuming they are characterised by having a v sin i > 200 km s⁻¹), which is very similar to the observed one. Furthermore, de Mink et al. provided some predictions regarding the expected type and percentage of binary products populating the tail of fast rotators in the v sin i distribution of O-type stars. This mainly includes mergers and mass gainers orbited by a hot stripped star or a compact degenerate object (neutron star or black hole).

In this context, we should also note that the merger scenario is considered as one of the mechanisms producing magnetic fields in massive stars (e.g. Ferrario et al. 2009; Schneider et al. 2016, 2019). In such magnetic massive stars, rotation is braked fast, so that mergers would lead to slow rotators. Since the debate remains open on the exact outcome of mergers, the merger scenario should still deserve particular attention when investigating the origin of fast rotators.

In addition, some of these fast-rotating O-type stars would be expected to be detected as runaway stars resulting from a disrupted binary after supernova explosion of the initially more massive component of the system (see e.g. Blaauw 1961; Gies & Bolton 1986; Walborn et al. 2014). The latter would imply a complementary or alternative explanation to the dynamical ejection scenario from a stellar cluster to the occurrence of runaway events in the massive star domain (Portegies Zwart et al. 1999).

The empirical confirmation of these theoretical scenarios has important consequences for several topics of modern astrophysics, especially those influenced by our specific knowledge about massive star formation, evolution, and feedback. In these respects, it is important to recall that stellar rotation is known to play an important role in the evolution of high-mass stars (e.g. Maeder & Meynet 2000; Ekström et al. 2012). It not only modifies the evolutionary paths followed by these stars, as well as their lifetimes and final fates, but it has also been proposed to induce the transport of core-processed elements to the stellar surfaces. There is clear evidence that a non-negligible percentage of massive O-type stars rotate at velocities fast enough to be affected by the aforementioned effects. Thus, we need to be careful while using the empirical information compiled for these fast-rotating O-type stars in order to constrain our theories of high-mass star formation and evolution. This is because the theories will be different if the stars have acquired their angular momentum during the star formation process or through any type of binary interaction. For example, if a large percentage O-type stars having a v sin i larger than ∼200 km s⁻¹ are post-interaction binaries, proposed by de Mink et al. (2013, 2014), using these stars to investigate the efficiency of rotational mixing in single star evolution models may lead to erroneous results and conclusions (e.g. Hunter et al. 2009; Cazorla et al. 2017a).

Another point worth mentioning is the exotic nature of the possible companions to fast rotators. Indeed, short-period binary systems that are passing through the common-envelope phase may produce binary black holes or neutron star systems that can be progenitors of gravitational wave events (e.g. Langer et al. 2020). These systems undergo a Roche-lobe overflow resulting in the donor becoming a Helium star while the accretor becomes a fast rotator of the OB-type. Taking into account subsequent evolution of such systems, BH/NS+OB configurations may arise. While some steps towards the detection of BH+OB systems have been taken (e.g. Villaseñor et al. 2021; Mahy et al. 2022; Banyard et al. 2023; Shenar et al. 2022a,b; Janssens et al. 2023), our programme sample is one of the best compilations that can be used in the search of such systems. Finally, we may mention in this context that cooler fast-rotating stars showing Be-type signatures have shown hints of a post-interaction nature (e.g. Bodensteiner et al. 2020; Wang et al. 2021; Klement et al. 2022).

As a continuation of the efforts devoted by the IACOB project (P.I. Simón-Díaz) to provide solid empirical foundations to our knowledge about the physical properties and evolution of massive OB-type stars (see e.g. Simón-Díaz & Herrero 2014; Simón-Díaz et al. 2017; Holgado et al. 2020, 2022), in this paper, we perform a complete and homogeneous search for empirical signatures of binarity in a statistically meaningful sample of several tens of fast-rotating Galactic O-type stars. Our ultimate goal is to evaluate the scenario proposed by de Mink et al. (2013, 2014) to explain the existence of a tail of fast rotators in this stellar domain. To this aim, we use as starting point the results presented in Holgado et al. (2022), where we performed a reassessment of the empirical rotational properties of Galactic massive O-type stars using the results from a detailed analysis of ground-based multi-epoch optical spectra obtained in the framework of the IACOB & OWN surveys (see Sect. 2). The spectroscopic observations considered in Holgado et al. (2022) can now be complemented with an extended multi-epoch spectroscopic dataset (including a minimum of 3–5 epochs for all stars in the investigated sample of fast-rotators) and a set of superb quality data about proper motions and photometric variability recently delivered by the Gaia and TESS space missions.

The paper is organised as follows. In Sect. 2, we introduce the programme sample and the compiled observational data set. Section 3 describes the analysis we have performed (including a radial velocity analysis and Gaia and TESS data) and the literature overview regarding the programme sample. In Sect. 4, we discuss and interpret the obtained results. Section 5 presents the main conclusions of the paper.

2. Sample definition and observations

2.1. Programme sample

To build the sample that forms the basis for the study detailed in the present paper, we benefitted from information presented in Holgado et al. (2022). There, a detailed investigation of the spin-rate properties was performed for a sample of 285 Galactic O-type stars identified as apparently single or single-line spectroscopic binaries (SB1), using high-quality optical spectroscopic data gathered by the IACOB and OWN surveys (last described in Simón-Díaz et al. 2020a and Barbá et al. 2017, respectively).

Holgado et al. (2022) started from an initial sample of 415 Galactic O-type stars and excluded 113 double-line spectroscopic binaries, plus another 17 peculiar stars (i.e. presenting signatures in their spectra that are typically associated with Oe, Wolf-Rayet, or magnetic stars). For the remaining sample of 285 stars, they obtained estimates for the projected rotational velocities and other stellar parameters that can routinely be procured by means of a quantitative spectroscopic analysis. To this aim, they applied for each star the methodology described in Holgado et al. (2018) to the best S/N spectrum (see also Holgado et al. 2020).

The estimated effective temperature (T_eff), the surface gravity (corrected from centrifugal forces, log g_c, see Herrero et al. 1992; Repolust et al. 2004), and v sin i allowed us to locate the sample in the spectroscopic HR diagram (sHRD, Langer & Kudritzki 2014), as well as to build the corresponding global v sin i distribution. As indicated in Sect. 1, and illustrated by the right panel of Fig. 1 (see also Holgado et al. 2022), this v sin i distribution is characterised by a dominant low velocity component and a tail of fast rotators extending up to ∼450 km s⁻¹.

Fig. 1. Location in the spectroscopic HR diagram of the 285 Galactic (likely-single and SB1) O-type stars investigated by Holgado et al. (2022; left) and associated v sin i distribution (right). In both cases, the sample of likely single and SB1 fast-rotating stars (with v sin i > 200 km s−1) are highlighted in red. Non-rotating evolutionary tracks from Ekström et al. (2012) are also included in left panel for reference.

Following de Mink et al. (2013), this tail of fast rotators – especially above 200 km s⁻¹ – is expected to be mostly populated by the evolved post-interaction binaries. To empirically evaluate this hypothesis, we defined as our initial programme sample of stars those targets identified by Holgado et al. (2022) as having a v sin i > 200 km s⁻¹ (see Table C.1). This corresponds to a total of 50 stars, distributed throughout the O-star domain in the sHR diagram as illustrated in the left panel of Fig. 1.

While this boundary in rotational velocity is somewhat arbitrary, as predicted by de Mink et al. (2013), it should allow us to minimise the presence of pre-interacting binary systems in the sample under study. Therefore, except for a small fraction of stars in the v sin i range between 200 and ∼300 km s⁻¹ (which could be short-period binary systems with individual components spun up by tides), the large majority of stars in our working sample should be (again in the context of de Mink’s scenario) mergers (i.e. genuinely single stars) or mass gainers in which the initially less massive star has now become an O-type star accompanied by a post-mass transfer object (i.e. a compact remnant or a stripped star) or a fast-rotating runaway star. We also note that while we have assumed the same boundary in equatorial velocity (v_eq) as in de Mink et al., we are actually using the quantity v sin i and not v_eq to build our working sample of fast rotators. Therefore, there could still be fast-rotating stars in Holgado’s sample that we are missing because their rotational axis has a low inclination angle. This effect will be taken into account in the discussion and interpretation of our results (see e.g. Sects. 4.2 and 4.6).

For the sake of completeness, while we were not able to obtain accurate individual v sin i measurements for the two components of those stars identified as double-line spectroscopic binaries (SB2), we carried out a careful inspection of all available spectra for the 113 SB2 systems in the initial sample to identify those in which at least one of the components could have a v sin i larger than 200 km s⁻¹. This information allowed us to have an estimate of the percentage of such systems with at least one of the two components being a fast rotator (see Sect. 4.4 and Table 2). The four identified systems fulfilling this criterion are quoted at the bottom of Table C.1.

2.2. Multi-epoch spectroscopic observations

As was the case of previous papers of the IACOB series dealing with O-type stars (see e.g. Holgado et al. 2018, 2020, 2022), the bulk of the spectroscopic observations used in this work comes from three high-resolution spectrographs: the FIES instrument (resolving power, R ≈ 46 000 or 25 000, Telting et al. 2014) attached to the 2.56-m Nordic Optical Telescope (NOT), the HERMES spectrograph (R ≈ 85 000, Raskin et al. 2011) attached to the 1.2-m Mercator telescope, and the FEROS spectrograph (R ≈ 48 000, Kaufer et al. 1997) presently installed at the MPG/ESO2.2-m telescope. The former two instruments are both located in the Roque de los Muchachos observatory (La Palma, Spain), while the latter has been operating in La Silla observatory (Chile) since 2002. Detailed information on the collected spectroscopic data set can be found in Table 1 of Holgado et al. (2018).

Our programme sample comprises Galactic O-type stars in both hemispheres covering a range in B magnitude from 2.5 and down to ∼11 mag. Those stars observable from the Canary Island observatory were observed as part of the IACOB survey (P.I. Simón-Díaz) during several observing runs allocated between 2008 and 2016. This survey initially included a minimum of three epochs per star; however, we also devoted several new additional campaigns since 2017 (P.I. Holgado) to increase the number of available epochs for the subsample of 23 fast rotators visible from the Canary Islands’ observatory. These extra observations, which have allowed us to cover a time-span of more than ten years with more intense time coverage between 2017 and 2018, do not only include FIES and HERMES spectra, but also up to ten additional epochs obtained with the SES high-resolution spectrograph (Strassmeier et al. 2004) attached to the STELLA1.2-m robotic telescope operating at Izaña observatory (Tenerife, Spain).

For those stars observable from La Silla, the multi-epoch spectroscopic dataset compiled for this work was mostly obtained in the framework of the OWN spectroscopic survey¹ (P.I. Barbá) between 2007 and 2017. In addition, these observations were complemented with a few spectra downloaded from the ESO-FEROS archive².

The full list of 54 fast-rotating stars comprising our programme sample³ is quoted in Table C.1, where we also indicate the corresponding spectral classifications (as provided in version 4.1 of the Galactic O-star catalog, GOSC, Maíz Apellániz et al. 2013), the number of FIES, HERMES, FEROS, and STELLA spectra initially available, and the total time-span covered by each set of multi-epoch spectra. Generally speaking, after discarding those spectra with low (less than 5) SNR in the 5875 Å region, we have been able to compile a minimum of 4 epochs for almost 90% of the stars in the sample (see Col. 8 of Table C.1), reaching up to 10 epochs in more than 50% of them. In addition, for some bright Northern targets, we could gather more than 25 epochs. Regarding the covered time-span, except for a few cases, we were able to reach a minimum of 3 years and up to 10–15 years for ∼20 stars (see last column of Table C.1).

Thanks to Gaia-EDR3, we also have information about parallaxes for all stars in our working sample (see also Sect. 2.3). By taking as a first approach the median of the geometric distances provided in Bailer-Jones et al. (2021), we see that 51 of the 54 stars are closer than 3 kpc, and only 6 of them have a RUWE value⁴ larger than 1.4 (see Table C.2).

As described in Sect. 3.1.1, having access to this information allowed us to obtain estimates for other stellar parameters (radii, luminosities) beyond what can be obtained by spectroscopic means, namely, projected rotational velocities, effective temperatures, surface gravities, or surface abundances of certain elements.

2.3. Other sources of empirical data considered in this work

In addition to our main multi-epoch spectroscopic observations described in previous section, we also benefitted from various other sources of empirical data including (a) photometric variability as provided by the light curves obtained by the TESS mission (Ricker et al. 2015), (b) the proper motions delivered by the Gaia mission (Gaia Collaboration 2021), and (c) information about possible close-by visual companions as provided by the Gaia mission, as well as other on-ground high spatial resolution surveys such as the AstraLux optical survey (Maíz Apellániz 2010), the Southern Massive Stars at High angular resolution survey (SMASH, Sana et al. 2014), and the Fine Guidance Sensor resolution survey (FGS, Aldoretta et al. 2015).

Lastly, we also searched for additional info in the literature about whether any of the programme sample of fast-rotating stars had been previously identified as a spectroscopic, eclipsing, and/or X-ray binary. In this regard, we paid special attention to some specific studies providing information about the orbital parameters of previously detected binaries (see Sect. 3.2).

3. Gathering empirical information

Table 1 summarises all the empirical information of interest compiled for this paper regarding our programme sample (excluding the four clearly detected SB2 systems), ordered by increasing rotational velocity. In addition to the spectral classification (Col. 2), we quote the projected rotational velocity (Col. 3), the peak-to-peak amplitude of radial velocity variability, as measured from all available spectra per star (Col. 4), and the assigned spectroscopic binarity status after a first visual inspection of the variability of the He Iλ5875 line-profile (Col. 6).

The above-mentioned information – mainly obtained from the spectroscopic data set – is also complemented with other information resulting from the analysis of the additional empirical data mentioned in Sect. 2.3. The latter includes the type of variability detected in the TESS light curves (Col. 7), a determination of whether some close-by visual companions have been detected within 1 and 2 arcmin, respectively (Col. 9), the runaway status resulting from the analysis of the proper motions provided by Gaia and other studies in the literature (Col. 10) and, lastly, the final binary status established for each individual target after revisiting the literature and a more detailed analysis of the available radial velocity curves (Col. 11).

We list the physical parameters of our sample, gathered from Holgado (2019) and associated papers, in Table C.2. We include estimates for the effective temperature (T_eff, Col. 2), the ℒ parameter (defined as /g_c, see Langer & Kudritzki 2014, Col. 3, where g_c is the surface gravity corrected from centrifugal forces), as well as the stellar luminosity and radius (Cols. 4 and 5, respectively), based on the individual Gaia-EDR3 distances (Col. 9) provided by Bailer-Jones et al. (2021). Below, we describe how each of these pieces of empirical information was obtained.

3.1. Empirical information extracted from the spectra

3.1.1. Spectroscopic and fundamental parameters

Throughout this paper, we consistently use the same set of spectroscopic parameters (basically T_eff and log ℒ) determined by Holgado et al. (2020) and later utilised by Holgado et al. (2022). The corresponding estimates and associated uncertainties are indicated in Cols. 2 and 3 of Table C.2. As commented in Sect. 2.1, these data allow us to locate our programme sample of stars in the sHRD (see red dots in the left panel of Fig. 1).

In addition, the above-mentioned (spectroscopic) parameters were also complemented with other fundamental parameters, such as the radii, luminosities, and spectroscopic masses, which could be derived by considering the extinction corrected V magnitudes provided by Maíz Apellániz & Barbá (2018b) and the Gaia-EDR3 distances quoted in Bailer-Jones et al. (2021). All these additional information can be found in Table C.2.

3.1.2. Projected rotational velocities

Regarding the projected rotational velocities (v sin i), instead of directly using the values obtained by Holgado et al. (2022), we decided to repeat the line-broadening analysis but this time using all available FIES, HERMES, and FEROS spectra per star. This way we wanted to investigate what is the associated uncertainty in the derived v sin i resulting from any potential source of spectroscopic variability affecting the line profiles. In addition, those cases in which an important scatter in the time-dependent v sin i measurements is present could also be an indication that the star is actually a double-line spectroscopic binary.

To this aim, we applied again the IACOB-BROAD tool (Simón-Díaz & Herrero 2014) to the same diagnostic line considered in Holgado et al. (2022). As commented in that paper, for the sake of homogeneity, we tried to use the O IIIIλ5591 line in all cases. However, it was only possible for ∼40% of the stars of the sample. For the rest of the stars, in which the O III line appear too weak and shallow due to the high v sin i, we needed to use either the He Iλ5015 line (∼37%) or the He IIλ5411 line (∼23%).

The results of this multi-epoch analysis are presented in Col. 3 of Table 1, where we indicate the mean values and associated standard deviations computed from the goodness-of-fit solutions provided by IACOB-BROAD. These new estimates are then compared with those obtained by Holgado et al. (2022) in Fig. 2. Generally speaking, there is a very good agreement between both determinations and, except for a few (6) stars, the standard deviation resulting from the multi-epoch analysis is not larger than 5–6%.

Fig. 2. Comparison of v sin i estimates obtained by Holgado et al. (2022) and this work. In the case of our measurements, we include as the error bars the standard deviation associated with the multi-epoch line-broadening analysis. Diagonal lines represent the 1:1 relation and the 10% deviation, respectively.

In fact, the obtained standard deviation of the v sin i measurements is better than 10% in all stars but three (HD 152200, BD+60°2522, and HD 5689), where it ranges between 15 and 20%. Interestingly, the TESS light curve of HD 152200 shows reflection modulation variability (see Sect. 3.3.1), likely produced by a deformation of the star, hence affecting the v sin i measurements depending on the considered orbital phase. From visual inspection of the variability of some of the line-profiles of BD+60°2522, we were not completely sure at first whether this star was a two-component spectroscopic binary (see Sect. 3.1.4), something which could explain the measured scatter in the v sin i estimates. Regarding HD 5689, we did not notice any peculiarities; however, we suggest using high deviation in v sin i as a hint for the possible identification of a faint companion or of peculiar photometric variability.

3.1.3. Radial velocities

Deriving radial velocities (RV) in stars with broad line profiles is a challenging task. Significant spectral broadening caused by fast stellar rotation affects the shape of the spectral line and simple Gaussian-Lorentzian fitting of such line profiles is not efficient. In addition, other effects such as line blending, the fact that rotational broadening makes the lines to become much shallower, or the occurrence of line profile variability within the broad line-profiles, hinder the application of such fitting technique. In this case, the use of a cross-correlation technique – being aware that it has some limitations as well, especially when applied to spectra with a limited signal-to-noise ratio – becomes a viable solution. While this technique allows the possibility to consider a spectral window including several lines, its application to a well isolated line is also a valid option.

For this work, we use the cross-correlation technique described in Zucker (2003). We refer the reader to that paper for a detailed description of how a RV estimate and its associated uncertainty can be obtained from the cross-correlation function (CCF) resulting from a specific observed spectrum and an adequate reference spectrum.

In particular, we decided to apply this technique to just one prominent and unblended diagnostic line which is present in all stars of our programme sample: He Iλ5875. This decision was taken after evaluating the possibility of using a larger set of diagnostic lines present in the optical spectra of O-type stars (6 in total). As illustrated in Fig. D.1, the main outcome of this exercise was that the attempt of using other lines beyond He Iλ5875 was neither improving the accuracy in the RV_PP estimations (see below for definition of this quantity) or importantly modifying the results obtained using just the most prominent (and always present) He I line.

Lastly, we are aware that also He Iλ5876 line can be affected by stellar wind in some specific cases, this effect should not significantly affect the shape of the CCF if the wind emission is not very pronounced as expected for our sample, which is mostly composed of late O-type stars. Indeed, one of the advantages of the cross-correlation approach is that the behavior of the CCF does not depend on the relative distribution of the input data, that is, if we cross-correlate spectral lines with different widths, it does not affect the resulting function. This is why we can securely cross-correlate one convoluted line profile template with the rest of the spectra even if the lines change width or shape.

Technically speaking, we started by fitting a rotationally convolved (synthetic) profile to the He Iλ5875 line of the first spectrum in the multi-epoch data set available for a given star (see e.g. the top panel in Fig. 3). Then, by taking this synthetic line profile as our reference, we cross-correlated the rest of available spectra in the same wavelength range, thus obtaining a CCF as depicted in bottom panel of Fig. 3.

Fig. 3. Example of an observed HeIλ5875 line-profile (top panel) and the best-fit synthetic line used as template for the computation of the cross correlation function (CCF, bottom panel). Following Zucker (2003), the radial velocity and its associated uncertainty can be computed from the maximum and the sharpness of the CCF.

As a result, we obtained RV measurements relative to the first spectrum and their associated uncertainties for all available spectra. Based on these values, we calculated the peak-to-peak amplitude of radial velocity variability (RV_PP) following the various steps described below.

First, to make sure that only reliable RV measurements are considered, we calculate the S/N of our diagnostic line for each spectrum in the time-series. This quantity is defined as:

(1)

where min(F₅₈₇₅) is the flux at the core of the line-profile and σ(F_{5855 − 5865}) is the standard deviation of the normalised flux between 5855 and 5865 Å, a spectral window which only includes continuum points. Basically, this quantity indicates how prominent the spectral line is with respect to the continuum noise.

After evaluating several options, we decided to exclude all RV estimates for spectra with a S/N_L < 5. Even with this conservative threshold, most of the initially available spectra endure to the end of this cleaning process, indicating the good quality of our spectroscopic data set. The estimates of S/N_L and overall S/N for all measured spectra as a function of uncertainty in RV_PP are presented in Fig. D.3.

Then, in order to compute a final estimate of RV_PP, we use the pair of RV measurements that maximises the quantity:

(2)

defined in Sana et al. (2013), where i and j refer to two different spectra from the time-series of a given star, and σ_i and σ_j are the individual uncertainties of the radial velocity measurements RV_i and RV_j, respectively. This approach helps to identify RV_PP from the most accurate measurements by taking into account the individual uncertainties of radial velocity measurements. Also, the uncertainty associated with the final estimate of RV_PP is computed as σ(RV_PP) = .

Both quantities are presented in Col. 4 of Table 1, while the number of spectra that are finally considered in computing them are quoted in Col. 7 of Table C.1. Generally speaking, the accuracy of our RV_PP estimates is better than 10 km s⁻¹ in most of the stars in our programme sample.

3.1.4. Visual inspection of line-profile variability

Orbital motion is not the only effect that can led to the detection of large radial velocity variations in O-type stars, especially in the case of fast rotators. As shown elsewhere (e.g. Fullerton et al. 1996; Aerts et al. 2014; Simón-Díaz et al. 2017), the line profiles of most O- and B-type stars are often subject to various types of variability due to, for instance, stellar pulsations or spots. In addition, in the case of the He Iλ5875 line a non-spherically symmetric wind could also produce some variability which could be erroneously interpreted as empirical evidence of the star being a spectroscopic binary (SB) if only the RV measurements are taken into account.

We illustrate this argument in Fig. 4, where we show the detected line-profile variability in six illustrative examples including two unambiguously identified SB1 systems among our working sample, plus another four cases in which, despite a relatively high RV_PP having been measured (reaching up to 48.5 km s⁻¹ in one of the cases⁵), the detected variation in RV is likely produced by intrinsic line-profile variability.

Fig. 4. Eight examples of the type of line-profile variability detected in our sample of 54 fast-rotating O-type stars. In each panel, we indicate the name of the star, its spectral classification, the number of available spectra, the measured v sin i and RVPP, and the visually assigned classification within the following categories: LPV, SB1, LPV/SB1?, LPV/SB2?, and clear SB2 systems. We also include in the various panels the (much narrower) Na Iλ5890 interstellar line that we use as a sanity check to make sure that the heliocentric correction has been adequately applied to each spectrum in the time-series. Stars are sorted by spectral classification.

Taking into account these ideas, we decided to complement the information about the measured RV_PP with the outcome from a visual inspection of the detected variability of the He Iλ5875 line profile in all stars in our programme sample, classifying each star as being part of one of the following subgroups: ‘LPV’, that is, line profile variables ‘SB1’, ‘LPV/SB1?’, and ‘LPV/SB2?’. This information is provided in ‘SB tag’ column of Table 1.

We used the LPV tag to identify stars which are likely single (LS) as we only detected small variations in the shape of the He Iλ5875 line-profiles. In an opposite situation, if a significant shift of the entire spectral line is detected visually, we classified the star as SB1. If the detected line-profile variability is not prominent, but there seems to be an overall shift of the line, we consider the star as an unclear SB1 and marked it as ‘LPV/SB1?’.

In addition to this, we identified five SB2 candidates for which we were not completely sure about their spectroscopic binary nature. Hence, these stars were provisionally labeled as ‘LPV/SB2?’. One interesting example of this type is HD 24912 (ξ Persei). In this case, the line-profile variability is most probably a consequence of the existence of co-rotating bright spots on the surface of the star, as detected by Ramiaramanantsoa et al. (2014), using photometric observations provided by the MOST satellite, and, hence, it is not an SB2 system. We excluded this star from the final statistics of ‘LPV/SB2?’. Further notes on the other four targets can be found in Appendix B, along with a final decision of their revised spectroscopic binary status.

When presenting the final statistics of detected spectroscopic binaries (Sect. 4.4 and Table 2) among our sample, the ‘surviving’ SB2 candidates will be added to the other four clearly detected SB2 systems⁶ included in the initial sample of 415 O-type stars for which we have found that at least one of the two components has a v sin i larger than 200 km s⁻¹. For reference purposes, we also show two illustrative examples of the detected line-profile variability in the case of these double-line spectroscopic binary systems (see Fig. 4).

In this regard, we also provide here some further comments about the strategy we followed to identify SB2 systems. Again, our main diagnostic line for visual detection of double line spectroscopic binaries has been the He Iλ5875 line. This is not only one of the stronger lines in the spectra of O-type stars, but is also a line which remains strong in the full B star domain. This characteristic makes the line perfect to detect any secondary component hidden in the spectrum, even in cases where the v sin i of this second component is large and, hence, the line is greatly diluted. As a consequence, we can state with a high degree of confidence that given the quality of our spectroscopic data set (in terms of resolving power and S/N), this first visual inspection allows us to detect all possible companions contributing down to ∼10–20% of the total flux of the system in the optical range, especially when we have enough epochs and the amplitude of RV of this (fainter) secondary component is larger than ∼70–80% of the v sin i of the primary.

Obviously, there will be certain cases in which this visual inspection will fail, specifically when the amplitude of RV variability of the fainter companion is less than ∼50% of the v sin i of the more luminous star. This case is expected to affect more importantly to stars with larger v sin i (i.e. the fast-rotator domain), thus leading to situations as those described above, where we are not sure if the detected line profile variability is due to any type of intrinsic variability in a single star or the presence of a companion (the LPV/SB2? case). To minimise the number of such cases, we specifically increase the number of available epochs in the sample of fast rotators and also explored more carefully other diagnostics lines which could help us to decide if we have a SB2 system or a LPV case. Also, we explored the possible detection of eclipses in the available TESS light curves (see Sect. 3.2.2), as well as the potential identification of close-by companions using high angular-resolution images (see Sect. 3.3.2) in order to complement the spectroscopic information with the aim of minimizing as much as possible the effect of observational biases in our detection of both single- and double-line spectroscopic binary systems.

Some additional notes about how the above-mentioned possible observational biases could be affecting the resulting statistics of detected spectroscopic binaries in stars in both the faster and slower rotating samples can be found in Sect. 4.4.

3.2. Empirical information extracted from Gaia and TESS data

The Gaia and TESS missions have provided a unique opportunity to have access to very valuable information of interest for our study in a homogeneous (and almost unbiased) way.

On the one hand, the proper motions delivered by Gaia-EDR3 (Gaia Collaboration 2021) can be used to identify runaway stars. Some of these runaways are expected to be produced by the dynamical ejection of the surviving companion in a high-mass binary system after a supernova explosion event. Following de Mink et al. (2013), an important fraction of the O-type stars with v sin i > 200 km s⁻¹ are the mass gainers in binary systems after Roche-lobe overflow of the initially more massive star. In this sense, it is interesting to know, not only which of the stars in our programme sample are identified as a runaway star, but also to compare the percentage of runaways detected in the slow- and fast-rotator samples of O-type stars investigated in Holgado et al. (2022).

On the other hand, the high quality light curves provided by the TESS mission for almost all our programme stars (with the caveats described in Sect. 3.2.2), allow us to search for signatures of hidden companions not detected through our multi-epoch spectroscopy. Also, although such study is out of the scope of this paper, a thorough investigation of the detected photometric variability (by means of standard asteroseismic data analysis techniques) can provide new insights about the evolutionary nature of fast rotators.

3.2.1. Detection of runaway star candidates among Galactic O-type stars using Gaia ED3 proper motions

It has been known since the 1950s that some OB stars move at high speeds through the Galaxy as a consequence of dynamical interactions between three or more bodies in stellar clusters or of supernova explosions in binary systems (Zwicky 1957; Blaauw 1961; Poveda et al. 1967). Some of the stars are ejected with velocities higher than 30 km s⁻¹. Those are called runaway stars (Hoogerwerf et al. 2001) and can be easily found with Gaia astrometry (Maíz Apellániz et al. 2018a). The slightly slower ones defined as walkaway stars are more common but more difficult to identify (Renzo et al. 2019).

In order to identify runaway stars one needs to calculate their 3D velocity with respect to their local standard of rest (LSR). Strictly speaking, one should differentiate between the ejection velocity and the current velocity due to the possible differences between them caused by the different locations in the Galaxy and the effect of the Galactic potential on the trajectory (see Maíz Apellániz et al. 2022 for examples). However, such differences are usually small (especially for recent ejections) and their calculation requires knowledge about the location of the ejection event, something that we do not currently have for most runaway candidates. For that reason, we consider only the current velocities here.

The 3D velocity is calculated from the 2D components of the tangential velocity and the radial velocity. To obtain the tangential velocity (v_t, lsr) of a star we need to know its distance and proper motion. In this domain, Gaia has opened up the door to a revolution in our knowledge. Radial velocities are a different story. Gaia will provide radial velocities for many stars using its radial velocity spectrometer, which operates in the calcium triplet window. Unfortunately, O-stars have few lines in that wavelength range and the most prominent ones there belong to the Paschen series, which are too broad to obtain precise radial velocities. Furthermore, O-stars suffer from their multiplicity that makes many of them spectroscopic binaries, hence requiring multiple epochs to determine their average radial velocities accurately. To complicate matters further, the spectral lines of O-stars are broad and affected by winds, pulsations, and other effects that lead to disagreements between measurements by different authors (Trigueros Páez et al. 2021).

Given the limitations described above, we estimated the number of runaways in the initial sample of 285 O-type stars not identified as SB2 stars by Holgado et al. (see Sect. 2.1) – splitting the information between the two samples of fast and slow rotators, respectively (see Fig. 1) – using only the information available from Gaia-EDR3 data, that is, the parallaxes and proper motions. We applied the astrometric calibration of Maíz Apellániz (2022), which includes a zero point for the parallaxes that depends on magnitude, color, and position in the sky and a correction on the parallax and proper motion uncertainties that depends on magnitude and the correction to proper motions of Cantat-Gaudin & Brandt (2021). Distances are calculated with the prior of Maíz Apellániz (2001, 2005) and the parameters of Maíz Apellániz et al. (2008), as those are the most appropriate for Galactic O-type stars. We use the Galactic rotation curve described in Maíz Apellániz et al. (2022) and the velocity of the Sun with respect to its LSR from Schönrich et al. (2010). Using those parameters, we calculated the velocity of the stars in the plane of the sky with respect to their LSR.

The uncertainties on the calculated velocities depend first on the parallax uncertainties: a large value implies a large uncertainty on the distance and from there on the tangential velocities. They also depend on the distances themselves, as for stars that are far away the subtraction of the LSR velocity may be biased by our assumed Galactic rotation curve. Hence, in order to avoid objects with large velocity uncertainties, we restricted our samples to those objects with relative distance uncertainties (σ_ϖ, ext/ϖ_c, see Maíz Apellániz et al. 2021 for the notation) lower than 1/3 and (average) distances smaller than 4 kpc. With those restrictions, the sample now contains 179 slow rotators (out of 235) and 41 fast rotators (out of 50).

Figure 5 compares the histograms of tangential velocities for the two samples. Overall, 14/41 (34.1%) of the fast rotators have tangential velocities above the runaway threshold of 30 km s⁻¹, while for slow rotators, only 35/179 (19.6%) are above the threshold. Those objects are secure runaways but it is possible that some objects have a radial velocity with respect to their LSR large enough to bring them across the threshold when combined with their tangential velocity. To estimate how many runaways we are missing, we can count how many have tangential velocities in the 20–30 km s⁻¹ range, as those are the most likely candidates to shift status. We label those as ‘possible runaways’ and their numbers are 7/41 (17.1%) for the fast rotators and 18/179 (10.1%) for the slow rotators.

Fig. 5. Histogram (in percentage of the sample) of tangential velocities with respect to their LSR for the slow- and fast-rotator samples. The two dotted vertical lines indicate the critical velocities vt, lsr of 20 km s−1 and 30 km s−1, respectively.

Jumping to a more detailed investigation of our working sample of fast rotators, Col. 10 of Table 1 (under the heading ‘runaway?’) quotes the 26 stars identified as runaways. We note that, in this case, we also take into account previous findings in the literature, not necessarily based on Gaia-EDR3 data (see e.g. Maíz Apellániz et al. 2018a, and references therein), but also reported by other surveys: HD 191423 (Li 2020), HD 117490 (Li & Howarth 2020), and HD 15642 (de Burgos et al. 2020). Interestingly, six out of the nine stars not fulfilling the distance criteria mentioned above are recovered as confirmed runaways, with a couple of them (BD+60°2522 and HD 41161) being tagged as runaways due to the existence of an associated bow shock, despite having tangential velocities in the 20–30 km s⁻¹ range (Green et al. 2019). These additions increases the number of ’bona fide’ runaways among the sample of fast rotators to 52%, namely, a value closer to that obtained from Gaia-EDR3 data, assuming a threshold in tangential velocity of 20 km s⁻¹ (instead of 30 km s⁻¹).

Some implications of these findings regarding the global sample, as well as further notes about the connection between the spectroscopic binary and runaway status of the sample of fast rotators can be found in Sects. 4.3 and 4.4.

3.2.2. Using the TESS light curves to identify eclipsing binaries and B-type hidden companions

The TESS mission (Ricker et al. 2015) is delivering an enormously rich amount of high quality data for the asteroseismic analysis of large samples of stars covering a broad range of masses and evolutionary stages. In addition, the light curves delivered by TESS serve, among other things, for the detection of eclipsing binary systems (not necessarily transiting exoplanets) and the identification of specific variability patterns associated with the presence of spots at the stellar surface, variable winds, disks, or magnetospheres.

Our interest in getting access to the TESS lightcurves for a good fraction of our sample of fast-rotating O-type stars was mainly twofold. On the one hand, we wanted to identify signatures of eclipses in the photometric data, which could be indicating the presence of a companion not necessarily detected via time-series spectroscopy. On the other hand, by analyzing the resulting periodograms, we would be able to detect frequency patterns associated with fainter, less massive companions (i.e. β Cep and SPB type pulsators, see Aerts et al. 2010 for an overview) whose identification would not be straightforward in the spectra.

An example of this latter situation was presented in Burssens et al. (2020, their Fig. 10), where it was shown that the periodogram obtained from the TESS lightcurve of the O7 V star HD 47839 – known to be ∼25 years period SB1 system (Gies et al. 1993, 1997). In addition, one of the longstanding standards for spectral classification clearly shows a high frequency peak (at ∼12.5 d⁻¹), which most likely corresponds to a faint B-type companion and not to the bright O7 V star.

With these two ideas in mind, we first searched and extracted the TESS 2-min short-cadence data from the Mikulski Archive for Space Telescopes (MAST⁷) for all stars in Table 1, if available (20 stars). We retrieve two light curves, the light curve extracted using simple aperture photometry (referred to as SAP) and the pre-conditioned light curve (Pre-search Data Conditioning Simple Aperture Photometry, PDCSAP). The latter has systematics removed that are common to all stars on the same CCD (Jenkins et al. 2016). Nonetheless, the SAP light curve may be preferable for certain stars as the pipeline is not optimised for OB stars.

By means of a visual comparison and based on predicted variability in the O-star regime, we selected the preferable light curve for each star. We additionally inspected the light curve aperture masks using the lightkurve software package (Lightkurve Collaboration 2018) to rule out any contamination by nearby sources, large sector-to-sector mask variations, or the presence of oversaturated pixels. Light curves for which this was the case were removed from further consideration (1 star).

For stars with no available 2-min short cadence data, we extracted 30-min long cadence data using the lightkurve software package. If available (21 additional stars), we performed simple aperture photometry using a watershed method. That is, we included pixels with a light contribution within 8σ in flux of the light contribution of the central pixel of the source of interest. The extracted light curve was then detrended using principal component analysis, following Garcia et al. (2022). Again, problematic light curves were removed from further consideration (HD 216532).

The extraction procedure yielded 39 TESS light curves, all of which show some form of variability. This includes stochastic low-frequency variability, coherent pulsation modes, and rotational modulation and eclipses, which are in line with general findings for the O-star regime (Pedersen et al. 2019; Burssens et al. 2020). We show three examples of typical pixel maps, light curves and periodograms in Fig. 6. A detailed discussion about the photometric variability for each interesting target is presented in Appendix A. All information regarding photometric variability is presented in Table 1 (column ‘Phot. var.’).

Fig. 6. Three examples of the typical TESS pixel maps, light curves, and periodograms for fast-rotating O-stars. In the left panels, we show the TESS pixel masks, with the pipeline mask in red and Gaia sources marked by white circles. The markers are scaled logarithmically with the Gaia magnitude of the source. In the middle panels, we show the one (or two) sector(s) TESS light curves extracted from the TESS pixel maps. In the right panels, we show the periodograms of the light curves. Significant frequencies, satisfying a signal to noise of 5 in a 1 d−1 window, are given by blue triangles. In the case of BD+36°4145 the highest amplitude frequency was not significant, and we therefore mark it by a green cross. The rotational modulation frequency range is given in gray. This range is estimated using the measured v sin i, an upper limit of veq = 450 km s−1 to account for the inclination and the stellar radius inferred from Gaia-EDR3 (cf. Table C.2).

3.3. Compiling extra information from the literature

In addition to our own spectroscopic, photometric and astrometric analysis, we also performed a careful search in the literature for extra relevant information about our programme sample of stars. In particular, we wanted to know if any of the stars in our sample had been previously identified as a spectroscopic binary, as well as to compile any type of orbital and dynamical information resulting from any existing (more detailed) study of specific targets. Also, we gathered published information about close-by companions as resulting from high angular resolution surveys. Lastly, we looked for papers investigating whether any X-ray emission and/or magnetic feature had been identified among our sample of fast rotators.

3.3.1. Spectroscopic binaries

An important fraction of the stars considered in this work have been previously studied elsewhere. For example, 25 of them were included in the investigation of chemical abundances in fast-rotating massive stars by Cazorla et al. (2017b). The authors also searched for spectroscopic signatures of binarity in their sample. In addition, Trigueros Páez et al. (2021), Mahy et al. (2022) have studied several targets in common with our sample. A detailed discussion about the stars in common and some further notes about how this cross-match between our results and those obtained by previous authors have helpped us to fine-tune and/or reinforce our classification of stars in the sample of fast rotators between SB1, LPV, LPV/SB1? and LPV/SB2? (see column ‘SB tag’ in Table 1) is presented in Appendix B.

3.3.2. Hunting for visual companions using high angular-resolution images

In order to check for any possible contamination of TESS data from other visually close stars to our programme targets, we searched for visual companions in different photometric surveys.

Taking into account that during the extraction of flux from the TESS full frame images, we applied a threshold mask, at the end we collected the flux from a different number of CCD pixels for each star (see Sect. 3.2.2). During the flux extraction, typically we chose the pixel with the greatest light contribution and select all pixels with light contributions within 8σ of flux. Thus we collected the flux, typically between four and five pixels; however, in some crowded areas, the final mask of pixels consists of a larger amount of pixels. The size of TESS pixel is 21 arcsec, then we should take into account any source contamination within 2 arcmin (a bit more than five pixels).

This search for any possible contaminants was initiated with Gaia-EDR3 data. We identified all stars within 2 arcmin that demonstrate a difference in G magnitude smaller than 3 mag. Then, we complemented the available information in Gaia by performing an additional search in the Washington Double Star Catalog (WDSC, Mason et al. 2001), which provides updated information about any identified companion in the literature. We also included results form specific high angular resolution surveys such as SMASH and Astralux, among others.

Generally speaking, we consider that if there is a nearby star with a difference in magnitude of less than 3 mag, the TESS lightcurves for a given target can be contaminated by the flux of a companion, which is not necessary gravitationally bound. This identification of the visual neighbours serves as a warning when extracting and interpreting the TESS photometric data. We note, however, that we also performed a thorough identification of potential contaminants which could be avoided when deciding on the final pixel mask used to extract the TESS light curve for each individual investigated target (see Sect. 3.2.2).

In the ‘Potential contaminants’ column of Table 1, we indicated the outcome of this search. Namely, we indicated number of companions found within 1 arcmin, and between 1 and 2 arcmin which are separated by the symbol ‘+’. More detailed information about the cases in which at least one visual component has been identified is summarised Table C.3, including their angular distances, difference in magnitude, and the corresponding reference from which the information has been extracted.

3.3.3. Incidence of X-ray emission

Close compact companions to massive stars may lead to the emission of hard and bright X-rays, as demonstrated in X-ray binaries (e.g. Reig 2011). Amongst our targets, 17 objects (including four SB1 systems) have been detected at X-ray wavelengths. For 14 targets, the log(L_X/L_BOL) estimates are available and have been published in various works. The corresponding flux ratio is listed in our Table 1 (see X-ray column). For the following targets, the measurements of X-ray fluxes are available in the literature, however, the log(L_X/L_BOL) estimates are not: HD 15137 (McSwain et al. 2010), HD 149452 (Fornasini et al. 2014) HD 46485 (Wang et al. 2007). Two additional ones, HD 76556 and HD 228841, are listed in the 4XMM source catalog. Their spectra were downloaded⁸ and analyzed within Xspec. All X-ray sources have a soft emission with log(L_X/L_BOL)∼ − 7, that is, their X-ray emission can fully be explained by the usual embedded wind-shocks of massive stars. There is therefore no indication for the presence of an accreting compact companion in any of our targets, nor of X-ray bright colliding winds arising in massive binary systems.

3.3.4. Incidence of magnetism

Another parameter that can be useful for a characterization of our sample is the presence of a magnetic field. The vast majority of O-type stars are known to be non-magnetic, but ∼7% of O-stars display strong, dipolar magnetic fields (Grunhut et al. 2017; Petit et al. 2019). Several theories have been proposed to explain this magnetism. In particular, Ferrario et al. (2009), Schneider et al. (2016) suggested that large-scale magnetism among massive stars originates in mergers. Because of the presence of a magnetic field in one component of Plaskett’s star, close binary interactions were also thought to be possible generators, but this assumption was discarded based on actual observations (Nazé et al. 2017).

Nevertheless, such channels (merging, binary interactions) may produce fast-rotating stars as an end product; thus, we decided to check available literature for the existence of magnetic field in our targets. Within the MiMeS survey (Grunhut et al. 2017; Petit et al. 2019), the Stokes V/I profiles have been modelled for seven stars from our sample: HD 210839, HD 36879, HD 24912, HD 192281, HD 203064, HD 46056A, and HD 149757 – and none of these stars were found to be magnetic. Given the small number of studied stars, we could not draw any conclusion regarding the correlation between the existence of a magnetic field with binary or runaway statuses of our fast rotators.

4. Results and discussion

4.1. Preliminary considerations

In this section, we use the empirical information compiled in Tables 1 and C.2 to evaluate the validity of the binary interaction scenario to explain the existence of a tail of fast rotators among Galactic O-type stars. To this aim, we compare some of the global properties of this working sample with those extracted from a complementary sample of stars with v sin i < 200 km s⁻¹.

Starting with the original sample of Galactic O-type stars investigated by Holgado et al. (2020; 2022, see also Sect. 2.1 and Fig. D.4) and guided by the objective of trying to minimise as much as possible the potential effect of observational biases when performing this comparison, we decided to exclude all those stars fainter than V = 10 mag and/or located at a distance from the Sun larger than 3 kpc. After considering this filter – which was not only applied to the initial sample of 285 LS and SB1 stars, but also to the 113 SB2 systems detected by Holgado et al. –, we ended up with 179 and 47 LS or SB1 stars⁹ having a v sin i below or above 200 km s⁻¹, respectively, plus 93 SB2 systems (four of them having at least one of the components with v sin i > 200 km s⁻¹). Hereinafter, we refer to them as the slow- and fast-rotating samples, respectively.

4.2. The RV_PP – v sin i diagram of Galactic O-type stars

Figure 7 depicts the distribution of the above-mentioned sample of 226 LS and SB1 stars in a RV_PP vs. v sin i diagram. Both the RV_PP and v sin i measurements for those targets in the slow-rotating sample are directly taken from Holgado (2019) and Holgado et al. (2022). In the case of the fast rotators (located to the right of the vertical grey shadowed band), we updated the estimates of these two quantities following the guidelines presented in Sects. 3.1.2 and 3.1.3. As a result, we are also able to indicate the associated uncertainties in both v sin i and RV_PP for the fast rotators.

Fig. 7. Distribution in the RVPP vs. v sin i diagram of the sample of 226 LS and SB1 O-type stars described in Sect. 4.1. The vertical gray area represents the artificial threshold between slow- and fast-rotating stars (180–200 km s−1). The horizontal gray area represents the RVPP threshold (20–30 km s−1) in the slow-rotating domain below which the identification of SB1 systems is more difficult due to the effect of intrinsic variability. We note that the small sample of 22 stars for which we only have one spectrum available (i.e. RVPP = 0) is also included for completeness. See Sect. 4.2 for the meaning of the inclined red and gray lines and other details. We also note that while we do not include results in this figure for the additional sample of 93 SB2 systems fulfilling the criteria indicated in Sect. 4.1, the corresponding number and percentage of stars of this type in both the fast- and slow-rotating domain can be found in Table 2.

Except for those stars with v sin i < 200 km s⁻¹ and RV_PP < 20 km s⁻¹ (see reasoning below), wed use different colored symbols to identify the spectroscopic and eclipsing binaries, the periodic photometric modulation variables (including ellipsoidal and reflection modulations, labelled as ‘periodic phot. mod.’), and those stars labelled as runaways. In the case of the fast-rotating sample, we take this information directly from Table 1. For the sake of the interest of the discussion, we have also performed a quick evaluation of the above-mentioned characteristics – following a similar strategy as in our main working sample – in the group of stars located in the upper left region of the diagram.

Red inclined lines depict the regions where RV_PP is 5%, 10%, 15% of the v sin i value in the fast-rotating domain. These lines have been drawn to illustrate how intrinsic variability can lead to values of RV_PP up to ∼15% of the corresponding v sin i. This spectroscopic variability does not only hamper the identification of low amplitude single-lined spectroscopic binaries among fast-rotating O-type stars¹⁰ – especially when only a few number of epochs is available – but can also lead to a spurious identification of SB1 systems among these stars (see e.g. the case of HD 203064 in Sect. 3.3). However, based on our analysis we can suggest to use the threshold of RV_PP > 0.15 v sin i to separate the clear binaries from apparently single stars in the fast-rotating domain.

In the same vein, and following the guidelines presented in Simón-Díaz et al. (2020b), we depict a horizontal grey shadowed band at RV_PP ∼ 20–30 km s⁻¹. This approximate threshold separates those targets below v sin i = 200 km s⁻¹ which can be clearly marked as SB1 from those whose detected line-profile variability could originate from intrinsic variability. In particular, we mark all those stars with RV_PP > 30 km s⁻¹ as SB1, but only those clearly detected as SB1 among the v sin i < 200 km s⁻¹ and 20 < RV_PP < 30 km s⁻¹ sample after performing a more thorough evaluation of their spectroscopic binary status. The remaining targets (with RV_PP < 20 km s⁻¹) are excluded from a similar study since most of them do not have enough number of epochs available for a reliable investigation of the spectroscopic variability.

Along this line of argument and taking into account the fact that the measured RV_PP is expected to importantly depend on the number of available spectra when this number is low, we highlight stars for which we count on less than three spectra using empty circles. As depicted in Fig. 7, thanks to the specific observational efforts motivated by this study (see Sect. 2.2), the new compiled spectra have made us possible to avoid the low number of epochs caveat for a large percentage of stars in the fast-rotating sample (even reaching more than 5 epochs in most of them). However, as mentioned above, this is not the case for the v sin i < 200 km s⁻¹ sample. Therefore, in our working sample of fast rotators, we are quite confident that the measured RV_PP is not going to importantly change by increasing the number of epochs and the percentage of identified SB1 systems will remain basically unaffected. On the contrary, in the low v sin i sample, those targets with less than five epochs (or even three) and RV_PP estimates below 20 km s⁻¹ could actually be unidentified SB1 systems. While this limitation will be taken into account for the discussion presented in Sect. 4.4 (where we compare the multiplicity statistics in the slow- and fast-rotating domains), it is interesting to also remark that there is a non-negligible number of stars (7) in the low v sin i sample with two to four epochs and a measured RV_PP larger than 30 km s⁻¹. Therefore, this fact seems to indicate that large amplitude SB1 systems are easily detected even with such a low number of epochs.

Grey inclined lines initiated in those data points corresponding to stars with RV_PP > 50 km s⁻¹ represent deprojected values of both velocities assuming an inclination angle ranging from 90° to 50° (and perfect alignment between the spin and orbit axis). Filled gray squares along those lines represent the values at i = 75° and 60°, respectively. This analysis actually shows us that the inclination effect cannot significantly affect the statistics of spectroscopic binaries in the fast-rotating domain.

Globally speaking, in addition to the nine previously mentioned fast-rotating SB1 stars, we identify up to 38 clear single-line spectroscopic binaries with RV_PP > 20 km s⁻¹ among the slow-rotating sample. Among them, and thanks to the availability of TESS data, we identify 3 eclipsing binaries (HD 36486, HD 152590, and BD+60°498), as well as another 4 stars showing clear signatures of ellipsoidal (or reflection) modulation in their light curves (HD226868 – aka Cyg X1 –, HD 12323, HD 53975, and HD 94024)¹¹. In addition, within the full sample of stars with v sin i < 200 km s⁻¹, there are 50 secure runaways (9 of them also tagged as SB1).

Inspection of Fig. 7 allows us to highlight several results of interest about the distribution of the global sample in the RV_PP–v sin i diagram. Firstly, all SB1 systems in the fast-rotating domain have RV_PP ≳ 30 km s⁻¹, and are mainly distributed among two main groups – with high (∼120–170 km s⁻¹) and low (∼30–60 km s⁻¹) RV_PP amplitudes, respectively (see also Fig. 8). Indeed, a similar distribution is found in the slow-rotating sample, where a clear gap in RV_PP is also detected. In addition, we confirm the result previously found by Holgado et al. (2022) that the percentage of SB1 stars with v sin i > 300 km s⁻¹ is very small, and basically zero for stars in the extreme tail of fast rotators. Further notes about the SB1 systems, along with an evaluation of the possible nature of the hidden secondary components (also using information from the literature and the TESS light curves) can be found in Sect. 4.6.

Fig. 8. Histograms (in percentage, with respect to the total number of stars in each subsample) of RVPP estimates for the fast- (top) and slow-rotating (bottom) O-type star samples. Number of SB1 systems and stars for which we have only one or two spectra available are also indicated in each panel in red and grey, respectively.

For completeness, we remind that, in addition to the mentioned difficulties to separate RV variations due to orbital motions in a binary system from intrinsic stellar variability in those stars with RV_PP < 20 km s⁻¹, there are 40 stars in the low v sin i sample for which we have less than three spectra (i.e. ∼15% of the sample). This clearly explain why we do not see a normal distribution in the low RV_PP domain of the corresponding histogram in Fig. 8 (contrarily to the case of the fast-rotating sample). Again, this will be taken into account when discussing the comparison of percentages of detected spectroscopic binaries in both samples.

Interestingly, all but two SB1 systems with RV_PP > 100 km s⁻¹ are identified as eclipsing binaries. One of the stars in this subsample not identified as EB is Cyg X-1 (HD 226868, O9.7 Iab p var, v sin i = 95 km s⁻¹), a well-known binary system (P ∼ 5.59 d) hosting an accreting stellar-mass black hole (see Caballero-Nieves et al. 2009, an references therein); the other one is HD 130298 (O6.5 III(n)(f), v sin i = 167 km s⁻¹, RV_PP = 143 km s⁻¹, P∼14.63 d) recently proposed by Mahy et al. (2022) to host a quiescent stellar-mass black hole.

If we now concentrate on the sample of SB1 stars with v sin i > 150 km s⁻¹ and RV_PP < 100 km s⁻¹, there is only one star detected as eclipsing binary (HD 46485). Taking into account that the higher the measured v sin i, the most likely the binary system is observed at a high inclination (i.e. a configuration which favors the presence of eclipses if the orbital and rotational axes are aligned, as often assumed), those fast-rotating SB1 stars for which TESS do not show any signature of eclipses are potential candidates to host a compact object, as will be further discussed in Sect. 4.6. Other SB1 systems of interest (in both the slow and fast-rotating samples) for which we should be able to provide further information about the evolutionary nature of the hidden companion (see Sect. 4.6) are those stars identified to show ellipsoidal variability in the TESS light curves.

Lastly, most of the LS stars in the fast-rotating domain are runaways (see also Table 2), and only one fast-rotating SB1 star is also identified as runaway. The higher fraction of runaway stars among fast rotators was already highlighted in Sect. 3.2.1, while further insights about this result are presented in the next sections.

4.3. Comments on the global statistics of detected runaways

In Sect. 4.4, we discuss in more detail our findings about the multiplicity and runaway incidence among fast-rotating O-type stars. We also evaluate to what extent the obtained empirical results can be used to confirm or reject the binary evolution scenario proposed to explain the existence of a tail of fast rotators. Prior to this, we consider it of interest to briefly discuss the global statistics of detected runaway candidates among Galactic O-type stars¹² presented in Sect. 3.2.1 (see also Fig. 5).

Two important results can be derived from these global statistics. The first one is that the fraction of runaways we find is very high. Combining the fast and slow samples, there are ∼22% certain runaways and an additional ∼11% possible runaways. Those numbers indicate that up to 1/3 of the population of O-stars in the solar neighborhood may be runaways. These numbers are in rough agreement with the runaway fraction of 27% derived by Tetzlaff et al. (2011), but we note that their sample and methods are quite different. Their sample is dominated by B stars of lower mass than our stars and they used Hipparcos astrometry of much lower precision than that of Gaia-EDR3, so they were only able to assign probabilities to each star. On the other hand, the fraction is significantly higher than the value obtained by Maíz Apellániz et al. (2018a), which was likely a consequence of the conservative nature of their methods.

The second result, of greater relevance for the present study, is that there is a significant difference between fast and slow rotators: for the first type we find that ∼35–50% are runaway stars, while for the second, we find a smaller number of 20–30%. This result confirms the finding of Maíz Apellániz et al. (2018a) that Galactic runaway O-stars rotate significantly faster on average than their non-runaway counterparts and is also in agreement with the recent study by Sana et al. (2022), who found that the runaway population of (presumed) single O-type stars in the 30 Doradus region of the Magellanic Cloud presents a statistically significant overabundance of rapidly-rotating stars, compared to its non-runaway population.

This non-negligible difference in the percentage of runaways between the slow- and fast-rotating samples of Galactic O-type stars is somewhat expected if we assume as valid the proposal that binary interaction plays a dominant role in populating the tail of fast rotators. In this case, an important fraction of fast-rotating runaways (if not all) would be originated by the disruption of a post-interaction binary after the first core-collapse in the system. Indeed, Renzo et al. (2019) predicted that ∼50% of a population of high-mass interacting binaries will become a disrupted binary in which the initially less massive star is still on the main sequence. This percentage has been calculated following the information presented in Fig. 4 of Renzo et al. (2019) study; namely, starting from 78% of binary systems which do not merge and considering that 86% of those are predicted to be disrupted, 75% of them including a high-mass main sequence object after core-collapse of the companion. In constrast, since this ejection mechanism is not expected to be operating so efficiently among the slow-rotating sample, the associated runaways would be more likely produced by dynamical ejections resulting from a multi-body interaction in a dense cluster.

There is, however, one important caveat that must be taken into account in the argumentation above. One of the main outcomes of the extensive numerical study of the evolution of massive binary systems performed by Renzo et al. (2019) is that, despite the large percentage of disrupted binaries resulting from the simulations, only a small fraction of them is predicted to acquire peculiar velocities above 20–30 km s⁻¹ (i.e. becoming a runaway from an empirical point of view). Therefore, if the estimations by Renzo et al. (2019) are correct, only a minor fraction of the detected runaways among the fast-rotating sample would come from the disruption of a binary.

Renzo et al. (2019) claim that this is a robust outcome of their simulations, also indicating that similar findings have been previously found by other authors (e.g. De Donder et al. 1997; Eldridge et al. 2011). However, along the next sections, we will provide some arguments supporting the idea that this theoretical result seems to be in tension with our empirical findings. In particular, if all detected runaways in our sample of Galactic O-type stars would have been produced by dynamical ejections, there would no reason for the significantly larger fraction of runaways found between the slow and fast-rotating samples. Indeed, we note that even if we consider the most extreme runaways (with v_t, lsr > 50 km s⁻¹), the fractions of stars with such tangential velocity are 19% in the fast-rotating domain and 8% in the slow-rotating domain (see Fig. 5). However, this is not the only argument and we present more details in the next section, where information about the detected binary status is also taken into account.

4.4. Multiplicity and runaway incidence amongst fast-rotating O-type stars

Table 2 summarises the number and relative percentage (with respect to the total number of stars in each group) of the various types of identified spectroscopic binaries within the slow- and fast-rotating samples. In addition, we indicate the number and percentage of runaways detected among the fast rotators¹³. Lastly, for reference purposes, we indicate that a preliminary estimation of the runaway fraction of stars in the slow-rotating SB1 and LS subsamples (based on the sample of stars considered in Sect. 3.2.1 and including ‘possible runaways’) results in a similar percentage of ∼30%, respectively.

Concentrating first on the fast-rotating sample, it can be noticed that the sample is mainly dominated by likely single stars (∼70%), with a considerable smaller contribution of the single-lined spectroscopic binaries (∼18%) and only ∼8% comprising the confirmed SB2 systems. These numbers assume that the two stars identified as “LPV/SB2?” are LS stars; however, if these were actually SB2 stars, the percentages of SB2 and LS stars would be slightly modified, but the main results remain valid.

In order to put these results in a wider context, we compare these numbers with the corresponding statistics for the stars with v sin i < 200 km s⁻¹. Before presenting our results and making comparisons with the fast-rotating sample, we must indicate that (as mentioned in Sect. 4.2) there is a sample of 22 slow-rotating stars (8.2% with respect to the total number of stars in the slow-rotating sample) for which we only have 1 spectrum available, plus another 18 (6.7%) with RV_PP measurements based on two spectra. Although these stars will be labelled as LS hereinafter, there is some probability that some of them will be identified as SB1 or even SB2 systems when considering more epochs. Therefore, the percentage of SB1 and LS stars presented in the left column of Table 2 must be considered as lower (SB1 and SB2) and upper limits (LS), respectively. However, these percentages are not expected to vary by more than ∼7% (i.e. half of the sample with 1–2 epochs). In addition, since we have not explored in detail the SB1 status of those stars with v sin i < 200 km s⁻¹ and RV_PP < 20 km s⁻¹, the percentages of SB1 and LS stars among the slow-rotating sample could again be somewhat larger (SB1) and smaller (LS), respectively, than those indicated in Table 2.

Taking all this into account, we can nevertheless extract interesting conclusions from inspection of the statistics presented in Table 2. The percentage of SB2 systems is smaller by a factor ∼4 in the fast rotating sample (with a difference between the relative percentages of the two samples of ∼25%). In contrast, the percentage of LS stars behaves in the opposite direction (with a minimum difference in percentage of ∼20% if we assume that no new SB1 systems in the low v sin i sample will be discovered whenever adding new epochs). Regarding SB1 systems, the percentage appears slightly larger among fast rotators. However, both percentages could become more similar if some SB1 stars are still hidden in our sample of stars with few epochs.

We can speculate that the observed fraction of SB2 systems in the slow-rotating domain is higher because these systems are pre-interaction binaries that have not yet evolved to a post-interaction binary system (including a fast-rotating O-type star orbiting a faint stripped star or a compact object), a merger, or a disrupted single runaway star after a supernova explosion event (e.g. de Mink et al. 2014).

Alternatively, the explanation could be related to the difficulty to detect SB2 systems among stars with broad line profiles, or the fact that the majority of SB2 fast-rotating systems are short-period contact binaries (see for example HD 100213, Penny et al. 2008) and this evolutionary phase is relatively short with respect to the other stages. However, we consider these two latter explanations as less probable than the abovementioned one. On the one hand, the quality and number of epochs of our compiled observations are good enough to have been able to detect the vast majority of SB2 systems (we note that only 2 out of 54 stars in our sample of fast rotators have survived as ‘LPV/SB2?’). On the other hand, we have found four (out of 51 stars, i.e. 8%) clear SB2 systems, all of them having v sin i in the range of what is expected in terms of spin-orbit synchronization due to tides.

Another result that seems to further support the binary evolution scenario is the following. Even if future analyses of high quality data with a more extended multi epoch data set could aid to detect faint companions among the SB1 samples, if we add together the percentages of SB1 and SB2 systems, we end up with 26% and 46% of secured spectroscopic binaries among the fast- and slow-rotating samples, respectively. Even if we add the two fast rotators labeled as LPV/SB2? as spectroscopic binaries, the percentage of this type of systems among fast rotators is still considerably smaller than in the case of stars with v sin i < 200 km s⁻¹. In addition, as commented, some of the stars in our sample with 1 or 2 spectra could be discovered as SB1 systems in future investigations, thus increasing the difference in percentages even further.

Another point of interest refers to the statistic of runaway stars. As commented prevously (in Sects. 3.2.1 and 4.3), 52% of the total sample of fast rotators are identified as runaways. Interestingly, all of them but HD 15137 (O9.5 II-IIIn, v sin i = 283 km s⁻¹, RV_PP = 44, SB1) are likely single stars, a result which, assuming the theoretical scenario proposed by de Mink et al. (2011), would support the dominance of the binary evolution scenario (after first core-collapse and supernova explosion) over the dynamical ejection (from a stellar cluster) channel (e.g. Perets & Šubr 2012). As we discuss earlier in this paper, a significant difference in the observed fraction of runaways in the slow- and fast-rotating domain supports the scenario of disrupted binaries. However, we should note that the mentioned theoretical predictions (including Renzo et al. 2019, i.e. simulations) are sensitive to many initial conditions, including the star formation history rate, etc. Thus, if we do not know at which age we are observing the population of stars (the case of our sample), we cannot directly predict the exact fraction of runaway stars at a given stellar age.

From a different perspective, if we concentrate on the (36) fast-rotating and likely single stars, we can highlight two main groups: (i) the runaways (64% of this subgroup versus 30% in the slow-rotating subgroup of these stars), as noted (see comments below), are expected to be mainly dominated by the products of a binary disruption event; and (ii) the rest, which could either be the end-products of a merger event or binary systems for which we are not able to detect any clear spectroscopic binarity signature of a companion.

In reference to the runaway stars in the sample, there is still a possibility that some of them are fast-rotating single stars ejected from their parental stellar cluster by means of one or several dynamical interaction events during the star formation process (Maíz Apellániz et al. 2022). However, to distinguish the dynamical ejection from the binary supernova scenario for a sample of our stars is a tricky task that will require further work. An example of a star with peculiar kinematic properties is HD 149757, the origin of which is still under debate. It has been suggested that this star is a post-interaction binary (Villamariz & Herrero 2005; Renzo & Götberg 2021). It was also suggested that this star has been ejected after the supernova explosion of its massive companion (Tetzlaff et al. 2010; Kirsten et al. 2015). Also, the runaway status of HD 93521 possibly originates from a close binary interaction, and in addition to this, it is a possible merger product (Gies et al. 2022). Interestingly, we detected a single runaway SB1 system (HD 15137, McSwain et al. 2010) amongst them: the only possibility to explain its existence would then be a dynamical ejection from a triple system in the past. In this context, it is interesting to note that HD 13268 and HD 15642 are runaways and members of the Perseus OB1 association (de Burgos et al. 2020). Thus, there is a possibility that they were dynamically ejected.

Regarding the LS fast rotators, the majority of stars in this group have RV_PP < 30 km s⁻¹, thus if any hidden secondary component would be present, it should have a relatively small mass compared to the O-star. Indeed, we recall that in this rapid rotation regime, we expect to see the stars with the rotational axis almost perpendicular to the line of view. Moreover, considering (zero order approximation¹⁴) that the orbit of fast-rotating star is expected to be perpendicular to the rotation axis, as usually assumed in case of short period or over-contact binaries (e.g. Fabry et al. 2022), this implies that the velocity variations due to orbital motions should be close to their maximum values and would have been detected if the secondary masses were high enough. Therefore, it is only possible to have some hidden low-mass companion amongst them. One example is HD 46485, for which we could not detect any clear spectroscopic signature of binary motion with the available data set (five spectra) but TESS data show eclipses. In addition to this, de Mink et al. (2014) simulated the appearance of various binaries in terms of the orbital velocity of the components. According to these simulations, the systems with K < 10 km s⁻¹ (RV_PP ≲ 20 km s⁻¹) could be either products of binary interaction (i.e. mergers or post-interaction binary systems) or effectively single stars. However, these simulations are valid for the general case of main-sequence objects without taking into account the rotation regime (presented in Sect. 4.7).

4.5. Distribution of the fast-rotating sample in the sHRD

The left panel of Fig. 9 depicts the location in the sHRD of our working sample of fast rotators, also including (for reference) the rest of the O-star sample investigated by Holgado et al. (2022), as well as the zero-age main sequence (ZAMS) line and the non-rotating evolutionary tracks computed by Ekström et al. (2012). Coloured symbols as in Fig. 7 are also used here to identify the various types of binaries as well as the runaway stars among the fast-rotating sample. For comparison, a similar sHRD diagram is depicted in the right panel of Fig. 9, this time including only the slow-rotating sample.

Fig. 9. Location of the investigated sample of fast-rotating stars on the spectroscopic HR diagram with the non-rotating solar metallicity evolutionary tracks and ZAMS from Ekström et al. (2012; left). SB1 systems, runaway stars, and eclipsing binaries are labeled by different symbols in the same way as in Fig. 7. Gray dots represent the sample of slow rotating O-type stars (from Holgado et al. 2022). The same diagram but only showing the slow-rotating sample of O-type stars (right).

The first interesting result (already pointed out by Holgado et al. 2022) is that all (9) fast-rotating SB1 systems are located below the 32 M_⊙ evolutionary track, where a mixture of presumably single stars identified both as runaways (8) as well as targets without peculiar kinematical properties with respect to their local environment (7) can be also found. Additionally, the higher mass sample of fast rotators is basically dominated by LS runaway stars, with only 4 of the 18 stars in this region being detected as non-runaway LS stars. Also remarkable is that no fast rotators are found in the more evolved region of the MS in the higher mass domain.

The situation is completely different for the case of SB1 systems in the slow-rotating sample, where the distribution is more homogeneous along the whole O star domain (although with a somewhat larger relative percentage above the 40 M_⊙ track). In addition, there is a clear concentration of runaway stars (mostly LS, but also some SB1) in the top right region of the diagram, where no fast-rotating stars are found. We should note that in a such comparison between fast- and slow-rotating domains we are dealing with target masses which are distributed uniformly; thus in our analysis and conclusions we do not have any bias toward specific mass domains. If we assume the scenario proposed by de Mink as valid, we could interpret the different distribution of runaways and SB1 systems in both panels in the following way. On the one hand, the fast-rotating higher mass runaways might be rejuvenated, spun-up gainers in a post-interaction binary system already disrupted after supernova explosion of the initially more massive companion. On the other hand, most of the more evolved, slow-rotating higher mass O-type runaways could be single stars which have been dynamically ejected from a cluster during the star formation process. Since there is no rejuvenation effect due to binary interaction, the probability of detecting more evolved O-stars among this latter sample is larger compared to the case of the fast-rotating sample. In the same vein, the SB1 systems also detected as runaways in the slow rotating sample could be interpreted as pre-interacting binaries which have been also ejected from their parental cloud during the star formation process – or as binary systems which previously were part of a triple system and in which one of the components exploded at some point as a supernova, hence producing the anomalous velocity (with respect to the local interstellar medium) of the detected SB1 system.

In summary, roughly speaking, we have detected two dominant sub-groups of fast-rotating stars in the sHRD, which likely represent two different stages of binary evolution: (i) the least luminous SB1 systems and (ii) presumably single luminous fast-rotating runaway stars coming from post-interaction systems. Such dominance of luminous runaway stars actually contradicts the simulations that are expected to have BHs as a companion of these systems, hence producing a smaller kick to the accretor (Sukhbold et al. 2016; Shenar et al. 2022a). To explain this distribution of runaways and SB1 systems will require detailed evolutionary modelling of each of these systems. In addition, these two main sub-groups are complemented with a third one, also mostly concentrated in the lower mass region of the O star domain, corresponding to the LS, non runaway stars, and representing ∼25% of the whole fast-rotating sample. This latter group could correspond to merger products or stars with relatively low mass (dim) post-interaction accretors which do not produce detectable RV variations nor eclipses in the TESS light curves. In this context, it is interesting to note that despite the superb quality of our spectroscopic dataset and the TESS photometry, the detection of small-RV signatures of orbital motion in a potential binary system and/or faint eclipses below RV_PP ∼ 20 – 30 km s⁻¹ and ΔT_p ∼ 5 – 20 mmag is hampered by the ubiquitous presence of intrinsic photometric and spectroscopic variability at the aforementioned level in the whole O-stars domain (see also Simón-Díaz et al. 2020b).

Interestingly, it is unlikely that some of the stars that are located on the upper part of the sHRD have some hidden SB1 components taking into account that almost all of them are classified as runaways. Thus, we do not expect to spot any other potential fast-rotating SB1 systems in that part of the sHRD, as a result, we do not have any observational bias on this sHRD distribution of fast rotators. As we have already pointed out, to explain the observed distribution of SB1 systems in the slow- and fast-rotating domain will require additional work in modelling the population of these systems. However, in any case it should serve as a constraint for further theoretical works that describe the evolution of binary fast-rotating systems.

4.6. Further notes about the properties of SB1 systems and the characterization of hidden companions

Table 3 provides a summary of all information of interest for the nine definitely confirmed fast-rotating SB1 systems¹⁵. In addition to the information obtained by ourselves for the purpose of this paper (top half of the table), we compile additional information about several orbital parameters (period, semi-amplitude of the RV curve, inclination, eccentricity, and mass function; namely, P, K₁, i, e, and f(m), respectively) in the middle part of the table. The information described above is also complemented with some estimates of the surface velocity the star would have if the spin and orbit were synchronised (v_sync), as well as the quantity v sin i/v_crit (were v_crit is the critical rotational velocity), and v_eq, namely, the estimated equatorial velocity of the star assuming the inclination angle quoted in the table. Lastly, at the bottom, we provide two estimations of the individual mass of the hidden companions by assuming that the mass of the primary is given by a) the spectroscopic mass resulting from the estimated surface gravity (corrected for centrifugal forces) and the radius derived taking into account the apparent magnitude of the star, its distance, and extinction, as well as b) an approximated evolutionary mass (25 ± 5 M_⊙), as obtained from the location of the fast-rotating SB1 sample in the sHRD (see left panel in Fig. 9).

Most of the information on the orbital parameters quoted in the table have been extracted from the recent study by Mahy et al. (2022), except for three stars; one of them (HD 152200) is a newly discovered SB1 system, while HD 46485 is confirmed ‘bona fide’ SB1 system and the third one (HD 165246) has been studied in detail by Johnston et al. (2021).

To first give a global overview of the properties of the sample, regarding the amplitude of RV variability and the detection of eclipses and/or ellipsoidal modulation in the TESS light curves we have: (i) three large-amplitude (RV_PP > 100 km s⁻¹) SB1 binaries, all of them being also eclipsing binaries (HD 37737, HD 165246, and HD 52533); (ii) four small amplitude (RV_PP < 50 km s⁻¹) SB1 binaries, including one eclipsing, non-runaway system (HD 46485), two non-runaway targets whose TESS light curves show periodic photometric modulations (HD 308813 and HD 152200), and one runaway star with no signatures of eclipses or ellipsoidal variability (HD 15137); (iii) two intermediate amplitude non-runaway SB1 binaries for which there is no TESS data available (preventing us from evaluating their EB/EV status, HD 163892 and HD 165174).

If we pay attention to the orbital periods of the systems, we find two main groups, one of them including six systems with orbital periods in the range P∼ 4.5 – 8 d (most of them also detected as eclipsing binaries and/or ellipsoidal variables), plus another one comprising three systems with a somewhat longer period (∼20 – 50 d). In this latter group, there is one eclipsing binary, one runaway star, and a third target tagged neither as runaway nor as eclipsing binary or ellipsoidal variable. Interestingly, there does not seem to be any correlation between RV_PP and P within stars in any of these two main groups.

The first conclusion we can extract by inspecting the information compiled in Table 3 (more specifically, by comparing the measured v sin i with the estimated equatorial velocities assuming spin-orbit synchronization) is that it is quite unlikely that the spectroscopically detected components of these nine fast-rotating SB1 systems have acquired their relatively large equatorial velocities through tidal effects, leaving again the binary interaction channel as the most probable one.

According to theory, the dim (hidden) companion of a post-interaction high-mass binary can be either a compact object (e.g. black holes and neutron stars, BH/NS; see Langer et al. 2020; Mahy et al. 2022) or a stripped helium star or subdwarf object (e.g. de Mink et al. 2013; Shao & Li 2014; Götberg et al. 2018).

Among the nine investigated SB1 systems, there are five targets for which the presence of a black hole or neutron star can be discarded. This includes the four stars for which two clear eclipses are detected in the TESS light curve, plus¹⁶ HD 163892, for which Mahy et al. (2022) claim to have identified spectroscopically the secondary by means of a spectral disentangling technique¹⁷. Among the remaining four stars, there are two targets for which TESS data show potential signatures of periodic photometric modulation (HD 308813 and HD 152200), a third one detected as a runaway SB1 system with an orbital period of ∼55 days and a quite large eccentric orbit (HD 15137), and a fourth target (HD 165174, P ∼24 d, e ∼0.16) not detected as runaway, but for which we do not have TESS data available to evaluate its photometric behaviour. All but HD 152200 have been investigated in detail by Mahy et al. (2022), who do not propose any of them as candidates to host a quiescent stellar mass black hole. Additionally, they show that in case the secondary would have a mass exceeding ∼3 M_⊙ (see last 2 rows of Table 3), it would have been detected in the spectrum using their disentangling technique.

HD 308813 displays clear phase-locked modulations in flux and radial velocities, however, they are difficult to interpret especially in view of the low number of spectroscopic observations. Nevertheless, no evidence indicate that this binary system hosts a black hole or neutron star. Further characterisations of this system will require more spectroscopic observations.

Another situation occurs in the case of the two newly confirmed SB1 systems HD 152200 and HD 46485. The shape of their light curve shows the typical behavior of a binary system in which half of the stellar surface of the cooler component is illuminated by the hotter star (reflection modulation; i.e. maximum and minimum flux at conjunction moment). This leaves us with HD 165174 as the only possible fast-rotating O-type star in our sample for which the presence of a quiescent compact object cannot be definitely ruled out (pending the acquisition of new TESS data or high-quality ground-based photometry and further ultraviolet observations).

Lastly and rather importantly, our only fast-rotating SB1 system also detected as a runaway (HD 15137) deserves further attention. This star was investigated in detail by McSwain et al. (2010), finding that it may contain an elusive neutron star, in the ejector regime or a quiescent black hole with conditions unfavorable for accretion at the time the XMM observations used in that study were obtained. Mahy et al. (2022) also included this star in their search for quiescent O+BH/NS systems, estimating the mass of the hidden companion to be in the 1–6 M_⊙ range (with a more probable mass of ∼2.5 M_⊙) and indicating that they could have identified this dimmer companion in the spectra if it would have a mass above ∼3 M_⊙ while the star is not actually detected. Our estimated mass of the secondary also points to this boundary limit which leaves the exact nature of the hidden companion still elusive.

As commented in Sect. 4.2, there are another four stars among the low v sin i sample that deserve further attention in this context. These refers to those SB1 stars with v sin i < 200 km s⁻¹ for which we have detected signatures of ellipsoidal (or reflection, depending on the considered orbital period of the system) variability ion the TESS light curves. One of them is HD 226868 (aka Cyg X-1), a well studied binary system comprising a late-O supergiant orbiting an accreting stellar-mass black hole (see Caballero-Nieves et al. 2009, an references therein). Another two are: HD 12323 (ON9.2 V, v sin i = 121 km s⁻¹, RV_PP = 58 km s⁻¹) and HD 94024 (O8 IV, v sin i = 162 km s⁻¹, RV_PP = 54 km s⁻¹). By combining the available RV and TESS data we can confirm that the period of these binary systems are 1.92 and 2.46 days, respectively, and that the folded light curves correspond to ellipsoidal modulation (as also pointed out by Mahy et al. 2022). Interestingly, these authors claimed as unlikely that any of these two systems host a stellar mass black hole by considering that the higher mass of the hidden companion is obtained when assuming spin-orbit synchronization. However, the situation changes if, instead of tidal synchronization and in view of the detected ellipsoidal modulation, we assume that these stars actually have a much faster equatorial velocity, reaching values on the order of v_eq/v_crit∼0.4 – 0.6, as in the case of the other SB1 system in our fast-rotating sample (see Table 3). In this case, the estimated masses of the hidden companion would be in the range ∼3–6 M_⊙ and 2.5–5 M_⊙, respectively; hence, leaving the door to be considered as binary systems hosting a stellar-mass black hole (or at least a neutron star) still open.

Regarding HD 53975 (O7.5 Vz, v sin i = 179 km s⁻¹, RV_PP = 47 km s⁻¹), the last star in this group, we do not have yet enough epochs to check if the orbital period is ∼12.5 or ∼6.25 days. Therefore, more RV measurements are needed to decide if the photometric variability detected in the TESS light curve is associated with ellipsoidal or reflection modulation, an evaluate the evolutionary nature of the hidden companion.

4.7. Comparison with theoretical predictions

As we have pointed out, based on the observational properties of the investigated sample of fast rotators, it is likely that the vast majority of them can be assumed to be post-interaction binary products. Using the information presented in Table 2, Fig. 10 provides a global overview of the main characteristics of our sample of 51 Galactic O-type stars with v sin i > 200 km s⁻¹. In particular, it highlights the relative percentage of presumably single stars and spectroscopic binaries, as well as of detected eclipsing binaries and runaways in each subsample. As we can see from Fig. 10, only 25% of all fast rotators are likely single stars – and not runaways. These targets could be mergers, following the theoretical predictions regarding the origin of rapid rotators of de Mink et al. (2013). The detailed characterization of fast rotators presented in this work is a crucial testbed for further theoretical studies of stellar evolution. The next step is therefore to compare model predictions with our observational results.

Fig. 10. Overall observed fractions of spectroscopic binaries, runaways (including walkaways) and eclipsing binaries for all working sample of fast rotators within 3 kpc distances (according to Table 2) shown on the left. Simulated fractions of spectroscopic binaries and runaway (incl. walkaway) stars over all population of O-type fast rotators (based on the BPASS models) shown on the right.

Recent simulations of OB binary system configurations by Langer et al. (2020) suggest that fast-rotating systems with a range of orbital velocities from K₁ ∼ 50 km s⁻¹ (i.e. RV_PP ∼ 100 km s⁻¹, case A – donor is on the main sequence) to K₁ ∼ 100 km s⁻¹ (i.e. RV_PP ∼ 200 km s⁻¹, case B – donor is evolving to the Red Giant phase) most likely have black holes as companions (see Fig. 6 in the cited paper). These simulations can be placed in terms of rotational velocity as well. Thus, theoretically, the OB+BH systems have a higher probability to be detected with equatorial rotational velocity (v_eq) ∼ 100 km s⁻¹, K₁ ∼100 km s⁻¹ (slow-rotating domain) and v_eq ∼500 km s⁻¹, K₁ ∼50 km s⁻¹ (fast-rotating domain, C. Schürmann, priv. comm.). As we discuss in Sect. 4.6, we did not detect any SB1 system with v sin i > 350 km s⁻¹. However, we did detect LS runaway stars with v sin i ∼400 km s⁻¹ that could be proven to be disrupted systems following a supernova explosion. From another aspect, according to Fig. 7, the majority of SB1 systems with possible or confirmed BH companions are located in the slow-rotating domain at RV_PP ∼ 100±50 km s⁻¹. These observational results should be taken into account in future theoretical works regarding the prediction of the existence of OB+BH systems.

Moreover, within the present work, we also aimed to evaluate the observational appearance of fast-rotating O-type stars based on recent theoretical predictions. To do so, we have used the synthetic stellar populations from the Binary Population And Spectral Synthesis v2.2.1 results (Eldridge et al. 2017; Stanway & Eldridge 2018) to estimate the expected properties of fast rotators in various O-star subpopulations, i.e. SB1, SB2 and singles. We use the fiducial BPASS population, this includes a mix of single stars and binary stars as described by Moe & Di Stefano (2017), although the majority of O-stars are in binary systems. We also use prescribed initial mass ratio and period distributions from the same source. The applied IMF is Kroupa et al. (1993) with an upper mass limit of 300 M_⊙.

BPASS stellar models calculate the detailed structure and evolution of each star in sequence using a custom version of the Cambridge STARS code. First, it evolves the more massive primary in detail and approximates the secondary evolution using the stellar evolution equations of Hurley et al. (2000). It then calculates the evolution of the secondary in detail either as a single star or a binary with a compact companion. Thus, in making synthetic O-star sub-populations we use the following constraints to put stellar models into each group. We classify an O-star as having T_eff ≥ 30 kK, a luminosity of log(L/L_⊙)≥4.5 and a surface hydrogen mass fraction above 0.2.

In what follows, we use specific definitions of SB1, SB2, and single stars. It assumes that the SBs are all post-interaction products. This is an extreme assumption but it provides important constrain on what binary interaction can produce and, thus, help to interpret the observational results.

Indeed, to accurately predict the number of SB2 stars expected from an observation survey is difficult. It is a complex interplay between the relative luminosity of the two stars but also whether the orbital period is short enough that radial velocity variations can be observed. In the clearest cases of similar luminosities and short periods, a SB2 identification is simple. Actually, we performed the simulations of the line-profile appearance by varying the flux contributions of the secondary component with different v sin i regimes of the primary component (see our simulations in Fig. D.5). It becomes more difficult as the luminosities become more different and the periods become longer (i.e. smaller orbital velocities). In some ways determining if a star is a SB2 or SB1 depends more on the luminosity, while the LS sample will have some wide SB1 binaries. However, the error introduced by not taking all these factors into consideration is on the order of a few % as the majority of our LS arise from unbound binaries and single stars. A true detailed reproduction of the observed population, taking into account all the complex selection effects, is beyond the scope of this paper. In a present work, for our BPASS modelling we used the following definitions. The SB2 O-stars are taken from the initial models when both stars in the binary are O-stars or when the luminosities of the two stars are both within an order of magnitude. The SB1 are all stars in binary systems that do not count as an SB2 system. This includes systems where the secondary star is an O-star and the primary is the post-binary interaction donor star or a compact remnant. It also includes systems with an O star as primary and a secondary being a lower mass main-sequence star with a difference in luminosity greater than an order of magnitude. The single O-stars are either O-stars that were originally single (a small number), post-interaction binaries after a merger event or the initially less massive companions in a binary system that become unbounded after the supernova explosion of the star that was the primary when the system originated.

Current BPASS models do not include a detailed model of rotation. However, rapidly rotating O-stars are identified by tagging stellar models where the accretor increases its mass by more than 5% of its initial mass or those O-stars that have merged with their companion. In both cases, such stars will be expected to be rapidly rotating and be far above the v sin i threshold assumed in our observations (de Mink et al. 2013). While this is an approximate treatment, it allows us to estimate a lower limit to the population of rapid rotating O-stars. For example, rotation will extend the life of the O-star and we did not include models in which there is a lesser degree of mass transfer, but where stars will still be boosted up the extent that they end up being observed as rapidly rotating. However, these are second-order effects and our predictions provide a robust estimate of the importance of mass-transfer as the source of rapid rotation. In varying the parameters, we have chosen to define the synthetic O-star populations, such as the magnitude difference for SB2 or SB1 systems, leading to an uncertainty in our predicted sub-populations on the order of a few percent of the population for each predicted value.

We used two metallicity models appropriate for solar metallicity: Z = 0.014 and 0.020. This allows for variations and uncertainties in the abundances of the observed stellar population. Amongst the modelled O-star populations, the fraction of fast rotators is ∼20–25% among, which we can emphasise the next sub-groups of targets depending on their binary status: (i) LS: ∼70% (among which ∼14% are runaways and ∼11% are walkaways); (ii) SB1: ∼15% (among which ∼9% are runaways, and ∼3% are walkaways); (iii) SB2: ∼15%.

The fraction of runaway stars was calculated by considering the impact of the first supernova. For single star runaways, these are assumed to be O-stars whose binaries have been unbound in the first supernova, while the SB1 runaways are all systems that have experienced a supernova, forming a neutron star and not a black hole.

To calculate the O-star velocities, we modelled the effect of the primary star’s supernova on the binary. We used the formalism of Tauris & Takens (1998) and Tauris et al. (1999), and assuming the neutron star receives a kick from the Hobbs et al. (2005) distribution. The O-star velocities are essentially the pre-supernova orbital velocity of the O-star in cases where stars are unbound. We estimate that the number of systems where a kick will be strong enough to accelerate the accretor up to 30 km s⁻¹ is only 3% (i.e. walkaways) and 9% is the fraction of runaways with velocity more than 30 km s⁻¹ of SB1 systems (or approximately 2% of the entire rapid rotator population). For the single stars, we estimate the number of walkaways to be 11% and 14% of runaways with a velocity of more than 30 km s⁻¹ (in total this is ∼18% of the entire population). As we may notice, according to our simulations and observations, the fast rotators will be more likely to be runaways than walkaways.

Schematically, these results are presented on the right panel of Fig. 10 with the same labels as on the left panel, to ease comparison with the observational results. These values are on the same order as the rates that we observationally inferred, especially if we will take into account that some of the ‘LPV/SB2?’ systems may be shown to be LS objects after all. This good agreement implies that the majority of fast-rotating stars appear to be LS (de Mink et al. 2011) and mass transfer is the primary physical process responsible for creating rapid rotation, as already suggested by de Mink et al. (2013). This is also in line with expectations from Haemmerlé et al. (2017) and Bodensteiner et al. (2020) regarding Be stars, that such a rapid rotation would be difficult to be already present at birth. We note that the prediction here of the number of mergers becoming rapid rotators (i.e. 52%) should be treated as an upper limit. A recent work by Schneider et al. (2019) suggests that merger products may be highly magnetic. This would then lead to a rapid spin down of such objects, suggesting some merger products may, in fact, be slow-rotating (Schneider et al. 2020).

BPASS provides further detail on the nature of SB1 systems. In these models, a SB1 is any binary system for which there is only one O-star in the binary. If both stars are main-sequence stars, this occurs when the primary is ≥0.5 dex more luminous than the secondary since the secondary is not an O-star in this case. In addition to the case above, as SB1 systems BPASS considers any O-type binaries with compact remnants (NS) as secondaries. Thus, BPASS yields: (i) the fraction of SB1 systems with a rapidly rotating star is about ∼6% of all O-type SB1 systems; (ii) for SB1 systems with a slow-rotating O-star, ∼92% are the pre-interacting stars of their systems, ∼5% are systems in which the O-star was initially the secondary star but its status changed due to mass transfer, and ∼3% are O-stars paired with a compact remnant; (iii) for SB1 systems with a fast-rotating O-star, mass transfer has occurred in all cases (from the definition of fast rotators in the models, see above). The fraction of O-star paired with a compact remnant is about ∼84% of the systems, while the fraction of the O-star paired with a post-mass transfer object (e.g. stripped star) is ∼16%.

The illustration of these results is presented in Fig. D.6. A more detailed investigation of the SB1 system’s companions requires the combination of photometric and radial velocity curve property studies in order to reconstruct the orbit of the system that could give a hint about which of the two predicted subgroups of companions we are dealing with.

Given that these BPASS models explain many other varied observations of massive stars and related transients (e.g. Eldridge et al. 2017; Stanway & Eldridge 2018; Massey et al. 2021; Briel et al. 2022), it is positive to see that the same fiducial population can provide estimates that are of the same order as the observed population. However, the agreement is not perfect and, thus, refinements to the models of massive O-stars are required. Especially it concerns the predicted fraction of runaway stars – the predicted runaway population includes runaways caused by supernova only – if the dynamical runaways were included the total runaway population could be significantly increased (see e.g. Eldridge et al. 2011).

5. Conclusions and future prospects

In the present work, we studied a statistically meaningful sample of several tens of Galactic, fast-rotating O-type stars. By analyzing their multiplicity and runaway status, we can draw the following main conclusions:

• The fraction of runaway stars among fast-rotating O-stars is ∼33–50%, which is significantly higher than for slow rotators (∼20–30%). Notably, we have detected several fast-rotating runaways with significant high tangential velocity, that is, v_t, lsr > 50 km s⁻¹.

• The fraction of SB1 systems in the fast-rotating domain is ∼18%. If we assume that they are all post-interaction binaries, the comparison with models indicates that they are products of mass-transfer with a secondary component which can be a compact remnant (black hole or neutron star) or a post-mass transfer object (stripped helium star or subdwarf object).

• The fraction of SB2 systems among O-type fast rotators is about ∼8–12%, which is significantly lower than in the slow-rotating domain ∼33%. Most likely, this statistic indicates how binary systems are evolving from slow-rotating pre-interacting systems to the post-interaction system with rapid rotation. In addition, to support this hypothesis, we detected that the fraction of SB1 systems in the fast-rotating domain (∼18%) is slightly higher than the fraction of these systems in the slow-rotating domain (∼13%).

• We found a couple of intriguing fast-rotating SB1 systems that will require further detailed investigation, namely: HD 152200 and HD 46485, which possibly host a post-interaction low-mass companion (compact remnant or stripped star or subdwarf object), as well as HD 308813, which will require further characterization.

Taking into account the fact that single fast-rotating runaway stars are post-mass transfer systems without a secondary component, whereas SB1 systems are actually post-interaction systems too, we can conclude that the majority of fast-rotating O-type acquired their fast rotation via mass-transfer (∼65%). In addition, some fast rotators are possibly passing by tidal synchronization scenario as the overcontact SB2 systems (∼10%). The rest of the stars could be considered as merger products (e.g. HD 93521), although their exact fraction remains elusive. Indeed, the issue of how to distinguish the population of mergers from effectively single fast rotators (if any) or from fast rotators with hidden companions remains an open question.

The last point that we did not cover in this work is the study of the chemical evolution of fast-rotating stars which will require further investigation. We did not perform a detailed investigation for all fast rotators in our sample, however, helium and nitrogen abundances are available from the literature (Cazorla et al. 2017a; Holgado 2019). Notably, Cazorla et al. (2017a) showed that it is possible to reproduce the stellar atmospheric parameters and abundances for half of their sample of fast-rotating stars by using single-star evolutionary models but others require binary evolution models. In order to perform a detailed investigation of the current evolutionary stages of fast-rotating stars from the chemical composition point-of-view, we need to derive the CNO abundances for the full sample of targets.

Acknowledgments

We thank the referee, Tomer Shenar, for an extremely constructive and respectful report which has helped us to improve the first version of the paper, and also provided interesting ideas for future development. N.B. acknowledges support from the postdoctoral programme (IPD-STEMA) of Liege University. N.B., S.S-D, G.H 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, PID2021-122397NB-C21, 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. This paper is 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 (partly) on data obtained with the STELLA robotic telescopes in Tenerife, an AIP facility jointly operated by AIP and IAC. 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. J. M. A., G. H. and M. P. G. acknowledge support from the Spanish Government Ministerio de Ciencia e Innovación through grant PGC2018-95 049-B-C22. Y.N. acknowledges support from the Fonds National de la Recherche Scientifique (Belgium), the European Space Agency (ESA), and the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Programme (contract linked to XMM-Newton).

References

Appendix A: Overview of TESS results

In this section, we discuss in detail photometric peculiarities of some fast rotators, especially the ones we presented in Fig. 6. In the top row we show the TESS data for the HD 46485 system, recently shown to be eclipsing by Burssens et al. (2020) using cycle 1 TESS data. The eclipsing signal is clearly visible in the light curve and we measure a frequency, f = 0.1453(4) d⁻¹, corresponding to the binary period P = 6.88(2) d. Because of the large TESS pixel size (21 arcsec), and therefore pixel masks, it is important to consider possible background contaminants. A detailed discussion about how we looked for background flux contaminants is presented in Sect. 3.3.2. In the case of HD 46485, possible contaminants in the pixel mask are much fainter (>3 mag) such that HD 46485 is likely the source of the signal. We identified three eclipsing binaries (EB), HD 37737, HD 52533, and HD 46485 (see Table 1, ’var tag’) for which we can confidently exclude a background origin. In addition, detailed investigation of TESS data of HD 152200 reveals the presence of a small-amplitude periodicity of ∼4.5 days. It could be linked to the presence of a compact companion (e.g. ellipsoidal or reflection modulations) or be linked to rotational features. Our further spectroscopic observations were able to clarify the nature of this system. Thus, we classify HD 152200 as a prominent photometric variable with reflection modulations. For some targets the TESS data is not available yet (marked as ’_ ’ in Table 1); however, we indicate if these targets have been identified as EB based on other photometric surveys (e.g. HD 165246; Johnston et al. 2021).

In addition, we looked for signals in the light curves that might reveal a pulsating hidden companion. We demonstrate this by means of a comparison in the second and third row of Fig. 6. The second row shows the O8.5 V(n) star, BD+36°4145, whose light curve and periodogram show low amplitude stochastic variability at low frequencies, typically found in O dwarfs. There are no significant frequencies in the typical SPB and β Cep frequency regimes. By comparison, in the third row, we show the ON9.5 IIn star HD 91651. In addition to low-frequency variability, we detect a signal with f = 7.3893(1) d⁻¹. Before we confirm the presence of a hidden companion, we first need to consider potential background contamination again. While there are about a dozen background sources in the TESS pipeline mask of HD 91651 (see Fig. 6) they are much fainter than the central source. Moreover, a previous study by Pigulski & Pojmański (2008) measured a similar frequency for HD 91651 using 525 ASAS-3 observations over ∼2800 days. The signal is therefore likely intrinsic to the HD 91651 system. However, since HD 91651 is a late O9.5 star, we also need to account for the possibility that the signal may be from the O-star itself and not from a hidden lower mass companion. Indeed, some late O-stars do show β Cep type pulsations, such as HD 46202 (Briquet et al. 2011). The situation for HD 91651 is therefore less clear than in the case of the O7 V star HD 47839 mentioned earlier. We find that all other measurements of frequencies above 4 – 5 d⁻¹ in stars in this sample also occur in O9 – 9.5 stars: HD 89137, HD 28446A, HD 117490, and HD 102415. In all these cases there is therefore the possibility that the signal comes from a lower mass B-type companion but we cannot fully confirm this given the available data. On the other hand, we did not find clear evidence of SPB type pulsations (1 – 3 d⁻¹) in any of the stars in the sample. This is due to the ubiquitous presence of the stochastic low-frequency variability detected in non-eclipsing light curves, making it difficult to disentangle stochastic and coherent modes.

Finally, HD 210839, HD 14434, HD 14442, HD 192281, BD+60°2522, HD 24912, HD 149757, and HD 93521 have been already extensively studied in terms of photometric variability to reveal the nature of non-radial pulsations in On-type stars (see Rauw et al. 2012; Rauw & Nazé 2021). Also, due to significant photometric variations, HD 165174 has been proposed as a candidate β Cep variable (Chini et al. 2012). We confirm the variability in those stars for which a TESS light curve is available.

Appendix B: Overview about individual targets regarding their spectroscopic binary status.

In this section, we discuss in detail the spectroscopic binary classification of individual targets and the literature overview regarding it.

Our first separation between SB1 and LPV stars (based on the measured RV_PP and the visual inspection of the line profiles) is in good agreement with the results presented in Cazorla et al., except for two cases: HD 203064 and HD 52266. While these authors claim that these two stars are SB1 systems with periods P ∼ 5.1 d and ∼ 75.8 d, respectively (see the appendix in the aforementioned article), we identified them as LPV/SB1? and LPV, respectively.

For the case of HD 203064, we count on 41 spectra covering a total time-span of more than 10 years (compared to eight spectra over 6 years in Cazorla et al. 2017b), and with some spectra distributed in blocks of several spectra per night during 4 nights. As illustrated in Fig. D.2, the detected variability is not compatible with a 5.1 d orbit, but more with the expected hourly variability resulting from stellar oscillations. Also, we do not find any clear peak at this period in the periodogram computed from the full dataset. Hence, we can exclude this star to be a SB1 system and modify its status to likely single (LS, see last column of Table 1). We note, in addition, that the same conclusion was reached by Trigueros Páez et al. (2021), see below.

Regarding HD 52266, we did not find a very convincing phase-folded RV curve using the period quoted in Cazorla et al. (2017b). In addition, by performing various periodogram techniques using all the existing radial velocity data for this target, we could not detect any prominent peak in the associated periodogram. Thus, we consider this target as a likely single star with significant radial velocity variations. We will consider it as LS for the purposes of the discussion presented in Sect. 4, but tag it as LPV/SB1? in the last column¹⁸ of Table 1.

The study of spectroscopic binarity among Northern Galactic O-type stars performed by Trigueros Páez et al. (2021) has five stars in common with our sample: HD 229232, HD 52533, HD 192281, HD 37737, and HD 15137. Among these targets, HD 15137 is quoted in that paper as a runaway SB1 system (P ∼ 55.4 d, see also McSwain et al. 2010). Our preliminary assessment of the detected line profile variability led us to propose it as LPV/SB1? despite measuring a RV_PP of 44 km s⁻¹. However, folding our RV measurements with a period of 55.4 d, we found a nice RV curve in agreement with the orbital solution proposed by Trigueros Páez et al. (2021), McSwain et al. (2010), and the recent study by Mahy et al. (2022, see also below). Therefore, we decided to modify its binary status to SB1.

Aiming at the identification of quiescent compact objects in massive Galactic single-lined spectroscopic binaries, Mahy et al. (2022) have investigated 32 O-type stars reported as SB1 in the literature. We have six stars in common, all of them also identified as SB1 by us. Interestingly, for two of these systems, HD 52533 and HD 163892, Mahy et al. could identify a secondary component using spectral disentangling. However, for consistency with our own classifications (since we do not detect any clear signature of the secondary component in our dataset) we keep these stars labelled as SB1, but mark them with a (†) in the last column of Table 1 to highlight the findings by these authors.

We have also two additional stars in common with Williams et al. (2013), namely HD 308813 and HD 229232, both identified as SB1 systems with periods of 6.34 and 6.2 d, respectively. The binary status of HD 229232 has been more recently revised by Trigueros Páez et al. (2021) to a single star status, which we confirm. Regarding HD 308813, we confirm its single-lined spectroscopic binary status (see also Mahy et al. 2022). In the same vein, HD 165246, which we identified as SB1, was also found by Johnston et al. (2021) as an SB1 eclipsing binary system¹⁹ with a period of 4.5927 d. After revisiting the literature, we can highlight that (to our best knowledge) we have identified one clear SB1 system (HD 152200) and another potential SB1 system (HD 46485) previously studied in Burssens et al. (2020).

HD 152200 was already highlighted by Pozo Nuñez et al. (2019) as a possible high-mass eclipsing binary; however, we note that given the photometric variability detected in the TESS light curve, it is more likely that the star is a variable with a reflection modulation (RM). Indeed, the RM nature of variability for this target is confirmed by our on-going FEROS observations of this target, whose preliminary radial velocity measurements seems to indicate that the orbital period is ∼4.5 d (instead of the 8.89 d tentatively suggested by Pozo Nuñez et al. 2019, when interpreting the photometric variability as two eclipses).

HD 46485 is one of those cases in which, despite we did not identified the target as clear SB1 from the compiled radial velocities (RV_PP = 29.5 km s⁻¹), the TESS light curve shows two clear eclipses (see top right panel of Fig. 6 and Burssens et al. 2020). Thus, we modified its spectroscopic binary classification to SB1. This star is presently being studied in more detail in a separated papers using an extended spectroscopic dataset Simón-Díaz et al., in (prep.) and Nazé et al. (2023).

Regarding the four double-line fast-rotating spectroscopic binaries we identified among the whole list of Galactic O-type stars investigated by Holgado (2019), two of them, HD 100213 – O8 V(n)z + B0 V(n) – and HDE 228854 – O6 IVn + O5 Vn, have been previously studied in detail by Penny et al. (2008) and Abdul-Masih et al. (2021), respectively. Both are (eclipsing) over contact binaries – with an orbital period of 1.387 and 1.886 d, respectively – in which the two components fill their Roche lobe. The third object, HD 175514, for which we measured v sin i = 288 and 108 km s⁻¹, respectively, is actually a triple system comprising an inner eclipsing binary (O5.5 V((f)) + B0.5: V with a period of ∼1.62 d; Mayer et al. 2005) and a third outer component – O7.5 IV((f)) – orbiting the other two with a period of at least 50 years (Maíz Apellániz et al. 2019). For the fourth one (HD 165921, O7 V(n)z + B 0:V:) we obtained v sin i = 224 and 151 km s⁻¹, respectively, and did not find any detailed study, nor information about whether it is an eclipsing system, in the literature. A preliminary analysis of our spectroscopic data set for this star²⁰ indicates that the period of this system is ∼1.7442 d, hence another newly discovered over contact binary. Example of line profile variability for some of the SB2 systems is presented in Fig. 4.

Lastly, among five targets that we initially labelled as LPV/SB2?, as already mentioned in Sect. 3.1.4, one of them (HD 24912) was detected by Ramiaramanantsoa et al. (2014) to present co-rotating bright spots, something that might explain the peculiar spectroscopic variability. Thus, we label it as LS in last column of Table 1.

Regarding the other four targets, BD+60°2522 is the main ionizing source of the Bubble Nebula and its classified as O6.5(n)fp, indicating that it has peculiar and variable line-profiles likely produced by rotational modulation of its strong and likely non-spherically symmetric wind (see study of Rauw et al. 2003; Nazé et al. 2021). It is also detected as having a high proper motion with respect to its local interstellar medium pointing towards a possible dynamical ejection after a supernova event. Therefore, we conclude that intrinsic variability is the most probable explanation for the spectroscopic features which could be interpreted as signatures of a secondary companion and, hence, also modify the spectroscopic binary status of this star to LS.

The same decision was taken for HD 175876 after checking a more extended spectroscopic data set (not included here) obtained in the final phases of development of this paper, and taking into account the fact that this star was also independently classified as presumably single by Cazorla et al. (2017b), using a different set of spectra. This decision is also supported by the study presented in Rauw & Nazé (2021), where the authors indicate that the peculiar behavior of the spectral line profiles with strong photometric variability (as the case of this two stars) could be just a consequence of the occurrence of non-radial pulsations.

In contrast, we decided to keep the other two targets (HD 124314 and HD 91651) as ’LPV/SB2?’. On the one hand, HD 124314 was found by Sana et al. (2014, see also Sect. 3.3.2 and Table C.3) to be composed by two objects²¹ separated by less than 2 mas and having a difference in magnitude of ∼1.4 mag. Also, Sota et al. (2014) marked this star as a likely SB2 from a preliminary inspection of some available high resolution spectra, and De Becker & Raucq (2013) identified this system as a colliding-wind binary. On the other hand, we refer the reader to Sect. 3.2.2, were a detailed discussion about the detected variability in the TESS light curve of HD 91651 is presented, indicating that the frequency peak located at ∼ 7.4 d⁻¹ could be produced by a lower mass β Cep companion instead of the star itself – although this statement cannot be definitely confirmed without a proper asteroseismic modelling of the stars and a much larger spectroscopic data set.

Appendix C: Supplementary tables

In this appendix, we present the number of spectra we got for each target, by listing the spectrographs we used as well (Table C.1). In Table C.2, we present the fundamental physical parameters of fast rotators based on Holgado (2019) and Holgado et al. (2020). Table C.3 presents the detailed information about visible components (if any) in a programme sample we found by analyzing Gaia data and WDSC. In Table C.4, we list some information of interest for the SB1 stars in the slow-rotating domain.

Appendix D: Supplementary figures

Figure D.3 presents the quality of RV_PP measurements as a function of v sin i and local S/N around HeIλ5875 line. In Fig. D.4, we present the distance distribution for all O-type stars in the slow- and fast-rotating domain. In Fig. D.6, we illustrate the possible type of second components based on the BPASS model simulations.

Fig. D.1. Example of line profile fitting of HD 228841 for six diagnostic lines of O-type stars (left panel). RVPP estimates based on all presented lines for six stars (sorted by v sin i values) selected as illustrative cases: BD+36°4145 [O8.5 V(n)], BD+60°134 [O5.5 V(n)((f))], HD 228841 [O6.5 Vn((f))], HD 117490 [ON9.5 IIInn], HD 102415 [ON9 IV:nn], and ALS12370 [O6.5 Vnn((f))] (right panel). Note how, for those lines that are too weak, the cross-correlation technique provides individual RV measurements with a larger uncertainty (e.g. the O IIIλ5591 and He Iλ5015 lines in HD 228841, see left panel), or no results at all. While the He Iλ5875 line is always available and provides quite accurate RVPP measurements, this is not the case for the other five diagnostic lines, for which the possibility to use them critically depends on the spectral type and v sin i combination.

Fig. D.2. Small portion of the radial velocity curve of HD 203064 illustrating that the detected variability is not compatible with a 5.1 d orbit, as previously proposed, but more likely is associated with line-profile variability produced by stellar oscillations, wind variability or spots on the surface.

Fig. D.3. Distribution of individual radial velocity errors vs. signal-to-noise of the HeIλ5875 line (S/NL, left panel) and overall signal-to-noise of the spectrum (S/N, right panel) respect to the v sin i (the colour bars on the right) for all analyzed spectra from the working sample. The gray region represents the limit of S/NL < 5.

Fig. D.4. Distance and V magnitude distribution of the 285 LS+SB1 stars and 113 SB2 systems comprising the sample of Galactic O-type stars investigated by Holgado et al. (2020, 2022) Distance estimates from Bailer-Jones et al. (2021). We highlight in grey the ranges in distance and V magnitude of the sample of stars which has been excluded for the discussion presented in Section 4.

Fig. D.5. Simulations of HeIλ5875 line profile variation based on the TLUSTY (Hubeny & Lanz 1995) stellar atmosphere models. We assume that we have two stars with Teff = 32.5 [kK] and logg = 4.0, for the primary we set up two regimes of v sin i, namely 200 and 300 km s−1, without orbital velocity K1 = 0 km s−1 (left and right panel, respectively). In the case of the secondary component, we assume it has ten random K2 values in the range from 0 to 50 km s−1 and constant v sin i of 100 km s−1. Then we modelled the composite spectrum with the different contributions of the secondary component by varying the flux ratio F2/Ftotal from 0.1 to 0.3.

Fig. D.6. Schematic illustration about possible nature of second components of O-type SB1 systems according to the BPASS model simulations. See Sect. 4.7 for details.

All Tables

All Figures

Fig. 1. Location in the spectroscopic HR diagram of the 285 Galactic (likely-single and SB1) O-type stars investigated by Holgado et al. (2022; left) and associated v sin i distribution (right). In both cases, the sample of likely single and SB1 fast-rotating stars (with v sin i > 200 km s−1) are highlighted in red. Non-rotating evolutionary tracks from Ekström et al. (2012) are also included in left panel for reference.

In the text

	Fig. 2. Comparison of v sin i estimates obtained by Holgado et al. (2022) and this work. In the case of our measurements, we include as the error bars the standard deviation associated with the multi-epoch line-broadening analysis. Diagonal lines represent the 1:1 relation and the 10% deviation, respectively.
In the text

	Fig. 3. Example of an observed HeIλ5875 line-profile (top panel) and the best-fit synthetic line used as template for the computation of the cross correlation function (CCF, bottom panel). Following Zucker (2003), the radial velocity and its associated uncertainty can be computed from the maximum and the sharpness of the CCF.
In the text

Fig. 4. Eight examples of the type of line-profile variability detected in our sample of 54 fast-rotating O-type stars. In each panel, we indicate the name of the star, its spectral classification, the number of available spectra, the measured v sin i and RVPP, and the visually assigned classification within the following categories: LPV, SB1, LPV/SB1?, LPV/SB2?, and clear SB2 systems. We also include in the various panels the (much narrower) Na Iλ5890 interstellar line that we use as a sanity check to make sure that the heliocentric correction has been adequately applied to each spectrum in the time-series. Stars are sorted by spectral classification.

In the text

	Fig. 5. Histogram (in percentage of the sample) of tangential velocities with respect to their LSR for the slow- and fast-rotator samples. The two dotted vertical lines indicate the critical velocities vt, lsr of 20 km s−1 and 30 km s−1, respectively.
In the text

Fig. 6. Three examples of the typical TESS pixel maps, light curves, and periodograms for fast-rotating O-stars. In the left panels, we show the TESS pixel masks, with the pipeline mask in red and Gaia sources marked by white circles. The markers are scaled logarithmically with the Gaia magnitude of the source. In the middle panels, we show the one (or two) sector(s) TESS light curves extracted from the TESS pixel maps. In the right panels, we show the periodograms of the light curves. Significant frequencies, satisfying a signal to noise of 5 in a 1 d−1 window, are given by blue triangles. In the case of BD+36°4145 the highest amplitude frequency was not significant, and we therefore mark it by a green cross. The rotational modulation frequency range is given in gray. This range is estimated using the measured v sin i, an upper limit of veq = 450 km s−1 to account for the inclination and the stellar radius inferred from Gaia-EDR3 (cf. Table C.2).

In the text

Fig. 7. Distribution in the RVPP vs. v sin i diagram of the sample of 226 LS and SB1 O-type stars described in Sect. 4.1. The vertical gray area represents the artificial threshold between slow- and fast-rotating stars (180–200 km s−1). The horizontal gray area represents the RVPP threshold (20–30 km s−1) in the slow-rotating domain below which the identification of SB1 systems is more difficult due to the effect of intrinsic variability. We note that the small sample of 22 stars for which we only have one spectrum available (i.e. RVPP = 0) is also included for completeness. See Sect. 4.2 for the meaning of the inclined red and gray lines and other details. We also note that while we do not include results in this figure for the additional sample of 93 SB2 systems fulfilling the criteria indicated in Sect. 4.1, the corresponding number and percentage of stars of this type in both the fast- and slow-rotating domain can be found in Table 2.

In the text

	Fig. 8. Histograms (in percentage, with respect to the total number of stars in each subsample) of RVPP estimates for the fast- (top) and slow-rotating (bottom) O-type star samples. Number of SB1 systems and stars for which we have only one or two spectra available are also indicated in each panel in red and grey, respectively.
In the text

Fig. 9. Location of the investigated sample of fast-rotating stars on the spectroscopic HR diagram with the non-rotating solar metallicity evolutionary tracks and ZAMS from Ekström et al. (2012; left). SB1 systems, runaway stars, and eclipsing binaries are labeled by different symbols in the same way as in Fig. 7. Gray dots represent the sample of slow rotating O-type stars (from Holgado et al. 2022). The same diagram but only showing the slow-rotating sample of O-type stars (right).

In the text

	Fig. 10. Overall observed fractions of spectroscopic binaries, runaways (including walkaways) and eclipsing binaries for all working sample of fast rotators within 3 kpc distances (according to Table 2) shown on the left. Simulated fractions of spectroscopic binaries and runaway (incl. walkaway) stars over all population of O-type fast rotators (based on the BPASS models) shown on the right.
In the text

Fig. D.1. Example of line profile fitting of HD 228841 for six diagnostic lines of O-type stars (left panel). RVPP estimates based on all presented lines for six stars (sorted by v sin i values) selected as illustrative cases: BD+36°4145 [O8.5 V(n)], BD+60°134 [O5.5 V(n)((f))], HD 228841 [O6.5 Vn((f))], HD 117490 [ON9.5 IIInn], HD 102415 [ON9 IV:nn], and ALS12370 [O6.5 Vnn((f))] (right panel). Note how, for those lines that are too weak, the cross-correlation technique provides individual RV measurements with a larger uncertainty (e.g. the O IIIλ5591 and He Iλ5015 lines in HD 228841, see left panel), or no results at all. While the He Iλ5875 line is always available and provides quite accurate RVPP measurements, this is not the case for the other five diagnostic lines, for which the possibility to use them critically depends on the spectral type and v sin i combination.

In the text

	Fig. D.2. Small portion of the radial velocity curve of HD 203064 illustrating that the detected variability is not compatible with a 5.1 d orbit, as previously proposed, but more likely is associated with line-profile variability produced by stellar oscillations, wind variability or spots on the surface.
In the text

	Fig. D.3. Distribution of individual radial velocity errors vs. signal-to-noise of the HeIλ5875 line (S/NL, left panel) and overall signal-to-noise of the spectrum (S/N, right panel) respect to the v sin i (the colour bars on the right) for all analyzed spectra from the working sample. The gray region represents the limit of S/NL < 5.
In the text

	Fig. D.4. Distance and V magnitude distribution of the 285 LS+SB1 stars and 113 SB2 systems comprising the sample of Galactic O-type stars investigated by Holgado et al. (2020, 2022) Distance estimates from Bailer-Jones et al. (2021). We highlight in grey the ranges in distance and V magnitude of the sample of stars which has been excluded for the discussion presented in Section 4.
In the text

Fig. D.5. Simulations of HeIλ5875 line profile variation based on the TLUSTY (Hubeny & Lanz 1995) stellar atmosphere models. We assume that we have two stars with Teff = 32.5 [kK] and logg = 4.0, for the primary we set up two regimes of v sin i, namely 200 and 300 km s−1, without orbital velocity K1 = 0 km s−1 (left and right panel, respectively). In the case of the secondary component, we assume it has ten random K2 values in the range from 0 to 50 km s−1 and constant v sin i of 100 km s−1. Then we modelled the composite spectrum with the different contributions of the secondary component by varying the flux ratio F2/Ftotal from 0.1 to 0.3.

In the text

	Fig. D.6. Schematic illustration about possible nature of second components of O-type SB1 systems according to the BPASS model simulations. See Sect. 4.7 for details.
In the text