Free Access
Issue
A&A
Volume 660, April 2022
Article Number A122
Number of page(s) 14
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202142418
Published online 25 April 2022

© ESO 2022

1. Introduction

The formation of massive stars remains an open question in astronomy today. For low-mass stars, the core accretion paradigm – starting with the collapse of pre-stellar cores into a single or binary protostar with the subsequent accretion of matter through a Keplerian disk (e.g., Shu et al. 1987Inutsuka 2012) – is expected to produce stars over timescales of 10−50 Myr. The problems in scaling up the low-mass star formation models to the high-mass regime indeed come from the short timescales involved (e.g., Zinnecker & Yorke 2007; Brott et al. 2011). Therefore, for massive stars, other mechanisms such as competitive accretion (Bonnell et al. 2001) and protostellar mergers (Bonnell et al. 1998; Bally & Zinnecker 2005; Moeckel & Clarke 2011) have been proposed. Except for collisional models, most high-mass star formation theories are in agreement with regard to the presence of dense and massive accretion disks in sustaining accretion in the presence of radiation pressure. These disks are presumably unstable to fragmentation and this can lead to the formation of stellar companions (Kratter & Murray-Clay 2011).

A clear observational evidence is that the multiplicity frequency is significantly higher among massive stars (Duchêne & Kraus 2013). Several spectroscopic surveys of OB stars both in our Galaxy, and in the LMC (see e.g., Chini et al. 2012; Kobulnicky et al. 2012; Sana et al. 2012, 2013; Sota et al. 2014) have shown that the binary (or multiple) frequency may be > 70% for binaries with physical orbital separation smaller than 1 AU. The occurrence of longer period binaries has been explored through speckle interferometry by Mason et al. (1998, 2009), adaptive optics (AO) by Turner et al. (2008) and Close et al. (2012), and lucky imaging by Maíz Apellániz (2010) and Peter et al. (2012). These studies also demonstrate the high incidence of binaries and multiples among longer-period systems. In order to fill the observational gap between classical imaging and spectroscopic surveys, the Southern MAssive Stars at High angular resolution survey (SMaSH+, Sana et al. 2014) has combined optical interferometry (VLTI/PIONIER) and aperture masking (NACO/SAM). The SMaSH+ survey is sensitive to mostly bright companions (ΔH < 4) between and for a large sample of O-type stars. The occurrence of fainter (ΔH < 8) companions at larger separations (up to 8″) was also probed in entire NACO field-of-view (FoV). The main conclusion of this work was that nearly all massive stars have at least one companion in the separation range covered by the observations and that over 60% are part of a higher order multiple.

It is thus critical to probe the frequency of even fainter (ΔH < 4) companions at these separations to determine the shape of the period and mass ratio distribution and to estimate the total binary frequency. The properties of the binary population may indeed serve as a useful diagnostics tool to discriminate between different formation models, particularly at separations that approximately correspond to the size of the accretion disk, where we expect to find the low-mass companions formed from the remnants of the fragmented disk.

In this regard, the first paper of the Carina High-contrast Imaging Project of massive Stars (CHIPS, Rainot et al. 2020) represents a proof of concept that the extreme AO, implemented at the VLT through the Spectro-Polarimetric High-contrast Exoplanet REsearch instrument (SPHERE, Beuzit et al. 2019), provides the necessary spatial resolution and dynamics to look for the faintest companions to nearby massive stars.

In this paper, we analyze the multiplicity properties of a small sample of 18 dwarf O-type stars from the Galactic field, loose associations, or denser clusters to have a first statistics on the occurrence of faint (ΔK < 12) companions in the to 6″ angular separation regime.

First, we describe the sample we considered (Sect. 2). The setup for the observations and the data reduction are presented in Sect. 3. The image post processing is described in Sect. 4. Then, in Sect. 5 we show our results and we discuss them in Sect. 6. Finally, we offer our conclusions in Sect. 7.

2. Sample

Our sample consists of 18 dwarf O stars with spectral types ranging from O3 V to O9.7 V. Dwarf O stars only represented 20% of the SMaSH+ sample due to the magnitude-limited (H < 7.5) quality of the survey. These represent, however, the most natural targets in the search for observational constraints on the outcome of the massive star formation processes because they are less evolved than giants and supergiants and because they are intrinsically less luminous, allowing us to probe the low-mass end of the companion mass function. Thus, in the framework of this small-scale project, we exclusively targeted a sample of dwarf O stars covering a range of environments (clusters, diffuse OB associations, field) and masses (from 15 to 60 M). The properties of the objects are described in Table 1. The ages of the targets are between and 1 to 5 Myr and the distances from ∼1 to 4 kpc.

Table 1.

O-type stars used in this work.

3. Observations and data reduction

The 18 O type stars were observed as part of a SPHERE Science Verification program and a standard service mode program in 2015. The high-contrast imaging instrument SPHERE is mounted at the Naysmith platform of Unit 3 telescope (UT3) at ESO’s VLT, and consists of an extreme adaptive optics system, coronagraphic masks, and three different sub-systems. The observations were carried out in the IRDIS and IFS extended mode (IRDIFS_EXT) by simultaneously using the Integral Field Spectrograph (IFS) and the Infra-Red Dual-beam Imaging and Spectroscopy (IRDIS) sub-systems (Galicher et al. 2018).

The IFS images have a pixel scale of 7.4 mas and with a total size of 290 × 290 pixels cover a FoV on the sky. IRDIS instead has a FoV of 12″ × 12″ with a pixel scale of 12.25 mas (i.e., 1024 × 1024 pixels in total). The IRDIFS_EXT mode allows us to combine the YJH band observations with IFS to dual K-band observations with IRDIS. Due to the size of the FoV and its spectroscopic capabilities, IFS enables the detection and characterizion of companions at close separations, whereas the larger FoV of IRDIS provides statistics for companions at larger separations and for the local field density of objects.

