© The Authors 2025

1. Introduction

Massive stars go through complex evolutionary scenarios that are significantly influenced by their high multiplicity fraction (Mason et al. 2009; Chini et al. 2012; Sana et al. 2014; Kobulnicky et al. 2014; Dunstall et al. 2015; Moe & Di Stefano 2017; Ramírez-Tannus et al. 2024) and frequent binary interactions (Sana et al. 2012; de Mink et al. 2014). Spectroscopic studies in the Galaxy and the Large Magellanic Cloud (LMC) consistently indicate that more than 50% of O- and B-type stars are part of binary or multiple systems within the orbital period range where interactions are expected (≲3500 d). Including longer orbital periods, the overall multiplicity fraction can exceed 90% for the O-type stars (Sota et al. 2014; Bordier et al. 2022). Sana et al. (2012) found an intrinsic multiplicity fraction (after correction for observational biases) of f_mult = 69 ± 9% for a sample of 71 O-type stars in young Galactic clusters. In the Cygnus OB2 association, Kobulnicky et al. (2014, and references therein) presented orbital solutions for 48 massive stars with spectral types ranging from O5 to B2, resulting in an intrinsic multiplicity fraction near 55% from a sample of 45 O-type and 83 B-type stars. More recently, Banyard et al. (2022) studied the B-type content of the young Galactic cluster NGC 6231, determining f_mult = 52 ± 8% from 80 stars covering the entire B-type spectral type range (B0-B9). While further multiplicity studies have been conducted in the Milky Way (MW), many have not been corrected for observational biases (e.g. Barbá et al. 2017; Berlanas et al. 2020; Ritchie et al. 2022), have targeted narrow spectral ranges (Abt et al. 1990), focused on peculiar types (e.g. chemically peculiar stars; Schöller et al. 2010), or addressed larger separations (Mason et al. 2009; Sana et al. 2014). These factors hinder direct comparison of multiplicity fractions due to the different observational strategies and biases.

In the 30 Doradus (30 Dor) region of the LMC (Z_LMC = 0.5 Z_⊙), the VLT-FLAMES Tarantula Survey (VFTS; Evans et al. 2011) performed multiplicity studies on 360 O-type (Sana et al. 2013) and 408 B-type stars (Dunstall et al. 2015, of which 361 were main-sequence (MS) and giant stars), resulting in binary fractions¹ of f_mult = 51 ± 4% and f_mult = 58 ± 11%, respectively. The consistency across these studies suggests that the multiplicity properties of massive stars may be independent of metallicity. Furthermore, the distribution of orbital periods is remarkably similar across most of these samples (Almeida et al. 2017; Villaseñor et al. 2021; Banyard et al. 2022), suggesting that O- and B-type stars in multiple systems are formed with similar efficiency and orbital configurations within the metallicity ranges studied so far, emphasising the need for multiplicity studies of massive stars in lower metallicity environments.

In the case of solar-type binaries, a strong anti-correlation between binary fraction and metallicity has been found (Badenes et al. 2018; Moe et al. 2019; Offner et al. 2023). Moe et al. (2019) computed intrinsic binary fractions for five different observed samples comprising spectroscopic and eclipsing binaries (EBs), and found a combined close binary fraction that drops from 53 ± 12% at [Fe/H] = −3.0 to 10 ± 3% at [Fe/H] = +0.5. This stands in stark contrast to the findings for OB-type stars (see also Moe & Di Stefano 2013) and also to solar-type stars in wider orbits (a ≳ 200 au; El-Badry & Rix 2019). Both Moe et al. (2019) and El-Badry & Rix (2019) agree that the formation of wide binaries is governed by turbulent core fragmentation which is independent of metallicity, but, in the case of close binaries, the main formation mechanism is disc fragmentation up to separations of a ∼ 200 au. Moe et al. (2019) suggest that the difference between the sharp increase in the binary fraction of solar-type stars towards low metallicities and the close-to-constant fraction observed in massive stars lies in the disc of massive protostars. At solar metallicity, these discs are already highly unstable and have a high probability of fragmentation, leaving little room to further increase multiplicity at reduced metallicities.

The findings in the MW and LMC suggest a universality in the multiplicity properties of massive stars, and motivate further investigation on whether trends persist at lower metallicity. Since the LMC metallicity is only half solar, it is essential to perform multiplicity studies of massive stars at even lower metallicity. In this context, the Small Magellanic Cloud (SMC) is a unique environment due to its proximity (62 kpc; Graczyk et al. 2020), low metallicity (Z_SMC = 0.2 Z_⊙) representative of the high-redshift universe (Nakajima et al. 2023), and high content of massive stars (Evans et al. 2004). Furthermore, theoretical predictions indicate that the winds of massive stars at SMC metallicity are strongly reduced (Vink et al. 2001; Björklund et al. 2021), with empirical evidence showing that they might be even weaker than predicted (Ramachandran et al. 2019; Rickard et al. 2022). If there is a connection between metallicity and the multiplicity properties of massive stars, the SMC offers a prime environment for testing the universality of these properties. For other low-metallicity galaxies in the Local Group, resolving a significant sample of massive stars with the current generation of telescopes is not possible due to their greater distances.

The Binarity at Low Metallicity (BLOeM; Shenar et al. 2024) campaign addresses this knowledge gap by focusing on the multiplicity properties of massive stars (M_ini ≳ 8 M_⊙) within the SMC. BLOeM has observed 929 stars in eight different regions of the SMC. The sample has been divided into subsamples based on their spectral types, emission properties, and luminosity classes, enabling an interpretation of the results within the context of each subsample. The five samples are: the O-type stars (Sana et al. 2025); the B-type supergiants (spectral types B0-B3, Britavskiy et al. 2025); the B-type dwarfs and giants (this work); the cooler supergiants (spectral types later than B3, Patrick et al. 2025); and the emission-line stars (Oe/Be, Bodensteiner et al. 2025).

In this study, we present the first detailed analysis of the multiplicity properties of a large sample of early B-type dwarf and giant stars at the low metallicity of the SMC. In Sect. 2, we introduce the sample of B-type stars and their main properties. Section 3 details our methodology to measure radial velocities (RVs) and our criteria to classify spectroscopic binaries. Our results on the multiplicity properties are presented in Sect. 4, and our search for periodicity is described in Sect. 5. We discuss our results in a broader context in Sect. 6, and finally provide our conclusions in Sect. 7.

2. The sample

The B-type dwarfs and giants sample comprises 309 stars with spectral types ranging from B0 to B2.5 and luminosity classes V to III². Spectral types for the full sample were presented in Shenar et al. (2024). Specifically, our sample consists of all early B-type “non-supergiant” objects, including MS (class V), subgiant (class IV), and relatively unevolved giant stars (class III). We expect our sample to primarily consists of MS stars undergoing core-hydrogen burning and stars that have just recently exhausted hydrogen in their cores. Bright giants and supergiants (classes II and I) form a separate early B-type supergiants subsample, which is the focus of Britavskiy et al. (2025). Hereafter, we refer to our sample as the III/V sample. This classification aligns with that used for the 30 Dor B-type stars, where Dunstall et al. (2015) defined “unevolved” stars as those with surface gravity log g > 3.3 dex and analysed their multiplicity properties separately from the supergiant sample. However, unlike the 30 Dor unevolved sample, our sample is mainly composed of giant stars (see Sect. 4.2). In the case of Be stars, we assume that most of them are binary interaction products (Bodensteiner et al. 2020b; Dallas et al. 2022; Dodd et al. 2024), therefore they will be the subject of a separate work (Bodensteiner et al. 2025).

Our limit on the spectral type is given by the BLOeM survey magnitude cut-off of G < 16.5 mag. A spectral type of around B2–B3 roughly corresponds to the mass limit above which we expect stars to be core-collapse-supernova (CCSN) progenitors (M ≈ 8 M_⊙) assuming single-star evolution. In this way, we ensure that the sample is representative of the population of CCSN progenitors that might become gravitational-wave sources.

BLOeM has currently obtained and reduced nine epochs from the first semester of observations using the FLAMES instrument (Fibre Large Array Multi Element Spectrograph; Pasquini et al. 2002) in GIRAFFE mode (R = 6200) on the UT2 unit telescope of the Very Large Telescope (VLT) at the European Southern Observatory (ESO) in Chile. For ten stars in our sample the number of epochs is below nine, fluctuating between three and seven epochs, and in five additional stars one of the epochs has too low a signal-to-noise ratio (S/N ≲ 15), thus only eight epochs were used. Additional epochs are expected by the end of 2024, with a total of 25 epochs planned by the conclusion of the campaign in September 2025. For a complete overview of the sample selection, data reduction, and spectral classification, see Shenar et al. (2024).

Each field’s central coordinates and the modified Julian dates (MJDs) per observation are provided by Shenar et al. (2024, Table A.1). Here, we present the total baseline of observations and minimum cadence for each field in Table 1. We have a mean observational baseline of 45 d and a minimum cadence of 1 to 2 d. This cadence allows us to constrain most short-period orbits, except for unevolved contact binaries with periods of less than P_orb ∼ 1 d. However, these are likely to be EBs and therefore detectable by the Optical Gravitational Lensing Experiment (OGLE; Pawlak et al. 2016) through their light curves (see Sect. 5.1).

The total baseline allows us to recover periods of up to about 45 d for most of the fields, and close to 60 d for Field 1. Given that about 75% of binaries with early-B primaries and orbital periods below 500 d are expected to have periods shorter than 40 d (Villaseñor et al. 2021), the cadence of the first nine epochs can effectively detect these short-period binaries while also offering a significant total baseline. However, while nine epochs provide sufficient data to robustly estimate the multiplicity fraction and identify interesting systems, they are at the threshold of what is necessary to constrain orbital periods. Accurate periods, orbital solutions, and the distributions of orbital properties for the full sample will be addressed once all 25 epochs are assembled.

Figure A.1 shows the distribution of the B-type III/V sample in each field observed by BLOeM. Field 6 contains the largest number of objects, 71 stars, followed by Field 8 with 53, whereas Field 2 is the least populated with 18 stars. Notably, Field 4 includes NGC 346, the brightest giant H II region in the SMC with a rich OB population (Massey et al. 1989; Evans et al. 2006), however, we have excluded stars in the centre of the region to avoid crowding. Other less studied clusters/H II regions are NGC 371/N76 close to the centre of Field 1 (Ripepi et al. 2014), NGC 465/N85, NGC 460/N84 and NGC 456/N83 in Field 6 (Dapergolas et al. 1991; Caplan et al. 1996) with most stars coming from NGC 465, and NGC 267/N22 is at the centre of Field 3, surrounded by several other H II regions (N25, N26, N21, N23; Testor et al. 2014), and NGC 261/N12A in the top right (Sano et al. 2019).

