© The Authors 2023

1. Introduction

The formation mechanism of massive stars is still open to debate in present-day astronomy (e.g. Tan et al. 2014). The reason for this is that direct observations of forming massive stars are challenging, due to the rareness of such stars and the large distances at which they are found, their short formation timescale, and the fact that their formation stages are shrouded in dense environments obscured by dust and gas (e.g. Zinnecker & Yorke 2007).

Various massive star formation theories have been proposed (e.g. Tan et al. 2014). For example, massive stars may form through the merging of lower-mass stars if the density of stars is high enough to enable stellar encounters (Bonnell et al. 1998). Most theories, however, are based on accretion, such as monolithic collapse (McKee & Tan 2003; Krumholz et al. 2009) and competitive accretion (Bonnell et al. 2001; Bonnell & Bate 2006). These models consider the option that massive stars are not just formed by the accretion of cloud material, but that the accretion is channeled through a massive disk, allowing for the radiation pressure to be overcome. When these disks become gravitationally unstable, disk fragmentation likely leads to the formation of (low-mass) companions (Kratter & Matzner 2006; Kratter et al. 2010; Oliva & Kuiper 2020).

One of the most remarkable properties of massive stars is that they form almost exclusively in multiple systems (Sana & Evans 2011; Sana et al. 2013, 2014). The multiplicity characteristics of massive stars thus hold crucial information to unravel the secrets of their formation, since proposed massive star formation theories should be able to reproduce the observed multiplicity characteristics.

Several studies have characterised the multiplicity properties of massive stars at various separations, using spectroscopy, interferometry, or classical imaging. Mason et al. (2009), Sana et al. (2012, 2013), and Kobulnicky et al. (2014) identified small separation binaries and found multiplicity fractions close to 70% for O-type stars. The Southern Massive Stars at High-angular resolution survey (SMaSH+, Sana et al. 2014) used optical interferometry from the Very Large Telescope (VLT) with the Precision Integrated-Optics Near-infrared Imaging ExpeRiment (PIONIER), Sparse Aperture Masking (SAM) with the Nasmyth Adaptive Optics System combined with the Near-Infrared Imager and Spectrograph (VLT/NACO), and classical imaging in the entire NACO field of view (FoV) to search for companions around a sample of O-type stars between 1 mas and 8″. They reached contrast magnitudes up to ΔH < 4 at close separations and up to ΔH < 8 at larger separations. More than 90% of the O-type stars in their sample have at least one companion and 60% are members of a multiple system containing three or more stars.

These previous studies have not been able to probe the low-mass end of the companion mass function, as this would require even larger contrasts. Nevertheless, it is crucial to explore the parameter space of the faintest companions in order to estimate the total multiplicity frequency and to establish the shape of the mass-ratio and separation distribution.

The ongoing Carina High-contrast Imaging Project of massive Stars (CHIPS, Rainot et al. 2020) proved that with a combination of high-contrast imaging and extreme adaptive optics (VLT/SPHERE, Spectro-Polarimetric and High-contrast Exoplanet REsearch instrument) faint companions (down to contrasts of ≈12 mag) around bright O-type stars at angular separations between 0$\stackrel{″}{.}$15 and 6″ can be detected. This separation regime approximately corresponds to the size of the accretion disk, where we expect to find low-mass companions formed through disk fragmentation. Therefore, high-contrast imaging is essential to investigate an unexplored domain of stellar companions and acquire a more complete overview of the multiplicity properties of massive stars.

In this paper we study 20 O- and early B-type stars in the Scorpius OB1 association, which is an active star forming region located at 1530 pc (Sana et al. 2005), through VLT/SPHERE observations, thus probing the region from 0$\stackrel{″}{.}$15 to 6″ at contrast magnitudes down to ≈12 mag. In this paper, we focus on the stellar-mass companions, that is companions with an estimated mass M > 0.08 M_⊙.

Compared to the Carina region, Scorpius OB1 has an older population of stars (≲4 Myr for Carina and ≈4 − 8 Myr for Scorpius Sung et al. 2013; Damiani et al. 2016). This provides us with an ideal laboratory to investigate whether and how the multiplicity characteristics of massive stars change with age, for instance due to stellar interactions, dynamical evolution in clusters or inward and outward migration of the companions during their lifetime. In addition, some wide binaries might be disrupted, leaving a smaller number of companions for older stars. Sco OB1 is also a much looser, less active association. The comparison of Carina with Sco OB1 will thus allow us to investigate whether the multiplicity properties in the ∼200 − 9000 AU range are universal or potentially depend on age or environment.

2. Sample and distance

The sample consists of 13 O-type and 7 B-type stars. The initial sample was assembled by selecting all O-type stars and B-supergiants within a 3°-diameter region around NGC 6231 under ESO programme ID 0103.C-0803(A) and 105.20QJ.001, yielding 39 objects. Seventeen of these were observed. In addition, we re-analysed three O-type dwarf stars from Reggiani et al. (2022) which are also in Sco OB1 (ESO programme ID 095.D-0495(A)), but were not re-observed, yielding a total observed sample of 20 objects.

We used the Gaia EDR3 geometric and photogeometric distance estimates to check for membership (Fig. 1, Bailer-Jones et al. 2021). The geometric distance is based on the parallax, while the photogeometric distance takes into account parallax, colour, and magnitude (Gaia EDR3, Gaia Collaboration 2016, 2021). The median geometric and photogeometric distances are 1538 pc and 1517 pc, respectively. These are in excellent agreement with the distance measurement of ${1528}_{-109}^{+117}$ to the eclipsing binary CPD −41°7721 which is part of NGC 6231, the most prominent young star cluster in Sco OB1 (Sana et al. 2005). Figure 1 shows that there are four outliers whose geo- or photogeometric distances are not consistent within 3σ with the distance measurement by Sana et al. (2005). Firstly, HD 152 685 (B1 Ib) has geometric and photogeometric distances of ${905}_{-25}^{+24}$ pc and ${900}_{-19}^{+20}$ pc, respectively. These measurements have small errorbars and agree well with each other, so we conclude that HD 152 685 is likely a foreground star that is not part of Sco OB1. Secondly, HD 152 623 has much larger geo- and photogeometric distances of ${5276}_{-2214}^{+3422}$ pc and ${5116}_{-2263}^{+2674}$ pc. However, the uncertainties are very large and the high RUWE value of 17.9 indicates a poor astrometric fit, leading to erroneous values (Gaia Collaboration 2016, 2021). The star has been found to be a binary by Reggiani et al. (2022), which could have affected the Gaia measurement. We consider it part of the cluster for the analysis presented in this paper. Similarly, ζ¹ Sco (RUWE = 0.82) and HD 151 804 (RUWE = 1.35) have geometric and/or photogeometric distance measurements outside of the 3σ interval around the Sana et al. (2005) measurement, yet their Gaia uncertainties are large enough to be compatible with the eclipsing binary measurement at 1σ and 2σ, respectively, so that we consider them cluster members too.

Fig. 1. Geometric and photogeometric distances for stars in the sample from Gaia EDR3 (Bailer-Jones et al. 2021). The dashed lines are the median geometric and photogeometric distances of the sample, while the dash-dotted line indicates where the geometric distance is equal to the photogeometric distance.

2.1. HR diagram