All observations were carried out in pupil-tracking mode to allow for image post-processing through angular differential imaging (ADI, Marois et al. 2006) techniques. However, most of the objects were not observed during meridian passage and only a limited parallactic angle variation was achieved.

For both IRDIS and IFS, the observing sequence was composed of three types of observations. Our science frames (OBJECT, O) were obtained by blocking the light coming from the bright central stars with SPHERE’s apodized Lyot coronagraphs. We also obtained CENTER (C) frames, which were acquired by applying a sinusoidal pattern to the deformable mirror to infer the position of the star behind the coronagraph. Finally, for spectro-photometric calibration, we took uncoronagraphic FLUX (F) images of the stellar point spread function (PSF) by offsetting the central star from the coronagraphic mask and we used a neutral density filter (ND2.0) to avoid any saturation of the detector. The same F-C-O sequence was repeated three times for each target. For HD 123 056 only IFS observations were obtained.

The choice of detector integration times (DITs) and number of DITs (NDIT) for our OBJECT and FLUX exposures for each target are presented in Tables 2 and 3 for IFS and IRDIS, respectively. In total, for every star, we obtained four-dimensional (4D) IFS and IRDIS data cubes. The IFS cubes are composed of 290 × 290 pixel images for each of the 39 wavelengths channels (from 0.9 to 1.6 μm) and each sky rotation. The IRDIS data cubes contain 1024 × 1024 pixel images for each one of the two wavelengths channels (K1 and K2) and each sky rotation. The summary of the observing conditions and total parallactic angle variation for each object and is also presented in Tables 2 and 3.

Table 2.

Observing setup and atmospheric conditions for FLUX (F) and OBJECT (O) IFS observations.

Table 3.

Observing setup and atmospheric conditions for FLUX (F) and OBJECT (O) IRDIS observations.

The data reduction of IRDIS and IFS images was carried out by the SPHERE Data centre (Delorme et al. 2017, DC) at the Institut de Planetologie et d’Astrophysique de Grenoble (IPAG)1. The SPHERE-DC applies a standard data reduction to the science and PSF frames by removing bad pixels, dark and flat frames and estimating the bias in each exposure. They calibrated the astrometry using the on-sky calibrations from Maire et al. (2016), with a true north correction value of 1.75 ± 0.08° and a plate scale of 7.46 ± 0.02 mas pix−1 and 12.255 ± 0.009 mas pix−1 for IFS and IRDIS, respectively.

As three uncoronagraphic PSF observations were taken during our observing sequence at each of the wavelength channels (2 for IRDIS and 39 for IFS), we computed the median of the three PSF frames for each wavelength to increase the signal-to-noise ratio (S/N). We measured the total flux of the central object on the median-combined images at each wavelength and its uncertainty by computing the standard deviation of the flux measured in the three PSF frames.

4. PCA image processing

To post-process the reduced data cubes we used the python open-source Vortex Imaging Processing package2 (VIP, Gomez Gonzalez et al. 2016), which was developed to analyze high-contrast imaging datasets for exoplanet detection. Based on the type of data and user choice, it also performs angular, reference, and spectral differential imaging (ADI, RDI, SDI, respectively), or simultaneous ADI+SDI, all based on a principal component analysis (PCA, Amara & Quanz 2012; Soummer et al. 2012) approach. In our study, for each object, we applied PCA/ADI separately to the two K1 and K2 IRDIS cubes and PCA/SDI on all IFS channels simultaneously to get the final post processed images. The resulting reduced PCA/ADI K1-band IRDIS frames for all our targets are shown in Fig. 1. Only the final PCA/SDI IFS images that present possible companion detections are presented in Fig. 1.

thumbnail Fig. 1.

Final PCA/SDI IFS images for the targets with a detected source. The position of the central star is indicated by a white cross.

5. Results

The visual inspection of the final IFS and IRDIS PCA images displayed in Figs. 13 reveals the presence of six possible companions in the IFS frames and many other point-like sources in the IRDIS FOVs. To evaluate which ones are true detections, we estimated their S/N values using the appropriate function implemented in VIP. This module computes the S/N at every pixel of an image by measuring the signal in a 1 FWHM diameter aperture and comparing it to the standard deviation of the other resolution elements in an annulus at the same radial distance from the center of the frame. It uses the approach described in Mawet et al. (2014) on small sample statistics, which applies a student t-test to determine the S/N and contrast in high contrast imaging observations. In our study, we adopted a S/N threshold of 5 to assess if a source on the image is a true detection.

thumbnail Fig. 2.

Final post-processed PCA/ADI IRDIS images for the first 9 targets.

thumbnail Fig. 3.

Final post-processed PCA/ADI IRDIS images for the last 8 targets.

5.1. Source characterization

Once we identified the true physical objects on our images, we extracted position and flux for each source. As described by Rainot et al. (2020), we adopted two different techniques based on the radial separations of the sources. For all IFS sources and all IRDIS source within 2″, we used a negative fake companion approach, implemented in VIP (see Sect. 5.1.1). For sources with separation beyond 2″, we adopted a PSF fitting routine (Bodensteiner et al. 2020), as these objects do not suffer from the central star PSF influence in the derotated and collapsed cube.

5.1.1. Sources within 2″

