© T. Shenar et al. 2022

1 Introduction

The number of LIGO/VIRGO gravitational-wave (GW) detections originating in coalescing neutron stars (NSs) and black holes (BHs) is rapidly growing. The number of merger events is approaching the 100 mark (Abbott et al. 2019; The LIGO Scientific Collaboration 2021), and it is expected to grow by orders of magnitude in the coming years. There is an overwhelming international effort aimed at constraining the pathways and rates of these merger events (see Mandel & Broekgaarden 2022, for a recent review). This endeavour is intricately tied with our understanding of the progenitors of BHs and NSs – massive stars – which are born with masses higher than 8 M_⊙.

A central problem in the field of massive stars concerns the impact and incidence of stellar multiplicity (e.g. Paczyński 1967; Vanbeveren et al. 1998; Langer 2012; Sana et al. 2012). A fundamental input for binary population syntheses in this context is the initial mass ratio distribution of the binary components in massive binaries. In principle, direct measurements of the mass ratios of spectroscopic binaries can only be obtained from double-lined spectroscopic binaries (SB2). However, SB2 binaries are, by definition, systems in which both stars can be observed, and they are hence biased towards mass ratios close to unity. Moreover, SB2 binaries are most easily spotted when the spectral lines are not blended. For this to happen, the binary components need to exhibit appreciable radial-velocity (RV) amplitudes, which again introduces a bias towards short periods and high-mass companions. Short-period (P ≲ 10 d) binaries specifically entail the danger of having undergone an interaction in the past, as testified by short-period Algol-like systems (e.g. Mahy et al. 2020a,b; Janssens et al. 2021; Sen et al. 2022). In contrast, single-lined spectroscopic binaries (SB1) are dominated by binaries with more extreme mass ratios and longer periods, where a past binary interaction can be excluded in the case of two non-degenerate components. They, therefore, provide a vital extension to SB2 binaries for a derivation of the natal mass ratio distributions.

Furthermore, SB1 binaries contain an additional population: that of post-interaction binaries (de Mink et al. 2014). Of special interest are OB+BH binaries. Such binaries represent an intermediate evolutionary phase that appears in the majority of evolutionary sequences leading to the formation of BH+BH mergers (Belczynski et al. 2002; Giacobbo & Mapelli 2018; Kruckow et al. 2018; Marchant et al. 2019). As such, OB+BH binaries provide crucial constraints on binary evolution, supernova physics, and the presence of BH kicks (e.g. Fryer et al. 2012; Vanbeveren et al. 2020; Sen et al. 2021). Langer et al. (2020) predicted that about 2 to 3% of massive stars should host a BH companion, implying that hundreds of such binaries are lurking in Galactic and extragalactic OB-star catalogues. And yet, such binaries are sparsely known, and the ones observed are primarily high-mass X-ray binaries (HMXBs; Corral-Santana et al. 2016). However, the X-ray irradiation, mass transfer, and tidal interaction in X-ray binaries makes their analysis and interpretation complex (e.g. Giménez-Garcia et al. 2016; Hainich et al. 2020; Ramachandran et al. 2022). The strong tidal circularisation in these binaries also blurs information stored in their orbit regarding the ejecta mass and kick of the BH progenitor. Identifying X-ray quiescent O+BH binaries is thus highly desirable. Indeed, the HMXB population is expected to represent just a small fraction of a much larger population of X-ray quiescent binaries hosting a BH (Langer et al. 2020; Sen et al. 2021).