The distribution of spectral types per field is shown in Fig. 1. The overall distribution is very similar to that found by Ramachandran et al. (2019) in the supergiant shell SMC-SGS1 in the Wing of the SMC, a region that was not covered by BLOeM. Our sample is dominated by spectral types B1.5 and B2, corresponding roughly to masses between 9 and 12 M_⊙. There is a good range of spectral types across different fields, but in Fields 7 and 8, the number of B1.5 and B2 seem to be overabundant with respect to the earliest stars (B0–B0.2), potentially suggesting an older population. We computed the ratio of latest-to-earliest type stars for these spectral types (i.e. N(B1.5–B2)/N(B0–B0.2)), finding ratios of 4.4 and 18.0 for Fields 7 and 8 respectively, and a mean of 2.9 for the other six fields. Performing a Z-test returns that these differences are significant, with a z-score of 5 (Field 7) and 49 (Field 8) and p-values of 5.7 × 10⁻⁷ (Field 7) and 0.0 (Field 8). This is consistent with Fields 7 and 8 being areas without prominent star-forming regions (see Fig. A.1 and details in the main text). In contrast, Field 3, which exhibits several H II regions, has the lowest ratio of latest-to-earliest stars with 1.5. We note, however, that an uncertainty of one spectral subtype should be expected. Such uncertainty could at least partially explain the lack of B0.2 and B0.5 stars. These are stars that show weak He II, and since He IIλ4686 (stronger than He IIλ4542) is not covered by the FLAMES LR02 setup, they are more difficult to classify. Given the small number of stars in some of the fields, this variation could also result from sampling effects or constraints in the observational design.

Fig. 1. Spectral types of the B-type III/V sample, colour coded by field. The number of stars per field is shown in parentheses.

Examining the star formation history of the SMC (Harris & Zaritsky 2004, see their Figures 3 and 6) we can see that the SMC has had a relatively constant star formation rate (SFR) from 2.5 Gyr until approximately 40 Myr ago, when the SFR started to decrease. All BLOeM fields have experienced star-forming episodes in the last 40–10 Myr according to the SFR maps of Harris & Zaritsky (2004). Fields 6, 7, 8 are located in regions that experienced a strong SFR 10 Myr ago, but not in the last 6–4 Myr, whereas Field 3 can be associated with strong star formation at 10 and 6 Myr, which is consistent with our suggestion that Field 3 might contain a slightly younger population from the spectral types distribution analysis. Other fields that show a high SFR in the last 4 Myr from the maps of Harris & Zaritsky (2004) are 1 and 4, containing NGC 371 and NGC 346 respectively. We note however that our sample excludes O-type and supergiant stars, so while the presence (or absence) of early B stars can give a qualitative sense of recent star formation, it does not capture the youngest or most evolved populations. Thus, care must be taken when associating these stars alone with specific star-formation epochs. A quantitative analysis of the BLOeM stars by Bestenlehner et al. (in prep.) finds median ages for Fields 1–8 in the range of 8–13 Myr, consistent with our identification of star-formation activity in the last 10 Myr, but reflecting a wider range of masses and evolutionary stages than those considered here.

3. Radial velocities

To measure RVs we have used the same tools as for the B-type stars in 30 Dor (Villaseñor et al. 2021), which are now part of the MINATO (Massive bINaries Analysis TOols, Villaseñor et al. 2023) python repository³, in a module called ravel. Briefly, it is a spectral line profile fitting tool, which uses Gaussian or Lorentzian (in the case of Balmer lines) profiles to fit a set of spectral lines. In the early B-type stars’ spectra, there is a rich set of spectral lines covered by the wavelength range of FLAMES; He Iλ4026, 4144, 4388, 4471, and the Balmer lines Hδ and Hγ. There are also a few metal lines that are usually used in the analysis of later-type stars such as Si IIλ4128-31 and Mg IIλ4481. However, these are very weak in the stars in our sample due to their spectral range and due to the low metallicity of the SMC, thus we have not included them in our RV measurements. For the earliest types, He IIλ4542 is commonly present, but it is weaker than the He I lines and therefore provides less accurate RVs. This is similar for the case of Si IIIλ4553, which is usually too shallow for precise RV measurements as a result of the SMC’s low metallicity. All mentioned lines are fitted, but ravel selects the ones used to compute a weighted mean RV for each epoch based on the errors of the Gaussian or Lorentzian fit using a median-absolute-deviation (MAD; Hampel 1974) test. For this reason, in most cases, only the He I and Balmer lines (Hγ, Hδ) are used, and Si II and Mg II lines did not pass the MAD test for any of the systems. An example fit for a typical SB1 system is shown in the top plot of Fig. 2. Two He I lines and one Balmer line are displayed to illustrate the use of Gaussian and Lorentzian profiles.

Fig. 2. Top: spectral-line fits for the SB1 system BLOeM 4-097. Three different spectral lines are shown for epochs closest to quadrature. Observed spectrum (black) is overlaid by best fit and 1σ confidence levels (red); dashed grey profiles indicate initial guesses. Bottom: Similar to above, but for SB2 system BLOeM 1-055. Two He I lines at three epochs (near quadrature and conjunction) are shown. The dashed red line marks the mean SMC velocity of 172 km s−1 (Evans & Howarth 2008), with a ±200 km s−1 deviation (dashed yellow lines).

3.1. Double-lined spectroscopic binaries

A significant improvement in our methodology, compared to the approach used for the 30 Dor population by Villaseñor et al. (2021), involves the fitting of double-lined spectroscopic binaries (SB2s). The previous line-by-line fitting method was effective for single-lined binaries (SB1s) but struggled with SB2s, particularly when the binary components had similar flux ratios. This often led to misidentifications and the rejection of affected epochs by the code, severely affecting the number of epochs used in the analysis.

To address these challenges, we have implemented a simultaneous fitting approach using probabilistic programming with NumPyro (Bingham et al. 2019; Phan et al. 2019). Similar to the technique used by Sana et al. (2013) and Almeida et al. (2017), our method simultaneously fits all selected spectral lines and epochs. The approach allows the RV of each component to vary across epochs; however, all lines from the same component share the RV at a given epoch, while maintaining consistent line profiles (amplitude and width) for each spectral line across epochs. By modelling the data in this simultaneous way, we improve our ability to distinguish and separate the spectral components, even when their flux contributions are comparable.

The probabilistic programming framework offers several advantages (e.g. Nightingale et al. 2021): (i) the direct incorporation of uncertainties in the model parameters, which enhances the robustness of our fits without introducing biases; (ii) comprehensive uncertainty quantification through full posterior distributions for all fitting parameters, enabling a more complete understanding of parameter uncertainties and their correlations; (iii) efficient data handling by facilitating the simultaneous modelling of multiple spectral lines and epochs, improving computational efficiency and the accuracy of parameter estimations.

These enhancements lead to more reliable identification of SB2 components, improvements in the determination of RV uncertainties, and reduced the need to reject epochs due to misidentifications. However, challenges remain, particularly in cases with limited epochs, contributions from third objects, or significant profile variability. Future work will aim to refine our models to better address these complexities. The bottom plot of Fig. 2 demonstrates the SB2 fitting method by showcasing a double-component fit for the He Iλ4144 and He Iλ4388 lines of star BLOeM 1-055. The central panel (epoch 6) displays a phase near conjunction, where the two components are fully blended. In contrast, epochs with higher RV variations (e.g. epochs 5 and 7) provide better constraints on the line profiles of the individual components, enabling a reliable fit even for highly blended epochs.

3.2. Binary criteria

After weak or absent lines are omitted and epochs with a very low S/N (S/N ≲ 15) are rejected (see details in Villaseñor et al. 2021), we measure the peak-to-peak RV difference (ΔRV), which is the maximum difference between any two of the available nine epochs. These measurements are used to classify stars as spectroscopic binaries based on the criteria used by Sana et al. (2013), Dunstall et al. (2015), Bodensteiner et al. (2021), Banyard et al. (2022). The first criterion is based on the peak-to-peak RV variability. Stars with ΔRV > C are classified as binary candidates, where C is a specified threshold. Dunstall et al. (2015) adopted a value of C = 16 km s⁻¹ for the B-type stars in the LMC, based on their analysis of the binary fraction versus ΔRV threshold, which showed a distinct break at 16 km s⁻¹ for the B supergiants in their sample. Given that this feature is absent from our results, we adopt a value of C = 20 km s⁻¹, consistent with other studies (Sana et al. 2013; Bodensteiner et al. 2021; Banyard et al. 2022).

Furthermore, it is well known that B-type stars often exhibit pulsations (Aerts et al. 2009; Bowman 2020). These pulsations can cause intrinsic RV variability with amplitudes peaking in the 10–20 km s⁻¹ range (Stankov & Handler 2005; Hey & Aerts 2024). By adopting C = 20 km s⁻¹, we aim to exclude most stars displaying intrinsic variability caused by pulsations, thereby reducing contamination in our binary classification.

The second criterion, introduced by Sana et al. (2013), classifies a system as a spectroscopic binary if, in addition to the first criterion, any pair of its RV measurements satisfies the following condition:

$$\begin{array}{c}\hfill {\sigma }_{\mathit{ij}}=\frac{|v_i-v_j|}{\sqrt{{\sigma }_i^2+{\sigma }_j^2}}>4,\end{array}$$(1)

where v_i and v_j are the individual RV measurements taken at two different epochs, and σ_i and σ_j are the uncertainties associated with each RV measurement. We define σ_detect = max(σ_ij), representing the detection significance, which quantifies how significantly different two RV measurements are relative to their combined uncertainties. Since the RV uncertainties are critical for the determination of σ_detect, we test their reliability in Appendix B. If either of the two criteria are not met, we classify the star as a RV variable (RV var) if σ_detect > 2, otherwise, we do not consider the variability significant and classify those stars as RV constant (RV cst). Stars classified as RV var inevitably include unidentified binaries in long-period orbits or with low-mass companions, in addition to stars showing intrinsic RV variability (such as pulsations or atmospheric variability), which may or may not be truly single stars. For this reason, we correct for the observational biases and compute the intrinsic binary fraction (see Sect. 4).

For the binary classification procedure, we have first identified SB2 systems by visual inspection, and flagged them as SB2s so that ravel computes their RVs adequately. With the RV measurements computed, the binary classification is done automatically according to the two criteria, but some classifications were reevaluated by visual inspection of the fits to spectral line profiles and refitted. In the same way, we have classified three systems as SB3: BLOeM 5-062, BLOeM 6-062, and BLOeM 7-057. We base these classifications on the presence of one or two narrow cores in the He I lines, with broad profiles on both sides of the narrow lines, or a clear SB2 that cannot be successfully fitted with two Gaussian profiles (see examples in Appendix C.1). Furthermore, we found potential signals of contribution from a third star in the spectra of 15 SB2 systems. These have been flagged in the comments of Table D.1 (see also notes in Appendix C.2).