We positioned the stars in a Hertzsprung-Russel (HR) diagram (Fig. 2). To do so, we estimated their effective temperature based on their spectral type (Martins et al. 2005; Arnó 2006) or from literature values (references in Sect. 4.2.4). The B and V band magnitudes are taken from the Tycho-2 catalogue (Høg et al. 2000), the bolometric corrections from Martins et al. (2005) for O-type stars and Straizys & Kuriliene (1981) for B-type stars, and (B − V)₀ is taken from Martins et al. (2005) for O-type stars and Straižys (1992) for B-type stars. We converted the Tycho-2 BT and VT magnitudes to Johnson B and V using the conversion formula from ESA (1997). The distance modulus was calculated assuming all stars except HD 152 685 are part of the cluster and located at 1.53 kpc. For HD 152 685 the Gaia photogeometric distance of 900 pc was used. We took into account the uncertainties on the B and V magnitudes and on the distance measurement. For the temperature, we assumed a 1σ uncertainty of 500 K. Comparison with evolutionary tracks and isochrones allowed us to obtain an estimate of their mass and age, assuming the stars have been evolving as single stars (Brott et al. 2011). Most ages range between 3 and 8 Myr, which is in agreement with previous age estimates (Sung et al. 2013; Damiani et al. 2016). Similarly, the initial masses range from 10 to 50 M_⊙, with a median value ∼25 M_⊙.

Fig. 2. HR diagram for all targets in the sample. Evolutionary tracks and isochrones are from Brott et al. (2011) at galactic metallicity and initial rotation of ∼110 km s−1. The red line indicates the 20 M⊙ evolutionary track. Open and closed symbols are stars below and above 20 M⊙, respectively, used for the low and high mass bin in the calculation of the companion fractions (Sect. 6.2.2).

3. Observations and data reduction

The observations were obtained with the Spectro-Polarimetric and High-contrast Exoplanet REsearch instrument (SPHERE, Beuzit et al. 2019). SPHERE is positioned at the Unit 3 (UT3) Nasmyth focus of ESO’s Very Large Telescope (VLT) and includes an adaptive optics (AO) system, coronagraphic masks, and various sub-systems, two of which were used for the purpose of this paper: the Integral Field Spectrograph (IFS, Claudi et al. 2008; Mesa et al. 2015) and the Infra-Red Dual-band Imager and Spectrograph (IRDIS, Dohlen et al. 2008; Vigan et al. 2010). The observations were executed in IRDIS and IFS extended mode (IRDIFS_EXT), simultaneously using IFS and IRDIS.

The IFS science frames contain 290 × 290 spaxels that with a plate scale of 0.0074″/pix probe an area of 1$\stackrel{″}{.}$73 × 1$\stackrel{″}{.}$73 on the sky. In total, there are 39 wavelength bands between 0.95 and 1.65 μm. The IRDIS observations consist of frames of 1024 × 1024 spaxels with a pixel size of 0.01225″, resulting in a FoV of 11″ × 12$\stackrel{″}{.}$5 on the sky. Each target was observed in two wavelength bands (K₁ and K₂) at 2.11 μm and 2.25 μm. The IFS observations are useful to detect and characterise companions at close to intermediate separations (less than ∼0.9″), while the IRDIS observations give information about companions at larger separations and the overall local source density in the field around the central object.

The observing sequence used for the data acquisition with SPHERE consisted of three different types of observations. Firstly, FLUX (F) observations were made, which are off-axis observations of the central star, allowing to obtain a point spread function that was later used to calibrate the companion spectrum. A neutral density filter (ND2.0) was used to avoid saturation. Secondly, CENTER (C) frames were obtained, which were used to centre the coronagraph on the central star. Thirdly, OBJECT (O) observations were taken with a coronagraph that blocks the light from the central star to increase the contrast at close separations. The observing setup, atmospheric conditions, and parallactic angle (PA) variation are shown in Table 1 for IFS and Table 2 for IRDIS observations. For the coherence time τ₀, we provide the average value during each observation.

The data reduction of IRDIS and IFS data was performed by the SPHERE Data Centre (Delorme et al. 2017; Galicher et al. 2018) at the Institut de Planetologie et d’Astrophysique de Grenoble (IPAG)¹. After wavelength calibration, the data are organised in a 4D cube. These four dimensions are the two spacial dimensions, wavelength channels, and sky rotation. The latter is obtained by setting the telescope to pupil tracking mode, meaning that the central star is in the centre of the image and the sky rotates around it, and by taking multiple images of the object during the observing sequence.

All frames were corrected for bad pixels, dark, and flat field. In addition, the SPHERE Data Centre carries out the astrometric calibration following Maire et al. (2016). They applied a true north correction of −1.75″ ± 0.08″ and used a plate scale of 12.255 ± 0.009 mas/pixel for IRDIS and 7.46 ± 0.02 mas/pixel for IFS (Delorme et al. 2017).

4. Data analysis

The post-processing of the data was carried out using the python Vortex Image Processing package VIP (Gomez Gonzalez et al. 2017)². VIP is developed to analyse high-contrast imaging data sets and allows us to perform angular (ADI, Marois et al. 2006) and spectral (SDI) differential imaging through a principal component analysis (PCA, Amara & Quanz 2012; Soummer et al. 2012) approach.

We analysed the IFS images through SDI, taking advantage of the 39 wavelength bands for IFS. In addition, the spectrum of sources in the IFS images was retrieved by performing ADI on a channel by channel base. For the IRDIS images, we used ADI.

4.1. Source detection

For IFS detections, we computed the signal-to-noise ratio (S/N) of each pixel in a post-processed frame with VIP (Gomez Gonzalez et al. 2017), where the noise is computed in an annulus around the central star at the radius of that pixel with width equal to the full width at half maximum (FWHM) of the central star. We considered a true signal detection when the S/N was above 5.

However, the S/N map becomes problematic when multiple stars are located at the same radial distance (±1 FWHM), since the noise measurement in a certain annulus is contaminated by the signal of all the sources found at the same radial distance. For instance, when a bright star is detected at the same radius as a fainter star, the noise will be overestimated, such that the S/N of the fainter star might appear below the detection threshold.

Hence, for IRDIS detections, we instead calculated the standardised trajectory intensity mean (STIM) map (Pairet et al. 2019), which is unaffected by sources at the same angular separation. However, selecting the detection threshold is less straightforward than for the S/N map. Pairet et al. (2019) propose to select the detection threshold by calculating the inverse STIM map (using opposite values of the parallactic angles), which has a similar temporal dependence of the residual speckles as the STIM map, but true signals will be averaged to negligible values. For different thresholds, they determine the number of pixels n_ϵ in the inverse STIM map that are above a certain threshold ϵ. We took the smallest ϵ, so that n_ϵ/n ≤ 10⁻⁵, with n the total number of pixels in the detection map. For some targets, the region close to the central star was polluted by residual star light. Therefore, we adopted separate thresholds for the region within ∼2″ and the region outside ∼2″ for these targets. The post-processed SPHERE/IFS PCA-SDI (for targets with companions in the IFS FoV) and SPHERE/IRDIS PCA-ADI images are presented in Figs. 3 and 4, respectively. In order to calibrate the separation from arcseconds to AU, a distance of 1.53 kpc was adopted (Sana et al. 2005), except for HD 152 685 for which the photogeometric distance of 0.9 kpc was used. The number of principal components in the IFS SDI images was chosen such that the companions have a high S/N and are clearly recognisable in the images.