A few reports of quiescent OB+BH systems have recently emerged (e.g. MWC 656: Casares et al. 2014; LB-1: Liu et al. 2019; HR 6819: Rivinius et al. 2020; NGC 1850 BH1: Saracino et al. 2022; NGC 2004#115: Lennon et al. 2022; HD 96670: Gomez & Grindlay 2021). However, most of these discoveries were soon thereafter challenged (e.g. Abdul-Masih et al. 2020; Shenar et al. 2020; Bodensteiner et al. 2020; El-Badry et al. 2022; El-Badry & Burdge 2022; Frost et al. 2022), and the presence of a BH in MWC 656 and HD 96670 still requires validation. In the upcoming year, Gaia is expected to uncover hundreds of OB+BH binaries in our Galaxy relying on accurate astrometry (Breivik et al. 2017; Mashian & Loeb 2017; Janssens et al. 2022). However, to uncover these populations at subsolar metallicities in neighbouring galaxies such as the Small and Large Magellanic Clouds (SMC, LMC), spectroscopy of SB1 binaries is the method of choice. Shenar et al. (2022) identified VFTS 243 in the LMC as an O+BH binary as part of the study put forth here. Recently, a Galactic O+BH system was also reported (HD 130298, Banyard et al. 2022).

Owing to its well constrained distance (d = 49.97 kpc, Pietrzynski et al. 2013), modest reddening (Maíz Apellániz et al. 2014), ongoing star formation, and sub-solar metallicity content (Z ≈ 0.5 Z_⊙), the Tarantula nebula (30 Doradus) in the LMC provides an ideal laboratory for investigating massive binaries. It hosts a rich population of massive OB-type stars (Walborn & Blades 1997; Walborn et al. 2014), including the most massive stars and binaries reported to date (de Koter et al. 1997; Crowther et al. 2010; Hainich et al. 2014; Shenar et al. 2017, 2021; Bestenlehner et al. 2022; Brands et al. 2022). The VLT-FLAMES Tarantula Survey (VFTS, Evans et al. 2011) provided a homogeneous, multi-epoch spectroscopic dataset aimed to establish the stellar and binary properties of the massive-star content in this region. Among the ≈800 targets included in the VFTS, 360 were classified as O-type (Evans et al. 2011; Walborn et al. 2014). Of those, 116 (32%) were found to exhibit large radial velocity (RV) variations (∆RV > 20 km s⁻¹), making them excellent binary candidates (Sana et al. 2013). The Tarantula Massive Binary Monitoring (TMBM, PI: Sana, programme IDs: 090.D-0323, 092.D-0136) campaign acquired 32 additional epochs for the majority of the O-type binary candidates, as well as a handful of evolved B-type stars allocated to spare fibres. In total, 93 binaries were studied by Almeida et al. (2017), resulting in orbital solutions for 31 SB2 and 51 SB1 systems. The SB2s were recently analysed by Mahy et al. (2020a,b).

This paper aims to characterise the hidden companions in a sample of O-type and evolved B-type SB1 binaries in the TMBM sample utilising multi-epoch spectroscopy (Sect. 2). Our goal is twofold: to probe the lower-end of the mass ratio distribution, and to hunt for OB+BH candidates. This is achieved through the implementation of spectral disentangling, described in Sect. 3 and Appendix A. The results are presented in Sect. 4 and discussed in Sect. 5. We provide our final conclusions regarding the mass-ratio distribution and OB+BH population in this sample in Sect. 6.

2 Sample and observations

Our sample comprises the 51 SB1 binaries identified by Almeida et al. (2017) in the TMBM sample. In Table 1, we provide the list of the 51 targets, along with their previous spectral types (SpT), new SpTs derived from the disentangled spectra (Sect. 3), and the orbital periods. We note that Almeida et al. (2017) defined SB1 binaries as systems for which RVs could not be retrieved for a secondary source, which could either mean that the secondary is too faint to be detected, or that severe line blending prevents a trustworthy measurement of the RVs. Hence, for a few targets in the sample, the original SpTs by Almeida et al. (2017) already suggest the presence of a non-degenerate secondary.

Our analysis relies on 32 epochs of optical spectra acquired with the Fibre Large Array Multi Element Spectrograph (FLAMES/GIRAFFE) mounted on UT2 of the Very Large Telescope (VLT). The data and their reduction are described in detail in Almeida et al. (2017); here, we only repeat the essentials. Data were acquired between Oct. 2012 and Mar. 2014 with the L427.2 (LR02) grating, which covers the wavelength range 3964–4567 Å at a resolving power of R = 6400 and rebinned sampling of 0.2 Å. For each target, 32 epochs were secured, with the exception of VFTS 318 and 802, for which only 18 epochs are available. For highly eccentric or long-period (P ≳1.5yr) binaries, we also use five additional epochs from the original VFTS campaign, described in Evans et al. (2011). Table 1 provides the total number of epochs and average signal-to-noise ratios (S/N), which are typically of the order of 50–100 per wavelength element.

In Table 1, we also provide the mass functions f for each systems. Their calculation is primarily based on the results obtained by Almeida et al. (2017), though we use updated orbital parameters obtained in our study when applicable. The estimated flux ratios of the primaries in the visual f₁/f_tot(V) are also given (see Sect. 3.3). Table 1 also lists the effective temperatures T_eff,1 and bolometric luminosities log L₁ of the primaries. Generally, the effective temperature is obtained using SpT calibrations by Doran et al. (2013). For log L₁, we use bolometric corrections from Doran et al. (2013), V-band magnitudes and B–V colours from Evans et al. (2011), and reddening obtained from the difference between the observed colours and SpT-dependent colour (typically B–V₀ ≈ −0.31 ± 0.01, Conti et al. 2008), assuming a relative extinction of R_V = 3.5. We also account for the flux ratio of the primary. We implicitly assume here that the observed colours are not affected by a potential secondary source, which is valid given the early spectral types of the primaries and their companions. For the few targets that were analysed as single stars in previous studies (Bestenlehner et al. 2014; McEvoy et al. 2015; Ramírez-Agudelo et al. 2017), we adopted previously published T_eff and log L unless the secondary is found to contribute more than 20% to the visual flux.

We compute two independent estimates for the mass of the primaries: once based on SpT calibrations (M_1,SpT) and once based on evolution models (M_1,ev). To estimate M_1,SpT, we calculate the mean of the evolutionary masses of all apparently-single stars in the VFTS sample (Schneider et al. 2018) sharing a similar spectral types (i.e. differing by not more than one subtype and luminosity class ordinate). For the evolutionary masses, we use T_eff and log L estimates as input in the BONnSAI tool¹ (Schneider et al. 2014) using tracks from Brott et al. (2011) and Köhler et al. (2015). As priors, we used the Salpter initial mass function and the initial rotational velocity distribution from Ramírez-Agudelo et al. (2013). While some discrepancies between M_1,SpT and (M_1,ev) are observed, they all agree within 1σ, and no obvious trends are present. This is not entirely surprising, since we are comparing evolutionary masses computed by two different means, which suffer similar systematics (unlike comparing spectroscopic or dynamical masses, e.g. Herrero et al. 1992). As M_1,ev suffers from more assumptions than M_1,ev, we adopt M₁ = M_1,SpT for the analysis, but using M_1,ev would not impact our conclusions. From the binary mass function and the estimates on the masses of the primaries, constraints on the mass of the secondary M₂ can be derived, which are also provided in Table 1. The derivation of these constraints is described in Sect. 4. We note that low-mass bloated stripped stars could be contaminating the sample (Irrgang et al. 2020; Shenar et al. 2020), and to rule this out, spectroscopic analyses will be needed. We discuss this in Sect. 5.2.

Finally, Table 1 provides constraints on intrinsic X-ray luminosities obtained with the Chandra Visionary programme T-ReX Crowther et al. (2022). Evidently, only four targets (VFTS 333, 404, 532, 603) have finite detections, all of which smaller than 10³³ erg s⁻¹ and lying below or close to the canonical value L_x/L ≈ −7 (Pallavicini et al. 1981; Sana et al. 2006). This stands in contrast to the only HMXB candidate in the 30 Dor region, VFTS 399 (log L_X = 34.7 [erg s⁻¹], Clark et al. 2015). The few X-ray detections are discussed in Appendix B.

In addition to the spectroscopy, we also make use of I-band photometry obtained between Oct. 2001 and Mar. 2020 with the Optical Gravitational Lensing Experiment data reductions III and IV (OGLE, Udalski 2003; Udalski et al. 2008, 2015).

3 Analysis

3.1 Spectral disentangling

To investigate whether non-degenerate secondaries are present in our sample, we utilise the method of spectral disentangling (Bagnuolo & Gies 1991; Simon & Sturm 1994; Hadrava 1995; Ilijic et al. 2004; Sablowski et al. 2019). By exploiting the information stored in the spectral-line variability of a spectral time series of a binary, disentangling allows one to separate the component spectra in a binary and derive the orbital parameters. There are different algorithms of spectral disentangling, operating in either wavelength space or Fourier space. Here, we utilise the iterative shift-and-add scheme (Marchenko et al. 1998; González & Levato 2006; Mahy et al. 2010; Shenar et al. 2017). This technique generally yields comparable results to Fourier disentangling (e.g. Shenar et al. 2020), and, in some cases, appears to be more robust (e.g. Bodensteiner et al. 2020). The method of spectral disentangling is not limited to a specific types of stars: it has been implemented on O-type stars (e.g. Martins et al. 2012; Mahy et al. 2020a), B and Be stars (Hensberge et al. 2000; Bisikalo et al. 2000; Saad et al. 2005), lower-mass stars (e.g. Maceroni et al. 2014), and even Wolf-Rayet stars (e.g. Shenar et al. 2017, 2018, 2019, 2021). The underlying assumption is that the variability of the spectral lines is dominated by Doppler shifts of one or more sources (see also Appendix A.)

The spectral separation procedure assumes knowledge of the RVs of both components, RV_i and RV_2,i, where i ∈ 1, N runs over all available epochs of observation. Equivalently, the orbital parameters can be assumed, from which the RVs may be computed. The procedure is iterative. Starting with an approximation for the secondary’s spectrum at the jth iteration, B_j, the approximation for the primary’s spectrum at the jth iteration, A_j is computed by shifting-and-subtracting B_j from all observations according to the available RVs, and co-adding the residual spectra in the frame of reference of the primary. The same procedure is then used to compute B_j+1 from A_j. The first iteration assumes a flat spectrum for the (fainter) secondary, unless otherwise stated, such that the first approximation for the primary is simply the co-added spectrum. The convergence rate of the separation algorithm and hence the number of iterations depends on the RV amplitudes and line profiles, and may range between a few to hundreds, depending on the line profiles, RV amplitudes, and desired level of accuracy.

In many applications, however, the RVs of one or both components are not constrained. In this case, one can perform the disentangling across the K₁, K₂ axes (grid disentangling), where K₁, K₂ are the RV semi-amplitudes of both components. In this context, specific spectral lines (or sets of lines) are disentangled. By computing χ²(K₁, K₂) from the residual spectra: (1)

one can minimise χ² to retrieve the K₁, K₂ values that best reproduce the data, and separate the composite spectra with these values. In Eq. (1), N_λ is the number of wavelength bins in the selected range, N_epochs is the number of epochs (cf. Table 1). A_i,k, B_i,k are the disentangled spectra obtained for K₁, K₂, evaluated at k, and shifted to the appropriate RVs of the ith epoch. O_i is the observed ith spectrum, σ_i is the S/N in the continuum region in the vicinity of the spectral line, and v = N_λ(N_epochs − 2) is the number of degrees of freedom, since each pixel in each of the two disentangled spectra is considered a free variable. Naturally, one may extend this to other orbital parameters, such as the eccentricity (e) or the orbital period (P), albeit this becomes computationally expensive for the shift-and-add algorithm. However, since these parameters are typically well constrained in SB1 binaries, they can be kept fixed in our study, unless clear deviations are seen that require a revision of the parameters.

It is well known that the shift-and-add technique often results in apparent “emission wings” in the disentangled spectra, which are not of astrophysical origin, but are intrinsic to the method (Quintero et al. 2020). Removing these features typically requires very large number of iterations, making the disentangling procedure inefficient. A workaround is to enforce the disentangled spectra to lie below the continuum for absorption-line spectra, unless indications exist that the spectrum contains emission (Shenar et al. 2019; Quintero et al. 2020). This assumption is adopted throughout the paper, unless otherwise stated.

The method is thoroughly tested by creating mock data of different synthetic binaries, discussed in detail in Appendix A. Generally, the method works well in terms of reproducing the input spectra and RV semi-amplitudes whenever the semi-amplitudes are larger than half the resolution element (K₁, K₂ ≳ 20 km s⁻¹). Fortunately, this holds for the majority of targets. Depending on the S/N, we can retrieve companion spectra down to flux contributions of f₂/f_tot ≈ 1–3%, noting that the limit strongly depends on the data quality and the rotational velocity of the companion. In principle, pushing to fainter companions could have been possible using the strong Balmer lines. However, the strong nebular contamination in the TMBM data often make Balmer lines unusable for this purpose (see Sect. 3.4).

Finally, we highlight a few facts that appear to hold generally: (i) it is easier to detect the presence of a companion than to determine K₂. Therefore, in some cases, we can confidently conclude the presence of a non-degenerate companion even if we cannot determine K₂. (ii) Due to the large number of degrees of freedom, errors on K₁, K₂ from disentangling are larger than those obtained from orbital fitting. (iii) The fainter the secondary, the less of an impact K₂ has on the disentangled spectra. In contrast, even small deviations in K₁ can result in spurious features in the disentangled spectrum of the secondary. This typically results in a secondary spectrum that mimics that of the primary, which can lead to an erroneous detection of a non-degenerate companion. We critically assess these facts throughout our analysis, and discuss this for the individual targets in Appendix B.

3.2 Identifying non-degenerate companions

In the absence of eclipses, the disentangled spectra can only be retrieved up to a scaling factor that depends on the light contribution of the two components to the total flux, l₁ = f₁/(f₁ + f₂) and l₂ = 1 − l₁ = f₂/(f + f₂). The fainter the companion is, the larger the amplification factor 1/l₂ is. Fortunately, the spectral typing does not depend on these scaling factors.

It is important to realise that the method of disentangling, by nature, always returns a disentangled spectrum for a secondary star. It is then up to the user to decide whether the spectrum is “flat”, or whether it is of stellar origin. Whether or not the companion is classified as “non-degenerate” is left as a subjective decision, based on a visual inspection of its spectrum after applying the scaling, as well as comparison to spectral models. We generally avoid relying on Balmer lines due to the strong nebular contamination (Sect. 3.4), unless stated otherwise in Appendix B. We give much weight to the presence of He I and He II lines (depending on the level of nebular contamination). In most cases, when a non-degenerate companion is reported, the stellar spectrum of the secondary is evident, especially upon magnification through re-scaling.

In quite a few cases, it is not clear whether the observed features are of stellar origin, or whether they are caused by nebular contamination (Sect. 3.4) or sources of intrinsic variability (Appendix A). In these cases, we mark the systems as “uncertain SB2” or “uncertain SB1” (SB2:, SB1:), depending on the “degree of belief” we attribute to the spectral features of the secondary, such as their line profiles and widths, their appearance with respect to those of the primary, and the overall level of nebular contamination. Our reasoning for each target is individually discussed in Appendix B. While automated criteria could be developed to identify and classify the binaries, the presence of nebular lines makes it difficult to establish suitable objective criteria. Given the limited sample size, we prefer to keep the classification subjective, since developing an automated algorithm in this context would anyhow boil down to a subjective choice of criteria. The disentangled spectra are presented in Appendix C and are available at the CDS for the independent inspection of the reader.

3.3 Spectral classification and light ratios

For the spectral classification, we use the Marxist Ghost Buster (MGB) code (Maíz Apellániz et al. 2012, Apellâniz 2015) on the disentangled spectra, which allows the user to interactively compare the observed spectra with the spectral type-luminosity class 2-D grid of spectral standards. We use the latest spectral grid derived from the Galactic O-Star Spectroscopic Survey (GOSSS Maíz Apellániz et al. 2011) that includes stars as late as A0 (Maíz Apellániz et al., in prep.). As the GOSSS grid has a spectral resolution 2500, we degrade our spectra to that value.

Spectral classification of hot stars in the blue-violet is usually done in the 3900–5100 Å range. However, here we rely primarily on the disentangled spectra, which only extends down to 4560 Å. The primary luminosity classification criteria and various qualifiers for most O stars involve He II λ4686 and neighbouring lines, which are only present in the non-disentangled spectra. The ‘z’ qualifier (He II λ4686 in absorption and stronger than He I λ4471 and He II λ4542) cannot be established as that depends sensitively on the ratio of these lines, and we therefore omit it from the classification scheme (Arias et al. 2016; Maíz Apellániz et al. 2016). As we are also missing the 4630–4660 Å region in the disentangled spectra, the f and c qualifiers (involving the presence of N and C emission lines in this region) cannot be determined (Walborn et al. 2010; Sota et al. 2011). In this case, we adopted these qualifiers from the original classifications of Walborn et al. (2014), implicitly assuming that the features seen in the non-disentangled spectra belong to the primary. We therefore note that the spectral types provided here reflect those that best match our data (disentangled and non-disentangled combined), and may be revised in the future, when disentangling of the full 3900–5100 Å range might become available.

To obtain estimates for the light ratios and the correct scaling of the disentangled spectra, we rely on calibrations between spectral types and equivalent widths of He II24200, 4542, and He I λ14026, 4144, 4388, and 4471, as long as the equivalent widths are significant and not strongly contaminated by nebular lines. We use equivalent-width formulae provided by Sana et al. (in prep.), which are reminiscent of existing quantitative classification schemes by Conti & Alschuler (1971) and Conti & Leep (1974), but using modern data. If two non-degenerate stars are retrieved from disentangling, we compute the light ratio as a weighted mean between the light ratios obtained from the primary’s scaling and secondary’s scaling. If the spectral subtype of the secondary cannot be determined, we adopt a light ratio of 5% for putative non-degenerate companions and 3% otherwise.

3.4 Impact of nebular lines

Only a limited number of spectral lines are available in the data with sufficient S/N to be useful for the disentangling: the Banner lines Hδ and Hγ, the He 1 λ4026, 4144, 4388, and 4471 lines, the He II λ4200, 4542 lines, and rarely metal lines such as Mg II λ4481 and Si III λ4553, due to the limited S/N. Of those, only the Hell lines are free of nebular contamination, but as we are mostly hunting for faint B-type companions, Hell lines are of limited use. While sky spectra were subtracted to reduce the nebular contamination, strong residuals remain, likely due to the variability of the emission across the Tarantula. This problem worsens by the fact that the nebular residuals are not constant, but vary from epoch to epoch in strength and shape in a non-trivial manner, presumably due to the varying weather and seeing conditions between the different epochs. At the moderate resolution of our data (∆λ ≈ 0.7 Å), the nebular lines cover a significant part of the stellar profiles, and cannot be removed by simple interpolation without significantly biasing the data. To partly account for the nebular contamination, we extended the shift-and-add algorithm to three components, where component C represents the nebular lines (see also Abdul-Masih et al. 2019). Component C is enforced to be static, lie above the continuum, and achieve non-vanishing values only in regions overlapping with the nebular emission. To account for the varying strengths of the nebular lines, we scale the component C with the equivalent widths of the emission cores of Hγ.

While our method for removing the residual nebular-line contamination works reasonably well (see examples in Sect. 3.5 and Appendix B), it does not fully remove the impact of nebular lines. This is because the nebular lines vary in shape and strength, and are at times saturated or over-subtracted. This is especially important when the nebular contamination is very strong (typically in Balmer lines and often in the He 1 λ4471 line), or when the companion is faint, such that the nebular contamination dominates over the secondary’s features. Generally, not accounting for nebular contamination leads to spurious features that are especially apparent in the amplified spectra of the faint secondaries. These features can appear confusingly similar to actual stellar features, though they often exhibit strong asymmetries. As we are not able to fully remove the impact of nebular lines, for some systems, our results are ambiguous in terms of the presence of a non-degenerate companion (see Sect. 3.2).

3.5 Examples: VFTS 350 (O+O) and VFTS 779 (B+BH?)

We illustrate the method by using VFTS 350 and 779 as examples. VFTS 350 was classified 08 V by Walborn et al. (2014), revised to 08.5 V+09.5 V here. Almeida et al. (2017) derived period of P = 69.57 d and an eccentricity of e = 0.351 for this system, which they classified SB1. Indeed, inspection of the dynamical spectrum of VFTS 350 (Fig. 1) does not readily reveal indications for a non-degenerate secondary. However, an implementation of grid disentangling of lines such as He 124388 (Fig.2) reveals a well constrained minimum of χ²(K₁, K₂). We perform a similar analysis for all strong helium lines, avoiding He I λ4471 and the Balmer lines due to strong nebular contamination. We obtain the final K₁, K₂ values by computing the weighted mean of all measurements, amounting to K₁ = 70.3 ± 3.7 km s⁻¹, and K₂ = 91 ± 20 km s⁻¹. We note that the K₁ value derived here exceeds the value of 60.2 km s⁻¹ derived by Almeida et al. (2017). This is a typical bias seen when spectral lines are blended by two components (Bodensteiner et al. 2021; Banyard et al. 2022).

Figure 3 shows a comparison between the observations and the disentangled spectra for the derived K₁, K₂ values at two epochs corresponding to the RV extremes. Upon inspection of the line variability in these two epochs, the presence of a non-degenerate companion becomes readily apparent. Evidently, the secondary exhibits broader lines, making it more difficult to spot. The full disentangled spectra, corrected for line dilution, are shown in Fig. 4. Relying on Sota et al. (2011), we classify the secondary as 09.5 V, and estimated its flux contribution at 25%.

We now repeat this exercise for VFTS 779, which was classified B1 II-Ib by Walborn et al. (2014) and revised to B1 II-III here. Almeida et al. (2017) determined it to be an SB 1 binary with P = 59.9 d and a near-circular orbit. A dynamical spectrum is shown in Fig. 5. Like for VFTS 350, a companion is not readily seen in the dynamical spectrum. In Fig. 6, we show the χ²(K₁,K₂) map obtained for the He I λ4338 line. While K₁ can be well constrained, K₂ is fully unconstrained. Indeed, when inspecting the spectral variability of this line at RV extremes (Fig. 7), a companion cannot be identified. The spectral appearance of the disentangled spectra is virtually independent of the value of K₂. As the spectral separation still requires a K₂ value, we adopt the arbitrary value of K₂ = 3 × K₁ = 93 km s⁻¹ (i.e. assuming the companion is three times less massive). In Fig. 8, the full disentangled spectra are shown. To enhance the spectral features of a putative secondary, we adopt an extreme light ratio of l₂ = f₂/f_tot = 0.03 in the visual. Evidently, the disentangled spectrum is virtually flat, with the exception of features originating in diffuse interstellar bands (DIBs) and nebular contamination.

VFTS 779 therefore does not have an identified companion. With a minimum secondary mass of M₂ > 3.9 ± 0.6 M_⊙ (Table 1), this makes VFTS 779 a prime OB+BH candidate. However, given the that the system can accommodate a secondary as light as ≈3 M_⊙, we cannot rule out that it is a faint non-degenerate companion (main sequence star or helium star), which would contribute less than 1% to the flux. Indeed, VFTS 243 is the only candidate for which non-degenerate companions can be ruled out (Shenar et al. 2022).

3.6 Light curve analysis

The OGLE light curves are used to search for photometric variability in the systems, with focus on eclipses or ellipsoidal variations. We analysed available light curves using the Period04 (Lenz & Breger 2005) tool, which performs a Lomb-Scargle analysis on the data. Peak frequencies are refined using through phase dispersion minimisation (pwkit pdm Stellingwerf 1978; Schwarzenberg-Czerny 1997), and the light curves are phased using these periods and inspected. If a notable signal is observed, we present the light curve and discuss it in Appendix B. While not fundamental for our study, the light curves help uncover several eclipsing binaries and high-order multiples. An in-depth analysis of the light curves is beyond the scope of our study.

4 Results

We perform grid disentangling (1D or 2D) for all targets in the sample. A detailed account of our analysis for each system is given in Appendix B, and the scaled disentangled spectra are presented in Appendix C. Table 2 provides the final orbital parameters, mass ratios q, estimated secondary masses M₂, and classifications. Figure 9 summarises our findings visually by showing, for each system, the 68% confidence interval on the mass of the secondary, M₂. The derivation of these mass ranges is described in Sect. 4.1.

Using spectral disentangling, we can recover non-degenerate companions for the majority of our sample. In total, among the original 51 SB 1 binaries, non-degenerate companions were found in 43 systems (84%), of which 28 are considered certain and 15 less certain. For the majority of these new SB2 systems, K₂ RV semi-amplitudes could be derived, albeit often with substantial associated uncertainties. The hidden companions are typically early B-type stars, with some exceptions. As Fig. 9 illustrates, almost all companions with minimum masses above ≈5 M_⊙ could be recovered through disentangling. In the remaining eight SB1 binaries (16%), non-degenerate companions could not be uncovered, including the confirmed O+BH binary VFTS 243 (Shenar et al. 2022) and two additional OB+BH candidates, VFTS 514 and VFTS 779. We discuss this sample in Sect. 5.2.

4.1 Derivation of mass ratios and mass ranges

For systems in which K₂ could be derived, the computation of q and M₂ along with their 68% confidence intervals follows directly from the inferred mass ratio (K₁/K₂) and the primary mass, for which we adopt M₁ = M_1,SpT with the associated errors. For systems for which K₂ could not be constrained (including all SB1 systems), the computation is more involved. First, we asses the possible range of orbital inclinations for these systems. We assume a maximum inclination of 90° for all targets. The minimum inclination i_crit is estimated by equating the measured equatorial rotation velocity v sin i/sin i_crit with the critical rotation velocity v_eq,crit of the luminous components (Mahy et al. 2022). For simplicity, we adopt a typical value of the critical equatorial velocity as 500 km s⁻¹ (e.g. Townsend et al. 2004), such that: (2)

While we could estimate the critical rotation individually, the errors in spectral-type calibration and disentangling procedure easily dominate over the variance of v_crit over the relevant spectral range. To estimate υ sin i, we measure the full-widths half-maxima (FWHM) of the He I λ4388 and He II λ4542 lines of all spectra through Gaussian fitting, and use υ sin i-FWHM calibrations provided by Ramírez-Agudelo et al. (2015) to obtain v sin i. We only accept measurements for which the peak of the Gaussian exceeds the S/N by 6σ (i.e. the half-maximum is at least 3σ below the continuum). We provide the computed v sin i and i_crit values in Appendix D. We note that the method neglects other sources of star-dependent broadening mechanisms (pressure broadening, micro- and macroturbulence), such that the v sin i values should only be considered approximate. The method generally results in poor lower bounds on the inclination of the order of 10°, and is therefore hardly constraining. We note that this method may not be valid for binaries containing degenerate secondaries due to kicks, which could cause the stellar spin to mis-align with the orbital spin. However, since the method only impacts the upper bound on the mass, it does not have important consequences in the context of identifying OB+BH candidates.

With the mass functions, estimated primary masses (M₁ = M_1,SpT), and inclination ranges, a probability density for q can be computed for systems without a K₂ measurement. To do so, we perform a Monte Carlo simulation. We ignore errors on the mass function, which are negligible in this context. We sample a Gaussian distribution for the mass of the primary using the values in Table 1. We sample the inclination in the range i ∈ [i_crit, 90°] (see Eq. (2)), following the distribution P(i) ∝ sin i, which results from the projection of randomly oriented orbital orientations. For given i and M₁, the mass of the secondary M₂ and the mass ratio q can then be computed. For each SB 1 system (or SB2 for which K₂ is not known), we thus obtain probability densities P(M₂) and P(q). Using these probability densities, we calculate the medians M₂ and q and the 68% confidence intervals. These values are given in Table 2, and shown in Fig. 9.

4.2 High-order multiples

By combining spectroscopic and photometric information, three systems are found to be higher-order multiples. Hence, Table 2 provides 56 entries in total instead of 51. VFTS 64 is a triple system consisting of an outer 07.5 II tertiary orbiting an eclipsing 06+06 contact binary (P_out = 903 d, P_in = 1.3 d). VFTS 120 is found to be a quadruple system comprising an 09.5 IV + B0 V binary with a spectroscopic period of 15.65 d and an eclipsing B+B binary with an photometric period of P = 17.41 d. Finally, VFTS 702 is found to be a quadruple system owing to the presence of two distinct eclipse periods (1.98 d and 2.93 d), but the quality of the spectra does not enable us to determine the spectroscopic properties of the system unambiguously. Since the outer periods of both quadruples could not be determined, Fig. 9 includes 54 entries in total.

Given that high-order multiples were only found owing to eclipses, it is likely that additional ones are present in the sample. The FLAMES and OGLE angular resolutions (0.6″ and 0.26″, respectively) at LMC distance amount to separations of the order of tens of thousands of au, and orbital periods extending up to log P ≲ 8 [d]. In this range, the high-order multiplicity fraction is expected to be larger than 50% (Moe & Di Stefano 2017). However, spectroscopy and photometry are not efficient techniques for establishing the multiplicity of stars in the range log P ≳ 3–4 [d], where instead interferometry or direct imaging are needed, both of which are challenging in the LMC. We therefore cannot rule out that higher-order multiples contaminate our results. However, bright companions that do not trace the orbital period tend to result in K₂ amplitudes of the order of 0 km s⁻¹ (see, e.g. VFTS 64, Appendix B), and such rare cases are discussed in Appendix B. Relatively faint companions would be interesting to find, but have no implication on our analysis. Hence, it is unlikely that higher-order multiplicity impacts our results.

5 Discussion

5.1 Mass ratio distribution

In this section, we combine all available constraints on mass ratios coming from the homogeneous VFTS and TMBM surveys to reassess the power index k of the mass ratio distribution, modelled in the form f(q) ∝ q^K (with 0 < q ≤ 1). As discussed in Sect. 1, while SB2 systems allow for direct measurement of q = K₁/K₂, they are biased towards short periods, high component masses, and mass ratios close to unity. Disentangling is helpful in this sense, since it extends the sensitivity to wider periods and more extreme mass ratios. Hence, the results become less bias-dominated when disentangling is used. However, even with disentangling, the bias is not fully removed, since a substantial fraction of the SB1 systems for which no companions were retrieved are likely to host faint non-degenerate companions with even more extreme mass ratios. Finally, given that the binary detection probability becomes low for extreme mass ratios, one can expect a population of binaries with two non-degenerate companions that were classified as presumably single stars in the original VFTS survey (Sana et al. 2013). These undetected binaries need to be accounted for if we are to obtain an unbiased estimate for κ.

In Fig. 10, we show a histogram of the mass ratios measured for SB2 systems in the TMBM sample (Almeida et al. 2017) and newly uncovered SB2 systems in this study. We omit VFTS 603 due to its uncertain nature (see Appendix B), and we omit VFTS 243 due to it being identified as an O+BH system (Shenar et al. 2022). We further omit the five SB2 systems classified as post-interaction binaries by Mahy et al. (2020a), as we are interested in the natal mass ratio distribution. In addition, we show the median q values for the SB 1 systems, derived as describe in Sect. 4.1. The final population we need to consider comprises non-degenerate binaries that were classified as single stars. To visualise this population in the histograms in Fig. 10, we create a mock population of such binaries by using the binary detection probability function provided by Sana et al. (2013) as a function of q, f_det(q)∙ We digitise the function, extend it down to q = 0, where per definition f_det(0) = 0, and fit a polynomial through the data points (see Appendix E). For every q bin, the number of undetected binaries is given by l/f_det(q) × M(q), where M(q) is the number of detected binaries in the respective bin.

The distributions match the expected behaviour of the different samples. The known SB2 systems are preferentially found in the range 0.6 < q < 1, while the new SB2 systems, obtained from disentangling of the SB 1 systems in this study, probe more uniformly the range 0.2 < q < 1. The “true SB1” binaries (which also include the potential SB2 systems for which K₂ could not be constrained) cluster around more extreme mass ratios of the order 0 < q < 0.5, which is reasonable since a companion was not detected in them. Finally, the undetected binaries have the highest relative contribution for the lowest q values, as could be anticipated.

Assuming a power-law in the form f(q) ∝ q^κ, we now estimate κ in a Bayesian approach by computing the posterior of k given the set of measurements {q_i} as: (3)

We assume a flat prior on κ, such that P(κ) amounts to a normalisation constant, as does the probability for observing the data P({q_i}). We account for the errors on q_i by multiplying the mass-ratio distribution q^κ with the probability distributions of each q_i, P(q_i). The population of undetected binaries is accounted for via multiplication with the binary detection probability, f_det(q)∙ The posterior P(κ | {q_i}) is thus given by: (4)

where i runs over the populations of known SB2, new SB2, and SB1 systems. For SB2 binaries, P(q_i) is modelled as a Gaussian using the values given in Table 2. For SB1 binaries, P(q_i) is directly obtained from our Monte Carlo simulations described above.

The lower panel of Fig. 10 shows the results obtained for P(κ|q_i). The result obtained for the entire sample is κ_all = 0.2 ± 0.2, consistent with a uniform distribution (κ = 0). However, we also perform the analysis for various period cuts between 1 and 1000 d to investigate a possible dependence of κ on the period interval, as advocated by Moe & Di Stefano (2017). While most period cuts yield consistent results in the range −0.2 < κ < 0.2, the group of short-period binaries with P < 10 d exhibits a distinctively different distribution, with implying that very short period binaries tend to have mass ratios closer to unity. A possible explanation or this result could be that many such short-period binaries at more extreme mass ratios (q_ini ≲ 0.5) have already interacted and merged (Soberman et al. 1997; Langer et al. 2020), and are therefore not included in our sample. However, it is also possible that the star-formation process at very short orbital separations tends to form equal-mass binaries. This fundamental question should be explored in detail in the future.

The unbiased κ values derived here are consistent with a uniform mass ratio distribution (κ = 0). This is also evident when comparing the observed cumulative distribution of the sample with that of power laws for different κ values (Fig. 11). A flat mass-ratio distribution is in contrast with the value of κ = −1.0 ± 0.4 provided by Sana et al. (2013) using the same sample. However, Sana et al. (2013) only had access to RV variability information. A glance at their figure 6 reveals that RV diagnostics are poorly sensitive to k. Hence, the discrepancy could be easily masked by statistical and systematic biases. On the other hand, a flat mass ratio distribution is more in line with the findings of Sana et al. (2012) for Galactic O-type stars (κ = −0.1 + 0.6), and also better matches results for Galactic O-types stars by Kobulnicky & Fryer (2007).

It is possible that a few of the SB 1 binaries considered here host neutron stars or black holes, and should have been omitted from the sample of SB1 binaries. In Sect. 5.2, we estimate that ≈4–5 OB+BH are expected to be present in our sample. Such systems have typical mass ratios of the order of 0.5–1 (Langer et al. 2020). Having omitted the confirmed O+BH binary VFTS 243 (Shenar et al. 2022), and given our sample size, the impact of a few O+BH binaries is likely negligible. While OB+NS binaries would strongly skew the distribution towards extreme mass ratios, such binaries are unlikely to be part of our sample. The RV amplitudes of OB+NS binaries are typically of the order of 10–20 km s⁻¹ or less. Hence, such binaries are very unlikely to have passed the threshold of ∆RV > 20 km s⁻¹ imposed by the VFTS survey to be identified as spectroscopic binaries.

5.2 The OB+BH population

Among our sample, only one system could be robustly shown to be an O+BH binary: VFTS 243 (Shenar et al. 2022). In principle, all systems for which a non-degenerate companion could not be found are potential OB+BH candidates. However, not at the same likelihood. The likelihood that a binary in the sample hosts a BH can be written as P(BH) = P(BH|co) · P(co), where P(co) is the probability that the companion is a compact object, and P(BH|co) is the probability that the companion is a BH assuming it is a compact object. The latter is relatively easy to compute. The mass of the most massive neutron star known slightly exceeds 2 M_⊙ (Fonseca et al. 2021), and neutron stars are not expected to exist above 3 M_⊙ given constraints from their equation of state (Oppenheimer & Volkoff 1939; Akmal et al. 1998; Lattimer & Prakash 2004; Oertel et al. 2017). While the bulk of known BHs in the Galactic neighbourhood exceed 5 M_θ (Özel et al. 2010; Farr et al. 2011; Corral-Santana et al. 2016), the existence of lower-mass BHs is still debated (Kreidberg et al. 2012; Wyrzykowski & Mandel 2020; Lam et al. 2022). Hence, to reasonable approximation, P(BH|co) ≈ P3 := P(M₂ > 3 M_⊙). In Sect. 4.1, we outlined our Monte Carlo method for establishing probability distributions for M₂ and q. We can use the probability densities P(M₂) to estimate P3.

Less trivial is the computation of P(co), which is the complementary probability of the companion being a non-degenerate companion (assuming it is a single star). In principle, one can perform simulations to narrow down the possible types of companions (e.g. faint main sequence stars, helium stars; see Shenar et al. 2022; Mahy et al. 2022). One would then need to combine this with assumptions regarding the likelihood of these companions, which involves the initial mass function and binary evolution. Given the non-trivial uncertainties and complications involved, we simply consider P3 as a proxy for P(BH). These probabilities are given in Table 2. For uncertain SB2 systems (SB2:), we provide the values in parentheses, since P(co) is systematically lower in this case.

Aside from the confirmed O+BH system VFTS 243, two SB 1 systems approach or exceed a probability of 90%: VFTS 514 and 779. We mark these systems as prime OB+BH candidates, noting that non-degenerate companions could only be fully excluded in VFTS 243. It is also possible that BHs hide among the other five SB1 binaries, but at a lower probability. The 15 uncertain SB2 or SB1 binaries (SBL, SB2:) could also host BHs – especially those with high P3 probabilities. More data will be needed to clarify this question (see below).

We now estimate the number of OB+BH binaries expected to lurk among the 51 SB1 systems considered here. According to Langer et al. (2020), about 2% of an unbiased population of massive OB-type binaries are expected to host BH companions². There are 360 O-type stars³ included in the VFTS survey, and hence we expect roughly seven O+BH binaries among them. With typical RV semi-amplitudes in excess of 50 km s⁻¹, these O+BH binaries are all expected to have been flagged as binaries in the VFTS survey. The TMBM follow-up targeted 69% of the flagged binaries (Almeida et al. 2017), hence there should be 4–5 O+BH binaries in the TMBM sample. Naturally, none of those would have been classified as SB2 by Almeida et al. (2017), and so we can conclude that these 4–5 O+BH binaries would be present in our SB1 sample (8–10%). With three promising OB+BH candidates (one of which is unambiguous) and 20 less likely candidates, our results seem qualitatively in line with predictions. However, we refrain from over-interpreting this result until the OB+BH nature of the remaining candidates can be verified.

It is interesting to note that the orbit of the confirmed O+BH binary VFTS 243 is near-circular, implying negligible kick or ejecta for the BH progenitor (Shenar et al. 2022). This is in contrast to the recent report of a BH in the eccentric (e ≈ 0.5) binary HD 130298, which otherwise shares a similar orbital period and spectral type to VFTS 243 (Mahy et al. 2022). The orbit of the OB+BH candidate VFTS 779 is near-circular, while VFTS 514 is eccentric. If all of these systems are truly O+BH binaries, this may give evidence that some BHs are born with strong kicks, while others are not, in line with recent predictions (e.g. Fryer et al. 2012).

Future high-resolution spectroscopy should significantly improve the detection thresholds, since it should allow for a better removal of the nebular lines and thus make the Balmer lines usable. The Balmer lines are key diagnostics to reject the presence of companions contributing as little as ≈0.5% to the flux. Moreover, high resolution is vital for binaries with RV amplitudes comparable or smaller to the resolution element of FLAMES (K₁ ≲ 40 km s⁻¹). Additionally, UV spectroscopy should help reject the presence of non-degenerate companions through dilution of P Cygni lines (e.g. Georgiev et al. 2011; Shenar et al. 2016).

6 Conclusions

In this study, we characterised the nature of the hidden companions in 51 SB1 O-type and evolved B-type binaries identified in the framework of the Tarantula Massive Binary Monitoring (TMBM) campaign (Almeida et al. 2017). We implemented the shift-and-add disentangling algorithm to separate the component spectra and establish the RV amplitudes of the two components, whenever possible. We also investigated OGLE light curves in search for eclipses or other orbital modulations. Below, we summarise our main results:

1. Out of the 51 SB1 binaries, 43 are found to host non-degenerate companions: 28 are considered certain, and 15 are considered less certain. We also find one triple (VFTS 64) and two quadruple (VFTS 120 and VFTS 702) systems, and eight eclipsing binaries. The remaining eight targets retain their SB1 classifications.

2. Combining constraints from our sample with previously known SB2 binaries (Almeida et al. 2017), and accounting for detection bias, we derive κ = 0.2 ± 0.2 for the exponent of the mass-ratio distribution f(q) o q^κ in the range 0.05 < q < 1 and 1 < P < 1000 d. This drops to κ = −0.2 ± 0.2 if binaries with periods shorter than 10 d (potential post-interaction products) are removed, and increases to κ = 1.2 ± 0.5 for binaries with P < 10 d, which we propose is a result of binary interaction among the tight binaries. Our results are therefore consistent with a flat natal mass-ratio distribution for O-type stars at LMC metallicity.

3. Through a probabilistic approach, we identify three O+BH candidates: VFTS 243, 514, 779, of which VFTS 243 is confirmed as an O+BH binary (Shenar et al. 2022). The other five SB1 binaries (VFTS 225, 631, 645, 743, 769) and 15 uncertain SB1/SB2 binaries (VFTS 73, 171, 184, 256, 332, 409, 479, 603, 619, 657, 736, 750, 812, 827, 829) could also be O+BH binaries, but at a lower probability.

We strongly encourage further investigations of the various important sub-samples uncovered in this study: the OB+BH candidates, the remaining SB1 binaries (including uncertain SB2), the eclipsing binaries, and the higher-order multiples. The VFTS and TMBM samples offer unparalleled laboratories to homogeneously study binary interaction and evolution at subsolar metallicity. Step by step, the nature of each and every star in this population is illuminated, providing indispensable constraints on evolution models of massive single, binary, and multiple stellar systems at subsolar metallicity.

Acknowledgments

This research has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement number 772225: MULTIPLES). TS acknowledges support from the European Union’s Horizon 2020 under the Marie Sklodowska-Curie grant agreement No 101024605. This work is based on observations collected at the European Southern Observatory under programme IDs 090.D-0323 and 092.D-0136. L.M. thanks the European Space Agency (ESA) and the Belgian Federal Science Policy Office (BELSPO) for their support in the framework of the PRODEX Programme. P.A.C. acknowledges support from the UK Science and Technology Facilities Council research grant ST/V000853/1. A.H. acknowledges support by the Spanish MCI through grant PGC-2018-0913741-B-C22 and the Severo Ochoa programme through CEX2019-000920-S. J.M.A. acknowledges support from the Spanish Government Ministerio de Ciencia, Innovaciόn y Universidades through grant PGC2018-095 049-B-C22. P.M. acknowledges support from the FWO junior postdoctoral fellowship No. 12ZY520N. This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 945806). S.T. acknowledges support fromthe Netherlands Research Council NWO (VENI 639.041.645, VIDI 203.061 grants). M.G. is supported by the EU Horizon 2020 research and innovation programme under grant agreement No 101004719.

