© The Authors 2026

1. Introduction

The evolution of single massive stars is regulated by their mass, rotation, and mass loss (e.g., Langer 2012; Ekström et al. 2012). Most massive stars are also part of multiple systems with at least one nearby companion (e.g., Mason et al. 2009; Sana et al. 2011). Because of this high multiplicity fraction, interactions through mass transfer or mergers are thought to fundamentally modify the evolution of high-mass stars (Sana et al. 2012). When this interaction happens, and what its outcome is, primarily depends on the initial properties of a given system. Therefore, a good understanding of the overall multiplicity properties of massive star populations is fundamental to predict their evolutionary life cycle and final outcome.

In the last decade, the characterization of the multiplicity properties of massive stars has mostly been performed through spectroscopy (e.g., Banyard et al. 2022; Dunstall et al. 2015; Kiminki & Smith 2018; Kobulnicky et al. 2014; Maíz Apellániz et al. 2019; Villaseñor et al. 2021, 2025) and high-angular-resolution imaging techniques (Maíz Apellániz et al. 2018; Rainot et al. 2022; Pauwels et al. 2024). If limited to these approaches, the view of massive star multiplicity would be biased toward tight (< 1 AU) and wide (> 10³ AU) systems (Mason et al. 1998). For nearby stars, interferometry offers a way to probe the intermediate separation range. However, the higher contrast (Δm > 2 mag) and closest angular separation (d < 75 milli-arcsec) regimes were challenging to probe with the first generations of optical interferometers (Sana & Evans 2011).

However, progress in the robustness and sensitivity of long-baseline interferometry in the 2000s offered new possibilities. The Southern MAssive Stars at High angular resolution (SMASH+) survey (Sana et al. 2014, henceforth Paper I) combined observations performed with the sparse aperture masking (SAM) mode (Lacour et al. 2011) of NACO at the Very Large Telescope (VLT) and with the four-beam combiner PIONIER (Le Bouquin et al. 2011) at the VLT Interferometer (VLTI). NACO probes separations in the range 30–250 mas, while PIONIER opens the 1–45 mas window. Additionally, NACO provides adaptive-optics- (AO-) corrected imaging at large working angles, which allows the detection of binaries at separations > 300 mas. SMASH+ targeted all nearby O stars in the southern hemisphere (δ < 0°) with H-band magnitudes brighter than m_H = 7.5. In total, 117 O-type stars were observed with PIONIER, and 162 O-type stars with NACO/SAM. Of these objects, 105 were observed with both instruments. The SMASH+ sample is summarized in Table 1.

The resolved pairs found by the SMASH+ survey are shown in Fig. 1, which also shows the separation ranges and median sensitivities of the instruments. The survey revealed a large number of companions, with a fraction of 0.53 stars having at least one resolved companion within 200 mas. Including known spectroscopic or eclipsing companions, the multiplicity fraction reaches 0.91 at separations smaller than 8″ and the average number of companions is 2.2 ± 0.3. For luminosity class (LC) V, the fraction of bound companions reaches 100% at 30 mas. The SMASH+ survey therefore demonstrated that massive stars form almost exclusively in multiple systems and that the majority of massive stars are in triples or higher-order multiples.

Fig. 1. Magnitude difference versus angular separation for the pairs detected in the SMASH+ survey. Circled objects are detected by both SAM and PIONIER. The solid lines indicate the median H-band sensitivity of the various instruments. Figure reproduced from Sana et al. (2014).

A decade later, the SMASH+ survey has remained a reference for the separation range beyond the spectroscopic regime, providing a database to constrain the distributions of wide binaries and higher-order systems (e.g., Moe & Di Stefano 2017; Perets 2025). For this study, we converted the observed angular separations and magnitude differences of individual companions detected in the SMASH+ survey to projected physical separations and mass ratios, respectively, to derive their distributions and investigate possible correlations between physical and orbital properties. The methods used are described in Sect. 2. Section 3 investigates the sensitivity domain of the SMASH+ survey while Sect. 4 presents and discusses orbital parameter distributions. We summarize our findings in Sect. 5.

2. Methods

In this section we describe the methods used to convert the observed angular separations and magnitude differences to physical units, i.e., projected separation expressed in astronomical units (AU) and mass ratios. The distance determination needed to convert the angular separations to physical separations is presented in Sect. 2.1, followed by the mass determinations based on the magnitudes in Sect. 2.2.

2.1. Distance determination

To convert the observed projected separations of the companions on the sky to a projected physical separation, the distances to the targets are needed. While Gaia distance estimates are valuable for many targets, for bright stars and for binary systems they can be unreliable, especially when dealing with binary systems in the separation range covered here. Specifically, the astrometric measurements of the photocenter can be impacted by the scanning direction of Gaia with respect to the instantaneous orientation of the binary system on the sky. The Gaia measurements for bright, hot stars might further be subject to a number of systematics (Pantaleoni González et al. 2025). Some of the stars in our sample are also simply not in Gaia DR3, or have extreme uncertainties. To provide a homogeneous distance determination of the whole sample of stars, we developed an alternative strategy based on an absolute magnitude calibration as a function of spectral type (SpT) and luminosity class (LC). We then used the stars that do have reliable Gaia distances to assess the accuracy of our calibration.

2.1.1. Spectrophotometric distances

Our calibration was derived from the absolute magnitudes of Martins & Plez (2006, their table 2). As magnitudes are only given for LCs V, III, and I, we first derived magnitudes for luminosity classes IV and II by linearly interpolating between the closest luminosity classes. We then derived a relation for the absolute magnitude as a function of spectral type for each of the luminosity classes, so that we were able to derive absolute magnitudes for those spectral types not listed in Martins & Plez (2006, e.g., O2, O9.7). The procedure and resulting relations are described in detail in Appendix A. The relations were derived for the K band, but they can also be applied for the H band using (H − K)₀ = −0.10, which is valid for all spectral subtypes and luminosity classes in the O-star regime (Martins & Plez 2006). We adopted an uncertainty of half a spectral subtype and half a luminosity class to estimate the uncertainty on the absolute magnitude.