Fig. 3. Post-processed SPHERE/IFS PCA-SDI images with one principal component for HD 152 623, 10 principal components for CPD −41°7721, HD 152 042, HD 152 200, and HD 152 385, and 40 principal components for HD 151 805. The central star is at the centre of the images (‘+’).

Fig. 4. Post-processed SPHERE/IRDIS PCA-ADI images with one principal component. The wavelength band of the images is K1.

Fig. 4. continued.

4.2. Source characterisation

4.2.1. Position and flux estimation for sources within 2″

The position and flux of a source were estimated by evaluating the flux wavelength by wavelength, so that we obtained a spectrum of 39 wavelength bands for IFS sources. Since these sources are also in the FoV of IRDIS, two extra wavelength bands were added to form a total of 41 channels. For sources that are outside of the IFS FoV but within the IRDIS FoV, we only have an estimate of their K₁- and K₂-band fluxes.

We characterised sources within 2″ from the central massive star with VIP using the negative fake companion technique (NEGFC, Marois et al. 2010; Lagrange et al. 2010), which injects a PSF template (obtained from FLUX observations of the central star) with negative flux into a post-processed cube and attempts to cancel out the companion as well as possible. Firstly, the coordinates of the companion were estimated in the post-processed frames by eye and a first guess for the flux was obtained by performing aperture photometry on a derotated and median combined frame. Secondly, these first guesses for the position and flux were passed to the algorithm that injects negative fake companions in this location, keeping the angular separation ρ and position angle θ (measured from north to east) fixed, but varying the flux. Finally, a Simplex Nelder-Mead optimisation was performed where ρ, θ, and the flux were allowed to vary.

The uncertainties were calculated through a Monte Carlo (MC) method for every wavelength by injecting fake companions with the same separation and flux in a pre-processed cube, while varying the position angle. Consequently, we post-processed the cube through PCA-ADI and we recovered the flux and position of the fake companion using the same method as described above. This process was repeated for ten different angles. Finally, the uncertainties on the separation, position angle, and flux were computed as the standard deviation of the recovered values with respect to the injected ones.

4.2.2. Position and flux estimation for sources beyond 2″

For sources beyond 2″ (in IRDIS data), the contribution of light from the central star and the effects of the coronagraph can be neglected as the background noise dominates. Therefore, it is not needed to perform time-consuming ADI and SDI techniques in order to characterise the companions. In addition, sources that are too close to the edges of the frames cannot be characterised through the negative fake companion technique, since it works in an annulus around the central star. Hence, for sources beyond 2″ we adopted a PSF fitting technique, as described in Bodensteiner et al. (2020) and already applied to SPHERE/IRDIS observations in Rainot et al. (2020, 2022) and Reggiani et al. (2022). The IRDIS images were first derotated and median combined, so that two images remained for K₁ and K₂. Subsequently, the FLUX observations for K₁ and K₂ were taken as PSF models and fitted to each of the sources to estimate their flux and position, as well as their uncertainties.

4.2.3. Total uncertainty on position angle and separation

In order to obtain the total error on the position, errors related to pre-processing were also taken into account. For the position angle, a true north error of 0.08″ was assumed (Maire et al. 2016).

The separation was first converted from pixel units to physical units by multiplying ρ by the platescale (0$\stackrel{″}{.}$0074/pix for IFS, 0$\stackrel{″}{.}$01225/pix for IRDIS). Therefore, we took into account the uncertainty on the platescale, for which we assumed a value of 0$\stackrel{″}{.}$00002/pix (Maire et al. 2016). In addition, we included a star centering error of 0$\stackrel{″}{.}$0012 and dithering error of 0$\stackrel{″}{.}$00074 (Zurlo et al. 2016). These errors were added quadratically to the errors on the measured positions of the detected companions.

4.2.4. Flux calibration of the companion spectrum

The contrast flux of the companion with respect to the central star was obtained by dividing the obtained spectrum for the companion by the spectrum for the central star from the FLUX observations. In order to convert this into physical units, the contrast flux was multiplied by a calibrated spectrum for the central star. In the ideal case, a flux calibrated, measured spectrum for the central star would be used. However, since these are not available for our targets, we calculated a model SED for the central star with the non-LTE atmosphere code FASTWIND (Puls et al. 2005; Rivero González et al. 2011). FASTWIND provides a model SED based on an effective temperature (T_eff), radius (R), surface gravity (log g), mass-loss rate (Ṁ), and terminal wind velocity (v_∞). We took the first three parameters from literature, while the mass-loss rate and terminal wind velocity were computed with the formulas proposed in Vink et al. (2001). For the effective temperature, mass, luminosity, log g, and radius, we preferably used measured parameters. If among the radius, temperature, and luminosity at least two parameters were known, the third parameter was estimated through $\frac{L}{L_{\odot }}={\left(\frac{R}{R_{\odot }}\right)}^2{\left(\frac{T_{\mathrm{eff}}}{T_{\mathrm{eff},\odot }}\right)}^4$. If among the mass M, the gravitational acceleration g, and radius R at least two parameters were know, the third parameter was estimated with $g=G\frac{M}{R^2}$, with G the gravitational constant. If the parameters were not found in literature, they were estimated from their assumed spectral type (Martins et al. 2005; Arnó 2006; Trundle et al. 2007). The adopted parameters are presented in Table 3.

ζ¹ Sco has a very large radius of 103 R_⊙ (Clark et al. 2012), which resulted in an unrealistically high mass-loss rate by adopting the Vink et al. (2001) formulas (∼10⁻⁴ M_⊙ yr⁻¹). Therefore, the mass-loss rate and terminal wind velocity were assumed from Clark et al. (2012), who determined the terminal velocity from UV spectra and presented a clumping-corrected mass-loss rate.

We verified our model flux by comparing them with the observed 2MASS K-band magnitude (Cutri et al. 2003) for each star. Therefore, we convolved the FASTWIND SED with the 2MASS K_s transmission curve to evaluate the K-band flux in the SED (Rodrigo & Solano 2020). We adopted a distance of 1.53 kpc (Sana et al. 2005) and corrected the observed magnitudes for extinction using the Fitzpatrick et al. (2019) extinction law with R_v = 3.1 and E(B − V) = 0.43 (Yalyalieva et al. 2020). The comparison between the K-band magnitude from the FASTWIND SED and the observed 2MASS K-band magnitude resulted in a scaling factor that was applied to the K₁ and K₂ fluxes obtained from the FASTWIND SED. This scaling to the observed K_s magnitude is the most important aspect in estimating the central source flux-calibrated SED. Other parameters such as effective temperature, surface gravity, and mass-loss rate play a minor role given the low spectral resolving power and the fact that near infrared spectra of OB stars are well within the Rayleigh-Jeans tail. Finally, the flux calibration requires the star to be put at an artificial reference distance from the star, for which we arbitrarily adopted 100 R_⊙, although the choice of this value does not affect in any way the results of the analysis.

4.3. Sensitivity limits