Appendix A

The applicability regime of the shift-and-add method

In the context of our study, disentangling comprises two main steps: the retrieval of the RV amplitudes (K₁, K₂), and the retrieval of the disentangled spectra. The robustness of the disentangling procedure depends on many factors, including the individual S/N, Doppler-space coverage and number of observations, line profiles, light ratios, and RV amplitudes.

The underlying assumption of spectral disentangling, as implemented here, is that each observed spectrum is the sum of stellar spectra (two in the case of a binary) that are Doppler-shifted with the respective RVs of the stars, contaminated only by Gaussian noise and static nebular lines that can vary in strength from epoch to epoch (Sect. 3.4). Hence, intrinsic non-Doppler variability can impact the results, especially if it is (quasi-)periodic. In the case of O-type stars, the most important source of contamination is likely stellar pulsations, which can result in apparent peak-to-peak RV variability of up to ≈ 20 km s−1 (Ritchie et al. 2009; Aerts et al. 2009; Sana et al. 2013). Since all binaries in our sample have K₁ > 10 km s−1 per definition, we expect intrinsic variability to play a secondary role for most targets. Moreover, we note that the procedure of disentangling tends to smear-out features that do not follow the orbital period. We therefore do not expect pulsations to impact the qualitative or quantitative results significantly. An exception is when the dominant pulsation period is mistaken for the orbital period, or when they are otherwise comparable. A notable example is VFTS 603, which has a derived RV semiamplitude of K₁ = 11.4 km s−1 and a period of P = 1.76 d. The low-amplitude motion, combined with the short period, raises doubts that pulsations cause the underlying variability. Disentangling does suggest the presence of two non-degenerate stars, but it is possible that the second source is a spurious result of non-Doppler variability (see Appendix B).