For all IFS detections and for all sources with angular separations smaller than 2″ from the central star in the final IRDIS images, we first measured the flux with aperture photometry at each wavelength. We used this initial guess as a starting point of a Simplex Nelder-Mead optimisation implemented in VIP. VIP estimates position and flux for each source by applying a NEGative Fake Companion technique (NEGFC), which consists in inserting negative artificial companions in each individual frame before running PCA, varying at the same time their brightness and location. The residuals are then calculated in the final PCA images and compared to the background noise of all resolution elements in an annulus at the same radial distance. The combination of brightness, separation, and position angle that minimizes the residuals is estimated through a Nelder-Mead minimization algorithm. The artificial companions are obtained from the uncoronagraphic images of the central star. As the Nelder-Mead optimization does not return the uncertainties on the estimates of the parameters, we implemented a set of Monte-Carlo simulations to compute the accuracy of our algorithm. For each wavelength, we inserted 25 artificial sources at the same radial distance and with the same flux as a given detection, but varying their position angles. We once again measured the flux and location of the injected fake sources using the negative fake companion algorithm and compared them to the true values. The standard deviation of the measurements of each parameter gives us an estimate of the 1σ error.

5.1.2. Sources beyond 2″

Beyond 2″, the contribution of the central star is negligible and it is the background noise that is dominant (see Sect. 5.6). Therefore, for point-like objects beyond this separation, the use of ADI and SDI techniques is not necessarily needed to derive precise astrometry and photometry. Following Rainot et al. (2020), we adopted a standard PSF-fitting technique as described in Bodensteiner et al. (2020), which is based on the photutils3 python package along with an effective PSF model developed by Anderson & King (2000).

We implemented this method on the derotated and collapsed images for both K1 and K2 and we adopted the IRDIS uncoronagraphic images in each band as accurate PSF models for the fit. The PSF is then fitted to each source individually to estimate the best positions and flux values in K1 and K2, together with their uncertainties. This technique is particularly useful for sources that are too close to the edges of the frames for the NEGFC technique to work.

5.1.3. Final error budget

The methods described in Sects. 5.1.1 and 5.1.2 only take into account the uncertainties related to the algorithm used. For the total errors on the photometry we also accounted for the flux variations of PSFs of the central stars described in Sect. 3. To estimate the total uncertainties on the separation and position angle of each source, we adopted the plate scale and astrometric calibration precision given by Maire et al. (2016) and the ESO SPHERE user manual. The final errors are obtained by a quadratic sum of the algorithm measurement errors (either from the Monte-Carlo simulations or PSF fitting), the star’s center position uncertainty (1.2 mas, from Zurlo et al. 2016), the plate scale precision of 0.02 mas pix−1 for IFS and of 0.021 mas pix−1 for IRDIS, the dithering procedure accuracy (0.74 mas, Zurlo et al. 2016), and the true north uncertainty (±0.08 deg). The summary of the properties and corresponding uncertainties of all the sources found around each star are available at the CDS.

5.2. Spurious association probabilities

To identify which of the detected sources could be bound companions, we estimated on a statistical base the probability of chance alignment association. Following Rainot et al. (2020), we defined the probability of spurious association (Pspur(ρi|Σ(Ki))) as the probability that at least one object is found by chance at an angular separation from the central star ρ smaller or equal than that of the ith companion (i.e., ρ ≤ ρi), given the local source density Σ of stars at least as bright as i (K ≤ Ki). To compute Pspur, we first evaluated the local field density Σ(Ki) = Nobj(K ≤ Ki)/(πr2) of objects at least as bright as the companion i in a πr2 surface. To do so, we used the Gaia (Gaia Collaboration 2018) DR2 catalog to evaluate the number and brightness of companions in a r = 2′ radius region around each target. To convert the Gaia magnitudes into K-band magnitudes, we used the color relations given by Evans et al. (2018). For each i source, we then used a Monte Carlo approach to generate 10 000 populations of Nobj(K ≤ Ki) stars uniformly distributed in πr2. The probability of spurious association is thus obtained as the fraction of Monte Carlo simulations in which at least one star is found at the separation ρ ≤ ρi.

All properties and probabilities for sources with Pspur < 5% are given in Table 5. For the objects with low spurious association probability (e.g., Pspur < 0.05), which are most likely bound companions, it is essential to obtain a confirmation of common proper motions and a characterization of orbital motion in the future to definitively demonstrate a true physical association.

5.3. Absolute flux values

Knowing the flux calibrated spectrum of the central star in the wavelength range covered by our SPHERE observations (Y to K) is required to compute the absolute fluxes of possible bound companions in the images. In fact, under the assumption that the same extinction affects both the primary and its companions, the unreddened primary spectral energy distribution allows us to derive the companions spectral energy distribution and thus characterize them through a comparison with atmosphere models (see Sect. 5.4). As such a spectrum is not easily available, we modelled the spectral energy distribution of the targets with likely bound sources (those with Pspur < 5%) with the non-local thermodynamic equilibrium (non-LTE) atmosphere code FASTWIND (Puls et al. 2005). When the central object is composed by a spectroscopic multiple system, each component is modelled separately at first, and then combined to obtain a unique spectrum. This is the case for HD 150 135, BD −13° 4929, and V3903 Sgr. The assumed input parameters for FASTWIND for each component of the central object are presented in Table 4. Since we only characterized sources with Pspur < 5%, Table 4 only includes the stars hosting likely bound companions. The parameters for the computation were based on the spectral type characterization found in the literature (see Table 2) and the observational O-star calibration tables from Martins et al. (2005). Parameters for the B star components were instead found in Trundle et al. (2007). As inputs for FASTWIND, we calculated the mass-loss rate () and terminal wind velocities (v) for each star following Vink et al. (2001). We remark that in this process and later on in Sect. 5.5, it is necessary to adopt a reference radius for the sphere at the surface of which the flux is computed. Without loss of generality, we arbitrarily adopted a value of 100 R, although this value has no physical meaning, or impact in our calculation.

Table 4.

Assumed stellar parameters for the FASTWIND calculation of the calibrated spectra.