The sensitivity limits of the observations are expressed by maximum contrast curves, which give the detectable contrast as a function of the angular distance ρ to the central star. These contrast curves were measured with VIP by injecting fake companions and estimating which companion contrast corresponds to a 5-sigma detection at each angular separation. Figure 5 shows the contrasts that are reached for each target in the K₁-band. The detectable contrast magnitudes at 2″ range from ΔK₁ ∼ 10 to ΔK₁ ∼ 15 with an average value of ΔK₁ ∼ 13. This corresponds to an absolute K-band magnitude between 7 and 10 mag.

Fig. 5. Contrast curves. Top: 5-σ contrast curve, which gives the K1-contrast as a function of angular separation. Full lines represent O-type stars and dashed lines B-type stars. Middle: absolute K-band magnitude that is reached as a function of separation (5-σ). Linestyles are the same as top figure. Bottom: contrast curves translated to mass ratios. Linestyles are the same as top and middle figures.

We translated the contrast curves to mass ratio by associating to each absolute K-band magnitude the mass corresponding to an age of 6 Myr in the BT-Settl models from Allard (2014). Subsequently, we divided the masses in each curve by the mass of the central star (Fig. 2). Beyond 1$\stackrel{″}{.}$5, a mass ratio as low as 2 × 10⁻³ is reached for all targets.

5. Results

5.1. Background contamination

A reliable estimate of the expected background contamination is crucial to check whether detected companions are likely the result of a spurious association due to line of sight alignment or whether they are more likely bound to the central star. Because this is a critical aspect in estimating the true number of companions, we used various approaches.

The Besançon model of stellar population synthesis of the galaxy (Czekaj et al. 2014) simulates a catalogue of galactic stars in different directions on the sky and in various photometric bands. We performed a Besançon simulation in the K-band in a 0.01 deg² solid angle centered on HD 152 385 and computed the cumulative stellar number density per magnitude bin up to a distance of 150 kpc (Fig. 6). This shows that a large number of background stars is expected to be present. According to the Besançon model, a total of 82 946 stars are counted in a 0.01 deg² FoV up to K-band magnitude 20. This results in an average of about 92 counts per IRDIS image.

Fig. 6. Cumulative stellar densities as a function of K-band magnitude estimated from Besançon compared with observations from different catalogues: Gaia, VVV, and 2MASS (see legend).

To obtain a model independent estimate of the source number densities, we also attempted to use various optical and near infrared catalogues. We compared the prediction by the Besançon model with Gaia DR2 photometry (Gaia Collaboration 2018; converted to K-band using the colour relations from Evans et al. 2018) and two infrared surveys: the Vista Variables in the Via Lactea (VVV, Minniti et al. 2010) and the Two Micron All Sky Survey (2MASS, Skrutskie et al. 2006). The 2MASS data agree well with Besançon up to a magnitude of ∼14, while the VVV data do so between 12 and 16 mag. The lack of VVV targets at the brighter end might result from saturation effects, while the deviation at the fainter end probably reflects the sensitivity limits of the 2MASS and VVV catalogues. The Gaia observations seem to lack K-band sources. We consider the option that the photometric relationships from Gaia with infrared photometric systems are not accurate for Sco OB1 up to such faint magnitudes, since they were calculated for stars with G < 13. Further evidence for this scenario is found in our Besançon G-band simulation, which we converted to K-band using the colour relations from Evans et al. (2018). Figure 6 shows that this simulation agrees very well with the Gaia data, but they are both distinct from the other catalogues and the Besançon K-band simulation.

We prefer to use predictions from observed data whenever possible, since they provide a better representation of the local source density, as the Besançon model does not include cluster members from Sco OB1 itself. Because the VVV survey lacks bright stars (K < 12) compared to 2MASS, we decided that the background contamination is best described by 2MASS for K-band magnitudes brighter than 14 and the Besançon model for fainter sources.

5.2. Probability of spurious association

The probability that a source is a background source instead of a true companion is given by the spurious association probability P_spur. For a specific source, it is given by the probability that an object at least as bright as that source is randomly found in a circular region around the central star with radius equal to the angular separation ρ of the source (Rainot et al. 2020). It is based on the overall source density in the local region around the central star, which we estimated from 2MASS K_s-band photometry for magnitudes brighter than K_s = 14 in a box with sides 0.1 deg (corresponding to an area of 0.01 deg²) and the Besançon model for fainter sources in a 0.01 deg² solid angle, both centered on the central star. We converted the measured contrast magnitudes (in K₁) of the candidate companions to K-band magnitudes by adding the K-band magnitude of the central star (Cutri et al. 2003).

For a given source with separation ρ, we first counted the number of objects n in the 2MASS photometry in a box with sides s = 0.1 deg (K < 14) and Besançon simulation in a solid angle of 0.01 deg² (K > 14) around the central star that are at least as bright as that source. We calculated the spurious association probabilities through a Monte Carlo approach by generating 1000 times a population with n objects uniformly distributed in a circular region with area s ⋅ s = 0.01 deg². Let N_det be the number of times at least one source was found in the area πρ² around the central star, the spurious association probability is then:

$$\begin{array}{c}\hfill P_{\mathrm{spur}}=\frac{N_{\det }}{1000}·\end{array}$$

All detected sources with P_spur ≤ 20% that are in the stellar mass regime are shown in Tables 4 and 5. It is important to note that this excludes the substellar-mass IFS companions of HD 151 805 (see Fig. 3), which will be discussed in an upcoming paper.

5.3. Physical parameters

The flux calibrated companion spectrum was compared with a grid of models to obtain estimates for the mass, temperature, luminosity, radius, log g, and age of the companions. We used (pre-)main sequence evolutionary tracks from Siess et al. (2000) (< 7 M_⊙) or Brott et al. (2011) (> 7 M_⊙) associated at each time step with PHOENIX (< 3500 K, Husser et al. 2013) and ATLAS9 (> 3500 K, Castelli & Kurucz 2003) LTE atmosphere models. In addition, we also used BT-Settl evolutionary tracks and atmosphere models that overlap in mass (0.08−1.4 M_⊙) and temperature (< 4000 K) with the Siess tracks and PHOENIX and ATLAS9 models (Allard 2014). All models were compared with the flux calibrated companion spectrum through a least squares fit, where we accepted all models with a χ²-value below the 3 sigma (99.7%) threshold.

The PHOENIX and ATLAS9 temperature grid has steps of 250 K for temperatures below 13 000 K and steps of 1000 K beyond 13 000 K. The BT-Settl temperature grid has steps of 100 K. Log g runs for both grids from 0 to 5.5 in steps of 0.125. We accepted best fit ages between 4 and 8 Myr, corresponding to the age estimations of Sco OB1 (Sung et al. 2013; Damiani et al. 2016). Finally, wavelengths of each model were rebinned onto the 39 IFS wavelength channels and the two IRDIS bands.

For sources in the IFS frames, there are 41 wavelength channels in total (39 IFS + 2 IRDIS), which allowed us to extract a low-resolution spectrum and compare it with the model grid to constrain the fundamental parameters of the companion. For sources in the IRDIS frames, the two wavelength bands did not allow us to get accurate parameter constraints, since there is no information on the shape of the spectrum, but it is still possible to estimate a first guess of the mass and age by comparing the K₁ and K₂ fluxes with models within the 4−8 Myr age range. When fitting IRDIS bands, we increased the uncertainties to 10% of the flux value to avoid that our uncertainties were smaller than the resolution of the model grid.

The parameter estimates of the companion depend on the ability of the scaled FASTWIND models to reproduce the real spectrum of the central star and on the reliability of the evolutionary models. We do not take into account these uncertainties in our further analysis.