Fully mapping the impact of all parameters and line-profile variability on the output of spectral disentangling is beyond the scope of our study, and would require a dedicated study with prespecified models of line variability. Critical for our study is to establish threshold light ratios below which the spectrum of a non-degenerate secondary would be discerned from the noise. To estimate it, and to test the applicability regime of our method, we run a set of simulations of synthetic noisy composite spectra that mimic the resolving power, S/N, and timing of the TMBM programme. Specifically, we consider a hypothetical binary with P = 100 d, a circular orbit, and a mass ratio of 0.5 (i.e. K₂ = 2 K₁). We use synthetic TLUSTY spectra (Lanz & Hubeny 2003, 2007; Hubeny & Lanz 1995) for the primary and secondary. For the primary, we adopt T₁ = 35 kK, log g = 4.0 [cgs], and v sin i = 200 km s−1. For the secondary, we assume T₂ = 26kK, log g = 4.0 [cgs], and v sin i = 100 km s−1. We mimic the data quality of VFTS 231, which has varying S/N for the individual spectra between ≈ 10 and 40, with an average of ≈ 30. VFTS 231 is among the 10% worst-quality targets in terms of S/N, and is chosen to ensure that our estimates are conservative.

We constructed 12 mock binaries in total, spanning the values K₁ = 200, 50, and 20 km s−1 (and hence K₂ = 400, 100, and 40 km s−1), and secondary light contributions of l₂ = f₂/f_tot(V) = 0.5,0.1,0.03, and 0 (the latter corresponding to a black hole secondary). We then apply the 2D shift-and-add disentangling technique to derive K₁, K₂ and disentangle the spectra using those values. The results are shown in Fig. A.1. Evidently, the retrieval procedure works well across the parameter space, although discrepancies become apparent for lower RV amplitudes and fainter secondaries. Figure. A.1 illustrates that, at the low data quality of the simulations, companions contributing as little as ≈ 5% can be retrieved. This improves to 1 – 3% for higher S/N data available for some of our targets.