4. Multiplicity fraction

Figure 3 shows the distribution of ΔRV, for SB1, SB2, SB3, RV var, and RV cst stars. We classified 46 objects as RV cst, of which 42 have ΔRV below 10 km s⁻¹. Among a total of 110 RV var stars, 13 show ΔRV higher than 20 km s⁻¹ but do not meet our second binary criterion. All these 13 systems have low S/N, and most of them exhibit line-profile variability that might be attributable to undetected companions or pulsations, which can mimic or blur binary signatures. We also see eight SB2 systems with ΔRV < 50 km s⁻¹, these systems have been clearly identified as SB2 systems through our visual inspection, but their RVs might be affected by the contribution of a third star. In four of these cases, the optically dominant star in the spectrum presents the largest RV variation, suggesting it is the less massive star in the system. We define ΔRV as that of the currently more massive star (that with the smaller RV variations in SB2/3 systems), but in cases where mass transfer has taken place, the dominant star in the spectrum might be the donor of an Algol binary (MS stars with very high luminosity to mass ratio; van Rensbergen et al. 2011; Sen et al. 2022) or a low-mass bloated stripped star (e.g. Shenar et al. 2020; Bodensteiner et al. 2020a; Villaseñor et al. 2023; Ramachandran et al. 2024), much less massive than the gainer.

Fig. 3. Maximum peak-to-peak RVs for our sample, colour-coded by binary status. We define RV cst stars as those with peak-to-peak RV variability not reaching 2σ significance, which we consider as constant or non-variable. Whereas stars that surpass the 2σ level but do not meet both multiplicity criteria are considered variable (RV var). ΔRV values are grouped in bins of 10 km s−1.

Given that our sample includes giant stars (luminosity class III) it is certainly possible that interactions might have already occurred. A lower surface gravity value, or analogously a larger radius, will affect the minimum orbital period before interaction (P_crit) with consequences for the distribution of ΔRV_max and the multiplicity fraction (Badenes et al. 2018). To quantify the effect of including moderately evolved stars, we compute P_crit (following Mazzola-Daher et al. 2022) for a range of primary masses (M₁ = 6 − 16 M_⊙), surface gravities (log(g/cm s⁻²) = 3.0 − 4.5), and mass ratios (q = 0.1 − 1), representative of our sample of B-type stars. From P_crit, we compute the allowed ΔRV_max values for a range of orbital inclinations between 20 and 90 deg. Figure 4 (top panel) shows how P_crit increases significantly towards lower log(g) values for a fixed mass ratio of q = 0.5, illustrating that systems with orbital periods below 5 d would have interacted once the primary reaches log(g)≈3.5. Using q = 1 would shift the critical periods roughly by -0.2 dex in log(g). In the bottom panel we show the corresponding maximum peak-to-peak RV amplitude (ΔRV_max = 2K₁) for a typical inclination of i = 60°. Even at moderate inclinations and q = 0.5, the minimum ΔRV_max remains around 100 km s⁻¹ for most of our parameter space, safely above the 20 km s⁻¹ threshold used in our binary criterion. Only in extreme cases of very low mass ratios (q ≲ 0.1) and low orbital inclinations (i < 40°) would ΔRV_max fall below 20 km s⁻¹ (dashed white contours). Therefore, even though the larger radii of giant stars shift the shortest detached periods to higher values, in most cases the resulting RV amplitudes remain well above our detection threshold.

Fig. 4. Top panel: Contour plot showing the critical period (Pcrit) as a function of primary mass (M1) and surface gravity (log(g)) for systems with i = 60° and q = 0.5. Dashed white contour lines highlight selected Pcrit values (1, 2, 5, and 10 days). Bottom Panel: As the top panel, this time displaying the corresponding change in RV (ΔRVmax) computed at Pcrit. Dashed white curves mark the contours where ΔRVmax = 20 km s−1 at different inclinations and q = 0.1.

4.1. Observed multiplicity fraction

According to the binary criteria, we have classified 153 stars as spectroscopic binaries, 91 SB1s, 59 SB2s, and 3 SB3s. To illustrate the distribution of our classifications based on the two criteria, we show a 2D diagram in Fig. 5. There are three SB2 systems that do not reach the 4σ significance threshold, these are BLOeM 2-023, BLOeM 3-020, and BLOeM 4-056. All these systems have been flagged as potential binary interaction products, and as discussed previously, the more massive star displays smaller RV variations and larger RV uncertainty. Because of their clear SB2 status we have included these systems in the determination of the multiplicity fraction.

Fig. 5. Spectroscopic binaries, RV var, and RV cst stars classified according to our binary criteria. Our adopted thresholds of σdetect > 4 and C > 20 km s−1 are shown with dashed grey lines. Green region highlights systems fulfilling both criteria.

Here we defined the multiplicity fraction as that of the B-type primaries in our sample. Due to the unknown nature of most of the companions, we do not include them in our binary statistics. The multiplicity fraction quoted henceforth is effectively a close binary fraction that takes into account binaries at separations a < 10 au, or log(P_orb/d) < 3.5, which roughly corresponds to the separations that will allow for binary interactions to happen. Consequently, this is the range of orbital periods for which observational biases are accounted for when determining the intrinsic binary fraction (see Sect. 4.3). The 153 binaries in our sample yield an observed fraction of $f_{\mathrm{mult}}^{\mathrm{obs}}=50\pm 3\%$, where the uncertainty represents the standard error of a proportion based on binomial statistics, corresponding to a 68% confidence interval. This is a significantly larger observed fraction than that found for the 30 Dor “unevolved” stars (25 ± 2%; Dunstall et al. 2015). However, the detection probabilities of both surveys need to be carefully taken into account for a meaningful comparison. For example, the nine epochs analysed in this work provide an advantage over the six LR02 VFTS observations considered by Dunstall et al. (2015). On the other hand, VFTS was much more sensitive to moderate-to-long period binaries due to its longer observational baseline of 10–12 months (and up to 22 months for 31 targets; Dunstall et al. 2015). We further discuss the differences between the samples of B-type stars and the implications for bias correction in Sect. 4.3.

It is, however, possible to compare the observed multiplicity of our sample with that of the BLOeM O-type stars. Sana et al. (2025) determined an observed binary fraction of $f_{\mathrm{mult}}^{\mathrm{obs}}=45\pm 4\%$ for the 139 O-type stars in the BLOeM sample, which is in agreement with our findings for the B-type stars. This is interesting because the B-type population may suffer from different biases compared to the O-type stars due to their fainter magnitudes and lower S/N spectra. These biases result in a detection probability that drops at shorter orbital periods and larger companion masses for B-type stars (see Sana et al. 2025, and Sect. 6.1). A clear example of this is 30 Dor, where the observed multiplicity fraction was found to be $f_{\mathrm{mult},\mathrm{O}}^{\mathrm{obs}}=35.0\pm 2.5\%$ for the O-type sample (Sana et al. 2013) and $f_{\mathrm{mult},\mathrm{B}}^{\mathrm{obs}}=25\pm 2\%$ for the B-type stars (Dunstall et al. 2015). Given their faster evolution and higher detectability as binaries, the O-type population may be more affected by binary interactions. Binary interactions can effectively reduce the number of observable binaries through mergers or envelope stripping. Envelope stripping through stable mass transfer can lead to large mass ratios, making it difficult to detect RV variations and the flux contribution from the companion in optical wavelengths. Two possible scenarios might explain the similar observed multiplicity fractions of O- and B-type BLOeM stars: (1) the initial multiplicity fraction of O-type stars is higher than that of B-type stars in the SMC but has decreased due to binary interactions; or (2) both populations have similar initial multiplicity fractions, but binary interactions are less frequent during the early phases of evolution at low metallicity. While the first option seems the most logical, there is also evidence that might support the latter possibility. As a result of the significantly reduced atmospheric opacity compared to stars in the LMC and MW, MS massive stars in the SMC are more compact, hotter, and more luminous (Maeder & Meynet 2000; Götberg et al. 2017), potentially delaying Roche-lobe overflow (Smith 2014). Based on evolutionary models, Klencki et al. (2020) found that at metallicities comparable to that of the SMC, the most common post-MS donors in massive binaries ( ∼16–50 M_⊙) are slowly expanding core-He burning (CHeB) stars. This means that interactions take place on average at later phases of evolution due to the smaller radii of low-metallicity massive stars. Even if this effect is not large enough to make a difference during the MS, it might play an important role for stars crossing the Hertzsprung gap, during which case-B mass transfer occurs. Although not the scope of this work, binary population synthesis modelling could help evaluate the significance of metallicity effects across different evolutionary phases and their impact on the multiplicity fraction of a population of massive stars, such as that of the BLOeM OB-type sample.

In Fig. 6 we show the dependence of our observed multiplicity fraction on the adopted value of the C cut-off. We do not observe a significant feature that might indicate a transition between variability caused by intrinsic features (e.g. winds or pulsations) and that caused by orbital motion. This is in contrast to what has been observed for the O-type stars and B-type supergiants in 30 Dor (Sana et al. 2013; Dunstall et al. 2015), where a distinct “kink” is visible at 16 km s⁻¹. It is not clear if this is a consequence of the less active atmospheres of the B-type stars at low metallicity or what the role of pulsations might be. However, since the winds of massive stars are much reduced at SMC metallicity (Vink et al. 2001; Björklund et al. 2021), we do not expect winds to play a major role in the RV variability of the B III/V stars (e.g. Ramachandran et al. 2019). Similarly, if pulsations were to strongly contribute to the RV variability, we would expect a feature as a kink close to 20 km s⁻¹ where the fraction of large RV amplitude variations drops significantly. We will come back to the role of pulsations in Sect. 6.1, where we discuss the impact of adopting different ΔRV threshold values.

Fig. 6. Observed multiplicity fraction as a function of the cut-off value C. The 20 km s−1 threshold used by previous studies (Sana et al. 2013; Bodensteiner et al. 2021; Banyard et al. 2022) is marked with the dashed line. The adoption of C = 20 km s−1 yields the observed binary fraction of $f_{\mathrm{mult}}^{\mathrm{obs}}=49.5\pm 2.8\%$. At C = 0, the observed multiplicity fraction is purely given by Equation (1).

4.2. Variations in the observed multiplicity fraction

We have computed the observed binary fraction for each of the eight FLAMES fields in the BLOeM survey, as shown in Table 2. We do not find a significant variation across fields, with most of them exhibiting multiplicity fractions close to 50% within 1σ. The exception is Field 2, which exhibits a multiplicity fraction of 74 ± 10%. Compared to the mean fraction of the other seven fields ($f_{\mathrm{mult}}^{\mathrm{obs}}=50\pm 2\%$), this deviation appears statistically significant (Z − score = 2.4 and p = 0.02). However, Field 2 has a relatively small sample size (N = 19) compared to the other fields. The high uncertainty increases the likelihood that the observed difference is due to statistical fluctuations inherent in small sample sizes. We have not identified any environmental factors unique to Field 2 that would result in a higher multiplicity fraction. Therefore, we interpret this result cautiously and attribute the higher fraction in Field 2 to small number statistics rather than a genuine variation in the underlying stellar population.