To test the above calibration, we applied it to both the primary and secondary in double-lined spectroscopic binary (SB2) systems that have been resolved by PIONIER, and for which the spectral type and luminosity class of the secondary is known. The resulting predicted magnitude difference, which only depends on the spectral classification of both components and is independent of extinction and apparent magnitudes, can then be compared to that observed by PIONIER, as shown in Fig. 2. In general, we find good agreement between the predicted and observed magnitude difference.

Fig. 2. Magnitude difference of SB2 systems resolved by PIONIER versus those obtained from the absolute magnitude calibration.

Before we can use the derived absolute magnitudes of the primaries to determine the distance, the apparent magnitudes need to be corrected for the contribution from nearby companions. The photometry listed in Table 1 of Paper I originates from the 2MASS survey, which has an aperture of 4″. We thus corrected these magnitudes for the observed magnitude difference of all resolved companions within 4″ (see Appendix B).

Apart from the resolved companions, the observed magnitudes may also contain a contribution from unresolved companions, i.e., the inner spectroscopic binaries. To correct the observed magnitude, we further used the calibration given in Appendix A to estimate the magnitude difference. Therefore, we could only do this for systems for which the spectral types of both components are known, i.e., double-lined spectroscopic binaries (SB2). This method cannot be applied to systems for which the companion spectral type could not be determined, mostly single-lined spectroscopic binaries (SB1), unless the spectroscopic binary is so wide that the secondary was resolved by PIONIER. This typically occurs for systems with orbital periods of a few months to a couple of years, and we could identify four cases in the SMASH+ database. For the remaining SB1 systems, however, the non-detection of the companion in the spectrum likely implies that it is faint compared to the primary star, and thus will only have a small contribution to the combined magnitude, and we thus neglected their contribution. These systems are marked in Table F.1 for reference.

As a last step, we corrected the magnitudes for extinction using the available VJHK-band photometry and the Fitzpatrick (1999) extinction curve (see Appendix C). We then calculated the distance to the system using the absolute H-band magnitude and the apparent H-band magnitude of the primary star corrected for extinction and the contribution from companions. The resulting distances are given in Table F.1.

2.1.2. Accuracy of the distance determination

To assess the accuracy of the derived distances, we compared the median distances of stars in clusters to literature values of the distance to those clusters. We applied this to each cluster that has at least four members in the SMASH+ sample. This comparison is shown in Fig. 3. We find that on average the derived distances are slightly larger than the literature values, with a mean offset of +186 pc for the 10 clusters.

Fig. 3. Derived distances for clusters with at least four members in the SMASH+ sample. Black stars indicate the distances for the individual stars and green stars the median. Blue stars indicate literature values for the cluster distances, and red stars the distances based on Tycho-2 (Kharchenko et al. 2005; Mel’Nik & Dambis 2009).

In Fig. 4, we also compare our individual distance with recently available distance estimates from the Alma Luminous Star III catalog (Pantaleoni González et al. 2025, ALS III). The latter relies on Gaia but applies a series of corrections appropriate for bright stars. We rejected stars with RUWE > 1.4 in this comparison as this is generally indicative of a poorer astrometric solution, as well as a small number of stars with ALSIII-Gaia distances over 4 kpc and large distance errors. While the dispersion is significant, most of it results from typical error sizes of ∼200 pc.

Fig. 4. Comparison of the spectro-photometric distances derived in this work to the Gaia distances from Pantaleoni González et al. (2025). The dashed line gives the 1:1 relation. The red symbols indicate the locations of the outliers discussed in the text while the blue line and associated shaded area are the systematic deviations presented in Eq. (1).

We observe a small but systematic offset of ⟨d_H − d_Gaia⟩ = 156 ± 22 pc. Applying a 3σ clipping to remove outliers (see below), this difference lowers to 111 ± 25 pc. While statistically significant, this offset is smaller than the mean uncertainties from either method. A closer look at the 9 outliers with > 3σ deviations between the two distance estimates provides insight into the physical cause of this disagreement.

Three objects that have an ALSIII-Gaia distance significantly larger than the spectro-photometric distance are all Oe stars (HD169515, O9.7Ibep; HD93190, O9.7:V:(n)e; LS4067A, O4.5Kfpe). These objects are too faint compared to the Martins’ calibrations that we adopted, which probably results from a lower surface brightness due to their rapid rotation or to partial extinction by their disks.

Six stars have an ALSIII-Gaia distance smaller than their spectro-photometric distance, i.e., they are brighter than expected from the calibration. Four of these six stars are supergiants (HD151018 O9Ib, HD154811, OC9.7Ib; HD156212, B0Iab; HD169582, O6Iaf) and one is a bright giant (HD171589 O7.5II(f)). The supergiants are known to display significant spread in their luminosities at a given SpT as objects of different masses evolve almost horizontally while nearing the end of the main sequence, so that supergiants of a given spectral subtype can have different luminosities.

The last system (HD148937, O6f?p) is a magnetic system where the magnetic star is believed to be a merger product (Frost et al. 2024). In addition, it has been suggested that some blue supergiants could also be the result of binary interaction, as, for example, suggested by the strong decrease in nearby companions for O9.7 supergiants in the VFTS survey (Sana et al. 2013).

Most of these outliers seem to have physical reasons to deviate from the calibrations. As an interesting corollary, the spectro-photometric calibration that we provide could be used in conjunction with independent distance estimates such as those from Gaia to single out candidate objects that have a different evolutionary history for further investigations.

In summary, both the cluster distances and the recent individual ALSIII-Gaia distances suggest that the spectro-photometric distances that we derived are slightly overestimated. For a mean distance of 1915 pc of the stars in the full sample, the systematic error corresponds to less than 10% (clusters) and by about 5% (ALSIII-Gaia). This suggests that the Martins calibrations are actually slightly too bright, which is an interesting finding per se. The systematics could be calibrated out by correcting the spectro-photometric distance estimate by the slope of the relation

$$\begin{array}{c}\hfill d_{\mathrm{H}}=-122.3\pm 64.7+(1.19\pm 0.04)d_{\mathit{Gaia}}.\end{array}$$(1)