Appendix B

Comments on individual targets

This appendix provides a detailed commentary on the individual targets and the nature of the unidentified secondaries. We also discuss the light curves of the targets (when available). Unless stated otherwise, cited periods and eccentricities are taken from Almeida et al. (2017).

VFTS 64, O8 II:(f) + (O7 V + O7 V), was reported to be a long-period, eccentric binary with P = 903 d and e = 0.53. Among our sample, it has the largest mass function, with f (m) = 10.7 ± 0.4 M_⊙. Spectral disentangling unambiguously shows the presence of a companion. However, disentangling of various He lines results in unrealistically small K₂ values, and the disentangled spectra of the secondary appear to centre around a different systemic velocity than that derived by Almeida et al. (2017). This suggests that this system may be a multiple.

Available OGLE photometry confirms the triple nature of the system. Figure B.1 shows the phased I-band light-curve of VFTS 64. A period of P = 1.327181 ± 0.000003 is found from the light-curve analysis. This period is clearly different than the one derived from RVs (903 d). The light curve suggests that the secondary is in fact an eclipsing contact system consisting of two nearly-identical stars, such that the system is a triple with an outer period of P_out = 903 d and an inner period of P_in = 1.33 d.

We extended the 2D disentangling technique to account for three components. Denoting with the subscript ‘2’ the contact binary as a whole, we now need to constrain the RV amplitudes K₃ and K₄ of the components of the contact binary. We fix the orbital parameters derived by Almeida et al. (2017) for the primary, including its RV amplitude of K₁ = 57.2km s−1. The period and time of primary eclipse P_in, T_0,in are fixed to those found from the light curve analysis (Fig. B.1). Given the properties of the light curve, we adopt e_in = 0 and ω_in = 90°, such that T_0,in corresponds to zero RV shift. Assuming that the inclinations of the inner and outer orbits are aligned given the relatively low separation between the binary and the tertiary, it can be shown from Newtonian relations that the RV amplitude of the contact binary itself, K₂, is obtained via (B.1)

We note that by requiring K₂ > 0 km s−1, Eq. (B.1) imposes the condition K₃ + K₄ ≳ 420 km s−1. We then measured K₃ and K₄ by disentangling the strong He II and He I λ4471 lines (e.g. Fig. B.2), whose weighted mean yields K₃ = 282 ± 36kms−1, K₄ = 323 ± 33 kms−1, yielding in turn K₂ = 41 ± 13. This unique system could be a prototypical progenitor of the triple WR system BAT99 126, which was also found to host a contact binary (Janssens et al. 2021).

VFTS 73, O9.5 IV + B: has a reported period of P = 151 d and an eccentricity of e = 0.20. The spectral variability does not readily suggest the presence of a non-degenerate companion (Fig. B.3). With a mass function of f(m) = 0.29 M_⊙ and an estimated primary mass of 19 M_⊙, the companion weighs at least 5.6 ± 0.8 M_⊙. Disentangling of the spectra does not yield a unique K₂ value. When adopting plausible values of the order of 70 km s⁻¹ very faint spectral features appear in the secondary’s disentangled spectra, primarily in Balmer lines and strong He I lines. However, they could be the result of nebular contamination. We tentatively classify this binary as SB2: with a B-type companion and an assumed a 5% light contribution, but we cannot rule out a BH. Analysis of the OGLE light curve did not reveal significant periods.

VFTS 86, 09.5 III + 08 IV:n has a reported period of P = 183 d and an eccentricity of e = 0.51. Inspection of the line profiles is suggestive of the presence of two non-degenerate stars in the system, both of which appear to contribute a significant amount to the visual flux (Fig. B.4). For this reason, we resort to 2D disentangling. We used the strong He I lines (which suffer very little nebular contamination) and the two He II lines. The lines yield broadly consistent results, albeit with large errors. The smallest errors are obtained from the He I/14388 line, yielding K₁ = 43 ± 8 km s−1 and K₂ = 48 ± 10 km s−1. With a derived light contribution of 44%, the secondary appears to be an 08 star, that is, of earlier spectral type than the primary. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 93, 08.5 V + B0.2: V has a reported period of P = 250 d and an eccentricity of e = 0.20. The primary moves with a very low RV amplitude of K₁ = llkms−1, making it difficult to properly disentangle the spectra. However, careful inspection of the line-profile variability of strong He I lines at extremes suggests the presence of a non-degenerate companion (Fig. B.5). We re-derived K\ by measuring the RVs of the He II lines, which are not contaminated by the secondary, and obtain K₁ = 12.8 + 1.1 km s−1. Fixing this value, the strong He I lines all result in minima of the order of K₂ = 10 – 30 km s−1. A weighted mean of all measurements yields K₂ = 20 ± 9 km s−1. The spectrum is dominated by He I absorption, and best corresponds to a B0.2 spectral type, with an estimated light contribution of 21 %. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 120, (09.2 IV + B) + (09.5 V + 09.5 V) was reported to be have a period of P= 15.7 d and an eccentricity of e = 0.28. Inspection of the He I lines reveals line-profile variations that are suggestive of a second non-degenerate star in the spectrum (Fig. B.6). For this reason, we implemented the 2D disentangling algorithm. However, the solution appeared to be strongly dependent on the lines used. Moreover, the He I λ4388 and He I λ4471 lines, which are the strongest and appear to be hardly contaminated by nebular emission, result in K₂ amplitudes of the order of 0 km s⁻¹. As in the case of VFTS 64, this could imply the presence of additional stellar components in the system.

Interestingly, Fourier analysis of the OGLE light curve shows a peak at 8.7 d and its harmonics. Folding the light curve with twice this value, that is, P = 17.41 d, reveals a clear signature of eclipses in the system (Fig. B.7). It can be fully ruled out that this period corresponds to the spectroscopic period. Given the similarity of the periods, an eclipsing tertiary can be ruled out. The only conceivable way to explain this is that VFTS 120 is a quadruple system. For this reason, we extended the shift-and-add algorithm to account for four components. To reduce the parameter space, we adopt for the eclipsing binary the photometric period, an eccentricity of 0 (as implied by the ≈ 0.5 phase difference in eclipses), and the time of periastron (defined here as time of primary eclipse). The four RV amplitudes K₁,K₂,K₃, K₄ are then scanned through a χ² minimisation. We find K₁ = 92 ± B km s−1, K₂ = 93 ± 17 km s−1, K₃ = 20 ± 7km s⁻¹, and K₄ = 21 ± 7 km s⁻¹. However, the results should be taken with caution, given the assumptions, limited resolution, and data quality available here. The outer period connecting the two binaries cannot be established here.

VFTS 171, 08.5 III:(f)) + B1.5: V, was reported to be a long-period binary with P = 677 d and e = 0.56. The low-amplitude motion of the primary (K₁ = 12 km s⁻¹) prevents us from obtaining unambiguous results for this system, and inspection of the line-profile variability is not readily suggestive of a non-degenerate secondary (Fig. B.8). Disentangling the spectra results in faint signatures in Balmer lines, He I lines, and Hell lines. However, we cannot derive K₂ from the relatively flat χ² map. We tentatively classify the companion as B1.5: V and estimate its light contribution at 8%, but given the low RV amplitude, we consider the SB2 nature uncertain. The OGLE light curve shows a variability pattern when phased at the orbital period (Fig. B.9), though it is not clear that this pattern is of really astrophysical origin.

VFTS 184, 06.5 Vn +OB: has a reported period of P = 32 d and eccentricity of e = 0.20. However, the orbital solution derived by Almeida et al. (2017) shows substantial scatter. The lack of line-profile variability at RV extremes does not readily suggest the presence of a non-degenerate companion (Fig. B.10). Given the scatter in the orbital solution, we attempted 2D disentangling across the K₁,K₂ axes, and probing all Hei and Hen lines, as well as Hδ, which is weakly contaminated by nebular lines. The results of the different lines are consistent, but the errors are very large, with final mean averages of K₁ = 29 ± 11kms−1 and K₂ = 61 + 29 km s⁻¹. Disentangling the spectrum for these RV amplitudes results in similar spectra for both components, which sheds doubt on the validity of the results, as this is often the case when the adopted K₁ is significantly different than the true one. When disentangling for K₁ = 12.1 km s⁻¹, the features in the spectrum of the secondary weaken, though they do not fully disappear. Our results are therefore uncertain. We tentatively classify the companion as OB:, and refrain from using our measurements for any further analysis. With a minimum mass of M₂ > 1.7 ± 0.4 M_⊙, the nature of the secondary is virtually unconstrained. Higher-resolution spectra will be necessary to disentangle this system. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 191, 09.2 V + 09.7: V: has a reported period of P = 359 d and an eccentricity of e = 0.22, although the orbital solution (especially the eccentricity) are uncertain. The spectral variability does not readily suggest the presence of a non-degenerate companion (Fig. B.11). However, grid disentangling consistently points towards the presence of a relatively bright companion entangled in the spectrum. We are not able to constrain a K₂ value using disentangling of the strong He I lines. However, adopting plausible K₂ values of the order of 2 – 3 K₁ ≈ 50 – 70 km s⁻¹ results in a secondary spectrum that appears to belong to a fainter, rapidly rotating O-type star, contributing roughly 13% to the total flux. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 201, 09.7 V + B1.5: V has a reported period of P = 15.3 d and eccentricity of e = 0.46. The variability seen in the He I lines (Fig. B.12) is suggestive of a second non-degenerate component in the system. Due to the lack of strong He II lines and the apparent contribution of the secondary, we resort to 2D disentangling. All He I lines are in broad agreement with each other, but show substantial errors on the RV amplitudes. Disentangling all strong He I simultaneously, we find K₁ = 79 ± 9 km s−1 and K₂ = 105 ± 26 km s−1. The disentangled spectrum of the secondary, which has an estimated light contribution of roughly 30%, shows strong He I lines and no He II lines, with only faint evidence for the Mg II/U4481.1,4481.3 doublet, and it is classified as B1.5: V. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 225, B0.7 III has a reported period of P = 8.2 d and a virtually circular orbit. Disentangling of the strong He I lines (Fig. B.13) does not reveal a significant minimum for K₂, and the disentangled spectrum of the secondary is virtually featureless in He I lines, and shows very little signatures in Balmer lines that appear to originate from inaccurate subtractions of the weak nebular lines. The system is thus a bona fide SB1 system. However, the minimum mass imposed by the orbital parameters and the estimated mass of the primary yield M₂ > 1.8 ± 0.3 M_⊙, such that the nature of the secondary remains unconstrained.