5.4. Notes on individual targets

CPD −41°7721 A is an O9.7 V:(n) type star (Maíz Apellániz et al. 2016). It is a visual binary with the B1.5 V type star CPD −41°7721 B at a separation of 5$\stackrel{″}{.}$8 (Sana et al. 2008). We also detect this companion at a separation of 5$\stackrel{″}{.}$63 with ΔK₁ = 1.43. From Gaia DR3 (Gaia Collaboration 2023, 2016), we find that this source has similar parallax and proper motion characteristics as the central star, likely indicating physical companionship (Table 6). In addition, we find one companion in the IFS FoV at $\rho =0\stackrel{″}{.}28$ for which we find a mass of 0.25 M_⊙ at an age of 4 Myr. The IFS spectrum with the best fit is shown in Fig. 7. Table 5 presents the best fit parameters. Finally, we detect 15 sources with IRDIS, among which three sources with P_spur ≤ 20% at separations of 1$\stackrel{″}{.}$09, 3$\stackrel{″}{.}$62, and 5$\stackrel{″}{.}$63 (mentioned above).

Fig. 7. IFS and IRDIS spectra for IFS companions of CPD −41°7721, HD 152 200, HD 152 042, HD 152 385, and HD 152 623 (blue). The best fit is given in orange. Dashed lines are used for wavelengths that we omitted from the fit, which are likely due to absorption in Earth’s atmosphere. The uncertainty plotted in light blue is the 3σ uncertainty.

HD 151 515 is an O7II-type star (Sota et al. 2014). It has been claimed to be a spectroscopic binary (Chini et al. 2012; Kervella et al. 2019). In addition, Turner et al. (2008) identified the presence of a visual companion at 6$\stackrel{″}{.}$2 from the Two Micron All Sky Survey (2MASS) (Cutri et al. 2003), which is just outside of our field of view. We detect 32 sources, among which one source with P_spur ≤ 20% at a separation of 1$\stackrel{″}{.}$3.

HD 151 804 is an O8Ia-type star (Sota et al. 2014) with a variable wind (Burssens et al. 2020, and references therein). Burssens et al. (2020) also report on stochastic low frequency variability. In addition, observational evidence for the existence of relativistic electrons has been found, indicating that the star is possibly a colliding wind binary (De Becker & Raucq 2013). We detect 19 sources, among which none with P_spur ≤ 20%.

HD 151 805 has spectral type B1Ib (Houk 1978). The multiplicity of this star in the range 0−150 mas has never been investigated. We detect 55 sources, among which three have P_spur ≤ 20%. Two of those are in the IFS FoV, but they have estimated masses in the substellar mass regime, so we do not discuss them in this paper. The third source is a bright companion at a separation of $5\stackrel{″}{.}9$ with an estimated mass between 1.7 and 2.2 M_⊙.

HD 152 003 is an O9.7 Iab type star. Sana et al. (2014) detected a companion through NACO/SAM observations at 38.54 mas. We detect 42 sources, but none of those have a spurious association probability below 20%.

HD 152 042 has spectral type B0.5 (Hamdy et al. 1993; Grillo et al. 1992; Garrison et al. 1977), but different luminosity classes have been proposed (III, Grillo et al. 1992; IV, Garrison et al. 1977). Our comparison with the observed 2MASS K_s-band magnitude shows that the best match is found for a spectral type B0.5IV. The multiplicity in the range 0−150 mas has never been investigated. We detect 45 sources, among which two have spurious association probabilities below 20%. One companion is in the IFS FoV at a separation of $0\stackrel{″}{.}53$. We could not find a solution for the fit within 3 sigma, so we increased the uncertainties until the smallest χ²-value was equal to the value corresponding to the 3 sigma threshold. This resulted in uncertainties of 7% of the flux value. We find a best fit at a mass of 0.25 M_⊙ (Table 5). The spectrum and best fit are shown in Fig. 7. In addition, there is a bright candidate companion at a separation of $3\stackrel{″}{.}37$ with P_spur = 6%. However, Gaia DR3 parallax and proper motion measurements do not seem to be compatible, suggesting a line of sight association (Table 6).

HD 152 147 has spectral type O9.7 Ib (Sota et al. 2014). The system is an SB1 according to Chini et al. (2012), which is confirmed by Burssens et al. (2020). Moreover, Sana et al. (2014) report on a companion at a separation of 0.77 mas with ΔH = 2.81. Though probable, it is unclear if Chini et al. (2012), Burssens et al. (2020), and Sana et al. (2014) detected the same companion, since Chini et al. (2012) and Burssens et al. (2020) do not provide information about the period of the binary system. We detect 32 sources, but none of those have a spurious association probability below 20%.

HD 152 200 is a O9.7 IV(n) type star (Sota et al. 2014). It is an eclipsing and spectroscopic binary with a debated period (4−9 days, Pozo Nuñez et al. 2019; Banyard et al. 2022). We detect 29 sources, among which one with P_spur ≤ 20% in the IFS FoV. The companion has a separation of $0\stackrel{″}{.}72$ and a mass of 0.1 M_⊙ (Table 5). The best fit is shown in Fig. 7.

HD 152 233 is the F component in a multiple system where the A component is HD 152 234 (Mason et al. 2001). It is a spectroscopic binary (Fa, Fb) with the companion resolved at a separation of 2.81 mas and with a contrast magnitude of ΔH = 1.96 (Sana et al. 2014). We find 48 sources in SPHERE, among which one has P_spur ≤ 20%. The source is a low-mass (0.08−0.20 M_⊙) candidate companion at a separation of $0\stackrel{″}{.}92$.

HD 152 248 is a known colliding-wind binary with a period of 6 days and spectral types O7 III + O7.5 III (Sana et al. 2001, 2004, 2008). Mason et al. (2009) reported on a companion at 52 mas with ΔV = 2. However, Sana et al. (2014) could not confirm this. In addition, Mason et al. (2001) detected two visual companions at $13\stackrel{″}{.}5$ and $24\stackrel{″}{.}0$ that we do not consider here. We detect 40 sources, among which one with a spurious association probability below 20% at a separation of $4\stackrel{″}{.}6$. The parallax and proper motion of this source are very similar to the central star, suggesting that it is a physical companion (Table 6).

HD 152 249 is an OC9 Iab-type star (Sota et al. 2014). Gosset et al. (2009) reported on evidence for line-profile variability in its spectrum. According to Burssens et al. (2020), the dominant type of variability is stochastic low frequency. We detect 54 sources, but none with P_spur ≤ 20%.

HD 152 385 is a B1.5V-type star (R. H. Barbá and J. Maíz Apellániz, priv. comm.). The multiplicity in the 0−150 mas separation range has never been studied. We identify 54 sources, among which three companions with P_spur ≤ 20%. One of those is in the IFS FoV at a separation of $0\stackrel{″}{.}36$. The best fit is found for a mass of 0.1 M_⊙ as shown in Fig. 7 and Table 5. The other two companions are found at separations of $1\stackrel{″}{.}4$, one of which is also detected by Gaia (Gaia Collaboration 2023, 2016). The proper motion and parallax of companion 5 are similar to those of the central object, revealing that it is likely a physical companion (Table 6).