However, the difference is small enough that it does not impact the empirical distributions of logarithmic (projected) separations (logSep) as the systematics only amount to ΔlogSep ≈ +0.02. As a consequence we proceed without applying the correction and use the spectro-photometric distances d_H to convert the angular separations ρ to projected physical separation as

$$\begin{array}{c}\hfill \left(\frac{\mathrm{Sep}}{\mathrm{AU}}\right)=\left(\frac{\rho }{\mathrm{mas}}\right)\left(\frac{d_H}{\mathrm{kpc}}\right).\end{array}$$(2)

The resulting projected separations are given in Table F.1.

2.2. Mass determination

In this section we provide an overview of the method to determine the masses of the primaries and of their companions. The details of the derivation, and the resulting calibrations, are given in Appendices D and E.

2.2.1. Primary masses

For the primary stars, we used the available spectral type and luminosity class to estimate the mass from the Martins et al. (2005) and Martins & Plez (2006) calibrations. First, H-band bolometric corrections were used to derive calibrations as a function of spectral type and luminosity class. Similar to the absolute magnitude calibration, we first derived H-band bolometric corrections for luminosity classes IV and II which are not in Martins & Plez (2006, see Fig. D.1). Then we derived relations for the bolometric correction for each of the luminosity classes as a function of spectral type (Fig. D.2). These relations were then used to convert the absolute H-band magnitude to the bolometric magnitude, and hence the absolute luminosities.

The masses of the primary stars can now be derived from the mass-luminosity relation. However, the exponent of the mass-luminosity relation changes both with luminosity, and with luminosity class (see Appendix E). For luminosity classes III and I, we used the Martins et al. (2005) calibration to derive the exponent of the mass-luminosity relation as a function of luminosity. For luminosity class V, we used the Brott et al. (2011) models evaluated at log g = 3.92, which is the value of the surface gravity g for all dwarf spectral subtypes in Martins et al. (2005). This allowed us to use the dwarf calibration at low luminosities, which could then also be applied to the companion stars. The latter are indeed expected to be predominantly dwarfs (see Sect. 2.2.2). For luminosity classes II and IV, we again interpolated between LC I–III and III–V, respectively. As above, we adopted half a luminosity class and half a spectral subtype to estimate uncertainties on the obtained masses. The resulting masses are given in Table F.1.

To assess the accuracy of the above approach, we compared the derived primary masses to literature values of SB systems with dynamical mass measurements (i.e., for systems with known inclination; see Table 2). This comparison is shown in Fig. 5. In general, our masses are in good agreement with the literature values. However, three stars (HD 115071, HD 100213, and HD 47129) deviate by more than 2σ. Interestingly, each of these stars is known to be either currently interacting or in a post-interaction state. HD 115071 is a semi-detached, over-luminous system, where mass-transfer has previously taken place (Penny et al. 2002). HD 100213 is a contact binary currently experiencing Roche-Lobe overflow (Penny et al. 2008). Lastly, HD 47129 (Plaskett’s star) is in a post-interaction state, and is a nitrogen-rich, evolved star (Linder et al. 2008) with a partially stripped donor and a magnetic accretor. That the relation that we derived from main-sequence calibrations does not work well for these interacting binaries is unsurprising, and we conclude that our calibration works well for pre-interaction systems, which is likely the case for the majority of the systems in our sample.

Fig. 5. Dynamical masses versus the masses derived in this work. Systems that are in a (post-)interaction state (see text) are indicated by red triangles. HD 135240 has an ambiguous luminosity class (III-V), the blue triangles indicate masses derived for LC III and V (20.3 and 26.7 M⊙, respectively). The dynamical mass of HD 159176 relies on the analysis of ellipsoidal variation and is quite uncertain.

2.2.2. Companion masses

As the luminosity classes of the companion stars are unknown, the above method cannot be directly applied to derive their masses. However, those companions that have a large magnitude difference with the primary should be considerably less massive, and therefore still be dwarfs regardless of the luminosity class of the primary. This is the case for the vast majority of our sample, and masses for these stars can directly be obtained from the calibration based on the Brott et al. (2011) models derived above.

However, for the brighter companions, the assumption of a dwarf luminosity may not be correct, in particular for systems with more evolved primaries. For example, for a system with an O9 I primary and a fainter companion with Δm = 1.0 mag, the assumption of a dwarf luminosity class for the companion would result in a mass ratio greater than unity. If this were true, the more massive companion should have evolved faster than the primary, and thus should have a more evolved luminosity class than the primary. Therefore, the assumption of a dwarf luminosity class does not hold. In these cases, we adopted the least evolved (i.e., the most massive for the given magnitude) luminosity class that gives a mass ratio smaller than unity, and adopted the resulting mass as an upper limit.

Figure 6 shows the distribution of adopted luminosity classes of the companions compared to the known luminosity classes of the primaries. The majority of companions have luminosity class V. Note that the distribution of primary luminosity classes is biased toward more evolved (visually brighter) stages as a result of the magnitude-limited approach of the survey (see Paper I), and hence is not expected to follow a canonical mass function.

Fig. 6. Histogram of the luminosity classes of the primaries and companions.

3. SMASH+ detection limits

The instrumental detection curves in the observational plane (ΔH versus ρ; e.g., Fig. 1) are well constrained and quite homogeneous throughout the survey (Paper I). However, each object has been observed only once, so that some binaries may have eluded detection due to unfavorable geometry or because they are (temporarily) within the inner working angle or beyond the outer angle limits of the instruments. In this section, we first investigate the relations between the instantaneous projected separation (Sep), the instantaneous true separation (r) and the semi-major axis (a) of the relative orbit. We then investigate the detection probability as a function of the orbital properties, accounting for the observational strategy of the survey.

3.1. Projection effects

The SMASH+ observations are snapshot measurements that provide the instantaneous angular separations ρ(t). The latter are easily related to the projected physical separation Sep(t) through Eq. (2). The relation between the instantaneous separation and the semi-major axis of the orbit is well known:

$$\begin{array}{c}\hfill \frac{r(t)}{a}=\frac{(1-e^2)}{1+ecos\nu (t)},\end{array}$$(3)

where e is the eccentricity and ν, the true anomaly. The ratio between the projected separation and the true separation is given by