The OGLE III data show substantial scatter and a long-term magnitude increase that is not supported by the OGLE IV data. We therefore remove them from the photometric analysis. Fourier analysis of the OGLE IV data reveals a peak at the orbital period. The phased light curve is shown in Fig. B.14.

VFTS 231, 09.7 V + B1.5: V has a reported period of P = 7.9 d and eccentricity of e = 0.41. Almeida et al. (2017) noticed the presence of a secondary star in the system, but could not derive its RVs. Indeed, comparing line profiles belonging to strong He I lines such as He 124026 (Fig. B.15) clearly shows the presence of a non-degenerate companion. We implemented 2D disentangling and find consistent results for all He I lines. A mean average of the measurements obtained for the strong He I lines (with exception of the He I λ4471 line, which suffers strong nebular contamination) yields K₁ = 99 ± 12 and K₂ = 130 ± 25 km s⁻¹. From the strengths of the spectral lines, we estimate a light contribution of 33% for the secondary. We classify it as B1.5:. Despite the short period of the system, we cannot find significant periodicities in the OGLE light curve.

VFTS 243, 07 V(n)((f)) + BH was shown to be a bona fide O+BH binary with a P = 10.4 d period and a near-circular orbit in Shenar et al. (2022), to which we refer for a detailed analysis of the spectra and light curve of this unique binary.

VFTS 256, 07.5 V: +OB: has a period of P = 246 d and an eccentricity of e = 0.63. Inspection of the line-profile variability at RV extremes does not readily reveal the presence of a non-degenerate companion (Fig. B.16). The low-amplitude motion (K₁ = 19.2 km s−1) and nebular contamination make disentangling of the system challenging. Given the potential contamination of the secondary, we performed 2D disentangling, but find a consistent K₁ value to that derived by Almeida et al. (2017), and therefore adopt their value. We cannot find a minimum in the χ²(K₂) map, and the disentangled spectra provide marginal evidence for the presence of a companion at faint He I lines features and Banner features, which may however be impacted by nebular contamination. Comparably weak He II features are seen, thought these features may be spurious due to an inaccurate K₁ value. We tentatively classify this as an SB2:. With a minimum mass of M₂ > 4.1 ± 0.5 M_⊙, the secondary could be a faint non-degenerate star or a BH. Analysis of the OGLE light curve does not reveal any periodicity at the orbital period, but does show periodicity at P = 341 d. This periodicity may well be related to the yearly cycle, and hence cannot be confirmed to represent an intrinsic period of the system.

VFTS 277, 09 V + B1.5: V has a period of P = 240 d and a high eccentricity of e = 0.93. The mass function and estimated mass of the primary imply a companion of at least ≈ 3 M_⊙. The spectral variability of He I lines weakly suggests the presence of a non-degenerate companion (Fig. B.17), though not conclusively.

Disentangling reveals significant He I absorption for the secondary but no He II absorption. However, K₂ cannot be constrained with the current data quality. Based on the appearance of the Hδ region and the Mg II, 124481.1, 4481.3 doublet, we tentatively classify the companion as B 1.5. However, as the spectral appearance of the secondary depends on the adopted K₂ value (here, we use K₁ = 2 × K₂ = 128 km s−1), the classification is not certain. Moreover, the presence of a faint, blue-shifted He II λ4200 absorption in the spectrum of the secondary implies that this system may be a higher-order multiple. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 314, 09.7 V(n) + B has a short period of P = 2.5 d and an eccentricity of e = 0.17. Almeida et al. (2017) noted the presence of a secondary in the individual spectra, but could not measure its RVs. Indeed, disentangling of the He I lines (Fig. B.18) results in a relatively well defined minimum around 200 km s−1 ; the most robust result is found from disentangling of the He I λ4026 line, which yields K₂ = 232 + 16 km s⁻¹. The corresponding disentangled spectrum shows a spectrum with clear He I absorption and no He II absorption, but a precise classification of the secondary is difficult due to the lack of metal lines. We estimate that the companion contributes 5% or less to the visual flux.

As could be anticipated from the short period, VFTS 314 is an eclipsing system. The light curve of VFTS 314 is shown in Fig. B.19. From analysis of the light curve, the period can be refined to P = 2.550786 ± 0.000005. Interestingly, the primary and secondary minima do not yield an identical period, suggesting that the system experiences apsidal motion.

VFTS 318, 09.5 V + 09.2 V is a nearly circular binary with a reported period of P = 14 d. Careful inspection of the line variability at RV extremes of He I lines such as He 124026 or He I λ4388 (Fig. B.20) suggests the presence of a second non-degenerate star in the binary. For this reason, we implemented 2D disentangling on various He I lines. A weighted mean of the RVs measured from the He I lines yields K₁ = 53 ± 12 Km s⁻¹ and K₂ = 93 ± 24 km s⁻¹. The Hell lines do not yield well-constrained measurements. The value of K\ is more than twice the amplitude derived by Almeida et al. (2017), a consequence of accounting for the secondary. The two stars a very similar spectral type, and the secondary is estimated to contribute 33% to the visual flux. It is therefore somewhat surprising that the secondary is found to have roughly half the mass of the primary, suggesting that the primary may be an evolved star. Analysis of the OGLE light curve did not yield significant periods.

VFTS 329, 09.5 V(n) + Bl: V: has a reported period of P = 7.0 d and eccentricity of e = 0.44. The spectral-line variability of He I lines suggests the presence of a fainter non-degenerate secondary in the spectrum, as is shown for the He I λ4338 line in Fig. B.21. Given the non-negligible contribution, we perform a 2D disentangling for the system. A mean average of the of strong He I, aside from the contaminated He I λ4471 line, yields K₁ = 93 ± 8, K₂ = 131 ± 24 km s⁻¹. The disentangled spectrum of the secondary, which has an estimated light contribution of 22%, matches a Bl V spectral type. Despite the relatively short period, analysis of the OGLE light curve did not yield significant periods.

VFTS 332, 09 III + 09.2 V is an eccentric long-period binary with P = 1025 d and e = 0.81. The coverage of the spectra in Doppler space is very poor, with only one epoch showing significant RV shifts compared to the others. Therefore, results from disentangling of this system should be taken with caution. The results point to the presence of a second, non-degenerate object in the system, although its presence is not readily seen from the spectral variability (Fig. B.22). Disentangling of all He lines results in a weighted mean of K₁ = 79 ± 10 km s⁻¹ and K₂ = 46 ± 13 km s−1. As in the case of VFTS 184, the similarity of the disentangled spectra suggests that the secondary spectrum may be contaminated by the primary due to a mismatch of the K₁ value. However, regardless of the K₁ value used for disentangling, a stellar spectrum is always obtained for the secondary, suggesting that it is likely real.

Taking our results at face value, it appears that the hidden secondary is the more massive star, with a spectral class of roughly 09. Despite of this, the disentangled spectra imply that it contributes only 24% to the total light, suggesting that the primary is perhaps more evolved. On the other hand, the spectral appearance of the secondary also corresponds to an evolved object. These results combined are hard to reconcile with each other, suggesting that the disentangled spectra suffer from cross-contamination between the components. High resolution spectroscopy, obtained during periastron passage, will be required to disentangle this object unambiguously. Analysis of the OGLE light curve did not yield significant periods.

VFTS 333, 09 ∏((f)) + 06.5 V: is an eccentric long-period binary with P = 980 d and e = 0.75,. Unlike VFTS 332, the spectra offer a good coverage in Doppler space. A non-degenerate companion can be seen at the RV extreme of strong He I lines (Fig. B.23). Disentangling of the strong Hell and He I lines yields values between 20–40 km s−1 for K₂, and suggest that the secondary is relatively bright. We therefore resort to 2D disentangling. In all cases, K₁ is found to be of the order of 40 km s⁻¹, with a weighted mean of K₁ =43.0 ± 2.5 and K₂ = 34 ± 14km s⁻¹. All measurements suggests that the object previously identified as the primary is in fact the less massive companion. We estimate it to contribute 42% to the flux based on the strengths of the spectral lines. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 350, 08.5 V + 09.5 V was thoroughly discussed in Sect. 3.5. Analysis of the OGLE light curve does not reveal any significant frequencies.

VFTS 386, 09 V(n) + Bl V: has a reported period of P = 20 d, an eccentricity of e = 0.25, and a low-amplitude motion of K₁ = 14.3 km s−1 for the primary. However, the orbital solution established by Almeida et al. (2017) exhibits substantial scatter, with a reduced χ² of 9.5. While the presence of a non-degenerate companion is not readily clear from the line profile variability, comparison of the variability in He I and He II lines reveals that this binary hosts two non-degenerate companions (Fig B.24). Given the low amplitude motion and the lack of significant variability in He I lines, we do not proceed with 2D disentangling. Instead, we re-derive the orbital solution of the system by focusing only on the Hell lines. The revised orbital solution yields a reduced χ² of 2.4 (Fig. B.25). The orbital parameters are then P = 20.45 ± 0.02 d, T₀ = 54816.5 ± 2.8 [JD - 2400000], K₁ = 30.8 ± 2.1 km s−1, and e = 0.15 ± 0.05. We then proceed with ID disentangling of the He I lines, fixing the orbital parameters derived above. The best constraints are obtained from the He I λ4388 (Fig. B.26) and He I λ4471 lines, which yield K₂ = 64 ± 12 km s⁻¹. The spectral type of the secondary corresponds to a B1 star, and we estimate its light contribution to be 16%. Analysis of the OGLE light curve did not yield significant periodicities.

VFTS 390, 05.5 V:((fc)) + 09.7: V: has a reported period of 22 d and an eccentricity of e = 0.49. Disentangling of the spectra proved to be very challenging due to very strong nebular-line contamination and the fact that the secondary appears to possess very weak lines (Fig. B.27). While disentangling of various He I lines supports K₂ values in the range 100 – 250 km s−1, the errors are very large. A weighted mean of the RVs obtained for He I λ4026 and He I λ4388 yields K₂ = 137 ± 34 km s−1, but the nebular emission and very low S/N may entail systemic errors. The secondary’s spectrum matches a 09.7 V star, and has an estimated light contribution of 15%.

Interestingly, the system shows a single eclipse (Fig. B.28), which substantiates beyond doubt the presence of a non-degenerate companion in the system, and allows us to slightly refine the ephemeris. From the orbital configuration, the observed eclipse occurs when the hotter, brighter primary eclipses the cooler, fainter secondary. Due to the orbital configuration, a second eclipse is not seen.

VFTS 404, 03.5: V:((fc)) + 05 V: has a reported period of P = 146 d and an eccentricity of e = 0.72. Almeida et al. (2017) noted the presence of a non-degenerate companion, but could not derive its RVs. Indeed, the secondary is clearly seen in both He I and He II lines (Fig. B.29). Given the significant contribution of the secondary to these lines and the fact that the components are constantly blended, it is very likely that the RVs measured by Almeida et al. (2017) for the primary were strongly impacted by the secondary. Therefore, we utilise 2D disentangling for this system. For this, we use the He II lines, which are free of nebular contamination and are stronger than the He I lines. A weighted mean yields K₁ = 98 ± 11 km s−1 and K₂ = 106 ± 11 km s−1. We note that the K₁ value derive here is roughly three times larger than derived by Almeida et al. (2017), which is a result of the secondary’s anti-phase motion that leads to an apparent lower amplitude of motion. The companion is estimated to contribute 42% to the visual flux, and is classified 05 V. Analysis of the OGLE light curve did not reveal significant frequencies. We note that VFTS 404 is one of the four targets for which X-rays are detected, and exhibits the highest X-ray luminosity in our sample (log L_X = 32.84 [erg s−1]). This luminosity lies slightly above the canonical value of L_X/L ≈ −7, perhaps due to the presence of colliding winds in the system.