HD 152 405 has spectral type O9.7II (Sota et al. 2014). Chini et al. (2012) states that the system is an SB1, which is confirmed by Sota et al. (2014) who identifies a period of 25.5 days. Sana et al. (2014) found a companion with a separation of 53.86 mas. We detect 13 sources, among which three with P_spur ≤ 20% at separations of $1\stackrel{″}{.}8$, $2\stackrel{″}{.}0$, and $2\stackrel{″}{.}2$.

HD 152 408 has spectral type O8 Ia (Sota et al. 2014). According to Chini et al. (2012) the system is an SB2, but this is not confirmed by Sota et al. (2014). However, observational evidence for the existence of relativistic electrons has been reported, indicating that it is possibly a colliding-wind binary (De Becker & Raucq 2013). In addition, HD 152 408 has a visual companion at $5\stackrel{″}{.}4$, which was labelled HD 152 408 B and classified as an F4 star by (Gahm et al. 1983). This source was also detected through NACO observations at a separation of $5\stackrel{″}{.}45$ with ΔK_s = 3.81, together with a source at $3\stackrel{″}{.}84$ with ΔK_s = 8.28 (Sana et al. 2014). We also detect these two companions in the IRDIS FoV, at a separation of $3\stackrel{″}{.}74$ with contrast magnitude ΔK₁ = 8.18 and at $5\stackrel{″}{.}34$ with ΔK₁ = 5.17. However, the brightest of these two sources also appears in the Gaia DR3 archive, with parallax and proper motion measurements that are quite different from the central star (Table 6). In particular, the larger parallax could indicate that it is actually a foreground star instead of a bound companion. In total, we identify 42 sources in the IRDIS FoV, including the two sources mentioned above. Among those, two additional stellar-mass sources have spurious association probabilities below 20% at separations of $1\stackrel{″}{.}68$ and $1\stackrel{″}{.}84$.

HD 152 424 has spectral type OC9.2 Ia (Sota et al. 2014). Chini et al. (2012) classifies the system as SB1. This is confirmed by Sota et al. (2014) who report a period of 133 days. In addition, Burssens et al. (2020) acknowledge the presence of a variable wind and stochastic low frequency variability. From the analysis of SPHERE data, we identify 58 sources, among which two with P_spur ≤ 20% at separations of $1\stackrel{″}{.}02$ and $1\stackrel{″}{.}28$ from the central star.

HD 152 623 is an O7 V(n)((f)) type star (Sota et al. 2014). It is a colliding wind binary (De Becker & Raucq 2013) with a period of 3.9 days (García & Mermilliod 2001). Sana et al. (2014) detected three companions at ρ = 28 mas, $0\stackrel{″}{.}25$, and $1\stackrel{″}{.}5$. We also detect the two furthest companions with ΔK₁ = 1.34 and 2.89, respectively. The companion at $0\stackrel{″}{.}25$ is detected with IFS with a separation of $0\stackrel{″}{.}25$ and we retrieve a mass of 30 M_⊙ (Table 5 and Fig. 7). In order to find a solution within 3 sigma, we had to increase the uncertainties on the flux to 14% of the flux value. The companion is also present in the Gaia DR3 archive, but both the central star and its companion have bad quality parallax and proper motion measurements that do not agree with the rest of Sco OB1 (Table 6). However, the RUWE is very large, so the astrometric measurements are not trustworthy. In total, we detect 17 sources, among which the two mentioned above and one additional companion with P_spur ≤ 20% at 4″ and with ΔK₁ = 6.04.

HD 152 756 is a B0 III-type star (Kobulnicky et al. 2019). The multiplicity in the 0−150 mas separation range has never been studied. We detect 60 sources, among which one source with P_spur ≤ 20% at a separation of $1\stackrel{″}{.}76$.

HD 154 313 is a B0 Iab-type star (Houk 1978). The multiplicity in the 0−150 mas range has never been studied. We detect 65 sources, among which two sources with P_spur ≤ 20% at separations of $1\stackrel{″}{.}66$ and $1\stackrel{″}{.}87$.

ζ¹ Sco has spectral type B1.5 Ia (Clark et al. 2012). It is a candidate Luminous Blue Variable (LBV) star (Clark et al. 2005). Mahy et al. (2022) detected a companion with Δmag = 6.3 at a separation of 11.54 mas from PIONIER observations. We identify ten sources in the SPHERE data, but none of those have P_spur ≤ 20%.

HD 152 685 has spectral type B1 Ib (Houk 1978). The multiplicity in the 0−150 mas range has never been investigated. We detect 59 sources, but none of those have P_spur ≤ 20%. It is important to note that we do not consider this star part of the Sco OB1 association.

6. Discussion

Figure 8 shows the absolute magnitude (ΔK₁ + K_s from Cutri et al. 2003, assuming a distance of 1.53 kpc Sana et al. 2005) of all sources as a function of angular separation, together with their spurious association probabilities. Assuming an age of 6 Myr for the cluster, the blue and red regions represent the brown dwarf and planetary mass regimes, based on the BT-settl models (Allard 2014). Our observations clearly probe the brown dwarf mass region, but not yet the planetary mass regime, and we detected many brown dwarf candidate companions. However, since the background contamination is significant, second epoch observations are highly necessary to confirm companionship of substellar sources by testing for common proper motion with the central star. In addition, different observations probe more or less deep into the brown dwarf regime (Fig. 5) so that detailed considerations of the limitations of each observation are needed to draw meaningful companion fractions in the brown dwarf regime. In this paper, we limit ourselves to a discussion of the possible stellar-mass companions.

Fig. 8. Absolute magnitude of all sources as a function of angular separation and their spurious association probabilities. Dots and triangles represent sources found in the neighbourhood of O- and B-type stars, respectively.

6.1. Distance

We assumed a distance of 1.53 kpc to Sco OB1 from Sana et al. (2005). The consequence of using a slightly larger distance would be that the magnitude-separation diagram is shifted upwards. However, since other distance estimates for NGC 6231 (and therefore Sco OB1) lie within 15% of the adopted distance (e.g. Kuhn et al. 2019), the effect would be small (∼0.3 mag).

6.2. SPHERE companion fraction in Sco OB1

The observed companion f_c is defined as the number of companions N_c in the sample divided by the sample size (N = 20) (Sana et al. 2014). We calculate the number of companions so that each companion counts as 1 − P_spur: $N_{\mathrm{c}}=\sum _{i=1}^N\sum _{j=1}^{N_{\mathrm{d}}}1-P_{\mathrm{spur}}^{i,j}$ where N_d is the number of detections for each target. The uncertainty is computed following Poisson statistics as in Sana et al. (2014):

$$\begin{array}{c}\hfill \sigma (f_{\mathrm{c}})=\sqrt{N_{\mathrm{c}}}/N.\end{array}$$

We do not include HD 152 685 since we do not consider it part of the Sco OB1 association.

6.2.1. O-type versus B-type