$$\begin{array}{c}\hfill \frac{\mathit{Sep}}{r}(t)={(1-{sin}^{2}i{sin}^2(\omega +\nu (t)))}^{1/2},\end{array}$$(4)

where i is the orbital inclination, and ω the argument of periastron passage. The relation to the semi-major axis of the relative orbit a is (e.g., Sana & Vrancken 2026)

$$\begin{array}{c}\hfill \frac{Sep(t)}{a}=\frac{1-e^2}{1+ecos\nu (t)}{(1-{sin}^{2}i{sin}^2(\omega +\nu (t)))}^{1/2}.\end{array}$$(5)

This results in a rather complex dependence of the measured projected separation on the orbital properties e, i, and ω. Figure 7 illustrates these various effects and shows the nontrivial mapping of a single snapshot measurement to the physical properties of the orbits.

Fig. 7. Ratio of the projected instantaneous separation Sep to the instantaneous separation r (top panel) and the semi-major axis a (bottom panel) as a function of the orbital phase for various orbital geometries (see legend).

Fortunately, the sample is large enough that we can hope to marginalize over the random 3D orientation of the orbits, defined by the set of Kepler’s angles (i, ω, Ω), and over the distribution of eccentricities. The argument of the ascending node Ω is critical for the orientation of the projected orbit on the plane of the sky but has no impact on the distribution of projected separations so that we ignore it. Figure 8 provides the probability density functions (PDFs) of Sep/a marginalized over the various elements defining the orbital configurations (i, ω, e).

Fig. 8. Probability density functions (PDFs) of Sep/a for various orbital configurations. The bottom panel gives the PDF marginalized over the 3D orientation and the uniform eccentricity distribution between e = 0.0 and 0.9. The bottom panel also provides the boundaries of the 50 and 68% high density intervals (HDIs).

It illustrates the critical impact of eccentricities, which widen the distribution, while large, more probable inclinations make the distribution asymmetric. Assuming a random orientation of the orbit in 3D space (f_{cos i} ∝ 𝒰(−1, 1); f_ω ∝ 𝒰(0, π)), and a uniform distribution of the eccentricities (f_e ∝ 𝒰(0, 0.9)), the bottom panel provides the overall PDF marginalized over all orbital configurations and summarizes the overall likelihood that the measured Sep is representative of the size of the orbit given by a. This PDF is the one that we can use to assess the overall impact of the geometry. Fortunately that distribution has a mode at Sep/a = 1, and a median Sep/a ≈ 0.9. Thus, on average, the measured projected separation Sep is a reasonable estimator of the semi-major axis a of the relative orbit, albeit with a broad dispersion (50%-HDI=[0.6:1.1]). As discussed in Sect. 4.1, we aim to constrain the distribution of semi-major axis in log a. As for the discussion of distance estimates in Sect. 2.1, and given that log a ∝ log Sep (Eq. 5), the small systematic error introduced by using Sep as a biased estimator of a remains small compared to the signal. Incidentally, it even almost counterbalanced the systematics on the distance quantified by Eq. (1) and we proceed without further corrections.

3.2. Detection maps

We used the Monte Carlo detection probability approach of Sana et al. (2013), adapted to relative astrometry methods as in Frost et al. (2025). Unlike Sana et al. (2013) and Frost et al. (2025) however, we did not draw the masses of the primary stars from a given initial mass function but we adopted the distance, primary mass, and H-band absolute magnitude of the stars in the SMASH+ survey as given in Table F.1. Indeed the SMASH+ primary masses do not follow a canonical mass function because of the overrepresentation of supergiants in the sample (see Sect. 2.2). We simulated 10 000 mock SMASH+ samples, adopting a uniform eccentricity distribution between e = 0.0 and 0.9, and random orbital orientations in 3D space. At each draw, the central stars were paired with a companion with a mass-ratio q = M₂/M₁ drawn uniformly between 0.01 and 1, and we computed the magnitude contrast and projected separation of a random observational epoch. These mock observables were then compared to the survey sensitivity curves presented in Fig. 1 to decide whether the system would be detected as a binary or not. Throughout this process, we assumed that companions are all dwarfs, which is valid for the majority of the sample.

Figure 9 shows the resulting detection probability maps. The sharp boundaries at orbital periods of about 1.5 months to 10⁵ yr (equiv. semi-major axes of ∼1 to 10⁴ AU) result from the inner and outer working angles of the survey (approx. 0.001 and 8″). Within these boundaries, the overall detection probability nears 100%. Similarly, the detection probability is uniformly excellent above the horizontal cut-off defined by the contrast limits of ΔH ≈ 4, 5, and 8 in the interferometric, coronagraphic, and adaptive-optics imaging regimes, respectively. These limits translate to mass-ratio limits of q_lim ≈ 0.2 and 0.005, respectively, with q_lim values slightly dependent on the luminosity class (Fig. D.3).

Fig. 9. Binary detection probabilities of the SMASH+ survey projected on the mass ratio versus the orbital period (left), versus the semi-major axis (middle), and the companion mass versus the semi-major axis (right) planes. The colored background and solid equi-probability curves are based on the sample of stars that have been observed both by PIONIER and by NACO. The dashed equi-probability curves (and cyan labels) show the detectability when restricting the sample to ΔH ≤ 4 (i.e., the “cleaned” sample). From top to bottom: the full sample, luminosity classes I and II, and luminosity classes III to V.

The contrast limits of the instruments that we used depend on magnitude difference and not the absolute brightness of the companion. As a consequence, the limiting companion mass that we are able to detect depends on the brightness of the central star. For this reason, we also separated the SMASH+ sample according to the luminosity class into roughly two samples of similar size. For each subsample, we recomputed the detection maps. The last column of Fig. 9 reveals that interferometric observations can detect a 4 M_⊙ companion near an LC V-III star but only a 7 M_⊙ companion near a LC II-I star. In the AO regime, SMASH+ is sensitive to solar-mass companions near LC I-II stars and to subsolar mass companions near LC III-V stars.