VFTS 409, 03.5: V:((f)) + B: has a reported period of 22 d and an eccentricity of e = 0.29. Unfortunately, the spectra are extremely contaminated by nebular lines that vary in strength, making only the He II lines available for a reliable analysis. Disentangling of the He II lines reveals no features. Disentangling of the He I lines reveals no clear minimum for K₂ is retrieved, with values in the range 20 – 300 km s−1 all yielding acceptable an χ². The disentangled spectrum shows faint signatures in He I lines (Fig. B.30), though nebular contamination questions the validity of these results. If these signatures are real, then the light contribution of the secondary is ≈ 0.05. We tentatively classify the secondary as B:, and the binary as SB2:. With a minimum mass of M₂ > 7.6 ± 1.1 M_⊙, the secondary could be an early B-type star, a helium star, or a BH. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 429, 07 V: + Bl: V: has a reported period of P = 30 d and an eccentricity of e = 0.56. This system is clearly SB2, as is readily seen from the spectral variability of strong He I lines. However, the spectra are strongly contaminated by nebular lines, making the derivation of K₂ challenging. The least contaminated line clearly showing both stars isolated is the He I λ4144 line (Fig. B.31). A weighted mean of the strong He I lines (with exception of He 124471 due to very strong nebular contamination) yields K₂ = 143 ± 27 km s⁻¹. Judging by its spectral appearance, we classify the secondary as Bl: V:, and estimate its light contribution in the visual at 10%.

Similar to VFTS 409, phasing the OGLE light curve of VFTS 429 with its spectroscopic ephemeris reveals a single eclipse (Fig. B.32). We used it to slightly refine the orbital period. Unlike VFTS 409 however, the eclipse in VFTS 429 occurs when the cooler B0 V secondary eclipses the hotter 07.5 V primary.

VFTS 440, 06: V:(f) + 08 V, is a long-period binary with P = 1019 d and e = 0.28. The low RV amplitude of the primary (K₁ = 11.6 km s−1) and the strong nebular contamination make the disentangling of this binary challenging. We first implemented the 2D disentangling of the He II lines to ensure that the low RV amplitude is not impacted by a secondary contributing significantly to the flux. We find K₁ that agree well with the value derived by Almeida et al. (2017), and we therefore fix it and continue with ID disentangling across the K₂ axis. While not readily apparent, careful inspection of the He II λ4542 line at RV extremes suggests the presence of a faint secondary (Fig. B.33). Disentangling of the He II lines results in poor constraints on K₂, since all values in the range [0,100] km s⁻¹result in acceptable χ² values (within 1er). However, disentangling suggests that the He I λ4009 line is present only in the secondary. We therefore use this line for disentangling as well. A weighted mean of the poorly constrained He II measurements and the He I λ4009 line yields K₂ = 29 ± 18 km s−1. The RV amplitude of the secondary is thus poorly constrained. We classify the companion as 08 V and estimate it contributes 16% to the total flux. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 441, 09.2 V + B0.5 V has a reported period and eccentricity of P = 6.9 d and e = 0.22. Inspection of the He I lines unambiguously reveals this system as an SB2 binary, albeit the strong nebular contamination of the He I lines makes the disentangling challenging. Because of the presence of a non-negligible secondary, we implement the 2D disentangling technique. The line least contaminated by nebular lines is He I λ4144, shown in Fig. B.34 at RV extremes. A weighted mean of all He I lines but He I λ4471 (due to dominating nebular contamination) yields K₁ = 73 ± 5 and K₂ = 108 ± 8 km s−1. The disentangled spectrum of the secondary, estimated to contribute 21% to the visual flux, shows no evidence for He II absorption. Its spectral appearance matches a B0.5 spectral type. Despite of the short orbital period, analysis of the OGLE light curve of VFTS 441 did not reveal significant periods.

VFTS 475, 09.7 V + B0 V has a reported period of P = 4.1 d and an eccentricity of e = 0.57. The spectral variability of He I lines implies the presence of a non-degenerate secondary in the system (Fig. B.35), and we therefore implement 2D disentangling. Plausible constraints are obtained with the He I λ4388 and He I λ4475 lines, albeit with substantial errors due to the low S/N of the data. A weighted mean of the measurements yields K₁ = 135 ± 33, K₂ = 169 ± 58 km s⁻¹. The He I λ4026 line results in low K₂ amplitudes, which do not agree with the other lines. The reason for this discrepancy is not clear, and could imply that the system is a higher-order multiple. The spectral appearance of the disentangled spectrum off the secondary matches a B0 V spectral type, and we estimate its light contribution at 35%. The OGLE light curve, phased with the orbital period, shows a clear periodic signature, potentially originating in ellipsoidal variations and irradiation effects (Fig. B.36).

VFTS 479, 04.5 V((fc)) + B: has a reported period of P = 14.7 d and an eccentricity of e = 0.31. The spectra are very noisy and contaminated by nebular emission, such that disentangling of the system is challenging (Fig. B.37). Disentangling of He I lines results in K₂ values of the order of 200 km s⁻¹ with very large errors. A weighted mean of the measurements obtained for the He I λ4026, He I λ4144, and He I λ4388 lines yields K₂ = 238 ± 96kms−1, but given the large error, we refrain from providing this value. The disentangled spectrum of the secondary exhibits He I features, but no Hen features. We therefore tentatively classify this system as SB2:, and the secondary as B:, although we cannot exclude that the nebular lines may bias our interpretation. We estimate the light contribution of the secondary at ≈ 5%. With a minimum mass of the secondary is likely a B0 star, though we cannot fully rule out a BH due to the nebular contamination. The OGLE light curve, phased with the orbital period, is shown in Fig. B.38, showing a very faint periodic signature.

VFTS 481, 08.5 V + 09.7: V: has a reported period of P = 142 d and a high eccentricity of e = 0.93. The presence of a non-degenerate companion is barely seen in the wings of strong He I lines such as He I λ4388 (Fig. B.39). Disentangling supports the presence of a non-degenerate companion. The K₂ value is difficult to constrain given the faintness of the companion. A weighted mean of the strong He I lines yields K₂ = 303 ± 50km s−1. The corresponding disentangled spectrum approx. matches an 09.7 V star, with an estimated light contribution of 10%. The OGLE light curve does not reveal any significant periods.

VFTS 514, 09.7 V, has a reported period of P = 185 d and eccentricity of e = 0.41. No clear minimum for K₂ is obtained, regardless of the disentangled line. The spectrum of the secondary appears to be featureless (Fig. B.40), with exception of the Balmer lines, which are strongly contaminated by nebular lines and are therefore not considered. We can rule out a companion contributing more than ≈ 5% to the visual light, corresponding roughly to B3 V. With a minimum mass of M₂ > 5.3 ± 1.8 M_⊙, we cannot rule out that the secondary is a faint non-degenerate star, but given the high BH probability (Table 2), we consider VFTS 514 an O+BH candidate. The OGLE light curve did not reveal any significant frequencies.

VFTS 532, 03.5: V:((f*)) + B III, has a reported period of P = 5.8 d and an eccentricity of e = 0.46. Given the very early spectral class of the primary, its orbital solution and RV amplitude of K₁ = 34.4 km s−1 can be well constrained using the He II lines, which are not impacted by nebular contamination. Similarly to VFTS 386, the He II lines exhibit clear Doppler motion, while the He I lines appear to be relatively static, which indicates that the secondary contributes mainly to He I lines (Fig. B.41). Because of the low S/N, spectral disentangling of various He I lines does not yield a clear minimum for K₂, and the results are generally line-dependent, with the K₂ values lying in the range 40 – 100 km s−1. A weighted mean of the measurements of the He I λ4026, He I λ4388, and He I λ4471 lines yields K₂ = 81 ± 32kms−1. The disentangled spectrum of the secondary matches best a B-type star, but the exact subtype is difficult to retrieve.

Analysis of the OGLE light curve yields a period of P = 5.796223 ± 0.000002 d, and reveals the unmistakable presence of two eclipses in the system (Fig. B.42). This further supports the fact that the system hosts two non-degenerate stars. The eclipse ratio implies that the secondary is much cooler than the primary, which agrees with the estimates spectral type. However, a B2 V star would contribute a mere 1% to the flux, and would exhibit RV amplitudes that are 4 – 5 times larger than derived here. Therefore, the secondary is likely a giant star, and we classify it as such. We note that VFTS 532 is one of the four objects in our sample for which X-rays were detected (Table 1).

VFTS 603, 04 III:(fc) +OB:, was reported as a short period binary of P = 1.76 with an eccentricity of e = 0.11. The primary has a very low reported RV amplitude of llkms−1, and the presence of a non-degenerate secondary is not readily evident (Fig. B.43). However, the orbital solution obtained by Almeida et al. (2017) has a fairly high reduced^² value of 6.9. Indeed, 1D disentangling of the strong He II lines results in a clear residual for the companion that appears distorted. We therefore resorted to 2D disentangling. We focus on the He II lines, including the He II λ4026 line, as the remaining lines are too weak and are dominated by nebular emission. However, we do not get consistent results from all lines. The K₁ amplitude is found consistently to be in the range 30 – 50 km s−1, but K₂ ranges from vanishingly small values for He II λ4542 to ≈ 50 km s−1 for He I λ4026. The latter is the only line that gives plausible results for the secondary, with K₁ =45 ± 8 and K₂ = 49 ± 7kms−1. These amplitudes are peculiarly low for such a short period, resulting in very small minimum masses of M₁ sin³ i = 0.08 M_⊙ and M₂ sin³ i = 0.07 M_⊙, and implying a very small inclination of i ≈ 6°. While the inclination is not likely (the probability of i ≤ 6.3° is 0.6%), there is an appreciable probability for this occurring ones in our sample of 51 stars (≈ 25%).

Given the short period and low RV amplitudes, it is possible that intrinsic pulsations led to a spurious classification of this object as a binary, and specifically an SB2. Hence, we do not readily adopt our solution in Tables 1 2. Higher-resolution data will be needed to verify the nature of this binary system.

Interestingly, despite of the short period, the OGLE light curve does not reveal any clear periodicities beyond integer numbers that are not of astrophysical origin. Moreover, this target is one of the four targets of our sample for which X-rays were detected, albeit at rather low luminosity (Table 1).

VFTS 613, 09 V + 07.5 V, was reported to have a period of P = 69 d and an eccentricity of e = 0.35. However, the RV measurements by Almeida et al. (2017) shows substantial scatter with respect to their derived orbital solution. Inspection of the spectral line variability suggests that the system is in fact SB2, consisting of two similar stars (Fig. B.44). We therefore performed 2D grid disentangling of various lines. The only lines to yield sensible constraints are the He I λ4388 and He II λ4542 lines, and a weighted mean of the measured RV amplitudes yields K₁ = 66 ± 15kms−1and K₂ = 96 + 39kms−1. The spectral appearance of the secondary suggests that, while less massive, it is hotter than the primary, with an estimated spectral type of 07 and visual light contribution of 49%. This may imply that the primary is an evolved giant or supergiant, though the spectra quality does not enable us to establish this empirically. There exists no OGLE data for this object.

VFTS 619, 08: V, was reported to have P = 14.5 d and e = 0.08. The spectra are strongly contaminated with nebular lines, which makes the disentangling challenging. We cannot find any clear minima in the χ² maps of K₂. The disentangled spectra of the secondary plausible K₂ values yields faint signatures for the secondary, but they are likely the result of the nebular lines (e.g. Fig. B.45). We therefore tentatively classify this system as SB1. With a minimum mass of 3.8 ± 0.4 M_⊙, the companion may be a B-type star, a helium star, or a BH. The OGLE light curve did not reveal any significant frequencies.

VFTS 631, 09.7 V, was reported to have a virtually circular orbit and an orbital period of P = 5.4 d. The spectral variability is not readily suggestive of the presence of a non-degenerate companion. Disentangling of the He I λ4388 line (Fig. B.46) suggests a χ² minimum at K₂ = 285 ± 33kms−1, but it is barely significant and therefore disregarded. Disentangling of the other He I lines does not result in a clearly defined minimum. Since the nebular contamination is relatively low, our conservative fiducial simulations suggest that we could probe down companions with brightness ratios down to ≈ 3%. Walborn et al. (2014) estimated M_V = −3.92mag for the system (dominated by the primary), such that our threshold would correspond to a companion brightness of M_V ≈ −0.1 mag, which matches roughly B6-8 V, or ≈ 3 – 5 M_⊙. With a minimum mass of 3.1 ± 0.4 M_⊙, the companion could therefore be a late type B or an early type A star, a helium star, or a BH. A Fourier analysis shows a peak at the orbital period, and the phased light curve shows a periodic pattern at a low amplitude when phased at the orbital period (Fig. B.47).