5.4. Spectral fitting

Similarly to what was done by Rainot et al. (2020) for QZ Car, we used the low-resolution IFS spectrum of all IFS detections to constrain their stellar parameters. To do so, we used both the ATLAS9 LTE atmosphere models (Castelli & Kurucz 2003), covering the 3500−50 000 K temperature range, and the LTE PHOENIX models (2300−12 000 K, Husser et al. 2013) for T < 3500 K. To each age value in the pre-main sequence (PMS) evolutionary tracks of Siess et al. (2000) below 7 M or of Brott et al. (2011) above, we associated an atmospheric model and we quantitatively compared it to the flux calibrated SED of each detection. For the comparison, we rebinned each model to the 39 IFS wavelength channels and the 2 IRDIS bands and we estimated the corresponding χ2 by taking into account the uncertainty on the flux calibrated spectrum.

Several combinations of stellar parameters are consistent with the observations. In the next section (Sect. 5.5), we summarize (object by object) the results of the SED fitting procedure and the best fit parameters corresponding to the 95% confidence interval.

Our IRDIS observations provide us with only two independent wavelength channels, K1 and K2, at 2.110 and 2.251 μm, respectively. This does not allow us to constrain the shape of the SED, however, it does offer the possibility to assess the object absolute K-band magnitude under the assumption that the IRDIS sources are located at the same distance as the central star. Therefore, for all IRDIS detections with Pspur < 5%, which are most likely to be bound companions, we compared the K1 and K2 absolute fluxes to the ATLAS9 and PHOENIX models. For many of them, we were able to find solutions that are in agreement with the central star age (as reported in Table 1). For those stars, the ranges of masses and ages that are consistent with the IRDIS observations are given in Table 5. For several objects, we could not find a good fit within the given age range, possibly indicating that they do not have physical companions in common. Finally, given our grid of models, we are also limited to PMS stars more massive than 0.1 M. Some of the very faint sources for which we could not find a proper fit could be shown to be objects at the stellar-substellar boundary that are not covered by our models.

Table 5.

Angular separations (ρ), position angles (PA), magnitude contrasts in the K1 and K2 bands (ΔK12), and spurious alignment probabilities (Pspur) for the sources with Pspur < 5% around our targets.

5.5. Summary

– HD 64 568 has been associated with an O3 V((f*))z SpT (Sota et al. 2014). It has recently been classified as runaway from the southern component of NGC 2467 by Maíz Apellániz et al. (2020, 2022). We did not find close companions in the IFS field nor in the IRDIS larger FOV. This supports the idea that all O stars are born in multiple systems.

– HD 93 128 (Trumpler 14 2) is part of Trumpler 14 and has been classified as O3.5 V((fc))z SpT (Sota et al. 2014). We do not detected any IFS companion. Two IRDIS sources are however found with Pspur ≤ 5% at separations of and and for which the fits of the K-band magnitudes is consistent with masses of 0.25−0.4 and 0.2−0.4 M, respectively, with ages between 0.3 and 1.5 Myr. A more detailed analysis of the system, including the reanalysis of the B visual companion (classified as B0.2 V by Maíz Apellániz et al. 2022) will be presented in Rainot et al. (2022).

– HD 155 913 is a classified as a possible runaway given that the proper motion points away from the young stellar cluster NGC 6822, which is located half a degree away from the star (Maíz Apellániz et al. 2018). HD 155 913 is a O4.5 Vn((f)) SpT according to the GOSS catalog (Sota et al. 2014). However, it is reported as SB2 in the OWN data (Barbá et al. 2010, 2017), so the width observed in the GOSS spectrograph could be due to unresolved orbital motion. A visual companion is reported by Aldoretta et al. (2015). No relevant sources are found in our IFS and IRDIS images.

– HDE 319 699 is an O5 V((fc)) type star (Sota et al. 2014) and it is part of NGC 6334 (Russeil et al. 2017). According to our IRDIS images, the star is part of a crowded field, but no close IFS companions or bright IRDIS sources are present in the data. Six sources within roughly 2″ are found with Pspur ≤ 5%.

– HD 124 314 is a O6 IV(n)((f)) Aa-Ab binary with a separation of 1.5 mas, separated by the O9.2 IV(n) visual components Ba-Bb by (Sota et al. 2014). A third visual component C was found by SMaSH+ at . No additional companions are seen in our IFS observations, but all visual companions are also detected in our IRDIS images with Pspur ≲ 1%. According to our fit, all visual components could be coeval with the central star system.

– HD 150 135: it has been recently classified as SB2 by Maíz Apellániz et al. (2020), composed by a O6.5 V((f))z primary and an O8 secondary component. SMaSH+ found it to be a Aa-Ab+B multiple system. The A-B components are separated by , whereas Aa-Ab by 0.95 mas. Although we do not detect any companion in the IFS field, the elongated IFS PSF suggests the close inner binary could be currently at a larger angular separation (with a tentative orientation of ∼45°). Our IRDIS images show the presence of the B components, which appears to be a possible double system with source 6 (Pspur ∼ 7%).

– HDE 326 775 (O6.5 V(n)((f))z, Maíz Apellániz et al. 2016) is part of RCW 113 HII region (Arias et al. 2016), situated in the southern part of the Sco OB1 association. We did not detect companions in the IFS data. Among the four sources with Pspur ≤ 5% in the IRDIS field, only two of them show K-band magnitudes that could be fitted with models in the age range of the central star.

– V3903 Sgr is a O7 V(n)z + B0: V binary system and it is part of the Sgr OB1 association (Sota et al. 2014). We found a most likely (spurious association probability of 1%) visual companion at a separation of and with a ΔK = 8.0 mag. Unfortunately, we cannot obtain an acceptable fit of the K-band magnitudes with an age that is consistent with the central star.