We have also investigated potential variations in the multiplicity fraction with spectral type, as a proxy for mass. In the histogram of Fig. 7 we show the number of stars per spectral type and luminosity class, with their corresponding observed multiplicity fractions annotated. To minimise the impact of small number statistics, spectral types have been grouped as follows: the B0 bin includes spectral types from B0 to B0.5, B0.7 is grouped with B1, and B1.5 objects remain ungrouped, as in Fig. 1. The few B2.5 stars have been incorporated into the B2 bin. The purple unfilled histogram represents the total number of stars in each spectral type bin, while $f_{\mathrm{bin}}^{\mathrm{Tot}}$ denotes the total multiplicity fraction within each bin, corresponding to the number of stars in the stacked histogram. Uncertainties in multiplicity fraction were computed using binomial statistics. Furthermore, the number of stars per luminosity class within each spectral type group is shown as stacked segments, each with its multiplicity fraction. It is not possible to draw strong conclusions from Fig. 7 due to the small number of stars with luminosity classes V and IV. For the multiplicity fraction of giants and total fraction, there are no significant trends with temperature, suggesting that within our spectral range, the multiplicity fraction remains constant.

Fig. 7. Histogram of our sample binned by spectral type. B0, B0.2, and B0.5 spectral types have been grouped in the B0 bin, as with B0.7 and B1 in the B1 bin, and B2 and B2.5 in the B2 bin, whereas B1.5 remains ungrouped. The unfilled histogram shows the total number of stars in each bin, whereas the filled histogram shows the number of identified binaries, stacked and colour-coded by luminosity class. Multiplicity fractions for each luminosity class have also been added as annotations, whereas the total multiplicity fraction is shown at the top.

4.3. Intrinsic multiplicity fraction

To compute the intrinsic multiplicity fraction of our sample, we used a similar methodology to that applied in previous studies of OB-type samples (Sana et al. 2013; Dunstall et al. 2015; Bodensteiner et al. 2021; Banyard et al. 2022) and other BLOeM samples (Sana et al. 2025; Patrick et al. 2025; Britavskiy et al. 2025).

We performed Monte Carlo (MC) simulations to estimate the detection probability of our observing campaign and correct for observational biases. We generated a synthetic population of stars where each star was randomly assigned to be either a single or a binary star based on a binomial distribution with a success probability f_bin, representing the intrinsic multiplicity fraction (Sana et al. 2013). For each binary system, we randomly drew orbital parameters, specifically orbital period (P), mass ratio (q = M₂/M₁), eccentricity (e), and primary mass (M₁), from assumed power-law distributions with specified indices and domains shown in Table 3.

Random orbital inclinations (i) and arguments of periastron (ω) were assumed to account for random orientations in space. Each synthetic star was matched to a real star in our sample, inheriting its time sampling (observation time) and RV errors. Using the assigned orbital parameters and observation times, we computed the expected RV measurements for the synthetic binaries, incorporating the real RV errors to simulate observational noise. For single stars, we assumed constant RVs with noise corresponding to the measurement errors.

We then applied our binary detection criteria described in Sect. 3.2 to the simulated observations to identify which synthetic binaries would be detected in our survey. By repeating this process 10 000 times, we obtained robust statistics and estimated the detection probability (p_det) as the fraction of synthetic binaries that met our detection criteria. This allowed us to compute the expected observed multiplicity fraction ($f_{\mathrm{obs}}^{\mathrm{simul}}$) given an intrinsic multiplicity fraction f_bin.

The simulated observed multiplicity fractions are compared with the actual observed multiplicity fraction from our data. By adjusting f_bin in our simulations to match the observed fraction, we estimated the intrinsic multiplicity fractions of our sample. This approach accounts for observational biases and detection limitations of the BLOeM campaign, providing a more accurate determination of the true multiplicity properties of our sample.

The detection probabilities for our sample of B-type binaries are shown in Fig. 8. As expected, we have a very high detection probability for short-period systems, $0.{99}_{-0.02}^{+0.01}$ for orbital periods up to 10 d. After 60–70 d the detection probability drops significantly, reaching ∼70% for P_orb = 100 d and ∼35% for P_orb = 200 d, with an integrated probability for orbital periods between 10 d and 1 yr of 0.74 ± 0.04. When considering the full possible orbital periods range (0 < log(P_orb/d) < 3.5) the detection probability of the BLOeM campaign for the B-type stars is $p_{\det }=0.{62}_{-0.06}^{+0.07}$. In the case of the mass ratio, the sharp drop is due to the minimum value considered (see Table 3). None of q, e, or M₁ have as much effect in the detection probability as P.

Fig. 8. Detection probability of B-type binary systems in the BLOeM campaign, as a function of the orbital properties P, q, e, and M1.

After accounting for the observational biases in the BLOeM campaign and given our detection probability, we found an intrinsic multiplicity fraction of f_mult = 80 ± 8%. We can compare this value to the intrinsic multiplicity fraction found for the B-type stars in 30 Dor of f_mult = 58 ± 11% (Dunstall et al. 2015). Performing a Z-test, we found these two multiplicity fractions to deviate by 3.2σ, which is statistically significant at the 1% α–level (p = 0.0015). The higher SMC fraction is even more relevant if we consider that the 30 Dor sample was dominated by luminosity class V stars (close to 60%) whereas our III/V sample contains a majority of giant stars, which as seen in Fig. 4, cannot avoid interactions below orbital periods of a few days. The difference is even larger (4.9σ) when compared to f_mult = 52 ± 8% of the galactic B-type stars in NGC 6231 (Banyard et al. 2022), and to f_mult ≈ 55% (3.2σ, assuming a mean uncertainty of 10) found for the OB-type stars in the Cygnus OB2 association (Kobulnicky et al. 2014). This suggests a trend with metallicity that will be explored in more detail in Sect. 6.1.

5. Determination of orbital periods

5.1. Eclipsing binaries

Some of our targets have been identified as EBs by OGLE (Pawlak et al. 2016). We found 41 EBs which are listed in Table D.2, including their photometric periods, I-band magnitudes, and light curve type. Their distribution of I-band magnitudes and orbital periods are shown in Fig. 9. The 41 EBs correspond to 27% of our binary sample. Among them, there are ten SB1 systems, 26 SB2s, all three SB3s, and two RV var stars. The latter two systems did not meet our binary criteria, but we found the correct orbital period for one (BLOeM 6-061), and the second system (BLOeM 7-054) has an orbital period just outside our frequency grid (1.099 d, see Sect. 5.2).

Fig. 9. Distribution of I-band magnitudes (top) and photometric orbital periods (bottom) of the systems in our sample flagged as EBs by OGLE (Pawlak et al. 2016).

Most of the EBs have been classified as detached or semi-detached by Pawlak et al. (2016), but there are two ellipsoidals and one contact system. Both ellipsoidal systems are SB2s that suffered from misidentification of their components. However, we still recovered half of the photometric orbital period in both cases, as is common for ellipsoidal variables. The only contact system, BLOeM 6-010, also suffered from component misidentification and it was flagged as a possible SB3, but we found the correct orbital period of 1.21 d. Despite finding the same period (or half in some cases) as the OGLE period for these systems, the significance in all three cases was too low to confidently list them as the correct periods (see details below).

We have done a preliminary inspection of the photometric data of the EBs mainly to find evidence and confirm third components. As mentioned in Sect. 3.2, we found clear signs of third companions in three systems that we have classified as SB3, but also evidence of potential contributions from further components in another 15 systems. Our photometric analysis indicates that there are at least eight systems where there is likely a third companion, which would represent a lower limit of 3% of triple systems in our sample. We provide comments on these and other systems in Appendix C.

5.2. Period search in the BLOeM data

We have applied the Lomb-Scargle (LS) periodogram (Lomb 1976; Scargle 1982), as implemented in Villaseñor et al. (2021), to our sample to search for periodicity in our RVs. Given the baselines for the different FLAMES fields, we could be sensitive to periods up to ∼40–60 d. The nine epochs of observations (less in some cases, see Sect. 2) allowed us to obtain significant signals for 45 systems, which corresponds to 29% of the binaries.

The significance of the signal has been defined with respect to the false-alarm-probability (FAP). All peaks above 1.1 d of period in the LS periodogram with a LS power (P_LS) above a false-alarm level (FAL) of 0.1% (i.e. FAP below 0.1%) are considered significant. It is possible that periods with a FAP above this threshold might be correct, but the signal is not strong enough due to the number of epochs, the intrinsic properties of the system (such as high eccentricity), or difficulties in the fitting procedure (e.g. component misidentification). To illustrate our criterion, we show in Figure 10 the P_LS of the highest peak in the LS periodogram (above 1.1 d) over the power at a FAL of 0.1% as a function of the orbital period. Systems are colour coded by binary classification (open circles), and we highlight systems with P_LS above a FAL of 1% with black filled circles. The green region encloses systems with orbital periods considered as significant. There are three systems with periods larger than 100 d; from the visual inspection of their RVs, they do present a variation consistent with long orbital periods, but we do not have enough coverage of the orbit to compute an accurate period.

Fig. 10. Significance of the LS periodogram peaks as a function of the recovered periods. We consider a period significant if its peak surpasses a FAL of 0.1%. Green rectangle shows the significant periods according to our criterion and the cadence and baseline of the BLOeM campaign.

5.3. Comparison with 30 Dor B-type binaries

5.3.1. Orbital periods

The B-type Binaries Characterisation programme (BBC; Villaseñor et al. 2021) studied 88 of the binary candidates identified by Dunstall et al. (2015) with 29 FLAMES observations over a baseline of 427 d. Villaseñor et al. (2021) found orbital solutions for 50 SB1 and 14 SB2 systems, and potential solutions for 20 more systems. We can compare the distributions of orbital periods for both populations, over a similar baseline. For this, we have selected nine epochs from the BBC campaign covering a baseline of ∼60 d, without including the low S/N epochs 10, 12, 13, 17 (see their Table 1 and data reduction section). We computed the LS periodogram for the corresponding RVs for all systems in the BBC sample.

The results can be seen in Fig. 11. Similarly to the periods for the BLOeM sample, we distinguish between those systems whose P_LS exceeds the 0.1% FAL (turquoise circles, 29 system), those between the 0.1 and 1% thresholds (orange rectangles, 13 systems), and those with P_LS below the 1% threshold (purple diamonds, 45 systems). We found excellent agreement for the first two groups. When limiting the results to a maximum period of 60 d, we found 93% agreement for systems with P_LS above the 0.1% threshold, and 90% agreement for those between 0.1–1%. For this reason we use the periods from both these groups in the comparison with the BLOeM sample, resulting in a sample of 34 systems.