Overall, our simulation results show that, within the sensitivity limits of the survey, the binary detection probabilities are very uniform. The one exception is the structure around P ≈ 10^3.5 yr (a ≈ 500 AU) that corresponds to the transition from NACO/SAM to NACO AO-imaging (Fig. 1). Given these results, the observed distributions obtained in the next sections can be considered to be directly representative of the true distributions, and do not require bias-corrections, as long as one remains in the high-detection-probability regions outlined in Fig. 9.

4. Results and discussion

By applying the methods presented in the previous section, the projected physical separations of all companions, and the masses of the primaries and of all their companions were determined. In this section we discuss these results.

4.1. Companion separations

The distribution of projected separations of all found companions within 8″ in the sample (Fig. 10) shows two separation ranges where the number of detected companions increases more sharply, at Sep ≈ 100 AU and ≈5000 AU. The first increase coincides with the switch between PIONIER and NACO/SAM at the average distance of the sample, and can hence be explained by the higher number of stars that have been observed with the latter instrument. The second increase occurs at the switch between the NACO/SAM and NACO/FOV modes and can be explained by the much deeper detection contrast of the latter mode (see Figs. 1 and 9). Many of the large number of faint objects at large separations are likely to be chance line-of-sight objects instead of true companions (see Paper I, Fig. 8).

Fig. 10. Cumulative distribution of the projected separations of all the detected companions within 8″ (blue), and of the detected companions with ΔH ≤ 4 mag, of systems that were observed by both PIONIER and SAM (green). Dashed lines indicate the drop in detector efficiency below ∼2 mas and at ∼300 mas at the average sample distance of 1915 pc.

To provide distributions that are less affected by selection biases, we defined a “cleaned” sample, which contains the 105 systems (297 companions, 288 within 8″) that have been observed by both instruments and for which we considered only detected companions brighter than a contrast limit of 4 magnitudes, resulting in a sample of 92 companions. Aside from two small separation ranges (below 2 mas (∼6 AU) and around 300 mas (∼650 AU)), the detection probability is indeed uniform within the adopted boundaries (Fig. 9). Furthermore, all companions in this sample are almost certainly physically bound instead of chance line-of-sight objects (see Paper I, Sect. 4.1).

Overall, the distribution of projected separations in the cleaned sample seems to follow an Oëpik law¹. There might be a small but significant deviation at ∼80 to 100 AU where we note a sharper increase in the number of detected companions. Either this is a remaining artifact of the transition between the PIONIER and SAM instruments, or it may be intrinsic. A larger sample would be desirable to confirm the presence of such a bump in the cumulative distribution function.

4.2. Masses

Figure 11 presents the cumulative distributions of the mass ratios, excluding 26 (7) companions in the full (cleaned) sample for which only upper limits on the mass ratios are derived. The overall distribution of mass ratios is skewed toward low-mass ratios, with almost 70% of the detected companions having q < 0.2 (Fig. 11, left). This abundance of low-mass ratio companions results from the large sensitivity of the NACO-AO observations at larger separations (Fig. 11, right). These companions disappear entirely when limiting ourselves to the inner 100 AU separation range as such low-mass-ratio companions are below the SMASH+ detection limit.

Fig. 11. Left: Cumulative distribution of the mass ratios of all the detected companions within 8″ (blue), and of the detected companions with ΔH ≤ 4 mag, of systems that were observed by both PIONIER and SAM (green). The red line shows a fit to the green distributions with its 95% confidence intervals. Middle: Same as on the left, but for separations < 100 AU. Right: Same as on the left, but for separations > 1000 AU.

Limiting ourselves to the cleaned sample, we used the Kuiper probability p_K as a goodness-of-fit metric to fit a power-law distribution f_q ∝ q^κ to the mass-ratio distribution. Outer envelopes are taken as the κ values for which p_K = 0.1. We obtained $\kappa \approx -0.7_{-0.6}^{+0.5}$, ${\kappa }_{<100}\approx -0.6_{-0.7}^{+0.9}$, and ${\kappa }_{>1000}=-1.6_{-1.2}^{+1.0}$ for the cleaned sample over the full separation, for Sep < 100 AU and Sep > 1000 AU, respectively. This quantifies that companions at larger separations have smaller mass ratios than if randomly drawn from a uniform distribution (see also discussion in Moe & Di Stefano 2017). One has to be careful to directly compare the best-fit exponent κ derived for Sep < 100 AU to the other two values as the lower boundary of the fitted range varies and it is known that the best-fit exponent value depends on the adopted range (Almeida et al. 2017). However, the differences here are large enough that we can conclude that the distribution is significantly more skewed at large separations compared to lower separations. To give an order of magnitude, the distributions in the inner (Sep < 100 AU) and outer (Sep > 1000 AU) ranges have medians at q ≈ 0.5 and 0.3, respectively.

Compared to mass-ratio distributions found for the tighter spectroscopic binaries (P ≲ 3000 d, a ≲ 10 AU) in the Milky Way (Sana et al. 2012, κ = −0.1 ± 0.6) and Large Magellanic Cloud (κ = +0.2 ± 0.2, Shenar et al. 2022), the index of the power-law distribution seems to decrease toward larger separations, a point already made by Moe & Di Stefano (2017). However, the κ value that we derived for the inner range ${\kappa }_{<100}=-0.6_{-0.7}^{+0.9}$ is much flatter than the value of κ_3 − 90 = −1.4 ± 0.4 derived by Moe & Di Stefano (2017) using a subsample of 21 SMASH+ targets. The reasons for such differences are unclear and may lie in different sample selection functions and survey sensitivity estimates combined with small sample statistics. One can argue that both values are within 2σ from one another, but the κ_{< 100} value that we derived would then also be compatible with the flat-mass ratio distribution of the tighter spectroscopic binaries. Determining the most appropriate mass-ratio distribution for intermediate-period binaries would be desirable as it impacts the outcome of case-B mass transfer, with potentially important differences for predictions of population synthesis computations.