– HD 152 623 is a known colliding-wind binary (De Becker & Raucq 2013), and a possible runaway (Maíz Apellániz et al. 2018, and references therein). The components B and C were detected by SMaSH+ at a separation of and , respectively. HD 152 623 B is detected in both IFS and IRDIS images and our best fit for the spectrum is obtained with a mass of 6 M at an age of 0.3 Myr (see Fig. 4). HD 152 623 C is detected with a ΔK = 6.13, and it is consistent with being a 3.5 M with an age of 1.5 Myr. Source 9 has also a probability of Pspur = 4%, suggesting it could be a bound companion to the system as well. It appears to be roughly coeval with other components and it can be fitted best with a 0.9 M star in the age range of 1.8−4.5 Myr

thumbnail Fig. 4.

Observed spectrum for HD 152 623 B.

– BD −13° 4929 it is part of NGC 6611 and it is found to be a SB3 by (Sana et al. 2009). The triple system is composed by an O8 V and two B0.5: V+B0.5: V SpT stars (Maíz Apellániz et al. 2022). We detected two companions in IFS, which are also confirmed by the IRDIS images. According to isochrones fitting, they are both estimated to be 0.2 M with a best age value of 1.8 and 1.9 Myr, respectively, in good agreement with recent age determinations of NGC 6611 (Bonatto et al. 2006; Martayan et al. 2008). The spectra of the companions are presented in Figs. 5 and 6. In IRDIS we also detect a fainter source at , which is not visible in IFS. The K-band fluxes are consistent with a 0.1 M star with an age in the range of 0.8−5 Myr.

thumbnail Fig. 5.

Observed spectrum for BD −13° 4929 C.

thumbnail Fig. 6.

Observed spectrum for BD −13° 4929 D.

– HD 101 191 (O8 V SpT, Sota et al. 2014) is a long-period SB1 system according to Sana et al. (2011) and Chini et al. (2012). No likely companion appears in IFS nor in IRDIS.

– HDE 323 016 is classified as a O8.5 V SpT and it is part of the S5 HII region according to Neckel & Chini (1981). Similarly to HDE 319 699, the IRDIS FoV reveals a densely populated region around the star, but none of the sources appears to have a Pspur ≤ 5%.

– HD 149 452 (O9 IVn) is a relatively isolated O type star (Sota et al. 2014). Sana et al. (2014) reported the detection of a source at from HD 149 452. Given the non-detection in their H band image, suggesting a strongly reddened object, they indicated the possibility of it being a background object. In our IRDIS observations, we found a ΔK = 3.9 mag source at , with a Pspur = 5%. Our K-band colors are consistent with 3−28 M object with an age of 2−5 Myr, depending on the evolutionary models used.

– HD 76 341 is classified as O9.2 IV by Sota et al. (2014). One companion has been detected by SMASH+ (Sana et al. 2014) at ρ = 169 mas and with a contrast of 3.7 mag in the H band. Sota et al. (2014) also observed variability in the spectrum of HD 76 341, indicating a possible hierarchical triple system. We re-detected the companion at 159 mas in both IFS and IRDIS observations and we characterized it as a 3.5 M object with an age of about 2.4 Myr. The observed spectrum of the companion is shown in Fig. 7.

thumbnail Fig. 7.

Observed spectrum for HD 76 341 B.

– BD −13° 4928 is a O9.5 V star according to Sana et al. (2009) and a fast rotator (J. Maíz Apellániz, priv. comm.). We did not find close companions in the IFS field. A source at a separation of is present in the IRDIS image. It has a ΔK magnitude of 7.9 and a probability of spurious association of < 2%. The best fit mass of this object is 0.1 M at 2.4 Myr.

– CPD −41° 7721 is part of NGC 6321, in the core of the Sco OB1 association, at about 1.5 kpc (Maíz Apellániz et al. 2022). According to Maíz Apellániz et al. (2016) it is a O9.7 V:(n) star. Together with the B1.5 V star CPD −41° 7721B, it forms a visual double star with the two components separated by . Besides confirming the already known visual companion with the IRDIS observations, we detected a new fainter and closer companion in IFS. We classified it as a 0.16 M star, with an age of 3.2 Myr. The spectrum of CPD −41° 7721 C is presented in Fig. 8. Another IRDIS source with a ΔK = 6.13 is detected at a separation of 1″ with a Pspur ∼ 1%, making it a likely bound companion. The best fit for this object is obtained with a best fit mass of 0.7 M and an age of 3.8 Myr.

thumbnail Fig. 8.

Observed spectrum for CPD −41° 7721 C.

– HD 123 056 is a field star and it is found to be a hierarchical triple system by Mayer et al. (2017). The PSF of IFS looks elongated, also indicating the multiple nature of the central object. No further IFS companions are found in the data, and unfortunately no IRDIS data were collected for this object.

– HD 152 200 is an O9.7 IV(n) star in NGC 6231 (Sota et al. 2014). It has been reported as an eclipsing binary with period of about 9 days by Pozo Nuñez et al. (2019). Close to the central star, at a separation of we found a 0.1 M star that is 5.2 Myr old. The spectrum of the companion is shown in Fig. 9.

thumbnail Fig. 9.

Observed spectrum for HD 152 200 B.

5.6. Sensitivity limits

We estimated the sensitivity we reached with our observations in terms of magnitude difference as a function of the angular separation to the central star. We used the VIP contrast curve function that calculates the contrast limits for a chosen σ level by injecting artificial companions (with scaled flux based on the unsaturated PSF of the central object) and calculated the noise and the algorithm throughput at different radial distances from the center. At close separations, we took into account the small sample statistics correction proposed in Mawet et al. (2014). The 5-σ sensitivity curves that we obtained for each target and for both IFS and IRDIS observations are presented in Figs. 10 and 11, respectively.