Fig. 11. Comparison of BBC orbital periods (Villaseñor et al. 2021) with periods derived from nine epochs. The dashed diagonal line represents equal periods. Data points are categorised based on their LS Power (PLS) relative to the power at a false alarm level (FAL). The shaded red area highlights period above 60 days.

In the case of the BLOeM periods, our test with the BBC sample reaffirms our confidence in the periods that met our significance criterion. However, we restrict the periods to those highlighted in Fig. 10, not including those with P_LS < 0.1% FAL. We do include the orbital periods from OGLE, for a total of 57 objects.

The cumulative distribution of orbital periods is shown in Fig. 12. For comparison, we also plotted the distribution of periods for all BBC systems with a period below 60 d (green triangles). When comparing both samples with nine epochs, there is only a marginal difference around 6 d. To assess the significance of this difference, we have computed the Kolmogorov-Smirnov (KS) test, obtaining a KS-statistic of D = 0.24 and a p-value of 0.16, which does not reach significance at the 10% level. However, when comparing with the full BBC sample up to periods of 60 d, we do find significant differences (D = 0.37, p = 0.0003). It is clear that both samples with nine epochs lack systems with periods larger than ∼15 d in contrast to the full BBC sample which approaches a uniform distribution in log P. Judging from Fig. 11, the significance of the periods found with the LS periodogram starts to decrease after 10 d, and from that point up to 60 d we see a majority of systems with P_LS < 1% FAL, even though we did find the correct periods for many of them.

Fig. 12. Cumulative distribution of the orbital periods of the BLOeM B-type binaries (blue circles) compared to those of the LMC sample from the BBC programme (Villaseñor et al. 2021) in orange rectangles (9 epochs, 34 systems) and green triangles (all epochs, 68 systems).

Once all the 25 epochs of BLOeM observations are assembled, we will be able to recover accurate orbital solutions for most of the binaries and compare with the BBC and other samples of B-type stars over the full range of orbital periods.

5.3.2. Distribution of ΔRV_max

We also examined the distribution of ΔRV_max values for the BBC systems with determined periods. Given the larger number of epochs (29) of the BBC programme and the reliable periods, their ΔRV_max values are robust. The distribution is shown in Fig. 13, where we have colour-coded the stacked histogram by the binary classifications used by Villaseñor et al. (2021). The top panel shows the ΔRV_max values computed from all the epochs with well determined RV measurements that were used in the computation of orbital solutions. In the bottom panel, we plot the distribution obtained by only considering nine epochs as in the case of the BLOeM sample. Comparing the two panels, we see that the distribution is shifted slightly towards lower values when only nine epochs are used, with long-period, eccentric binaries being the most affected. Indeed, the three SB1 systems that fall below 20 km s⁻¹ all have P_orb > 150 d. Notably, 34 of the BBC SB1 systems have ΔRV_max values below 80 km s⁻¹, using a large binary threshold such as this (see discussion in Sect. 6.2) would miss 64% of the SB1 systems. On the other side, a threshold of 20 km s⁻¹ would still skip some of the confirmed long-period binaries in BBC when only using nine epochs, whereas only three systems without orbital solutions remain above the threshold. This evidence also suggests that 20 km s⁻¹ is an appropriate limit to separate binarity from intrinsic variability.

Fig. 13. Distributions of ΔRVmax for the BBC sample. Top panel: ΔRVmax derived from epochs used for period determination in BBC (∼25 per star). The different binary classifications from Villaseñor et al. (2021) are shown in the stacked histograms. SB1* systems were defined as those with less robust orbital periods, whereas those without orbital solutions are shown in white. Bottom panel: same as above but using only nine randomly selected epochs per star. For comparison, we include the distribution of BLOeM SB1 systems in the dashed red histogram. The dashed vertical lines at 20 and 16 km s−1 highlight threshold values discussed in Sect. 3.2.

6. Discussion

6.1. Trend with metallicity

As seen in Sect. 4.3, the large multiplicity fraction we have found for the BLOeM B-type stars suggests a trend with metallicity. To test this dependence, we performed a simple linear regression f_mult = b + m log(Z) (as the one applied by Sana et al. 2025 to the O-type stars) to the results from Kobulnicky et al. (2014), Dunstall et al. (2015), Banyard et al. (2022), and our result, where we have adopted as uncertainty for the former the mean of the errors of the other studies, and for simplicity, in the case of asymmetric errors, the mean between them. The fit is shown in Fig. 14, with its 68% confidence interval in shaded blue. We found a slope of m = −0.38 ± 0.14, which is 2.6σ away from zero.

Fig. 14. Dependence of the multiplicity fraction of different mass regimes on metallicity. The B-type studies at different metallicities are represented by blue circles; ZSMC = 0.2 Z⊙ (this work), ZLMC = 0.5 Z⊙ (Dunstall et al. 2015), ZMW = Z⊙ (Kobulnicky et al. 2014; Banyard et al. 2022). The best fit to the data with the 68% confidence interval is represented by the blue line and shaded area. Red squares show the sample of G- and K-type dwarfs from Badenes et al. (2018), whereas red diamonds correspond to the combined data of five solar-type samples presented by Moe et al. (2019), see details in main text. The dashed purple line is the linear regression found by Sana et al. (2025) for the BLOeM O-type stars (purple triangles), included for comparison.

To illustrate and compare with the case of solar-type stars, we include in the plot a sample of G- and K-type IV/V stars observed by APOGEE and first presented by Badenes et al. (2018). Close binary fraction values as a function of metallicity were obtained from Moe et al. (2019, their figure 11), with the corresponding fit shown in red. We also include three additional values (red diamonds, not included in the linear regression) at [Fe/H] = −1.0, −0.2, and +0.5. These points represent a weighted moving average from five different solar-type samples, including the Badenes et al. (2018) sample, analysed by Moe et al. (2019, see their Figure 18). The averaged values are in good agreement with the former measurements, but we only fit the APOGEE sample since it is much larger (∼20 000 objects) and presents the lowest errors. For the details of the solar-type samples we refer the reader to Moe et al. (2019), but we note that all samples were corrected for incompleteness up to orbital periods of log(P/d) = 4. We have computed the linear fit to the APOGEE data and found a slope of −0.21 ± 0.03. The two slope values for solar and B-type stars are not significantly different from each other when compared with a t-test (p = 0.36), suggesting a similar trend of the multiplicity fraction with metallicity. Figure 14 also shows the linear fit found by Sana et al. (2025) for the O-type stars (dashed purple line) with its corresponding MW (Sana et al. 2012), LMC (Sana et al. 2013; Almeida et al. 2017), and SMC (Sana et al. 2025) multiplicity fractions. No clear relation was found for the O-type stars with metallicity by Sana et al. (2025)⁴.

To assess the significance of the trend between metallicity and multiplicity fraction for the B-type stars, we conducted a t-test under the null hypothesis that the slope of the linear regression is consistent with zero. The analysis yielded a p-value of p = 0.12, meaning that we cannot reject the null hypothesis at the conventional significance level of 0.05. This indicates that the observed trend is not statistically significant. The lack of significance may be attributed to the small sample size, which limits the statistical power of the test and increases the uncertainty in the slope estimate. With only four observations, the degrees of freedom are low, resulting in a t-distribution with heavy tails. This makes it more challenging to achieve statistical significance, even when the estimated slope suggests a potential trend. As an example, we repeated the test on the APOGEE and the solar-type averaged samples, with the latter returning p = 0.13 from the three measurements, whereas the APOGEE sample with seven measurements (and smaller errors) returned p = 0.008. Future studies at different metallicities are necessary to more reliably determine if the trend with metallicity for the B-type stars is a real feature.

The increased multiplicity fraction found for the SMC B-type stars in this work, contrasts with the case of EBs. Moe & Di Stefano (2013) found multiplicity properties of EBs with early B-type primaries to be invariant with metallicity in the range −0.7 < log(Z/Z_⊙) < 0.0. Recently, Menon et al. (2024) estimated the multiplicity fraction across EBs in the LMC and SMC to also be roughly constant. However, with EBs it is only possible to probe the close-binary regime due to the limited range of orbital inclinations that would lead to eclipses, which decreases with separation. Indeed, the sample from Moe & Di Stefano (2013) had orbital periods up to 20 d. This can also be seen in Fig. 9 where we can see only two EBs in our sample with periods above 10 d, and with most of the EBs having periods of less than 4 d. To reconcile the multiplicity fraction of short-period EBs with our increased multiplicity fraction at low Z, some type of variation in the formation mechanism of massive binaries would be needed.

Close massive binaries (a < 10 au) are believed to form through disc fragmentation followed by inward orbital migration via circumbinary accretion (Tokovinin & Moe 2020). Those in the closest orbits (P_orb < 10 d) require extremely early fragmentation, coupled with subsequent envelope mass accretion to harden the orbits to such short periods, a process that remains largely independent of metallicity. This would explain why EBs, which predominantly probe these very close binaries, exhibit multiplicity properties that are invariant with metallicity. This is similar to what we see for the more massive O-type stars. In order to form an O-type star, a high accretion rate and large disc surface density are needed, which do not depend on metallicity. This presents a natural explanation for why the multiplicity fraction of O-type stars remains high at different metallicities.

However, for B-type binaries at wider separations (0.2 < a < 10 au), metallicity might have played a more significant role. Metal-poor accretion discs, such as those in the SMC, are less optically thick and can cool more efficiently. This enhanced cooling and reduced opacity facilitate gravitational instabilities, leading to increased disc fragmentation. Therefore, for stars where fragmentation occurs at later times, there is not enough remaining disc mass or accretion to dissipate the orbital energy effectively and decrease separation beyond a < 0.2 au. This mechanism would potentially result in a larger multiplicity fraction in early B-type binaries with separations between 0.2 and 10 au at low metallicity.

6.2. Considering a higher RV threshold

The discussion so far has assumed that the ΔRV threshold of C = 20 km s⁻¹ has no effect on our binary fraction. Indeed, by correcting by the observational biases, we are taking into account the binaries that did not pass our detection criteria. The recovered fraction of binaries should be independent of the choice of the threshold, otherwise it is possible that we are heavily contaminated by single stars with intrinsic variability. To test this hypothesis, we have repeated our analysis adopting values of C = 50 km s⁻¹ and C = 80 km s⁻¹ and found an anti-correlation with the multiplicity fraction. For C = 50 km s⁻¹, the intrinsic multiplicity fraction is reduced to $f_{\mathrm{mult}}={70}_{-9}^{+12}\%$, whereas for C = 80 km s⁻¹ it further drops to $f_{\mathrm{mult}}={56}_{-8}^{+11}\%$.