Finally, we compare the distribution of observed masses of the primaries and the observed companions in Fig. 12. Also shown is a Salpeter initial mass function (IMF). As the primaries are all O-type stars, we only do this for masses greater than 16 M_⊙. This reveals a lack of primary stars with low masses (M ≤ 30 M_⊙) compared to the observed companions. This is a direct result of the magnitude-limited approach of the survey, causing an overabundance of (more massive) supergiants compared to dwarfs as these can be seen up to larger distances. Interestingly, the companions seem to follow the Salpeter IMF across most of the observed mass range, in agreement with the findings of GRAVITY Collaboration (2018) from a sample of 22 OB stars in the Orion Nebula.

Fig. 12. Cumulative distributions of the masses of all primaries (blue) and companions (red) with masses 16 M⊙ ≤ M ≤ 50 M⊙. Also shown is the Salpeter mass function in this mass range (green).

4.3. Separations versus mass ratios

Combining the results of Sects. 4.1 and 4.2 we can now show the results of Fig. 1 in physical units, i.e., projected separation versus mass ratio. This is done in Figs. 13 and 14, where the latter figure focuses on the cleaned sample limited to ΔH ≤ 4. While we note the absence of twin binaries (q > 0.95), the most striking result is probably the lack of high-mass-ratio companions at projected separations greater than 100 AU, i.e., in the top right area of the figure. At these separations, no companions with mass ratios greater than 0.8 are detected except for those with upper limits on their mass. This mass-ratio limit decreases toward larger separations, with no companions with mass ratios greater than ∼0.7 detected above 1000 AU, and no companions with mass ratios greater than ∼0.5 detected above 10⁴ AU. We tentatively derived an upper envelope delimiting the avoidance zone in the upper-right corner from the populated area below the curve:

$$\begin{array}{c}\hfill q_{\mathrm{crit}}\approx 1-0.2({log}_{10}\frac{\mathit{Sep}}{\mathit{AU}}-2)\mathrm{for}Sep>100\mathit{AU}.\end{array}$$(6)

Fig. 13. Projected separation versus the mass ratio for all the detected companions within 8″. Also indicated are the median detector sensitivities in the H band for an average sample distance of 1915 pc and assuming a dwarf primary.

Fig. 14. Same as Fig. 13 but for detected companions with ΔH ≤ 4 mag, for systems that were observed by both PIONIER and NACO.

It is beyond the scope of this paper to interpret these trends. Ultimately, the physical reasons for these trends are rooted in a complex set of factors. It is typically assumed, however, that components of binaries wider than a few hundred AU form nearly independently of one another (e.g., White & Ghez 2001; Moe & Di Stefano 2017; El-Badry et al. 2019; El-Badry & Rix 2019). Components that assembled at closer separations probably form in a highly correlated way. A decline in the maximum mass ratio q_crit for massive stars is also reported by Moe & Di Stefano (2017, Fig. 4), though the onset of the decline occurs at a much shorter separation (at about an orbital period of 10^1.5 d).

5. Summary and conclusions

In this paper we converted the SMASH+ observational measurements of Paper I into physical units, i.e., observed angular separations on the sky to projected physical separations, and observed H-band magnitude differences to mass ratios. To do this, we derived calibration relations based on physical properties of O-type stars presented in Martins et al. (2005) and Martins & Plez (2006), as well as the evolutionary models of Brott et al. (2011) for lower-mass stars.

The accuracy of the resulting distances was assessed by comparing them to known cluster distances, and to up-to-date ALSIII-Gaia distances. Overall, we find good agreement, with the spectrophotometric distances slightly overestimating the Gaia distances by just 5%, except for Oe stars, for a small subset of supergiants, and for binary evolution products. For the derived masses, primary stars with known dynamical masses were used to assess the accuracy of the relations used. We find excellent agreement, except for known (post-)interaction systems. Although this was not the original aim, our calibrations can be used with independent mass and/or distance estimates to identify peculiar stars whose physical parameters are potentially affected by binary evolution.

The derived projected separations follow an Oëpik law, i.e., a flat logarithmic distribution. The distribution of mass ratios follows a power-law distribution f_q ∝ q^κ, with κ changing with the separation range considered. For companions within 100 AU, the value that we derived ($\kappa \approx -0.6_{-0.7}^{+0.9}$) lies between the uniform mass-ratio distributions derived for spectroscopic binaries, and the steeper κ = −1.4 ± 0.4 value reported by Moe & Di Stefano (2017) based on the analysis of a subsample of the SMASH+ data. The value that we derived is compatible within errors with both regimes and a larger data set is probably needed to reach a firm conclusion. The distribution of companion masses seems to be well reproduced by a Salpeter initial mass function. Finally, we find an absence of massive companions at separations larger than 100 AU.

Data availability

Full Table F.1 is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/711/A236.

Acknowledgments