thumbnail Fig. 10.

Sensitivity of our IFS observations expressed as magnitude difference vs. angular separation to the central star. The curves for each objects correspond to the 5-σ contrast.

thumbnail Fig. 11.

5-σ contrast curves for our observations in K1 (left) and K2 (right) IRDIS bands.

6. Discussion

All the sources discovered in this work are presented in the Δmag vs. separation plot shown in Fig. 12. Despite the modest number of stars observed in this study compared to the SMASH+ survey, the IFS discoveries are clearly populating a region of the parameter space that has never been reached before.

thumbnail Fig. 12.

Δmag vs. separation diagram. The location of the newly discovered sources (black and grey dots), as well as the median IRDIS and IFS contrast limits are compared with the outcome of previous surveys (Sana et al. 2014). Black, grey, and light-grey dots corresponds to sources with Pspur ≤ 5%, 5%< Pspur ≤ 20%, and Pspur > 20%, respectively.

Among them, four are newly detected sources and have estimated masses below 0.2 M. Ranging from mass ratios of q = 0.002 (0.004) to 0.01, these estimates makes them the lowest mass-ratio companions ever discovered around O-type stars. These objects are located between 400 and 1500 AU, based on the distance of the central stars. Such separations are in agreement with recent observations of fragments and substructures in Keplerian disks around (proto-)O stars (e.g., Beuther et al. 2017; Ilee et al. 2018; Maud et al. 2019). Whether these substructures and fragments will survive, end up as companions at such large separations, or migrate inwards to form spectroscopic binary systems is still an open question that current hydrodynamics simulations are working to address (Oliva & Kuiper 2020).

Following the definitions given in Sana et al. (2014), we can calculate the fraction of companions as the mean number of companions per central star, that is, the ratio of the total number of companions to the sample size. The error on the fraction can be estimated with Poisson statistics (see Eq. (9) in Sana et al. 2014). Assuming that all sources detected with angular separations of less than are physically bound companions, the observed (uncorrected for bias) fraction of companions for O-type stars between 150 and 900 mas (based on the effective size of the IFS fov) is 0.39 ± 0.15. If we take into account the spurious association probability for sources with Pspur ≤ 5% in the larger IRDIS field of view (FOV), this fraction increases to 1.6 ± 0.3 in the separation range from to 6″. In order to compare our results with those of the SMASH+ survey, we need to restrict the comparison over the delta-magnitude and separation ranges covered by both studies. For angular separations between and , and contrasts < 4 mag, we observed a fraction of companions of 0.12 ± 0.07, whereas in the angular separation range of and Δmag < 8, we obtained 0.76 ± 0.21. Due to a brightness limitation, the SMASH+ survey only observed a small fraction of O-type dwarfs. For this subset of objects (50 in total), over the same ranges, they observed a companion fraction of 0.24 ± 0.07 between and and of 0.94 ± 0.14 between . Both fractions are consistent with our findings within the errors. The correction for observational and selection biases goes beyond the scope of the present work and it will be presented in an upcoming paper.

If we consider as bound objects all sources with Pspur ≤ 5%, the mean number of companions combining spectroscopic and eclipsing as well as visual multiples is 2.3. This number also includes known runaway stars. This implies not only that most massive stars are in multiple systems but also that triple or higher-order systems are more common than simple binaries. This outcome is in agreement with the results from the MONOS (Multiplicity Of Northern O-type Spectroscopic systems) project (Maíz Apellániz et al. 2019).

Finally, concerning the influence of the environment density on the companion fraction, we do not see any strong correlation on the total number of companions in the separation range, when comparing stars in denser cluster (e.g., Trumpler 14), OB associations, or rather isolated objects (see Table 6). Nevertheless, we note that the previously reported runaway stars HD 64 568 and HD 155 913 are confirmed to be single, according to our spurious association probability criterion, as well.

Table 6.

Summary of detected companions (Pspur ≤ 5%) in IFS and IRDIS for all objects.

7. Conclusions

In this work we used VLT/SPHERE in IRDIFS_EXT mode to simultaneously carry out observations with the IFS and IRDIS subsystems and characterize the multiplicity properties of a sample of 18 O-type stars from stellar clusters and loose associations between and 6″. We summarize the main results of our study below.

  1. Despite the small size of the sample, compared to previous high angular resolution surveys (e.g., SMaSH+, Sana et al. 2014), we added a considerable number of companions in the angular separation range. By reaching ΔH = 12, we also opened up a completely new region of the parameter space, with the possibility of exploring the existence of sub-solar companions around massive O-type stars.

  2. We found and characterized seven (five of which are previously unknown) companions within from the central star. The five newly discovered companions have estimated masses below 0.25 M, making them the highest mass-ratio binaries or multiple systems known thus far.

  3. In addition to the close stellar companions, we detected several other sources in the larger IRDIS FoV with Pspur < 5%. Although we expect many of them to be physically bound companions, only future proper motion observations will enable us to confirm their companionship.

  4. If we assume that all sources with angular separations below are physically bound companions, and by taking into account the spurious association probability for those with Pspur ≤ 5% from to 6″, the observed (uncorrected for bias) fraction of companions for O-type stars is 0.39 ± 0.15 from to and 1.6 ± 0.3 in the separation range from to 6″. These fractions clearly support the idea that massive stars form almost exclusively in multiple systems, with preference for triples or higher-order multiples.

  5. The results of this study demonstrate that probing extreme contrasts as allowed by large AO-assisted coronagraphic surveys is fundamental to fully constrain the multiplicity properties of massive star companions in regions of the parameter space that remained unexplored so far and to characterize the low-mass end of the mass function.