The lower intrinsic multiplicity fraction has clear consequences for the interpretation of our analysis. Perhaps the most relevant one is the effect on the metallicity trend. This lower value argues in favour of a constant multiplicity fraction across the metallicities of the MW, LMC, and SMC, very similar to the one found for the BLOeM O-type stars. The metallicity of the SMC is representative of galaxies at redshifts z = 3 (Sommariva et al. 2012) and it can be found at redshifts as high as z = 10 (Nakajima et al. 2023), suggesting an universality in the formation process of O- and B-type binaries. However, we note that the RV threshold here is only applied to the BLOeM sample, it would be interesting to see if the effect on the multiplicity fraction is also observed for the other samples of OB-type stars, but this is beyond the scope of our work.

However, it is crucial to understand where the dependence on the ΔRV threshold comes from if we want to draw any conclusions. One possible explanation for the variation of the intrinsic binary fraction is that there are other sources of variability in play contaminating our binary sample, different from binary motion. Indeed, our sample intersects the predicted instability region of β Cephei pulsating stars, which span the mass range of about 8 to 25  M_⊙, extending from the MS to the giant phase (Bowman 2020). However, the bulk of RV amplitudes caused by pulsations concentrates in the 10–20 km s⁻¹ range, with a tail that can go up to 40–50 km s⁻¹ (Stankov & Handler 2005; Hey & Aerts 2024), and periods roughly between 2 and 12 h (Aerts et al. 2010).

Even if we are heavily contaminated by β Cephei pulsators, it is hard to explain the drop in multiplicity fraction when increasing C from 50 to 80 km s⁻¹. We do not expect pulsation to play a significance role at that scale of variability. Furthermore, the mechanism causing pulsations is traditionally explained as variations in the star’s opacity driven by elements of the iron group (Aerts et al. 2010). Because the presence of significant iron is needed for the opacity enhancement, pulsating massive stars are rare in low metallicity environments, as in the SMC (Bowman 2020). Therefore, we would at least expect the role of pulsators to be smaller in comparison to the MW and LMC.

We are not aware of another source of RV variability that could play such a significant role at these high ΔRV values (≥50 km s⁻¹). With the full number of epochs available, robust orbital periods will allow us to more effectively disentangle the binary and pulsation contribution, since we are not sensitive to the periods produced by pulsations in β Cephei stars. Moreover, time series photometry may help elucidate (high-amplitude) pulsators among the BLOeM sample and the wider population of SMC massive stars (see Bowman et al. 2024).

An alternative explanation is that our assumed distributions of orbital parameters are incorrect. Our results are particularly sensitive to the orbital period distribution, for which we have assumed an Öpik law (see Table 3). As discussed in Sect. 5, the fraction of binary system with periods between ∼3–10 d seems to be larger for the BLOeM sample than for the BBC sample of binaries. Although our analysis of the orbital period distribution is inconclusive about the significance of the difference in the distribution of both samples according to the KS test, there are other pieces of supporting evidence that are worth examining. First, the fraction of SB2 systems is much larger among the BLOeM B-type binaries in comparison to the 30 Dor counterpart, with 39% of the binaries (41% if we include the SB3 systems) versus only 17% in the BBC sample. There are several factors influencing the detection of SB2 systems, some are intrinsic to the system such as orbital period, mass ratio, and rotational broadening, but external factors like S/N can also play a role. Villaseñor et al. (2021) estimated a mean S/N for the BBC sample of 33. In contrast, the BLOeM sample has a mean S/N of 55, which is significantly higher, probably responding to the higher fraction of giant stars. Moreover, the fraction of SB2s is similar to that of the VFTS O-type binaries (Almeida et al. 2017), which presented similar S/N to our sample. In conclusion, the fraction of SB2 systems, although large, is not impossible to explain.

Eclipsing binaries are also abundant among the BLOeM B-type III/V sample, reaching 27% of the binaries, compared to only 12% found by the BBC programme (Villaseñor et al. 2021). This could also be seen as evidence of an overabundant fraction of short-period systems. However, contrarily to BBC, our BLOeM sample is dominated by giants which are roughly a factor two larger than MS stars. Since the probability of eclipses scales with radius, giant primaries will have a higher probability of being seen as EBs. Only accounting for MS stars, we find that similar fractions of the BBC and BLOeM binaries are MS EBs (14% and 11% respectively).

We further test the possibility of deviating from an Öpik law with Markov chain Monte Carlo (MCMC) simulations in Sect. 6.4.

6.3. Binary fraction from the distribution of ΔRV_max

Even though we cannot explain the decrease in the intrinsic multiplicity fraction with the assumed RV threshold, we can explore other methods to obtain the multiplicity fraction. The distribution of ΔRV_max has been used to determine the binary fraction of various stellar populations, including white dwarfs (Badenes & Maoz 2012; Maoz et al. 2012) and low-mass stars (Badenes et al. 2018; Mazzola et al. 2020). These studies show that the distribution is characterised by a “core” of low ΔRV_max values, dominated by single stars and RV uncertainties, and a “tail” of higher values associated with close binaries. By reproducing these signatures with simulations, one can infer the intrinsic binary fraction of a sample without relying on explicit binary criteria.

We adopt an MC method similar to that described in Sect. 4.3, sampling the same physical and orbital parameters but extending the orbital period range down to log P = −0.15 to include sub-day periods observed by OGLE. We also extend the mass-ratio distribution to q = 0.01 to probe lower-mass companions. We assumed circular orbits (e = 0) for systems with orbital periods below 2 d (Moe & Di Stefano 2017), whereas for longer orbits we sample e from 0 to e_max (as given by Badenes et al. (2018)). Values of the distribution indices remain as in Table 3. To reduce stochastic noise, we generate 10 000 stars with a binary fraction given by f_bin and use a MCMC approach to fit the observed ΔRV_max distribution by varying f_bin. We scaled each simulated ΔRV_max distribution by the ratio of observed to simulated stars, and then compare with the real data using a Poisson-based log-likelihood (Cash 1979) in log-space intervals:

$$\begin{array}{c}\hfill log\mathcal{L}={\displaystyle \sum _i}[n_i^{\mathrm{obs}}ln(n_i^{\mathrm{mock}})-n_i^{\mathrm{mock}}],\end{array}$$(2)

where n_i^obs is the number of observed stars in bin i and n_i^mock is the corresponding (scaled) count of simulated stars.

The posterior distribution for f_bin is shown in Fig. 15, our simulations clearly favour a high binary fraction of f_mult = 79 ± 5% in good agreement with our value of f_mult80 ± 8% reported in Sect. 4.3. Figure 16 shows the observed ΔRV_max distribution (black histogram), it can be seen that there is no evident separation between a core and tail component, possibly due to the uncertainty in our RV measurements or the scarcity of very wide binaries (log(P/d) > 4) in comparison to the distribution of orbital periods of sun-like stars which peaks roughly at 10⁵ d (see Raghavan et al. 2010; Badenes et al. 2018). We include the distribution of ΔRV_max values for four simulated samples of 10 000 stars computed with an intrinsic binary fraction of 0.0, 0.5, 0.79, and 1.0 to illustrate the change in the distribution with f_bin. A lower binary fraction naturally increases the number of stars with low-amplitude RV variations, whereas a higher binary fraction adds more highly RV-variable systems. Notably, the transition value between single- and binary-dominated bins in the simulations occurs near 20 km s⁻¹, supporting the choice of this value as a sensible threshold.

Fig. 15. Posterior distribution for the binary fraction, fbin, from our MCMC simulations. Dashed grey lines indicate the 16th and 84th percentile.

Fig. 16. Distribution of ΔRVmax values for the BLOeM sample (black) and example simulated distributions for different binary fraction values, including our result of fbin = 0.79. The dashed grey line indicates the 20 km s−1 threshold.

Another point worth mentioning is the inclusion of slightly evolved stars in our sample. As shown in Sect. 4 and Fig. 4, for typical q = 0.5 and i = 60° detached systems with evolved primaries (log(g) = 3.5) can exhibit peak-to-peak RVs about 100 km s⁻¹ lower than those with log(g) = 4.0. These findings illustrate how stellar evolution can effectively shift some of the high-amplitude binaries towards lower values in the observed distribution of ΔRV_max. Because our simulations do not track orbital evolution with mass transfer or stellar ageing, these post-interaction systems might escape detection. Future BLOeM papers will explore these evolutionary effects in detail. In conclusion, our multiplicity fraction can thus be interpreted as a lower limit due to binary evolution likely removing part of the high-ΔRV_max tail from the distribution.

It is also interesting to examine the distribution of ΔRV_max values for f_bin = 0 (blue histogram in Fig. 16); there is a fraction of simulated single stars (9%) with values above 20 km s⁻¹. The RV values are simulated by adding a Gaussian noise with a standard deviation corresponding to the uncertainty of the RV value at each epoch of a real, randomly assigned, BLOeM star. One single bad epoch (e.g. large uncertainty due to low S/N or a bad multi-component fit) can produce a large ΔRV_max value. Figure 17 shows the distribution of the simulated sample of single stars in the binary criteria diagram. Out of the 10 000 stars, only four pass both binary criteria, meaning only 0.04% false-positive detections. This exercise not only confirms the high multiplicity fraction of early B-type dwarfs and giants in the SMC, but also further validates the choice of the binary criteria for our sample.

Fig. 17. As Fig. 5 for the simulated sample of single stars (fbin = 0), colour coded by number of stars. Top and right panels show the 1D distribution of each of the binary criteria.

6.4. Variations in the orbital period distribution

As mentioned in Sect. 6.2, the distribution of orbital periods plays a significant role in the determination of the binary fraction. As a final test, we treat the exponent of the orbital period distribution π as a free variable in our MCMC simulation. The resulting posterior distribution can be seen in Fig. 18. Fitting both π and f_bin simultaneously yielded $f_{\mathrm{bin}}=0.{85}_{-0.09}^{+0.07}$ and π = 0.16 ± 0.15. Although this represents a slight preference for longer periods (π > 0) the result is essentially consistent with π = 0. We find a significant correlation between both parameters, but overall, varying π does not drastically alter the binary fraction, suggesting that while the precise shape of the period distribution is not tightly constrained by the current data, the inferred binary fraction remains robustly high. Once the full 25 epochs of observations from the survey become available, it will be possible to put much stronger constraints on the distribution of orbital periods.

Fig. 18. Corner plot illustrating the posterior distribution for the binary fraction (fbin) and period distribution exponent (π). The central (off-diagonal) panel displays the 2D joint posterior density, with contours indicating (2D) 1σ, 2σ, and 3σ confidence levels. The two marginal panels (diagonal) show the 1D distributions for fbin and π, with vertical lines marking the 16%, 50%, and 84% quantiles.

7. Conclusions

In this work, we focused on the B-type sample of dwarf and giant stars from the BLOeM campaign in the SMC. We analysed 309 stars by measuring their RV variations to identify close binaries and constrain the multiplicity fraction. Our main conclusions are:

1. High multiplicity fraction at low metallicity: We detect 153 spectroscopic binaries (91 SB1, 59 SB2, 3 SB3) based on significant RV variability. For this, we consider as binaries those systems which simultaneously fulfil having peak-to-peak RV variations above a threshold of C = 20 km s⁻¹ and an RV significance above 4σ, consistent with previous studies of massive stars in the MW and LMC. This corresponds to an observed multiplicity fraction of $f_{\mathrm{mult}}^{\mathrm{obs}}=50\pm 3\%$, exceeding reported values for B-type stars in the MW and LMC, and suggesting a potentially larger multiplicity fraction at lower metallicities.

2. Orbital periods of the close binaries: Our multi-epoch strategy is well suited for identifying close systems with periods comparable to the total baseline of the observations. We obtained reliable orbital periods for close to 30% of the detected binaries. However, we cannot confirm any significant differences with short-period binaries at higher metallicities. The full set of 25 BLOeM observations will enable robust orbital solutions, providing a comprehensive understanding of the distribution of orbital parameters of our sample.

3. Intrinsic multiplicity fraction from independent methods: After correcting for observational biases with MC simulations, we infer an intrinsic close-binary fraction of f_mult = 80 ± 8%. Experimenting with a higher RV threshold reduces the multiplicity fraction down to ∼55% for C = 80 km s⁻¹. However, an independent MCMC analysis, which does not rely on binary criteria but considers the full ΔRV_max distribution, returns a binary fraction of f_mult = 79 ± 5%, confirming the high multiplicity. Letting the distribution of orbital periods vary further reaffirms our result, leading to $f_{\mathrm{mult}}={85}_{-9}^{+7}\%$, although the orbital period distribution cannot be well constrained with our data.

4. Potential anti-correlation with metallicity: The resulting multiplicity fraction is significantly higher in comparison to similar samples in the LMC and MW. This is consistent with an anti-correlation between metallicity and close-binary fraction among B-type III/V stars. While the fraction of O-type binaries seems to remain constant across different metallicities, B-type stars might be affected by how early or late disc fragmentation occurs during massive star formation. B-type binaries forming during extremely early fragmentation would have enough time and disc material to spiral inwards to very short orbits (P_orb < 10 d). At later times however, low-metallicity discs have had more time to cool down due to the reduced opacity, increasing fragmentation and therefore enhancing multiplicity. However, with less accreting material, the star is not able to migrate to such low separations. This implies that the increased multiplicity fraction at low metallicity predominantly involves binaries with orbital periods above ∼10 d, while the fraction of shorter-period binaries remains invariant with metallicity, consistent with previous studies of EBs.

5. Implications for the early Universe: The SMC’s metallicity is similar to that of many high-redshift, star-forming galaxies. A higher fraction of close binaries at low metallicity would translate to more frequent interactions, potentially increasing the number of exotic transients and envelope-stripped stars, relevant for cosmic reionisation. Our results also provide important input for population synthesis studies at low metallicity, and are critical for the population of gravitational-wave progenitors.

Data availability

Tables D.1 and D.2 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/698/A41

Acknowledgments

We thank the referee for his thoughtful and constructive comments that have greatly clarified our discussion and strengthened our conclusions. We acknowledge the support of the Data Science Group at the Max Planck Institute for Astronomy (MPIA) for their useful suggestions in the improvement of ravel. The authors acknowledge support from the European Research Council for the ERC Advanced Grant 101054731. DMB gratefully acknowledges funding from UK Research and Innovation (UKRI) in the form of a Frontier Research grant under the UK government’s ERC Horizon Europe funding guarantee (SYMPHONY; grant number: EP/Y031059/1), and a Royal Society University Research Fellowship (URF; grant number: URF/R1/231631). IM acknowledges support from the Australian Research Council (ARC) Centre of Excellence for Gravitational Wave Discovery (OzGrav; project number CE230100016). KS is funded by the National Science Center (NCN), Poland, under grant number OPUS 2021/41/B/ST9/00757. TS acknowledges support by the Israel Science Foundation (ISF) under grant number 0603225041 LRP acknowledges support by grants PID2019-105552RB-C41 and PID2022-137779OB-C41 funded by MCIN/AEI/10.13039/501100011033 by “ERDF A way of making Europe”. LRP acknowledges support from grant PID2022-140483NB-C22 funded by MCIN/AEI/10.13039/501100011033. DP acknowledges financial support from the FWO junior postdoctoral fellowship No. 1256225N. MP is supported by the BEKKER fellowship BPN/BEK/2022/1/00106 from the Polish National Agency for Academic Exchange and the Royal Physiographic Society in Lund through the Märta and Erik Holmbergs Endowment.

Note added in proof. The authors regret the omission of Dr Artemio Herrero and gratefully acknowledge his substantive contribution during the peer-review process.

References

Appendix A: Additional figures

A.1. Spatial distribution of the B III/V sample

Fig. A.1. Spatial distribution of the B-type dwarfs across the eight fields observed by the BLOeM campaign. The bright star-forming region NGC 346 is featured in Field 4, however, it has not been extensively covered by the BLOeM campaign (see main text for details).

A.2. Hertzsprung–Russell diagram

Figure A.2 shows the Hertzsprung-Russell diagram (HRD) of the full BLOeM sample, colour-coded by their different spectral type regimes. Each of the sub-samples has been investigated independently: the O-type stars by Sana et al. (2025), the B-type dwarfs and giants in this work, the early B-type supergiants (BSGs) by Britavskiy et al. (2025), the Oe/Be stars by Bodensteiner et al. (2025), and the late BAF supergiants by

Fig. A.2. Hertzsprung-Russell diagram of the full BLOeM sample. Stars are colour-coded by spectral type, as indicated in the legend, and Teff and log(L) values are from Bestenlehner et al. (in prep.). The black curves show evolutionary tracks for different initial masses at SMC metallicity from Schootemeijer et al. (2019) with an overshooting prescription from Hastings et al. (2021). Lines of constant radius are marked in grey, and isochrones of different ages are shown with dashed brown curves. The black cross on the left illustrates typical uncertainties in Teff and log(L).

Patrick et al. (2025). All these works have investigated the multiplicity fraction of their respective samples, and there is a clear trend with evolutionary phase. The two samples with the less evolved content, non-supergiant O- and B-type stars, have the highest multiplicity fractions, with ${70}_{-6}^{+11}\%$ found by Sana et al. (2025) and 80 ± 8% found in this work. The fraction rapidly decreases for early BSGs, reaching 40 ± 4%. Furthermore, Britavskiy et al. (2025) reports a drop in the observed binary fraction at spectral types later than B2. Looking at the HRD of Fig. A.2, we can see that the cooler BSGs have log(L/L_⊙) values above 4.5 dex. Therefore, it is likely that the drop in the binary fraction is related to the evolutionary phase of these stars; namely, the large radii of the BSGs have induced binary interactions, effectively removing the close binaries. This is supported by the study of Bodensteiner et al. (2025), who found a binary fraction of 18-32% among the emission-line stars. The stark contrast in the multiplicity properties of Oe/Be stars and OB-type stars strongly suggests that the former are mostly binary interaction products, as proposed by Bodensteiner et al. (2025). In many cases, the initially more massive star has lost a large fraction of its envelope mass and produces RV variations on the now more massive companion that are too small to detect, leading to a lower observed binary fraction. The more evolved stars in the BLOeM sample, the BAF supergiants, present the lowest binary fraction of all, with less than 18% of the late BSGs and $8_{-7}^{+9}\%$ of the A- and F-type stars showing binary signatures (Patrick et al. 2025). Although the intrinsic binary fractions for the supergiant samples are highly uncertain due to evolutionary aspects and the unclear role of intrinsic variability, the trend is clear. Binary evolution strongly impacts the evolutions of massive stars, as is evidenced by the large drop in the binary fraction from the MS to the late supergiant phases; if the intrinsic binary fraction of AF supergiants is really below 20%, then close to 60% of massive stars must have gone through envelope stripping or mergers before reaching the AF supergiant phase, and now appear as single stars.

Appendix B: Radial velocity uncertainties

Due to the importance of the RV uncertainties for the identification of binaries and ultimately for the intrinsic binary fraction, we have investigated the reliability of the RV uncertainties obtained form our spectral line profile fitting approach. In order to compare with a different and independent technique, we have performed cross-correlation in a similar way to that implemented by Mahy et al. (2022) and also applied to other BLOeM samples Bodensteiner et al. (2025), Patrick et al. (2025), Britavskiy et al. (2025)). Due to the difficulties of cross correlation handling SB2 systems, we exclude the SB2 and SB3 systems from the comparison.

The individual RV comparison is shown in the histogram of Fig. B.1. In both cases there is a tail after 7 km s⁻¹ which goes to almost zero after 10 km s⁻¹. The cross correlation uncertainties have a primary peak between 1-1.5 km s⁻¹, and a smaller secondary peak close to 4 km s⁻¹, whereas the uncertainties of the profile fitting technique peak at 2 km s⁻¹. However, both distributions have nearly identical medians close to 2.45 km s⁻¹. This simple comparison shows that there is no underestimation of our RV uncertainties that could lead to an overestimation of σ_d, and consequently to a larger multiplicity fraction.

Fig. B.1. Comparison of the uncertainties on the individual RV measurements provided by our profile fitting technique and those from cross correlation.

In Fig. B.2 we show the dependence of our RV uncertainty on the S/N of the individual observations. We can clearly see how the uncertainty anti-correlates with S/N. Spectra with S/N > 70 consistently have RV uncertainties of less than 3 km s⁻¹. Spectra with S/N < 50 have a larger spread on uncertainties, but most of them have uncertainties below 5 km s⁻¹.

Fig. B.2. Radial velocity uncertainty of our profile fitting approach as a function of S/N of the individual epochs.

Appendix C: Notes on individual systems

C.1. Confirmed SB3 systems

BLOeM 5-062: (OGLE-SMC-ECL-5096) Characterised by a narrow absorption feature moving with a large RV. As shown in Fig. C.1, the core of the line is not well fitted with a Gaussian profile, and in epochs 5 and 6 there are two other absorption features at both sides of the narrow absorption. Our inspection of the photometric data shows that this is likely a doubly eclipsing system, with a potential outer period of 77.2 d, close to the spectroscopic period (80 d). The OGLE (inner) photometric period is 5.2 d.

Fig. C.1. Gaussian fitting of the spectral lines He Iλ4388 for the SB3 system BLOeM 5-062. The dashed red line shows the mean SMC velocity of 172 km s−1 (Evans & Howarth 2008), and the dashed orange lines show differences of ±200 km s−1 with respect to the mean velocity. Epoch number increases from left to right and from top to bottom.