VFTS 645, 09.5 V, has a reported period of 12.5 d and eccentricity of e = 0.23. The disentangling of various strong He I lines such as He I λ4026 (Fig. B.48) did not result in a well defined minimum for K₂, and the disentangled spectra of the secondary for various K₂ values does not portray any features beyond Banner lines, which may be the result of nebular contamination. We therefore classify this system as SB1. Walborn et al. (2014) estimated a visual absolute magnitude of M_V = −3.86mag for the system, which is dominated by the primary. Adopting our conservative threshold estimate of 3% for the light contribution for an unseen companion, we can rule out dwarfs with spectral types as late as B6-8 V. With a minimum mass of M₂ ≳ 2.6 ± 0.5 M_⊙, the companion could be an A or B-type dwarf, a helium star, a NS or a BH. The OGLE light curve does not show any significant periods.

VFTS 657, 07 II:((f)) + OB:, has a reported period of P = 63 d and eccentricity of e = 0.48. The strong nebular line contamination makes it difficult to determine whether a non-degenerate companion is causing the observed variability (Fig. B.49), and also hampers our ability to disentangle the system robustly. Nonetheless, the He I λ4388 and He I λ4144 lines yield a consistent minimum for K₂, and the weighted mean of the measurements yields K₂ = 169 ± 38 km s−1. The spectrum of the companion best matches an early B-type star, but its line coincide with nebular lines, making the SB2 nature of the system questionable. With a minimum mass of 9.2 ± 2.7 M_⊙, the companion is either a non-degenerate early-type star or a BH. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 702, 08 V(n) + OB + (OB + OB), was reported to be a virtually circular short-period binary with P = 1.98 d. The spectra suffer from low S/N and very strong nebular contamination, hampering our ability to disentangle them robustly. However, Fourier analysis of the OGLE light curve reveals the presence of two eclipsing binaries in the system. The period of the first one closely matches the spectroscopic period derived by Almeida et al. (2017), refined to P = 1.981546 ± 0.000030 d (Fig. B.50). The second period is P = 2.932818 ± 0.000025 d. In Fig. B.50, we also show the light curve phased with this period, after subtracting the 1.98 d signal. Given the short periods, it is clear that the system must be a doubly-eclipsing quadruple system. From the comparable strengths of the eclipses, it is likely that all components are OB-type stars. Examination of available HST images does not reveal any nearby source within 0.1” (roughly 5 000 au), suggesting that the two binaries are gravitationally bound. Unfortunately, the data are not sufficient for disentangling the system, let alone with four components.

VFTS 733, 07.5 V + Bl II, was classified as an 09.7p binary with a period of P = 5.9 d and a virtually-circular orbit. Inspection of the line profiles readily reveals that this system is SB2. In fact, the He II lines belong almost solely to hotter primary, while the He I and a multitude of metal lines belong to the cooler secondary. Therefore, deriving the nature and orbital solution of this system does not require spectral disentangling, as the RVs can be well measured from isolated lines. Here, we use the He II λ4542 and the Si III λ4553 for the hotter primary and cooler secondary, respectively (Fig. B.51). The phased RV measurements and orbital solution are shown in Fig. B.52. We find K₁ = 39.4 ± 1.9 km s−1 and K₂ = 102.2 ± 0.9 km s−1 for the primary and secondary, respectively. The less massive secondary appears to contribute roughly half of the flux in the visual, suggesting that it is evolved (luminosity class III-I). This is supported by the spectral appearance of the secondary. The low minimum masses of M₁ sin³ = 1.26 M_⊙ and M₂ sin³ i = 0.49 M_⊙ imply a low inclination of i ≈ 20°. Eclipses are therefore not to be expected. The OGLE light curve shows a periodic pattern when phased at the orbital period (Fig. B.53).

VFTS 736, 09.5 V + B:, has a reported period of P = 68.8 d and eccentricity of e = 0.09. Disentangling of various He I lines does not yield a consistent solution for K₂, and the spectral variability is not readily suggestive of a non-degenerate companion (Fig. B.54). The disentangled spectrum of the secondary for plausible K₂ shows very faint He I signatures for strong lines, but these may be the result of nebular contamination. As M₂ > 3.9 ± 0.5 M_⊙, the companion could either be an B- or A-type star, a helium star, or a black hole. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 743, 09.5 V((n)), has a reported period of P = 14.9 d and a virtually circular orbit. The spectral variability and disentangled spectrum of the secondary for various K₂ values do not support the presence of a non-degenerate companion in the system (Fig. B.55). Some faint He I features, as well as Balmer features, are seen in the disentangled spectrum of the secondary, but their appearance suggests that they originate primarily in nebular contamination, and we therefore classify the system as SB1. Given the minimum mass of M₂ > 2.1 ± 0.3 M_⊙, the nature of the secondary is virtually unconstrained. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 750, O9.5 V + B:, is reported to be an eccentric long period binary with P = 417 d and e = 0.78. The spectral variability is not readily suggestive of the presence of a non-degenerate companion (Fig. B.56). When disentangling for plausible K2 values in the range 2–5 K₁, faint features of He I can be seen. Since the nebular contamination is modest in this star, it is likely that these features are real. However, due to the low-amplitude motion of the primary, the relatively poor Doppler coverage, and the faintness of the potential companion, we cannot retrieve an unambiguous value for K₂. We therefore tentatively classify the companion as B: and the system as SB2:. With a minimum mass of M₂ > 5.6 ± 1.2 M_⊙, the companion could be a non-degenerate object or a BH. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 769, O9.7 V, has a short period of P = 2.4d and a virtually circular orbit. The spectral variability of strong He I lines such as He I λ4388 (Fig. B.57) is not readily suggestive of the presence of a non-degenerate companion. The disentangled spectra of the secondary for various K₂ values are suggestive of Balmer absorption and hints of He I lines, but we cannot exclude that they are the result of nebular contamination. With a minimum mass of M₂ > 1.9 ± 0.2 M_⊙, the nature of the secondary is unconstrained. Analysis of the OGLE light curve reveals a significant peak at the orbital period of P = 2.365644 d and its harmonics. The phased light curve is shown in Fig. B.58.

VFTS 779, B1 III, has a reported period of P = 59.9 d and a virtually circular orbit. The spectral variability of the He I does not readily suggest the presence of a non-degenerate secondary (e.g. Fig. B.59). Disentangling of various He I lines results in nearly-flat ;χ²(K₂) maps, though an insignificant tendency towards values in the range K₂ ≈ 2 – 4 × K₂ is noted. The disentangled spectrum of the secondary for such K₂ values shows very faint He I absorption in some of the lines, though it is not clear whether this is due to low levels of nebular contamination. Given the relatively high quality of the data in this case, we can exclude companions contributing more than 2 – 3% to the total light, which would correspond roughly to a B3. We classify this system as SB1. As M₂ > 3.9 ± 0.6 M_⊙, the companion is either a mid- to late-type B dwarf, a He star, or a BH. Given the BH probability (Table 2), we consider this system a good B+BH candidate. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 802, O7 V: + O8 Vn, has a reported period of P = 182 d and an eccentricity of e = 0.60. The presence of a non-degenerate companion is readily seen in a multitude of He lines as a shallow, broad feature (Fig. B.60). Given the strong contribution of the secondary, we perform 2D Disentangling of the strong He I and He II lines. A mean average yields: K₁ = 70.4 ± 1.9 km s⁻¹, K₂ = 76 ± 14 km s⁻¹. The RV amplitude of the secondary is more difficult to constrain due to the substantial rotational broadening of its spectrum. The companion is estimated to contribute 36% of the visual flux, and is classified as O8. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 810, O9.7 V + B1 V, has a reported period of P = 15.7 d with an eccentricity of e = 0.68. The presence of a non-degenerate companion was already recognised by (Almeida et al. 2017), who however could not measure its RVs. The lines of the secondary become deblended in certain epochs, making its presence readily evident (Fig. B.61). 2D Disentangling of the strong He I lines yield a weighted mean of K₁ = 108 ± 9 km s⁻¹ and K₂ = 162 ± 27 km s⁻¹. The companion is estimated to contribute 29% of the visual flux, and is classified as B1. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 812, O4 V((fc)), has a reported period of P = 17.3 d and eccentricity of e = 0.62. No sign for a companion is seen in He II lines (e.g. Fig. B.62), and He I lines are contaminated by nebular lines, hampering our ability to use them robustly. No clear value for K₂ is retrieved from disentangling. The disentangled spectra of the secondary obtained for plausible K2 values show Balmer lines and some hints of He I lines, but not sufficiently clear to be able to classify this system SB2 unambiguously. With M₂ ≥ 5.1 ± 0.7 M_⊙, the companion may be a late-type OB star or a BH. Analysis of the OGLE light curve did not reveal significant frequencies.

VFTS 827, B1.5 III, has a reported period and eccentricity of P = 43 d and e = 0.24. The spectral line variability is not suggestive of the presence of a second non-degenerate companion (e.g. Fig. B.63). However, disentangling of the strong He I λ4471 favours K₂ values in the range K₂ < 70 km s−1. Disentangling for various K₂ values does not result a clear stellar spectrum for the secondary. However, strong Balmer lines are seen in the disentangled spectrum of the secondary. While the nebular contamination in this case is very low, the narrow widths of the Balmer lines in the disentangled spectrum of the secondary raise the suspicion that they originate in nebular contamination. We therefore classify the system as SB1:. With a minimum mass of M₂ > 2.5 ± 0.5 M_θ, the companion might also a late B-type star or an A-type star, and we cannot fully rule out that it is a NS or a BH. No OGLE light curve is available for this object.

VFTS 829, B1.5 III, has a reported period and eccentricity of P = 203 d and e = 0.27. The spectral line variability is not suggestive of the presence of a second non-degenerate companion (e.g. Fig. B.64), and disentangling of various He I results in relatively flat;χ²(K₂) maps. The disentangled spectrum of the secondary for K₂ values in the range 2 - 10 x K₁ is fully free of features with the exception of weak Balmer absorption that could be the result of nebular-line contamination. Given the low levels of nebular-line contamination and the relatively high S/Ns of the spectra, we classify this system tentatively as SB1:. However, the low minimum mass of M₂ > 2.0 ± 0.5 M_θ leaves the nature of the secondary virtually unconstrained.

VFTS 887, O9.7: V: + O9.5: V, has a reported period of P = 2.7 d and a virtually circular orbit. The spectral line variability (Fig. B.65) clearly reveals the presence of two components in the spectrum. We therefore utilise 2D disentangling on various He I and He II lines, which yield consistent results. A weighted mean yields K₁ = 105 ± 7 km s−1 and K₂ = 98 ± 12km s−1, implying that the object identified as the primary by Almeida et al. (2017) is the slightly less massive component. The disentangled spectra correspond to O9.7 V and O9.5 V spectral types, respectively, and are estimated to contribute 45% to the visual flux. Despite the short period, the OGLE light curve does not show a clear period pattern (Fig. B.66).

Appendix C Disentangled spectra

This section and its Figs. C.1–C.50 provides the disentangled spectra for each of our targets.

Appendix D Rotation velocities and minimum inclinations

Table D.1 provides v sin i values measured from the disentangled spectra using FWHM calibrations provided by Ramírez-Agudelo et al. (2015), as well as minimum inclinations using Eq. (2). We note that this relies on the assumption that the orbital axis is aligned with the rotational axis of the stars, which may not be valid for binaries with compact companions. Values should be considered as accurate to within 50 – 100 km s−1.

Appendix E Binary detection probability

We evaluated the numerical binary detection probability, f_det(q), provided by Sana et al. (2013), at ten individual points, and extended the function to f_det(0) = 0. We find that a sixth order polynomial is needed to adequately fit the data and exhibit a smooth behaviour in the range 0 ≤ q ≤ 1: (E.1)

where a = − 17.178762, b = 62.839749, c= −92.161434, d=69.536844, e=−28.963718, and f=6.7669361. The fit reaches an accuracy of 0.5% or better in the range 0 ≤ q ≤ 1 (Fig. E.1).

References

---

All Tables

All Figures