This project has received funding from the European Research Council under European Union’s Horizon 2020 research programme (MULTIPLES, No 772225), and from the KU Leuven Research Council (grant METH/24/012: SOUL and C16/17/007: MAESTRO). FT acknowledges support by grant PID2022-137779OB-C41, funded by the Spanish Ministry of Science, Innovation and Universities/State Agency of Research MICIU/AEI/10.13039/501100011033. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This work made use of Astropy (http://www.astropy.org): a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022).

References

Appendix A: Absolute magnitude calibration

To obtain the absolute K- and H-band magnitudes that are needed for the distance calculations, we used the Martins & Plez (2006) calibration. As Martins & Plez (2006) only gives values for luminosity classes V, III, and I, the first step was to obtain values for luminosity classes II and IV. To do this, we fitted a second degree polynomial of the form

$$\begin{array}{c}\hfill M=a\times {\mathrm{LC}}^2+b\times \mathrm{LC}+c\end{array}$$(A.1)

to the magnitudes of each of the spectral subtypes listed in Martins & Plez (2006). Here, LC is an integer in the range 1 to 5, representing the luminosity classes I to V.

Table A.1 lists the obtained values of a, b, and c for each of the spectral subtypes listed in Martins & Plez (2006). Figure A.1 shows the corresponding LC versus M_K relations.

Fig. A.1. LC versus MK calibration for each of the spectral subtypes listed in Martins & Plez (2006).

We used an identical approach to obtain relations for the absolute magnitude as a function of the spectral subtype (SpT). This was done to obtain absolute magnitudes for the spectral subtypes not listed in Martins & Plez (2006), that is, O2, O2.5, O3.5, O9.2, and O9.7. We fitted a second-degree polynomial of the form

$$\begin{array}{c}\hfill M=a\times {\mathrm{SpT}}^2+b\times \mathrm{SpT}+c\end{array}$$(A.2)

to the absolute magnitudes of each of the luminosity classes. For LC V, III, and I we used the Martins & Plez (2006) values, and for LC IV and II the values obtained from the LC versus M_K calibration presented above. This resulted in five SpT versus M_K relations, which are presented in Fig. A.2. The corresponding coefficients of the polynomials are given in Table A.2.

Fig. A.2. SpT versus MK calibration for each of the luminosity classes.

The obtained relations also provide the absolute H-band magnitudes using (H − K)₀ = −0.10, which is valid for all spectral subtypes and luminosity classes in the O-star regime (Martins & Plez 2006).

Appendix B: Magnitude correction

Before the distance can be derived from the H-band magnitude, this magnitude first needs to be corrected for the contribution to the 2MASS magnitude from companions within the 4″ aperture. The next subsection briefly describes this process for companions that are interferometrically resolved, and thus have a measured magnitude difference. Section B.2 describes the correction for the unresolved companions (i.e., the SBs).

B.1. Correction for resolved companions

We corrected for all observed companions that are within 4″ of the primary star, and for which either Δm_Ks or Δm_H is available. As we derived the distance from the H-band magnitude, we assume that Δm_Ks = Δm_H, which is valid when the companion is in the Rayleigh-Jeans domain in the H and K bands. While this assumption may not hold for the faintest companions, these stars have such a large magnitude difference from the primary that their contribution to the observed magnitude is negligible.

For a single companion, the corrected magnitude (m₁) of the central star is given by

$$\begin{array}{c}\hfill m_1=m_{\mathrm{obs}}+2.5log(1+{10}^{-0.4\mathrm{\Delta }m}),\end{array}$$(B.1)

where m_obs is the observed magnitude, and Δm the magnitude difference between the central star and the companion. If multiple companions are present within the aperture, the correction was iterated for each of the components.

B.2. Correction for unresolved companions (SBs)

Apart from the resolved companions, the observed magnitude also contains a contribution from the unresolved companions, i.e., the spectroscopic binaries. To correct the magnitude for this, we used the calibration given in Appendix A to estimate the magnitude difference. Therefore, we could only do this for systems for which the spectral types of both components are known. Thus, we could not correct for SB1 systems, apart from four SB1 companions that have been resolved by PIONIER. However, that the companion is not seen in the spectrum implies that it is likely faint compared to the primary star, and thus will only have a small contribution to the combined magnitude, i.e., the derived distance. These systems have been marked in Table F.1.

Appendix C: Extinction correction

To estimate the extinction toward each of the stars, we used the available VJHK-band photometry from 2MASS. Using Eqs. A3, A4, and A5 from Fitzpatrick (1999) we derived the following relations for E(B − V):

$$\begin{array}{c}\hfill E(B-V)=1.39\times E(V-J)/(R+0.02)\end{array}$$(C.1a)

$$\begin{array}{c}\hfill E(B-V)=1.19\times E(V-H)/(R-0.04)\end{array}$$(C.1b)

$$\begin{array}{c}\hfill E(B-V)=1.12\times E(V-K)/(R-0.02)\end{array}$$(C.1c)

To use Eq. C.3, we applied the average offset of 0.04 mag between the 2MASS K_s band and the K band, i.e., K_s = K + 0.04. We assumed the average value of the total-to-selective extinction of R = 3.1 for all systems. Figure C.1 shows the derived values of E(B − V) for each of the systems. Here, the given value is the mean E(B − V) derived from the three colors, and the error bars correspond to the standard deviation.

Fig. C.1. Derived E(B − V) for all stars as a function of the derived distance.

The resulting values of E(B − V) were then used to derive A_H and A_K using the model curve values of (Fitzpatrick 1999, their Table 2), which results in

$$\begin{array}{c}\hfill A_{\mathrm{K}}=0.36\times E(B-V)\end{array}$$(C.2a)

$$\begin{array}{c}\hfill A_{\mathrm{H}}=0.53\times E(B-V)\end{array}$$(C.2b)

Appendix D: Bolometric corrections

To convert the absolute magnitudes to bolometric magnitudes and hence the luminosity, we again used the Martins & Plez (2006) calibration to derive relations for the bolometric correction in the H band (BC_H). As for the absolute magnitude calibration, we first fitted BC_H as a function of the luminosity class using

$$\begin{array}{c}\hfill B{C}_H=a\times {\mathrm{LC}}^2+b\times \mathrm{LC}+c.\end{array}$$(D.1)

The results are shown in Fig. D.1 and the derived coefficients given in Table D.1. Using these results, we derived BC_H for luminosity classes IV and II, which are not given in Martins & Plez (2006). For each luminosity class we then fitted BC_H as a function of spectral subtype using

$$\begin{array}{c}\hfill B{C}_H=a\times {\mathrm{SpT}}^2+b\times \mathrm{SpT}+c.\end{array}$$(D.2)

Fig. D.1. LC versus BCH calibration for each of the spectral subtypes listed in Martins & Plez (2006).

The resulting relations are shown in Fig. D.2 and the coefficients given in Table D.2.

Fig. D.2. SpT versus BCH for each of the luminosity classes.

Finally, to also be able to derive BC_H for the companion stars, we fitted a second-degree polynomial for BC_H as a function of the absolute H-band magnitude M_H for each of the luminosity classes:

$$\begin{array}{c}\hfill B{C}_H=a\times M_H^2+b\times M_H+c.\end{array}$$(D.3)

The resulting coefficients are given in Table D.3 and the relations are plotted in Fig. D.3.

Fig. D.3. MH versus BCH for each luminosity class.

Appendix E: Mass-luminosity relation

To estimate the masses of the primaries and companions we used the mass-luminosity relation,

$$\begin{array}{c}\hfill \frac{L}{L_{\odot }}={\left(\frac{M}{M_{\odot }}\right)}^x.\end{array}$$(E.1)

The value of the exponent x varies with the luminosity and evolutionary state, which is shown in Fig. E.1 using the Brott et al. (2011) evolutionary models and the Martins et al. (2005) parameters. The exponent could be fitted well by a fourth-degree polynomial for the Brott et al. (2011) models, i.e.,

$$\begin{array}{c}\hfill x=a\times {\left(\frac{L}{L_{\odot }}\right)}^4+b\times {\left(\frac{L}{L_{\odot }}\right)}^3+c\times {\left(\frac{L}{L_{\odot }}\right)}^2+d\times \left(\frac{L}{L_{\odot }}\right)+e.\end{array}$$(E.2)

Fig. E.1. Exponents of the mass-luminosity relations derived from Brott et al. (2011) and Martins et al. (2005).

Fig. E.2. Derived masses from the adopted mass-luminosity relations for the primary stars (star symbols) compared to Martins et al. (2005) values (triangles) for each luminosity class.

The Martins et al. (2005) relations are nearly linear, and were fitted using a second-degree polynomial (i.e., a = b = 0). The coefficients of Eq. E.2 for both the Brott et al. (2011) models and the Martins et al. (2005) parameters are given in Table E.1.

Appendix F: Derived distances, separations, and masses

The distances, separations, and masses derived for each of the program stars and detected companions are given in Table F.1. The full table is available on CDS.

All Tables

All Figures

	Fig. 1. Magnitude difference versus angular separation for the pairs detected in the SMASH+ survey. Circled objects are detected by both SAM and PIONIER. The solid lines indicate the median H-band sensitivity of the various instruments. Figure reproduced from Sana et al. (2014).
In the text

	Fig. 2. Magnitude difference of SB2 systems resolved by PIONIER versus those obtained from the absolute magnitude calibration.
In the text

	Fig. 3. Derived distances for clusters with at least four members in the SMASH+ sample. Black stars indicate the distances for the individual stars and green stars the median. Blue stars indicate literature values for the cluster distances, and red stars the distances based on Tycho-2 (Kharchenko et al. 2005; Mel’Nik & Dambis 2009).
In the text

	Fig. 4. Comparison of the spectro-photometric distances derived in this work to the Gaia distances from Pantaleoni González et al. (2025). The dashed line gives the 1:1 relation. The red symbols indicate the locations of the outliers discussed in the text while the blue line and associated shaded area are the systematic deviations presented in Eq. (1).
In the text

Fig. 5. Dynamical masses versus the masses derived in this work. Systems that are in a (post-)interaction state (see text) are indicated by red triangles. HD 135240 has an ambiguous luminosity class (III-V), the blue triangles indicate masses derived for LC III and V (20.3 and 26.7 M⊙, respectively). The dynamical mass of HD 159176 relies on the analysis of ellipsoidal variation and is quite uncertain.

In the text

	Fig. 6. Histogram of the luminosity classes of the primaries and companions.
In the text

	Fig. 7. Ratio of the projected instantaneous separation Sep to the instantaneous separation r (top panel) and the semi-major axis a (bottom panel) as a function of the orbital phase for various orbital geometries (see legend).
In the text

	Fig. 8. Probability density functions (PDFs) of Sep/a for various orbital configurations. The bottom panel gives the PDF marginalized over the 3D orientation and the uniform eccentricity distribution between e = 0.0 and 0.9. The bottom panel also provides the boundaries of the 50 and 68% high density intervals (HDIs).
In the text

Fig. 9. Binary detection probabilities of the SMASH+ survey projected on the mass ratio versus the orbital period (left), versus the semi-major axis (middle), and the companion mass versus the semi-major axis (right) planes. The colored background and solid equi-probability curves are based on the sample of stars that have been observed both by PIONIER and by NACO. The dashed equi-probability curves (and cyan labels) show the detectability when restricting the sample to ΔH ≤ 4 (i.e., the “cleaned” sample). From top to bottom: the full sample, luminosity classes I and II, and luminosity classes III to V.

In the text

	Fig. 10. Cumulative distribution of the projected separations of all the detected companions within 8″ (blue), and of the detected companions with ΔH ≤ 4 mag, of systems that were observed by both PIONIER and SAM (green). Dashed lines indicate the drop in detector efficiency below ∼2 mas and at ∼300 mas at the average sample distance of 1915 pc.
In the text

Fig. 11. Left: Cumulative distribution of the mass ratios of all the detected companions within 8″ (blue), and of the detected companions with ΔH ≤ 4 mag, of systems that were observed by both PIONIER and SAM (green). The red line shows a fit to the green distributions with its 95% confidence intervals. Middle: Same as on the left, but for separations < 100 AU. Right: Same as on the left, but for separations > 1000 AU.

In the text

	Fig. 12. Cumulative distributions of the masses of all primaries (blue) and companions (red) with masses 16 M⊙ ≤ M ≤ 50 M⊙. Also shown is the Salpeter mass function in this mass range (green).
In the text

	Fig. 13. Projected separation versus the mass ratio for all the detected companions within 8″. Also indicated are the median detector sensitivities in the H band for an average sample distance of 1915 pc and assuming a dwarf primary.
In the text

	Fig. 14. Same as Fig. 13 but for detected companions with ΔH ≤ 4 mag, for systems that were observed by both PIONIER and NACO.
In the text

	Fig. A.1. LC versus MK calibration for each of the spectral subtypes listed in Martins & Plez (2006).
In the text

	Fig. A.2. SpT versus MK calibration for each of the luminosity classes.
In the text

	Fig. C.1. Derived E(B − V) for all stars as a function of the derived distance.
In the text

	Fig. D.1. LC versus BCH calibration for each of the spectral subtypes listed in Martins & Plez (2006).
In the text

	Fig. D.2. SpT versus BCH for each of the luminosity classes.
In the text

	Fig. D.3. MH versus BCH for each luminosity class.
In the text

	Fig. E.1. Exponents of the mass-luminosity relations derived from Brott et al. (2011) and Martins et al. (2005).
In the text

	Fig. E.2. Derived masses from the adopted mass-luminosity relations for the primary stars (star symbols) compared to Martins et al. (2005) values (triangles) for each luminosity class.
In the text