BLOeM 6-062: (OGLE-SMC-ECL-5853) Similarly to BLOeM 5-062, it shows an absorption feature with a narrow core, but in this case the core seems stationary or moving at a much lower velocity as it can be seen in Fig. C.2, suggesting it could belong to the component in the wider orbit. The amplitude of the narrow feature varies drastically among epochs, a sign of the movement of an inner pair. The two other absorption features are more clearly visible in epoch 4. The OGLE period is 1.13 d, whereas the spectroscopic period is 1.98 d, which we did not consider significant, and it is clearly affected by the contribution of the third component.

Fig. C.2. As Fig. C.1, this time displaying line ionHeiλ4471 for BLOeM 6-062.

BLOeM 7-057: (OGLE-SMC-ECL-4024) As shown in Fig. C.3, it clearly displays a narrower component, apparently stationary, and a broader one that is difficult to identify as two components. However, the profile variability in the broader component is evident, and consistent with large RV amplitudes of an inner pair, specially in epochs 2, 4, and 7. The photometric data suggest this is possibly a doubly eclipsing system, and a possible secondary period of 2.57 d (half the low-significance spectroscopic period).

C.2. Other eclipsing binaries

BLOeM 1-037: (OGLE-SMC-ECL-4510) Wrong period in the OGLE catalogue (1/3 harmonic), actual period of 9.975 d (twice the spectroscopic period).

BLOeM 2-051: (OGLE-SMC-ECL-2251) Eccentric orbit with apsidal motion, possibly due to a third companion. However, no evidence of a double eclipse was found.

BLOeM 3-050: (OGLE-SMC-ECL-6645) Possibly doubly eclipsing. No signs of outer period though.

BLOeM 3-072: (OGLE-SMC-ECL-1475) Possibly a doubly eclipsing system, possible secondary period of 10.03 d. This system was flagged as possible SB3.

BLOeM 4-095: (OGLE-SMC-ECL-4255) Probably doubly eclipsing. Secondary OGLE period of 2.72 d, also identified spectroscopically. This system was flagged as possible SB3.

BLOeM 5-002: (OGLE-SMC-ECL-4756) Probably doubly eclipsing. Secondary OGLE period of 2.06 d (spectroscopic period).

BLOeM 5-030: (OGLE-SMC-ECL-7623) Spectroscopic period is half of the photometric period.

BLOeM 6-115: (OGLE-SMC-ECL-5933) Spectroscopic period is twice the photometric period.

BLOeM 7-032: (OGLE-SMC-ECL-3848) Probably doubly eclipsing. Secondary OGLE period of 17.58 d (slightly longer than the spectroscopic period of 16.66 d). This system was flagged as possible SB3.

BLOeM 7-060: (OGLE-SMC-ECL-4028) Spectroscopic period is half of the photometric period.

BLOeM 8-029: (OGLE-SMC-ECL-2377) Spectroscopic period is twice the photometric period.

Fig. C.3. As Fig. C.1 for BLOeM 7-057.

BLOeM 8-105: (OGLE-SMC-ECL-2967) Spectroscopic period is a fourth of the photometric period. Eccentric orbit.

All Tables

All Figures

	Fig. 1. Spectral types of the B-type III/V sample, colour coded by field. The number of stars per field is shown in parentheses.
In the text

Fig. 2. Top: spectral-line fits for the SB1 system BLOeM 4-097. Three different spectral lines are shown for epochs closest to quadrature. Observed spectrum (black) is overlaid by best fit and 1σ confidence levels (red); dashed grey profiles indicate initial guesses. Bottom: Similar to above, but for SB2 system BLOeM 1-055. Two He I lines at three epochs (near quadrature and conjunction) are shown. The dashed red line marks the mean SMC velocity of 172 km s−1 (Evans & Howarth 2008), with a ±200 km s−1 deviation (dashed yellow lines).

In the text

	Fig. 3. Maximum peak-to-peak RVs for our sample, colour-coded by binary status. We define RV cst stars as those with peak-to-peak RV variability not reaching 2σ significance, which we consider as constant or non-variable. Whereas stars that surpass the 2σ level but do not meet both multiplicity criteria are considered variable (RV var). ΔRV values are grouped in bins of 10 km s−1.
In the text

Fig. 4. Top panel: Contour plot showing the critical period (Pcrit) as a function of primary mass (M1) and surface gravity (log(g)) for systems with i = 60° and q = 0.5. Dashed white contour lines highlight selected Pcrit values (1, 2, 5, and 10 days). Bottom Panel: As the top panel, this time displaying the corresponding change in RV (ΔRVmax) computed at Pcrit. Dashed white curves mark the contours where ΔRVmax = 20 km s−1 at different inclinations and q = 0.1.

In the text

	Fig. 5. Spectroscopic binaries, RV var, and RV cst stars classified according to our binary criteria. Our adopted thresholds of σdetect > 4 and C > 20 km s−1 are shown with dashed grey lines. Green region highlights systems fulfilling both criteria.
In the text

Fig. 6. Observed multiplicity fraction as a function of the cut-off value C. The 20 km s−1 threshold used by previous studies (Sana et al. 2013; Bodensteiner et al. 2021; Banyard et al. 2022) is marked with the dashed line. The adoption of C = 20 km s−1 yields the observed binary fraction of $f_{\mathrm{mult}}^{\mathrm{obs}}=49.5\pm 2.8\%$. At C = 0, the observed multiplicity fraction is purely given by Equation (1).

In the text

Fig. 7. Histogram of our sample binned by spectral type. B0, B0.2, and B0.5 spectral types have been grouped in the B0 bin, as with B0.7 and B1 in the B1 bin, and B2 and B2.5 in the B2 bin, whereas B1.5 remains ungrouped. The unfilled histogram shows the total number of stars in each bin, whereas the filled histogram shows the number of identified binaries, stacked and colour-coded by luminosity class. Multiplicity fractions for each luminosity class have also been added as annotations, whereas the total multiplicity fraction is shown at the top.

In the text

	Fig. 8. Detection probability of B-type binary systems in the BLOeM campaign, as a function of the orbital properties P, q, e, and M1.
In the text

	Fig. 9. Distribution of I-band magnitudes (top) and photometric orbital periods (bottom) of the systems in our sample flagged as EBs by OGLE (Pawlak et al. 2016).
In the text

	Fig. 10. Significance of the LS periodogram peaks as a function of the recovered periods. We consider a period significant if its peak surpasses a FAL of 0.1%. Green rectangle shows the significant periods according to our criterion and the cadence and baseline of the BLOeM campaign.
In the text

	Fig. 11. Comparison of BBC orbital periods (Villaseñor et al. 2021) with periods derived from nine epochs. The dashed diagonal line represents equal periods. Data points are categorised based on their LS Power (PLS) relative to the power at a false alarm level (FAL). The shaded red area highlights period above 60 days.
In the text

	Fig. 12. Cumulative distribution of the orbital periods of the BLOeM B-type binaries (blue circles) compared to those of the LMC sample from the BBC programme (Villaseñor et al. 2021) in orange rectangles (9 epochs, 34 systems) and green triangles (all epochs, 68 systems).
In the text

Fig. 13. Distributions of ΔRVmax for the BBC sample. Top panel: ΔRVmax derived from epochs used for period determination in BBC (∼25 per star). The different binary classifications from Villaseñor et al. (2021) are shown in the stacked histograms. SB1* systems were defined as those with less robust orbital periods, whereas those without orbital solutions are shown in white. Bottom panel: same as above but using only nine randomly selected epochs per star. For comparison, we include the distribution of BLOeM SB1 systems in the dashed red histogram. The dashed vertical lines at 20 and 16 km s−1 highlight threshold values discussed in Sect. 3.2.

In the text

Fig. 14. Dependence of the multiplicity fraction of different mass regimes on metallicity. The B-type studies at different metallicities are represented by blue circles; ZSMC = 0.2 Z⊙ (this work), ZLMC = 0.5 Z⊙ (Dunstall et al. 2015), ZMW = Z⊙ (Kobulnicky et al. 2014; Banyard et al. 2022). The best fit to the data with the 68% confidence interval is represented by the blue line and shaded area. Red squares show the sample of G- and K-type dwarfs from Badenes et al. (2018), whereas red diamonds correspond to the combined data of five solar-type samples presented by Moe et al. (2019), see details in main text. The dashed purple line is the linear regression found by Sana et al. (2025) for the BLOeM O-type stars (purple triangles), included for comparison.

In the text

	Fig. 15. Posterior distribution for the binary fraction, fbin, from our MCMC simulations. Dashed grey lines indicate the 16th and 84th percentile.
In the text

	Fig. 16. Distribution of ΔRVmax values for the BLOeM sample (black) and example simulated distributions for different binary fraction values, including our result of fbin = 0.79. The dashed grey line indicates the 20 km s−1 threshold.
In the text

	Fig. 17. As Fig. 5 for the simulated sample of single stars (fbin = 0), colour coded by number of stars. Top and right panels show the 1D distribution of each of the binary criteria.
In the text

Fig. 18. Corner plot illustrating the posterior distribution for the binary fraction (fbin) and period distribution exponent (π). The central (off-diagonal) panel displays the 2D joint posterior density, with contours indicating (2D) 1σ, 2σ, and 3σ confidence levels. The two marginal panels (diagonal) show the 1D distributions for fbin and π, with vertical lines marking the 16%, 50%, and 84% quantiles.

In the text

	Fig. A.1. Spatial distribution of the B-type dwarfs across the eight fields observed by the BLOeM campaign. The bright star-forming region NGC 346 is featured in Field 4, however, it has not been extensively covered by the BLOeM campaign (see main text for details).
In the text

Fig. A.2. Hertzsprung-Russell diagram of the full BLOeM sample. Stars are colour-coded by spectral type, as indicated in the legend, and Teff and log(L) values are from Bestenlehner et al. (in prep.). The black curves show evolutionary tracks for different initial masses at SMC metallicity from Schootemeijer et al. (2019) with an overshooting prescription from Hastings et al. (2021). Lines of constant radius are marked in grey, and isochrones of different ages are shown with dashed brown curves. The black cross on the left illustrates typical uncertainties in Teff and log(L).

In the text

	Fig. B.1. Comparison of the uncertainties on the individual RV measurements provided by our profile fitting technique and those from cross correlation.
In the text

	Fig. B.2. Radial velocity uncertainty of our profile fitting approach as a function of S/N of the individual epochs.
In the text

	Fig. C.1. Gaussian fitting of the spectral lines He Iλ4388 for the SB3 system BLOeM 5-062. The dashed red line shows the mean SMC velocity of 172 km s−1 (Evans & Howarth 2008), and the dashed orange lines show differences of ±200 km s−1 with respect to the mean velocity. Epoch number increases from left to right and from top to bottom.
In the text

	Fig. C.2. As Fig. C.1, this time displaying line ionHeiλ4471 for BLOeM 6-062.
In the text

	Fig. C.3. As Fig. C.1 for BLOeM 7-057.
In the text