Acknowledgments

The authors thank the referee, Jesús Maíz Apellániz, for their detailed and careful comments that helped to improve this manuscript. This work is based on observations collected at the European Southern Observatory under programs ID 60.A-9369(A) and 095.D-0495(A). We thank the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice) and Observatoire de Paris/LESIA (Paris) and supported by a grant from Labex OSUG at 2020 (Investissements d’avenir a ANR10 LABX56). We especially thank P. Delorme, E. Lagadec and J. Milli (SPHERE Data Centre) for their help during the data reduction process. We acknowledge support from the FWO 1280121N grant and the FWO-Odysseus program under project G0F8H6N. This project has further received funding from the European Research Council under European Union’s Horizon 2020 research program (grant agreement No. 772225, MULTIPLES) and under the European Union’s Seventh Framework Program (ERC Grant Agreement n. 337569).

References

  1. Aldoretta, E. J., Caballero-Nieves, S. M., Gies, D. R., et al. 2015, AJ, 149, 26 [Google Scholar]
  2. Amara, A., & Quanz, S. P. 2012, MNRAS, 427, 948 [Google Scholar]
  3. Anderson, J., & King, I. R. 2000, PASP, 112, 1360 [NASA ADS] [CrossRef] [Google Scholar]
  4. Arias, J. I., Walborn, N. R., Simón Díaz, S., et al. 2016, AJ, 152, 31 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bally, J., & Zinnecker, H. 2005, AJ, 129, 2281 [NASA ADS] [CrossRef] [Google Scholar]
  6. Barbá, R. H., Gamen, R., Arias, J. I., et al. 2010, Rev. Mex. Astron. Astrofis. Conf. Ser., 38, 30 [Google Scholar]
  7. Barbá, R. H., Gamen, R., Arias, J. I., & Morrell, N. I. 2017, in The Lives and Death-Throes of Massive Stars, eds. J. J. Eldridge, J. C. Bray, L. A. S. McClelland, & L. Xiao, 329, 89 [Google Scholar]
  8. Baume, G., Carraro, G., Comeron, F., & de Elía, G. C. 2011, A&A, 531, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Beuther, H., Walsh, A. J., Johnston, K. G., et al. 2017, A&A, 603, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Bodensteiner, J., Sana, H., Mahy, L., et al. 2020, A&A, 634, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Bonatto, C., Santos, J. F. C., Jr., & Bica, E. 2006, A&A, 445, 567 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Bonnell, I. A., Bate, M. R., & Zinnecker, H. 1998, MNRAS, 298, 93 [Google Scholar]
  14. Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 2001, MNRAS, 323, 785 [Google Scholar]
  15. Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Castelli, F., & Kurucz, R. L. 2003, in Modelling of Stellar Atmospheres, eds. N. Piskunov, W. W. Weiss, & D. F. Gray, 210, A20 [Google Scholar]
  17. Chini, R., Hoffmeister, V. H., Nasseri, A., Stahl, O., & Zinnecker, H. 2012, MNRAS, 424, 1925 [Google Scholar]
  18. Close, L. M., Puglisi, A., Males, J. R., et al. 2012, ApJ, 749, 180 [NASA ADS] [CrossRef] [Google Scholar]
  19. De Becker, M., & Raucq, F. 2013, A&A, 558, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Delorme, P., Meunier, N., Albert, D., et al. 2017, in SF2A-2017: Proceedings of the Annual Meeting of the French Society of Astronomy and Astrophysics, eds. C. Reylé, P. Di Matteo, F. Herpin, et al., 347 [Google Scholar]
  21. Duchêne, G., & Kraus, A. 2013, ARA&A, 51, 269 [Google Scholar]
  22. Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Galicher, R., Boccaletti, A., Mesa, D., et al. 2018, A&A, 615, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Gomez Gonzalez, C. A., Wertz, O., Christiaens, V., Absil, O., & Mawet, D. 2016, Astrophysics Source Code Library [record ascl:1603.003] [Google Scholar]
  26. Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Ilee, J. D., Cyganowski, C. J., Brogan, C. L., et al. 2018, ApJ, 869, L24 [Google Scholar]
  28. Inutsuka, S.-I. 2012, Progr. Theor. Exp. Phys., 2012, 01A307 [CrossRef] [Google Scholar]
  29. Kobulnicky, H. A., Smullen, R. A., Kiminki, D. C., et al. 2012, ApJ, 756, 50 [NASA ADS] [CrossRef] [Google Scholar]
  30. Kratter, K. M., & Murray-Clay, R. A. 2011, ApJ, 740, 1 [Google Scholar]
  31. Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016, SPIE Conf. Ser., 9908, 990834 [Google Scholar]
  32. Maíz Apellániz, J. 2010, A&A, 518, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Maíz Apellániz, J., Sota, A., Arias, J. I., et al. 2016, ApJS, 224, 4 [CrossRef] [Google Scholar]
  34. Maíz Apellániz, J., Pantaleoni González, M., Barbá, R. H., et al. 2018, A&A, 616, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Maíz Apellániz, J., Trigueros Páez, E., Negueruela, I., et al. 2019, A&A, 626, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Maíz Apellániz, J., Crespo Bellido, P., Barbá, R. H., Fernández Aranda, R., & Sota, A. 2020, A&A, 643, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Maíz Apellániz, J., Barbá, R. H., Fernández Aranda, R., et al. 2022, A&A, 657, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556 [Google Scholar]
  39. Martayan, C., Floquet, M., Hubert, A. M., et al. 2008, A&A, 489, 459 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Martins, F., Schaerer, D., & Hillier, D. J. 2005, A&A, 436, 1049 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Mason, B. D., Henry, T. J., Hartkopf, W. I., ten Brummelaar, T., & Soderblom, D. R. 1998, AJ, 116, 2975 [NASA ADS] [CrossRef] [Google Scholar]
  42. Mason, B. D., Hartkopf, W. I., Gies, D. R., Henry, T. J., & Helsel, J. W. 2009, AJ, 137, 3358 [Google Scholar]
  43. Maud, L. T., Cesaroni, R., Kumar, M. S. N., et al. 2019, A&A, 627, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Mawet, D., Milli, J., Wahhaj, Z., et al. 2014, ApJ, 792, 97 [Google Scholar]
  45. Mayer, P., Harmanec, P., Chini, R., et al. 2017, A&A, 600, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Moeckel, N., & Clarke, C. J. 2011, MNRAS, 410, 2799 [NASA ADS] [CrossRef] [Google Scholar]
  47. Neckel, T., & Chini, R. 1981, A&AS, 45, 451 [NASA ADS] [Google Scholar]
  48. Oliva, G. A., & Kuiper, R. 2020, A&A, 644, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Pantaleoni González, M., Maíz Apellániz, J., Barbá, R. H., & Reed, B. C. 2021, MNRAS, 504, 2968 [CrossRef] [Google Scholar]
  50. Peter, D., Feldt, M., Henning, T., & Hormuth, F. 2012, A&A, 538, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Pozo Nuñez, F., Chini, R., Barr Domínguez, A., et al. 2019, MNRAS, 490, 5147 [CrossRef] [Google Scholar]
  52. Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Rainot, A., Reggiani, M., Sana, H., et al. 2020, A&A, 640, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Rainot, A., Reggiani, M., Sana, H., Bodensteiner, J., & Absil, O. 2022, A&A, 658, A198 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Russeil, D., Adami, C., Bouret, J. C., et al. 2017, A&A, 607, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Sana, H., Gosset, E., & Evans, C. J. 2009, MNRAS, 400, 1479 [NASA ADS] [CrossRef] [Google Scholar]
  57. Sana, H., Momany, Y., Gieles, M., et al. 2010, A&A, 515, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Sana, H., James, G., & Gosset, E. 2011, MNRAS, 416, 817 [Google Scholar]
  59. Sana, H., de Koter, A., Garcia, M., et al. 2012, The Messenger, 148, 33 [NASA ADS] [Google Scholar]
  60. Sana, H., de Koter, A., de Mink, S. E., et al. 2013, A&A, 550, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Sana, H., Le Bouquin, J. B., Lacour, S., et al. 2014, ApJS, 215, 15 [Google Scholar]
  62. Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23 [Google Scholar]
  63. Shull, J. M., & Danforth, C. W. 2019, ApJ, 882, 180 [NASA ADS] [CrossRef] [Google Scholar]
  64. Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593 [Google Scholar]
  65. Sota, A., Maíz Apellániz, J., Morrell, N. I., et al. 2014, ApJS, 211, 10 [Google Scholar]
  66. Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28 [Google Scholar]
  67. Trundle, C., Dufton, P. L., Hunter, I., et al. 2007, A&A, 471, 625 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Turner, N. H., ten Brummelaar, T. A., Roberts, L. C., et al. 2008, AJ, 136, 554 [NASA ADS] [CrossRef] [Google Scholar]
  69. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Zinnecker, H., & Yorke, H. W. 2007, ARA&A, 45, 481 [Google Scholar]
  71. Zurlo, A., Vigan, A., Galicher, R., et al. 2016, A&A, 587, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Tables