We adopt an absolute magnitude cut of K = 6.5 mag (average K-band magnitude for a 0.08 M_⊙ object, see Fig. 8) in order to include only stellar companions. Such a cut further provides the advantage that the surroundings of all stars in our sample are homogeneously surveyed independently of the magnitude of the central star (Fig. 5, middle panel). In the region between $0\stackrel{″}{.}15$ and 6″ (225−9000 AU), the observed fraction of stellar companions for the 19 O- and early B-type stars in the sample is 2.9 ± 0.4. For O- and B-type stars separately, the companion fractions are 2.3 ± 0.4 and 4.1 ± 0.8, respectively. The observed B-type star companion fraction is thus significantly larger than the companion fraction for O-type stars in our sample, at least in the 0.15−6″ separation range. It is, however, important to note that our Sco OB1 sample is quite small (only 6 B-type and 13 O-type stars), so we are dealing with low number statistics. While these fractions are not bias corrected, Fig. 5 shows that the limit of K = 6.5 is reached around all stars except for ζ¹ Sco and HD 151 804 after 1″. Within 0.5″ we are insensitive to magnitudes fainter than K = 6 for most O-type stars, corresponding to a mass of ∼0.1 M_⊙ at an age of 6 Myr. Therefore, we possibly miss objects at the stellar-substellar mass boundary in the field of view of IFS.

6.2.2. Distinction based on mass

Some of the B-type stars in our sample are actually evolved O-type stars (e.g. ζ¹ Sco, Fig. 2), so making a distinction based on spectral type can lead to misleading results. Therefore, we divide the sample in a lower mass (< 20 M_⊙) and higher mass (> 20 M_⊙) part (indicated with open and closed symbols in Fig. 2), where there are seven stars in the lower mass sample and eleven stars in the higher mass sample. The companion fractions are 3.9 ± 0.7 and 2.3 ± 0.4 for the lower and higher mass sample, respectively, so we find that lower mass stars have more companions in the separation range covered by SPHERE. However, the contrast within 1″ is significantly worse for the more massive star sample than for the lower mass sample, with all of the lower mass stars already reaching an absolute magnitude of K = 6.5 at 0.3″, while for the higher mass sample this magnitude is reached for most stars at 0.8″. Indeed, we detect only one companion around members of Sco OB1 in the higher mass sample and it is very bright, while there are five companions at the stellar-substellar mass boundary in the lower mass sample. Excluding the companions in the IFS field of view, the IRDIS companion fractions are 3.3 ± 0.7 for the lower mass sample and 2.2 ± 0.4 for the high mass sample. Thus, the same trend still applies. This result is identical to the conclusion of Frost et al. (in prep.), who found that lower mass stars have more companions than higher mass stars based on an interferometric survey of ∼50 stars.

6.2.3. Mass ratio based cut

Instead of posing a magnitude cut to limit the discussion to stellar-mass companions only, we can also make a cut based on mass ratio. We define a mass-ratio limit as the stellar-substellar boundary mass of 0.08 M_⊙ divided by the lowest mass primary in our sample (10 M_⊙) in order to only include stellar-mass companions for every target, which is equal to a mass ratio of 0.008. We find companion fractions of 1.2 ± 0.3 and 3.6 ± 0.8 for O- and B-type stars, respectively. As expected, the companion fractions are slightly lower and the difference between the higher B-star and lower O-star companion fractions has become even more prominent.

6.2.4. Supergiants versus dwarfs and (sub)giants

Finally, we compare the multiplicity properties with respect to the luminosity class of the central star. The first sample contains luminosity classes Ia, Iab, Ib, and II (11 stars) and the second sample III, IV, and V (8 stars). We find companion fractions of 0.00 ± 0.10 (IFS) and 2.7 ± 0.5 (IRDIS) for the first sample and companion fractions of 0.6 ± 0.3 (IFS) and 2.5 ± 0.6 (IRDIS) for the second sample. The IRDIS companion fractions are the same in both samples, while IFS companions are only found in the sample of dwarfs and (sub)giants. The absence of companions in the 225−1300 AU separation range for more evolved stars can indeed be tangible, but it can also result from the difference in contrasts that are probed within the inner arcsecond around a supergiant and a dwarf (Fig. 5).

6.3. Total companion and multiplicity fraction

Table 7 summarises the number of companions for each object found through different observational techniques. The table includes all SPHERE candidate companions with a probability of 80% or larger to be bound to the central star, based on the probability criterion explained in Sect. 5.2. Ten out of 13 O-type stars have at least one previously reported companion in the 0−150 mas separation range. Only one out of six B-type stars has been reported to have such a companion and this star (ζ¹ Sco) is likely an evolved O-type star (Fig. 2). However, the multiplicity of the five other B-type stars has never been studied in the 0−150 mas separation range. The total number of companions per star ranges between zero and six. Overall, the B-type stars in the sample have more companions in the separation regime covered by SPHERE, while the O-type stars have more (detected) companions in the spectroscopic regime. Assuming all spectroscopic, interferometric, and visual companions from Table 7 are truly bound, as expected from separation and contrast limits of these techniques, and adopting the previously calculated companion fractions for SPHERE companions, the total companion fractions are 3.2 ± 0.5 and 4.3 ± 0.8 for O- and B-type stars, respectively. Thus, the total companion fraction for O-type stars is still smaller than for B-type stars. Given that the multiplicity of the O-type stars in our sample has been studied more extensively than the B-type star multiplicity (e.g. Sana et al. 2006, 2014; Sota et al. 2014), the lack of spectroscopic and interferometric companions around B-type stars is likely a result of the lack of previous observations for these stars. This suggests that the average number of companions around B-type stars is truly larger than that of O-stars in our sample. Larger samples would allow us to confirm this.

The total multiplicity fraction (0 − 8″) is 0.89 ± 0.07, where the error gives the 1σ uncertainty. The latter has been computed using binomial statistics and includes the effect of the sample size (Sana et al. 2014). We conclude that almost all of the stars in our sample are in binaries or multiple systems.

We note that there may be a difference between close companions, discovered primarily through spectroscopy, and distant companions, discovered through imaging (Offner et al. 2023). Formation channels for close companions suggest that they are formed through disk or core fragmentation, followed by migration. On the other hand, distant companions are more prone to capture.

6.4. Comparison with O-type stars in Tr14

Rainot et al. (2022) investigated the multiplicity of seven O-type stars in the dense and young (< 1 Myr, Sana et al. 2010) cluster Tr14 through SPHERE observations. They found a companion fraction of 2.0 ± 0.5, which is comparable to the companion fraction we derived for O-type stars in Sco OB1. This does not necessarily mean that the age of the cluster has no effect on the multiplicity properties, since Tr14 is a much denser cluster than Sco OB1. It is expected that the density of the stellar environment also has an effect on the multiplicity (Offner et al. 2023), which might compensate for a possible age effect. Investigation of clusters with similar ages but with different stellar densities is needed to answer this question.

6.5. Comparison with O-type dwarf stars in the field and loose associations