Table 1.

O-type stars used in this work.

Table 2.

Observing setup and atmospheric conditions for FLUX (F) and OBJECT (O) IFS observations.

Table 3.

Observing setup and atmospheric conditions for FLUX (F) and OBJECT (O) IRDIS observations.

Table 4.

Assumed stellar parameters for the FASTWIND calculation of the calibrated spectra.

Table 5.

Angular separations (ρ), position angles (PA), magnitude contrasts in the K1 and K2 bands (ΔK12), and spurious alignment probabilities (Pspur) for the sources with Pspur < 5% around our targets.

Table 6.

Summary of detected companions (Pspur ≤ 5%) in IFS and IRDIS for all objects.

All Figures

thumbnail Fig. 1.

Final PCA/SDI IFS images for the targets with a detected source. The position of the central star is indicated by a white cross.

In the text
thumbnail Fig. 2.

Final post-processed PCA/ADI IRDIS images for the first 9 targets.

In the text
thumbnail Fig. 3.

Final post-processed PCA/ADI IRDIS images for the last 8 targets.

In the text
thumbnail Fig. 4.

Observed spectrum for HD 152 623 B.

In the text
thumbnail Fig. 5.

Observed spectrum for BD −13° 4929 C.

In the text
thumbnail Fig. 6.

Observed spectrum for BD −13° 4929 D.

In the text
thumbnail Fig. 7.

Observed spectrum for HD 76 341 B.

In the text
thumbnail Fig. 8.

Observed spectrum for CPD −41° 7721 C.

In the text
thumbnail Fig. 9.

Observed spectrum for HD 152 200 B.

In the text
thumbnail Fig. 10.

Sensitivity of our IFS observations expressed as magnitude difference vs. angular separation to the central star. The curves for each objects correspond to the 5-σ contrast.

In the text
thumbnail Fig. 11.

5-σ contrast curves for our observations in K1 (left) and K2 (right) IRDIS bands.

In the text
thumbnail Fig. 12.

Δmag vs. separation diagram. The location of the newly discovered sources (black and grey dots), as well as the median IRDIS and IFS contrast limits are compared with the outcome of previous surveys (Sana et al. 2014). Black, grey, and light-grey dots corresponds to sources with Pspur ≤ 5%, 5%< Pspur ≤ 20%, and Pspur > 20%, respectively.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.