Reggiani et al. (2022) characterised the multiplicity of 18 O-type dwarfs (class IV or V) in the field, diffuse OB associations, and cluster environments ranging from 1 to 5 Myr with SPHERE. The sample also includes three stars from Sco OB1 (CPD −41°7721, HD 152 200, and HD 152 623), which we included in our sample and re-analysed (scaling the FASTWIND flux to the observed K-band magnitude (Sect. 4.2.4) and using 2MASS+Besançon instead of Gaia to calculate the spurious association probabilities (Sect. 5.1). Assuming all IFS companions to be bound, Reggiani et al. (2022) find a companion fraction of 0.39 ± 0.15. Excluding the three stars from Sco OB1, the companion fraction is 0.27 ± 0.13. This is comparable to the IFS companion fraction for O-type stars in the Sco OB1 sample, which is 0.23 ± 0.13 (again, assuming all IFS companions to be bound).

Nevertheless, this comparison suffers from multiple biases. Firstly, the ages that Reggiani et al. (2022) probe are younger than the Sco OB1 association. Therefore, dynamical interactions might have influenced the companion characteristics in Sco OB1. Secondly, Reggiani et al. (2022) observed targets that are at distances between 1 and 4 kpc from us. Compared to Sco OB1 at 1.5 kpc, a target at 3 kpc will miss out on the inner 225−400 AU, but its IFS FoV also probes a physical separation region that is twice as large (up to 2500 AU compared to 1250 AU for Sco OB1). Finally, another difference is that Reggiani et al. (2022) observed only O-type dwarf stars, while our sample contains mostly O-type giants and supergiants. These biases complicate the comparison, as it is difficult to isolate the effect of the density of the stellar environment from other influencing factors.

6.6. Mass and separation distribution

Figures 9 and 10 show the mass distribution of candidate companions for different separation bins and the separation distribution for different mass bins. These histograms are background corrected, assuming that the 4−5″ bin is representative of background contamination and is subtracted from the other bins. The vast majority of sources in this range indeed have P_spur > 0.9 (Fig. 8). Although the errorbars are large, both figures seem to indicate an overdensity of sources in the 1−2″ separation region. This is true for all sources with masses > 0.1 M_⊙, meaning that even low-mass sources could be associated with the central star at separations between 1″ and 2″. The companion mass distribution seems to be more skewed towards low masses than a flat companion mass-ratio distribution. We note that the lower density of sources in the 5−6″ regime (Fig. 10) is due to the fact that the observations are not complete in the 5−6″ region (IRDIS FoV is ${11}^{″}\times 12\stackrel{″}{.}5$).

Fig. 9. Stellar surface density (Σc) of sources per mass bin for different separation regimes (left) and Σc corrected for background (Σbkg) (right). The background is estimated from the source density in the 4 − 5″ separation region (assuming this regime is dominated by background sources).

Fig. 10. Stellar surface density (Σc) of sources per separation bin for different masses (left) and Σc corrected for background (Σbkg) (right). The background is estimated from the source density in the 4 − 5″ separation region (assuming this regime is dominated by background sources).

7. Conclusions

Using SPHERE observations to characterise the multiplicity properties of thirteen O- and six early B-type stars in the Scorpius OB1 association, we found sources with masses as low as 0.1 M_⊙ at separations between $0\stackrel{″}{.}15$ and 6″. Despite the small sample size compared to previous multiplicity studies (e.g. Mason et al. 2001, 2009; Sana et al. 2014), a significant number of new companions have been discovered, suggesting an overall multiplicity fraction of massive stars close to 100%. Five of these companions are located closer than $0\stackrel{″}{.}85$ (in the SPHERE/IFS FoV) and four of them have estimated masses below 0.5 M_⊙, making them some of the lowest mass ratio binaries discovered so far (q ∼ 0.01).

We find a higher observed companion fraction for B-type stars in the SPHERE range than for O-type stars. Although a detailed analysis of the biases involved and larger sample sizes are necessary to make strong conclusions, the total companion fraction including previously detected companions outside of the SPHERE range is still larger for the B-type stars than for the O-type stars, suggesting that the average number of companions is indeed larger for B-type stars in our sample. In addition, the same trend applies if we divide the sample in a higher mass (> 20 M_⊙) and a lower mass (< 20 M_⊙) part. Although this is in contrast with previous multiplicity studies of massive stars, where the multiplicity has been found to increase with spectral type (Duchêne & Kraus 2013), the conclusion that lower mass stars have more companions than more massive stars is identical to that of Frost et al. (in prep.) based on an interferometric survey of ∼50 stars.

In conclusion, with SPHERE we are exploring an uncharted area in the magnitude versus separation diagram of companions. This is crucial to achieve a complete overview of the multiplicity properties of massive stars and ultimately improve our understanding of massive star formation.

Acknowledgments

The authors thank the referee for their thoughtful comments and suggestions, which certainly improved the quality of this manuscript. This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon/CRAL, and is supported by a grant from Labex OSUG@2020 (Investissements d’avenir – ANR10 LABX56). Additionally, this publication has made use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. Moreover, this work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Furthermore, we are grateful for the support received from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement number 772225: MULTIPLES). Finally, T.P. and M.R. thank the Research Foundation – Flanders (FWO) for the PhD fellowship 1164522N and the postdoctoral fellowship 1280121N, respectively.

References

All Tables

All Figures

	Fig. 1. Geometric and photogeometric distances for stars in the sample from Gaia EDR3 (Bailer-Jones et al. 2021). The dashed lines are the median geometric and photogeometric distances of the sample, while the dash-dotted line indicates where the geometric distance is equal to the photogeometric distance.
In the text

	Fig. 2. HR diagram for all targets in the sample. Evolutionary tracks and isochrones are from Brott et al. (2011) at galactic metallicity and initial rotation of ∼110 km s−1. The red line indicates the 20 M⊙ evolutionary track. Open and closed symbols are stars below and above 20 M⊙, respectively, used for the low and high mass bin in the calculation of the companion fractions (Sect. 6.2.2).
In the text

	Fig. 3. Post-processed SPHERE/IFS PCA-SDI images with one principal component for HD 152 623, 10 principal components for CPD −41°7721, HD 152 042, HD 152 200, and HD 152 385, and 40 principal components for HD 151 805. The central star is at the centre of the images (‘+’).
In the text

	Fig. 4. Post-processed SPHERE/IRDIS PCA-ADI images with one principal component. The wavelength band of the images is K1.
In the text

	Fig. 4. continued.
In the text

Fig. 5. Contrast curves. Top: 5-σ contrast curve, which gives the K1-contrast as a function of angular separation. Full lines represent O-type stars and dashed lines B-type stars. Middle: absolute K-band magnitude that is reached as a function of separation (5-σ). Linestyles are the same as top figure. Bottom: contrast curves translated to mass ratios. Linestyles are the same as top and middle figures.

In the text

	Fig. 6. Cumulative stellar densities as a function of K-band magnitude estimated from Besançon compared with observations from different catalogues: Gaia, VVV, and 2MASS (see legend).
In the text

	Fig. 7. IFS and IRDIS spectra for IFS companions of CPD −41°7721, HD 152 200, HD 152 042, HD 152 385, and HD 152 623 (blue). The best fit is given in orange. Dashed lines are used for wavelengths that we omitted from the fit, which are likely due to absorption in Earth’s atmosphere. The uncertainty plotted in light blue is the 3σ uncertainty.
In the text

	Fig. 8. Absolute magnitude of all sources as a function of angular separation and their spurious association probabilities. Dots and triangles represent sources found in the neighbourhood of O- and B-type stars, respectively.
In the text

	Fig. 9. Stellar surface density (Σc) of sources per mass bin for different separation regimes (left) and Σc corrected for background (Σbkg) (right). The background is estimated from the source density in the 4 − 5″ separation region (assuming this regime is dominated by background sources).
In the text

	Fig. 10. Stellar surface density (Σc) of sources per separation bin for different masses (left) and Σc corrected for background (Σbkg) (right). The background is estimated from the source density in the 4 − 5″ separation region (assuming this regime is dominated by background sources).
In the text