© The Authors 2023

1 Introduction

Massive stars are key components of the evolution of their host galaxy. They are the main producers of heavy elements, and the momentum and kinetic energy involved in their death have an influence on a large part of their galaxy (Zinnecker & Yorke 2007). They are also the progenitors of the compact objects that, when they merge, produce gravitational wave bursts that we can currently detect (Abbott et al. 2016).

However, their short lifetime (a few million years) and their rapid formation process (10⁵ yr) make the observation of their early ages difficult (Tan et al. 2014). Indeed, their lifetime makes them rare, and so to observe a large number of massive stars (>100), one needs to look for them at significant distances, typically 1–3 kpc. In addition, the majority of young massive stars are still embedded in a cloud of gas and dust when they finish their formation process (Zinnecker 2006), making the observation of this formation step even harder. In consequence, the formation process of massive stars is still actively discussed.

Historically, the standard models of star formation could not explain the formation of stars with masses significantly higher than about 10 M_⊙. The main difficulty is overcoming the radiation barrier emitted by the protostar as soon as it starts burning nuclear fuel. So, specific formation models need to explain how massive stars can form. There are currently three main scenarios: (1) the core accretion (Terquem 2001; Yorke 2002), which uses a massive accretion disk to accrete more matter; (2) the competitive accretion (Larson 1978; Zinnecker 1982; Bonnell & Bate 2002; Bonnell et al. 2003) for which close protostar cores use the combined gravitational potential to attract matter from further away than with each individual gravitational potential; and (3) the collision (Binney et al. 1988; Bonnell & Bate 2005; Dale & Davies 2006) in which intermediate-mass stars collide to merge and form a more massive star. The outcomes of these three different formation models predict different multiplicity parameters. Hence, the study of the multiplicity of massive stars, after their formation process, should provide relevant constraints on these formation models.

Another motivation for the investigation of the multiplicity of massive stars, especially in the range of periods considered in this study, is the investigation of physical processes driven by the colliding winds. In particular, such systems are well suited for particle acceleration, hence the class of particle-accelerating colliding-wind binaries (PACWBs; De Becker & Raucq 2013; De Becker et al. 2017). Such systems, mainly revealed by synchrotron radio emission, are likely contributors to the population of lower energy Galactic cosmic rays. Appropriate knowledge of their orbit is required to interpret their behavior and model their shock physics.

To probe the full range of orbital separations of systems situated at a typical distance of 2 kpc, one needs to use different observational techniques. The close companions (up to 0.5 millisec of arc, mas, Mahy et al. 2020) can be probed by photometry (eclipsing binaries) or spectroscopy (radial velocity, Sana et al. 2012, 2013; Barbá et al. 2017), while wide companions (separation >50 mas) can be probed with techniques such as adaptive-optics-assisted direct imaging, aperture masking, speckle imaging, and coronography (Turner et al. 2008; Mason et al. 2009; Reggiani et al. 2022). For separations between 1 and 50 mas, the only technique we can use is optical long baseline interferometry (OLBI; see Fig. 1 in Sana 2017). But until recently, this technique was limited by its sensitivity and could only be applied to a modest sample of massive stars.

The advent of the PIONIER (Precision Integrated-Optics Near-infrared Imaging ExpeRiment) instrument at the very large telescope interferometer (VLTI; Le Bouquin et al. 2011) enabled the southern massive stars at high angular resolution (SMASH+), the first systematic interferometric large survey on massive stars in the Southern Hemisphere (Sana et al. 2014), probing the missing range of separation to have a complete statistical study of the multiplicity of massive stars. With this survey, Sana et al. (2014) observed 96 southern O-type star systems, nearly reaching the 100 targets required to get a statistical error <5% over the entire range of the multiplicity fraction. The observational constraints brought by the SMASH+ survey, especially the abundance of companions with a separation smaller than 100 AU, are qualitatively in overall agreement with the core accretion model leading to disk fragmentation. However, the statistic on subgroups of stars, such as masses or evolutionary stage, is too low to obtain a robust conclusion. Therefore, we aim to perform a similarly large survey (120 objects) in the Northern Hemisphere to double the total statistic.

To observe more than 100 O-type stars in the Northern Hemisphere, one needs to reach a limiting magnitude in the J-band (the spectral band available in the GOSC catalog, Maíz Apellániz et al. 2013) of J = 7.5, as shown in Fig. 1. In this figure, the limit of declination >–20^º corresponds to the limit of observability of the CHARA (center for high angular resolution astronomy) Array. Thanks to the recent implementation of the MIRC-X (Michigan InfraRed Combiner-eXeter) instrument (Kraus et al. 2018; Anugu et al. 2020; Lanthermann et al. 2019) at the CHARA Array (ten Brummelaar et al. 2016), located at Mount Wilson observatory, USA, this magnitude is now reachable in the H-band with the OLBI technique in the Northern Hemisphere. As O-type stars are hot stars and under typical reddening conditions, their magnitudes remain relatively the same in all the infrared spectral bands (Martins & Plez 2006), meaning that the limit on the magnitude of 7.5 in the J-band required to observe >100 O-type objects remains valid in the H-band.

In this paper, we present the results of the pilot survey performed on 29 systems, with the goal to demonstrate the feasibility of a large survey of more than 100 O-type stars with the CHARA/MIRC-X instrument. We present the observations in Sect. 2, with the definition of the sample, the description of the observation campaign, and the data reduction process. We then describe the data analysis that we performed in Sect. 3. Section 4 presents our results. We perform a statistical analysis of the results in Sect. 5. We finally discuss the results in Sect. 6 and conclude in Sect. 7.

Fig. 1 Cumulative histogram of the number of O-type stars as a function of their magnitude in the J-band. In blue, for the whole sky. In orange, for a declination greater than –20º to be reachable with the CHARA Array. Data from the Galactic O-star catalog (GOSC; Maíz Apellániz et al. 2013).

2 Observations

2.1 Observation sample

We built our sample using the Galactic O-star catalog (GOSC; Maíz Apellániz et al. 2013). We selected every O-star with a declination Dec > –20^º and with a magnitude in the J-band registered in GOSC of J < 7.0 mag. We chose the J-band criterion for several reasons. The first one is that the magnitude criteria available in GOSC are either in the B-band or the J-band. As we used the MIRC-X instrument, working in the H-band, the J-band is the nearest one of the two available. We adopted the threshold value of 7.0 or brighter because, during this pilot study, the official limiting magnitude offered by MIRC-X was H = 6.5 mag. As only the Rayleigh-Jeans tail of the O-type stars is observed in the J- and H-bands, their magnitudes are comparable in both bands. We took an extra 0.5 magnitude as a margin to be sure that our sample contains all targets observable with MIRC-X. We then looked for the H-band magnitude of the selected systems in the 2MASS All-Sky Catalog of Point Sources (Cutri et al. 2003) and performed the last selection on stars with a magnitude H < 6.5 to comply with the MIRC-X limiting magnitude. Our input sample is therefore magnitude limited. We note that some of the targets in our sample overlap with targets already observed by SMASH+, which can be used to validate our results.

2.2 Observation campaign

The MIRC-X beam combiner operates in the J- and H-band. The observations presented in this paper are performed only in the H-band (1.65 μm) because the J-band mode was still experimental during the pilot survey and only uses four telescopes. The six telescopes provide sufficient coverage of the uv-plane to constrain the multiplicity of a star in a single snapshot. This can also be done for data that combine only five telescopes when the conditions (weather, technical, operational) would not permit a six-telescope observation. The MIRC-X combiner allows different spectral resolutions. Our data have been taken with the PRISM-50 configuration, allowing a spectral resolving power of R ~ 50. This configuration was chosen to optimize the sensitivity of the beam combiner and because it brings the Outer Working Angle (OWA) of MIRC-X to:

(1)

where λ is the central wavelength which in the H-band is equal to 1.6 × 10^–6 m, and B is the length of the baseline, equal to 330 m for the longest baseline at CHARA. This OWA allowed us to fill the gap in angular separation that other techniques cannot reach. We note that detection of companions with a separation larger than the OWA is still possible, but the flux ratio will be biased by the bandwidth smearing, hence, hampering accurate measurement of the contrast (Hummel et al. 2016).

We note that the Inner Working Angle is the angular resolution of the instrument, meaning that a binary separated by less than this angle would not be resolved, and is defined by:

(2)

This angular resolution is about a factor of two to three smaller than the SMASH+ survey owing to the larger baseline B of CHARA compared to the VLTI.

For the calibration strategy, we alternated a calibrator with a science target. The calibrator was chosen with the tool SEARCH-CAL (Chelli et al. 2016), developed by the JMMC¹. The selection criterion was the calibrators needed to be at most 1.5 magnitudes brighter and 0.5 magnitudes fainter than the science target it would calibrate, and situated at a maximum angular distance of 3 degrees on the sky. Most of our calibration stars are of spectral type KIII, for which a sufficiently accurate diameter (a few percent) can be estimated from the apparent photometry (Chelli et al. 2016).

The observations with MIRC-X have been carried out during five runs spread over three observation semesters from 2018 to 2019. The observation time was obtained through NOAO² (now called NOIRLab) community access time (program IDs: 2018A-M12/NOAO2, 2018B-M17/NOAO4, 2019A-M16/NOAO4; PI: C. Lanthermann). We also used one night that studied O-type stars during which O-type stars were observed as backup targets because the original program could not be executed, in December 2017 (PI: S. Kraus).

During this campaign, we could obtain good quality data on 29 O-type stars, listed in Table 1 with information on the spectral type, distance, and the number of detected companions in various separation ranges, from the literature as well as those detected by this study. Figure 2 shows the histogram of O-type star systems for which we obtained good-quality data as a function of their magnitude in the H-band. As the official magnitude limit in typical conditions is H = 6.5 mag, the bulk of observed objects is around this limit. But as we advanced in time, the improved knowledge of the MIRC-X instrument acquired during commissioning allowed us to push for fainter targets, beyond the initial magnitude limit of the instrument. Ultimately, observing an O-type star with a magnitude up to H = 8.1 mag was possible in excellent seeing conditions, as well as a couple of systems around H = 7.5 mag in normal seeing conditions. These latter observations are important to demonstrate that one can observe a large number of O-type stars with MIRC-X. Indeed, as shown in Fig. 1, we can observe up to 120 systems if a magnitude limit of H = 7.5 mag can be reached.

Usually, to calibrate a science target with the OLBI technique, we use “CAL1-SCI-CAL2” sequences, meaning that we observe a first calibrator (CAL1), then the science target (SCI), and finally a second calibrator (CAL2). To choose CAL1 and CAL2, we take unresolved targets within a reasonable distance and a similar apparent brightness to the science target, as explained above. It appears that two consecutive science targets observed in this program are close enough in the sky and in magnitude that CAL2 of one science target can be used as CAL1 of the next science target. Therefore, we chose to follow an observing sequence such as “CAL1-SCI1-CAL2-SCI2…”, using each best-suited calibrator for each science target. This allows us to reduce the error due to calibration. When the same calibrator was best suited for two science targets observed one after the other(ex: SCI1 and SCI2), we preferred to use anothercalibrator for the second science target instead of using the same calibrator for both science targets (CAL1-SCI1-CAL1-SCI2). This reduces the risk that the results might be affected by a bad calibrator. A list of the calibrators observed but discarded because suspected to be actually multiple systems can be found in Appendix A.

2.3 Data reduction

We used the MIRC-X data reduction pipeline³ (Anugu et al. 2020). We set the maximum integration time (max-integrationtime keyword in the pipeline) for one calibrated file to 220 s. The max-integration-time set the maximum time one reduced file will cover, with one file giving one constituent of the calibrated parameters. This means that one file will give the average calibrated parameters over the given maximum integration time. This time allows us to bin our data recording sequences of about 10 min into three files with a similar signal-to-noise ratio (S/N). This is needed because the change of the uv-coordinates due to the rotation of the Earth affects the observed squared visibility (V2). V2 can indeed change significantly in a timescale of 300 s for binaries at the limit of the OWA (see Fig. 3 for an example).

To get rid of some outlier points due to bad S/N data, we applied a threshold on the S/N of four instead of the default value of 3. This limit allows us to select only data with high enough quality, while at the same time preserving as much data as possible. To increase the quality of the reduced data, one would like to coherently sum the data over as long a time as possible. However, the atmospheric conditions limit the amount of time one can coherently integrate, as the phase induced by the atmosphere blurs the fringes, therefore reducing the quality of the data. The optimum time interval over which we can sum data for V2 and closure phase (CP) can however be different. To find this optimum value, we reduced the data for different values of the coherent integration time (ncoherent keyword in the reduction pipeline), which is the number of frames we add together to increase photon S/N, a frame being one recorded image of the detector. The pipeline produces a measurement of the calibrated parameters for each of these coherently added frames, and it is these coherently added frames that are incoherently averaged in the final calibrated file. Then, for each baseline and each target, we plotted the mean value over the files of the V2 S/N and the CP error as a function of the coherent integration time. We finally merged the V2 of the files with the higher S/N on V2 and the CP of the files with the lower error on CP, for most of the baselines or triplets. For the V2 we chose the same coherent integration time throughout the night. For the CP we took the best coherent integration time for each target. The number of coherently added frames for every target is specified in Table B.1, each frame integrating the flux for about 3 ms.

The data reduction pipeline provides two options to unbias the CP. We chose the option that computes an estimation of the bias as it does not require tuning any extra parameters (in opposition to the second method). This unbiasing method is explained in Appendix B of Anugu et al. (2020).

Fig. 2 Histogram of the H-band magnitude of our observed sample in green, the original sample planned for this pilot survey in orange, and the planned sample for the large survey in blue.

Fig. 3 Example of change in the uv-plane modeled with the ASPRO2 software. Left: global view of the uv-plane. Different colors are for different pairs of telescopes (baselines). Right: zoom on the change in the uv-plane for the fastest moving baseline. Each straight black line indicates the uv-coordinate for this baseline for an instant snapshot, and each snapshot is separated by 5 min. The regularly spaced, orange lines display the V2 values modeled for a binary with a separation of 50 mas.

3 Data analysis

To analyze the reduced data, and look for companions in the interferometric signal, we used the CANDID⁴ software (Gallenne et al. 2015). CANDID is a tool developed to look for binarity signals in interferometric data as well as to determine the position and flux ratio of the companion(s). For nondetec-tions, CANDID also provides upper limits on the contrast of potential companions. More information about the algorithms and methods can be found in Gallenne et al. (2015). To streamline the analysis of the data given the relatively large size of our sample, we established an automated procedure to analyze all the data consistently and uniformly.

This procedure is summarized here. First, for each observation of a star, we input all the reduced data to CANDID. We then perform a first search for a companion, fixing the maximum separation to 50 mas, which is approximately our OWA (see Eq. (1)). We fixed the step of the search grid to 1.0 mas. This step size is fine enough to find the global minimum in the range of separation we are looking for, while it is large enough to preserve reasonable computation time (approximately an hour per observation). If a companion is detected (nσ > 5, with nσ being the significance of the binary model compared to a uniform disk model fitting the data), we use the bootstrap function of CANDID around the position of the found companion. This gives us errors on the position and flux ratio of the companion that are more realistic than those computed by the initial grid search method. Then, if a first companion is detected, we analytically remove the signal of this first companion from the interferometric data, and we perform a search for a second companion, using the same parameters as for the first search. We note that the detection of a second companion is made with an indirect method assuming the signal of the first companion is perfectly analytically removed, with no residual. Therefore, these detections should be taken with caution and would need further observations to confirm them. Finally, we determine the upper limiting ∆magnitude of a companion detectable in the data. If no companion is detected, we compute this limit on the reduced data directly. It gives us the limit in ∆magnitude for which we would detect a companion. If a first companion is detected, we perform this computation on the data with the signal of the first companion analytically removed. It then gives us the limit in ∆magnitude for which we would detect a second companion. We note that we can only remove the signal of one companion, so, the results of this detection limit might be biased if the signal of a second companion is present in the data.

Given the distance of these systems, even optical interfer-ometry with 330 m baselines is not resolving the diameter of the stars in those systems. Hence, we consider them as a point source in CANDID, reducing the number of free parameters in the search for companions.

4 Results

A summary of the results is shown in Table 2 for the first companion search and Table 3 for the search of a second companion, with parameters that characterize the systems.

4.1 New detections

We here summarize the newly detected companions.

Cyg OB2-5 A / BD+40 4220A / Schulte 5. This system includes a short 6.6 d period binary, along with other components on wider orbits (Rauw et al. 1999, 2019). Two distant companions were already known at separations of 0.93″ and 5.55″ (Maíz Apellániz 2010; Caballero-Nieves et al. 2014, 2020), which is largely outside the OWA of our observations. Cyg OB2-5 A is a known PACWB, with the nonthermal radio emission mainly associated with a wide orbit with a period of about 6.7 yr (Kennedy et al. 2010). We detect a companion in both observations of June 2018 and June 2019, at a mean separation of 12.18±0.35 mas and a mean ∆H = 1.63 ± 0.3 mag, and a second companion in the observation of June 2019 at a separation of 6.51±0.26 mas and a ∆H = 4.13 ± 0.01 mag. The computed detection limit for the search of a second detection for the observation of June 2018 is 4 magnitudes fainter than the primary component, while the detected second companion in June 2019 is 4.13 magnitudes fainter. This could explain why we did not detect it in the June 2018 data, but further observations are necessary to confirm this second companion. This is the first direct detection for both companions. The wind-wind collision between the combined wind of the inner system and the one of the newly detected companions is likely the cause of the synchrotron radio emission. Long-term monitoring should allow us to determine which is associated with the 7-yr period.

Cyg OB2-10/ 2MASS J20334610+4133010 / Schulte 10. Three companions have already been observed previously at separations of 0.21″, 0.74″ and 4.16″ with ∆K = 2.80 ± 0.78, 5.24 ± 0.05, and 6.03 ± 0.07 mag (Caballero-Nieves et al. 2020). In the observation of June 2018, we detected a companion at a separation of 7.35±0.20 mas and ∆H = 2.45 ± 0.03 mag, but we do not detect it again in the observation of June 2019. More observations could be useful to confirm this new companion detection.

HD 17505. This object is known as a hierarchical triple system (Hillwig et al. 2006; Sota et al. 2014). The system includes a close binary with an orbital period of 8.571 days, separated from the primary component by 2.161″ (Maíz Apellániz et al. 2019). We detect an additional companion at a separation of 15.43±0.02 mas from the primary component and ∆H = 0.35 ± 0.04 mag in two different observations separated by 1 day. This companion is detected for the first time, making this system a quadruple one.

HD 19820. This system is a known spectroscopic binary (SB) with a period of 3.36 days (Hill et al. 1994). We detect a companion at a mean separation of 13.87±0.03 mas and ∆H = 2.57 ± 0.01 mag in two different observations separated by 1 day, and another companion at a separation of 6.96±0.11 mas and ∆H = 4.16 ± 0.01 mag in only one of these observations. The detection of this second companion is made with an indirect method assuming the signal of the first companion is perfectly analytically removed, with no residual. Therefore, this detection should be taken with caution and will necessitate further observations to be confirmed.

HD 24431. This system has a companion situated at 0.72″ (Mason et al. 1998; Turner et al. 2008). We detect a new companion with a separation of 9.74±0.02 mas and ∆H = 1.37 ± 0.01 mag.

HD 28446. This system is supposed to be an SB2 (double line spectroscopic binary) with a period of 1.31 days in Mayer et al. (1994) but has not been confirmed since. A third component is known at a separation of 10″ (Eggleton & Tokovinin 2008). We detect a fourth companion in both our observations of February 2018 and June 2018 with a mean separation of 26.16+0.13 mas and a mean ∆H = 1.35 ± 0.08 mag.

HD 34078. This star is a known runaway (Hoogerwerf et al. 2001). Candidate companions have been detected by direct imaging at a separation of 8.4″ (Mason et al. 1998) and 0.35″ (Turner et al. 2008), but both those detections are suspected to be field stars observed in the line of sight of this star. We detected a new companion in two different observations, at a separation of 6.85±0.07 mas and ∆H = 2.76 ± 0.02 mag in December 2017 and a separation of 1.74±0.20 mas and ∆H = 3.29 ± 0.03 mag in September 2018.

HD 36861. This star (λ Ori A) is known to be variable (Fullerton et al. 1996), and part of a wide binary (HD 36861J). The components A and B of the system λ Ori may not be physically bound (Lindroos 1985; Mason et al. 1998). λ Ori A had no detected companion yet (Mason et al. 2009), but we detect a companion in our observation of February 2018, with a separation of 10.13±0.05 mas and a ∆H = 3.30 ± 0.02 mag. This companion has not been detected again in our observation of June 2018, despite the computed limiting magnitude of detection being ∆H = 5.03 mag. The reason for this nondetection on the second observation is unknown. A new observation of this system could confirm this new companion.

HD 47129. This is Plaskett’s star. It is a known SB2 with a period of 14.4 days (Linder et al. 2008). There are also two known visual companions at 0.78″ and 1.12″ (Turner et al. 2008). The SMASH+ survey resolved a faint companion at 36 mas with ∆H ~ 4.0 mag, with the NACO/SAM instrument, but it was too faint to be confirmed by PIONIER and the uncertainties on the separation found by NACO/SAM were large. We detect a companion at a separation of 32.29±0.06 mas and ∆H = 4.6 ± 0.01 mag. This detected companion is a confirmation of the candidate by NACO/SAM, providing compelling evidence for the existence of this companion. Interestingly enough, this long-period companion is compliant with the likely membership of this system to the category of PACWBs, based on radio results published by Kurapati et al. (2017).

HD 207198. This system has a known companion at a separation of 17.64″ (Mason et al. 2004). We detect for this system a companion at a separation of 41.07±0.04 mas and ∆H = 4.68 ± 0.01 mag. This companion is detected for the first time. We note that this detection is close to the OWA.

HD 229196. We detect for the first time in this system a companion, situated at a separation of 5.88±0.02 mas and amean ∆H = 2.80 ± 0.04 mag in two different observations separated by 1 day.

4.2 Already detected

We here summarize the redetection of companions that were already known from other techniques or previous optical inter-ferometric observations.

Cyg OB2-9/HIP 101419. This PACWB is a known very excentric SB2 (Blomme et al. 2013; Maíz Apellániz et al. 2019; Caballero-Nieves et al. 2020) with a period estimated at 2.4 yr and an excentricity of 0.713 (Nazé et al. 2010, Nazé et al. 2012) with another companion at a separation of 21” (Maíz Apellániz 2010). We detect a companion at a mean separation of 0.77±0.01 mas and a mean ∆H = 0.42 ± 0.04 mag in two observations separated by two days. With the orbital parameters in Blomme et al. (2013) we derive a minimum projected separation of ~2.04 AU with the known SB2 companion. Taking the distance of the system and the angular separation we detect our companion, we compute a projected separation of ~1 .38 AU. Taking the uncertainties into account, our detection is probably the SB2 component observed close to its periastron. Further observations could confirm it.

HD 47839/15 Mon. This system has a known companion at ≃ 0.1″ with a difference of magnitude of 1.6 in the visible (Hutter et al. 2021) and a third component at a wider (3″) separation (Mason et al. 1998; Sana et al. 2014). We detect a companion with a separation of at least 49.19±0.32 mas and ∆H = 1.81 ± 0.01 mag. This detection is at the limit of the OWA, which means that it is probable that the companion detected is further out, as discussed in Sect. 6.1. This companion is probably the known companion around 100 mas, as the differences in magnitude are compatible and as the separation close to the OWA cannot rule out this known companion.

HD 167971. This system is a known hierarchical triple system (Leitherer et al. 1987; De Becker et al. 2012; Le Bouquin et al. 2017; Sanchez-Bermudez et al. 2019), that also turns out to be the brightest synchrotron-emitting O-type PACWB in the catalog (De Becker & Raucq 2013). The central binary has an orbital period of 3.3 days and the third component is orbiting the inner binary on a timescale of 21.4 yr. We detect a companion at 19.89±0.01 mas and ∆H = 0.61 ± 0.01 mag. The separation of the detected companion is compatible with the one of the outer component of the system measured in De Becker et al. (2012) and Le Bouquin et al. (2017).

HD 193322. This is a complex multiple system. The A component consists of a single star Aa orbiting around a 312-day binary Ab in 35 yr (ten Brummelaar et al. 2011). Three other components are also known, with a separation of 2.6″ for the closest (Turner et al. 2008). We detect a companion at a separation of at least 47.33 mas and maximum mag in our three observations of this system in June 2018, in September 2018, and in June 2019. These detections are at the limit of the OWA, which means that it is possible that the companion is somewhat further out, as discussed in Sect. 6.1. The detected companion is most likely the already known pair Aa and Ab due to its separation being compatible with the known pair reported in ten Brummelaar et al. (2011).

HD 202214. One close companion is already known from spectroscopy, with an orbital period of 81.30 days (Mante 2002). The system also has two wider companions, one at a separation of 0.071″ (Mante 2002) and a ∆V = 0.6 (Mason et al. 2009) mag and one at 1.0″ with a ∆V = 0.3 (Mason et al. 2009) mag. We detect a companion at a separation of 47.27±0.05 mas and mag. This detection is at the limit of the OWA, which means that it is probable that the detected companion is actually further out, as discussed in Sect. 6.1. We, therefore, cannot rule out that the detected companion is the known one at 71 mas. The difference between our magnitude difference and the one from Mason et al. (2009) could then be explained by the separation potentially being larger than our OWA, introducing bias in our detection, and by the fact that our observing setup is not optimized for companions outside of the OWA. The derived contrast might therefore be systematically biased.

HD 206267. This system is a high order multiple system (Maíz Apellániz et al. 2019; Maíz Apellániz & Barbá 2020). The central component (AaAb) is composed of an SB2 (Aa) with a period of 3.71 days (Raucq et al. 2018) and another companion (Ab) separated from Aa by 0.1″ and a difference in magnitude of 1.63±0.3 at a wavelength of λ = 0.91 μm (SDSS z filter⁵; Maíz Apellániz et al. 2020). A third component (B) is situated at a separation of 1.7″ with a difference of 5.72±0.13 in magnitude at λ = 0.91 μm. Two other companions (C and D) are situated within 25″. We detect a companion at a separation of at least 49.52±O.22 mas with a difference of 1.64+0.02 mag in the H-band. This detection is at the limit of the OWA, which means that it is probable that this detected companion is further out, as discussed in Sect. 6.1. We, therefore, cannot rule out that this detection can be the known companion at 0.1"; the magnitude differences are compatible one with each other.

4.3 No detection

We do not detect interferometric companions around HD 30614 (SB1 system with a period of 3.68 days; Zeinalov & Musaev 1986), HD 47432 (one known companion at 0.78″; Turner et al. 2008), HD 188001 (single runaway star; Trigueros Páez et al. 2021), HD 195592 (spectroscopic binary with a period of a few days; De Becker et al. 2010), HD 201345 (one known wide companion at 7.38″; Turner et al. 2008), HD 206183 (two known wide companions at separations of 6.6″ and 11.6″; Mason et al. 1998), HD 209975 (four known wide companions at 3.79″, 4.14″, 19.8″, and 60.4″; Mason et al. 1998; Turner et al. 2008), HD 217086 (2 known companions at 2.8″ and 3.1″; Mason et al. 1998; Turner et al. 2008), and HD 228779.

4.4 Candidates

Some of our observations resulted in candidates with only a marginally significant detection criterion (3 < nσ < 5). Those systems are HD 45314 (Oe star Rauw et al. 2015), HD 210809, and HD 210839. More observations will be needed to validate or reject the presence of a companion. The position of those candidates can be found in Table 2.

5 Statistical analysis

Among the 29 systems observed with good-quality data, we confirm the detection of 19 companions for 17 multiple systems (see Sects. 4.1 and 4.2). Out of these 19 companions, 13 are detected for the first time. This gives us a multiplicity fraction ƒ_m = 17/29 = 0.59 + 0.09, and a companion fraction of ƒ_c = 19/29 = 0.66 + 0.13 in the range of separations to which we are sensitive. The uncertainty for the multiplicity fraction is obtained with a binomial uncertainty and the companion fraction uncertainty is obtained with a Poisson uncertainty.

We note that 38% of our sample corresponds actually to multiple systems with at least three components and that in 28% of the sample, the interferometric companion constitutes the detection of the outer orbit in a hierarchical triple system. This proportion of hierarchical triple systems should allow us to study the Kozai-Lidov effect (Naoz 2016), by studying their orbits.

Fig. 4 Histogram of detected companions as a function of the estimated mass ratio q, when taking into account only the first detected companions (blue), and when adding the second companions detected (orange).

5.1 Comparison with SMASH+

The multiplicity fraction in this study is marginally higher than the one in SMASH+ (0.41 ± 0.05). This difference could be explained by the fact that our preliminary magnitude limit (H < 6.5 mag) was brighter than the one of SMASH+ (H < 7.5 mag). This would push the observational bias to observe a larger fraction of multiple systems.

By taking into account all known companions, summarized in Table 1, the total multiplicity fraction (0.86 ± 0.07) and the total companion fraction (2.10 ± 0.27), are consistent with the one of SMASH+ (respectively 0.91 ± 0.03 and 2.1 ± 0.2). We also find that in our sample, 38 ± 9% of the systems contain spectro-scopic binaries, which is marginally lower than in the SMASH+ sample (49 ± 5%).

5.2 Estimated mass ratio distribution We define the mass ratio q as:

(3)

where M_comp is the mass of the detected companion and M_primary is the mass of the primary component of the system. This mass ratio can be estimated using the flux ratio between the two components. In this study, we use the same relation used in Le Bouquin et al. (2017) and derived from Martins et al. (2005), which gives a good approximation of the mass ratio for main-sequence stars:

(4)

where ƒ_H is the flux ratio in the H-band, given by the CANDID analysis. We note that if the central component is an unresolved binary, this method of computing the mass ratio is using the combined flux of both components, hence the estimated mass ratio will be biased.

Figure 4 displays the distribution of the estimated mass ratio of the detected companions. In this figure, we can see that the distribution seems to be bi-modal, with a lack of companions between q = 0.4 and 0.6, and favoring a lower mass ratio. The bi-modality could be explained by the small statistic we are using here, with only 17 companions. However, we ran a Kuiper test to compare the estimated mass ratio distribution with a uniform distribution. The result of this test is a value of D = 0.41, and a probability of obtaining the value D from a uniform distribution of 2.7%. From this low probability, we can conclude that the actual mass ratio distribution is most probably not uniform.

The distribution favoring low mass ratios goes against the observational bias of the survey being magnitude limited. This bias should favor the inclusion of binaries with bright companions, therefore, with a flux ratio, hence a mass ratio, close to 1. So, this tendency seems to come from the intrinsic O-type stars' mass ratio distribution. The results of the survey of the massive stars in the Orion region (GRAVITY Collaboration 2018) show a similar trend, with a mass ratio distribution following a power law ∝ q^α with α = 1.7.

Fig. 5 Histogram of detected companions as a function of the projected separation in Astronomical Unit (AU). The distribution in blue takes into account only the first companions. The distribution in orange considers all companions.

5.3 Projected separation distribution

In the absence of any estimate of the inclination of the orbit, the absolute physical separation cannot be determined. Rather, we determined the projected separation of the detected companions, using the distances published by Bailer-Jones et al. (2021). Figure 5 shows the distribution of the separation of the companions detected in this study. The distribution seems to favor the middle and high part of the probed separation, from 10 to 100 AU. We note that the companions detected close to the OWA may actually be located further out, meaning that the real distribution could have a tail at larger separations. This figure only shows the companions that were (re)detected in this study.

5.4 Estimated mass ratio as a function of the projected separation

Figure 6 shows the estimated mass ratio q as a function of the projected separation of the companions detected in this study. Error bars on separation ending with an arrow pointing right are to distinguish the companions detected close to the OWA, meaning that the physical separation on the plot is only a lower limit. The black points are for the companions detected with other techniques for which the mass ratio and physical separation are known or can be determined from the literature.

We notice the higher number of companions with an estimated mass ratio q < 0.6 as described in Sect. 5.2. However, we cannot discern any correlations in this plot between the estimated mass ratio q and the physical separation for the companions detected in this study as the companions seem relatively homogeneously distributed in the figure.

Furthermore, taking into consideration the companions previously detected by other techniques, a correlation seems to appear, with q being anticorrelated with the separation. This anticorrelation would be in favor of the competitive accretion formation process (GRAVITY Collaboration 2018). But one should note that this result could be due to the difference of biases in the detection of companions by the different techniques. A further study taking into account the different biases is necessary to bring a stronger conclusion.

We caution that the plot is affected by some factors that deserve a few comments. First of all, we certainly need a greater statistical sample if we want to be sure that there is or is not a correlation between these two parameters for massive star systems. For subsequent developments, we stress that we should ideally base our discussion on absolute physical separations, and not on projected ones, which requires knowledge about the inclination. This can be obtained through a suitable interferometric follow-up to derive the relative astrometric orbit of the systems. Finally, the best test for a potential (lack of) correlation with the mass ratio should rely on an estimate of the semi-major axis, and not on the measured separation at one specific epoch. The available measurements have been obtained at any orbital phase, and in the case of a significantly eccentric orbit, the measured separation is not necessarily a good proxy for the semi-major axis.

Fig. 6 Estimated mass ratio as a function of the projected separation in AU of the detected companions. We note that some error bars are hidden by the size of the markers. Error bars on separation ending with an arrow pointing right are to distinguish the companions detected close to the OWA. The black points are for known companions detected with other observational techniques.

6 Discussion

6.1 Outer working angle

As we limited our search for companions to OWA, the results for the companions found close to this limit (separation > 40 mas, HD 47839, HD 193322, HD 202214, HD 207198, and HD 206267) might correspond to a local minimum of χ², while the global minimum may not be probed in our search. The detection of those companions is still valid, as the interferometric signal still favors the presence of a companion compared to a uniform disk, but the found separation should be considered as a lower limit. Further observations with a wider OWA could confirm the real separation of those companions, but this requires a higher spectral resolution, which is doable with MIRC-X, however, it would reduce the sensitivity of the instrument. It can be achieved only for the brightest systems. Alternative methods such as sparse aperture masking may thus be preferred.

6.2 Candidates for orbital parameter measurements

While the statistical analysis will require the large survey data, the currently detected companions can already give us good candidates for follow-up orbital parameter measurements. To be a good candidate for follow-up with interferometry, we take here a limit on the orbital period of a maximum of 10 yr, which is a reasonable time scale for follow-up observations.

To obtain an estimation of the orbital period of the detected companions, we use the Kepler’s third law:

(5)

where P is the period in years, a is the semi-major axis of the orbital ellipse in AU, and M₁ and M₂ are the masses of the two objects orbiting each other, in solar mass.

For our estimation of P, we approximate the semi-major axis a with the projected separation that we computed earlier. This will lead to an overestimate of a, hence P, for eccentric systems as companions spend a larger fraction of the orbit at separations >a (Kepler’s law). For the mass of the two objects, we use the m_evol in Table 4 of Weidner & Vink (2010), which gives a theoretical mass of O-type stars as a function of their spectral type. We used the spectral types of our objects listed in Table 1 to get an estimation of the mass for the central object, and we used the estimated mass ratio q for the detected companion mass.

With the flux ratio, separation, and estimated mass ratio, we could estimate also the shift of the photo-center of our binary systems. This estimated shift could give us candidates of systems for which we could determine the orbital parameters with Gaia, as discussed in Le Bouquin et al. (2017). From equation five in that paper we can determine that:

(6)

For our estimation of μ, we used the separation of our detected companions as an estimation of a, the estimated mass ratio q, and the detected flux ratio ƒ.

Figure 7 displays the estimated photo-center shift μ in mas as a function of the estimated orbital period in years. The detection limit on μ detectable by Gaia is set to 0.1 mas or higher. This value comes from the fact that Gaia needs a shift of at least three times the Gaia accuracy of a transit, which is 0.034 mas on bright targets (G < 12, see Sect. 4 of Perryman et al. 2014), which is the case for all our targets. The conservative limit in the period for which Gaia can determine the orbital parameters is set to 6 yr. This limit is the one for which Gaia will be able to determine orbital parameters from a detectable astrometric shift of the photo-center in almost any case. But this limit is for the nominal mission period of 5 yr. Now Gaia has observed for 7.8 yr and the anticipated mission lifetime is now at least 10 yr. So we estimate the expected limit on Fig. 7 at 10 yr, which corresponds also to the period we set for the interferometric survey.

First, we see that eight companions have an estimated orbital period of less than 10 yr, which makes them suitable to be followed up by interferometry directly. Then, we see that Gaia should be able to obtain the orbital parameters of seven of our systems, which have an estimated orbital period of approximately 6 yr or less and an estimated photo-center shift of more than 0.1 mas. For the periods larger than 10 yr, the orbital parameters will not be well constrained by Gaia, but combined with the data from interferometry, it should be possible to constrain orbital parameters well. Furthermore, the orbital parameters with Gaia are interesting because, in combination with the information from interferometry, one could measure the individual masses of each component of a multiple system. Currently interferometry needs to be combined with spectroscopy to measure individual masses (e.g., Le Bouquin et al. 2017; Many et al. 2018; Fabry et al. 2021; Sanchez-Bermudez et al. 2022). The technique with Gaia would not need spectroscopy because it will provide the distance, which with the orbit size from interferometry would provide the total mass of the system, and by combining the orbit size and the orbit size of the photo-center provided by Gaia, one would obtain the mass ratio of the components of the system. Combining both total mass and mass ratio would provide the individual masses of these components. Gaia would also probe a different range of separation than the spectroscopy. These two techniques would then be complementary. In addition, one would be able to compare the results of the two methods and remove the potential biases that each of these methods could present, as the range of masses estimated by spectroscopy alone shows discrepancies with the measured masses by combining spectroscopy and interferometry (Le Bouquin et al. 2017).

Fig. 7 Estimated photo-center shift μ in mas as a function of the estimated orbital period in years. The red lines show the limit of capability of Gaia to measure the orbital parameters of a system, with the hatched area being not accessible by Gaia.

7 Conclusion

From the results of this pilot study, we can conclude that a large survey to study the multiplicity of northern O-type stars can be performed with the instrument MIRC-X at the CHARA Array. Indeed, we demonstrated that we can constrain the multiplicity of O-star systems with a magnitude in the H-band around 7.5 mag, in good atmospheric conditions. This magnitude allows the observation of more than 120 northern O-type stars. From the experience gained with this pilot survey, we can observe six to eight science targets per night with normal conditions. Taking into account an average loss of 25% of the night due to bad atmospheric conditions weather or technical issues, we estimate that the large program will require approximately 25 nights to be completed.

This study also detected 19 companions in 17 different systems, including 13 companions detected for the first time, notably the companion responsible for the nonthermal emission in Cyg OB2-5 A, and the confirmation of the candidate companion of HD 47129 previously suggested by SMASH+. The preliminary statistical study gives us a multiplicity fraction f_m = 17/29 = 0.59 ± 0.09, and a companion fraction of ƒ_c = 19/29 = 0.66 ± 0.13. Those results are consistent with the results of the southern large survey already performed, SMASH+.

We also demonstrated that a number of the detected systems are suitable for follow-up studies for orbital parameter measurement, either by interferometry (8) and/or with Gaia (7). The results obtained in this study are promising in terms of scientific returns of a more ambitious project focusing on a large sample, and involving repeated observations spread over several years.

Acknowledgments

Acknowledgements.

We would like to thank Carine Babusiaux, from IPAG, for her help in determining the capability of Gaia to help in the constraint of orbital parameters and Laetitia Rodet in her explanation of the relevance of the study of the Kozai-Lidov cycle. We also thank the referee for their constructive comments that resulted in a better paper. L.M. and E.G. acknowledge the European Space Agency (ESA) and the Belgium Federal Science Policy Office (BEL-SPO) for their support in the framework of the prodex programme. This work is based upon observations obtained with the Georgia State University Center for High Angular Resolution Astronomy Array at Mount Wilson Observatory. The CHARA Array is supported by the National Science Foundation under Grant No. AST-1636624 and AST-2034336. Institutional support has been provided from the GSU College of Arts and Sciences and the GSU Office of the Vice President for Research and Economic Development. Time at the CHARA Array was granted through the NOIRLab community access program (NOIRLab 2018A-0158, 2018B-0123, 2019A-0080; PI: C. Lanthermann). S.K., N.A., and C.L.D. acknowledge support from an ERC Starting Grant (Grant Agreement No. 639889), ERC Consolidator Grant (Grant Agreement ID 101003096), and STFC Consolidated Grant (ST/V000721/1). A.L. received funding from STFC studentship No. 630008203. J.D.M. acknowledges funding for the development of MIRC-X (NASA-XRP NNX16AD43G, NSF-AST 1909165) and MYSTIC (NSF-ATI 1506540, NSF-AST 1909165). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. DLV-772225: MULTIPLES). This research has made use of the Jean-Marie Mariotti Center Aspro service (available at http://www.jmmc.fr/aspro). This work was supported by the Programme National de Physique Stellaire (PNPS) of CNRS/INSU co-funded by CEA and CNES. This work has been partially supported by the LabEx FOCUS ANR-11-LABX-0013. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

Appendix A Calibrators

In this section, we list the stars that we used as calibrators in data reduction in Table A.1. In addition, we observed eight targets to use them as calibrators but ended up showing resolved features such as nonconstant closure phases, or obvious multiplicity signal. These targets are 2MASS J21472749+5746367, HD 47415, HD 167882, HD 194479, HD 237112, HD 239737, HD 281311, and SAO 49700: they will be reported and added into the bad calibrator list⁶ of the JMMC.

Appendix B Coherent time used in the data reduction

In this section, we display the number of frames that we coherently add to optimize the quality of the data in Table B.1.

In Figure B.1, we show an example of S/N of V2 as a function of the number of frames added coherently, for HD 210809, for the date 2019-06-26.

In Figure B.2, we show an example of the error on the CP as a function of the number of frames added coherently, for HD 210809, for the date 2019-06-26.

Fig. B.1 V2 S/N as a function of the number of frames added coherently, for HD 210809, for the different baselines.

Fig. B.2 CP error as a function of the number of frames added coherently, for HD 210809, for the triplets of telescopes.

References

All Tables

All Figures

	Fig. 1 Cumulative histogram of the number of O-type stars as a function of their magnitude in the J-band. In blue, for the whole sky. In orange, for a declination greater than –20º to be reachable with the CHARA Array. Data from the Galactic O-star catalog (GOSC; Maíz Apellániz et al. 2013).
In the text

	Fig. 2 Histogram of the H-band magnitude of our observed sample in green, the original sample planned for this pilot survey in orange, and the planned sample for the large survey in blue.
In the text

Fig. 3 Example of change in the uv-plane modeled with the ASPRO2 software. Left: global view of the uv-plane. Different colors are for different pairs of telescopes (baselines). Right: zoom on the change in the uv-plane for the fastest moving baseline. Each straight black line indicates the uv-coordinate for this baseline for an instant snapshot, and each snapshot is separated by 5 min. The regularly spaced, orange lines display the V2 values modeled for a binary with a separation of 50 mas.

In the text

	Fig. 4 Histogram of detected companions as a function of the estimated mass ratio q, when taking into account only the first detected companions (blue), and when adding the second companions detected (orange).
In the text

	Fig. 5 Histogram of detected companions as a function of the projected separation in Astronomical Unit (AU). The distribution in blue takes into account only the first companions. The distribution in orange considers all companions.
In the text

	Fig. 6 Estimated mass ratio as a function of the projected separation in AU of the detected companions. We note that some error bars are hidden by the size of the markers. Error bars on separation ending with an arrow pointing right are to distinguish the companions detected close to the OWA. The black points are for known companions detected with other observational techniques.
In the text

	Fig. 7 Estimated photo-center shift μ in mas as a function of the estimated orbital period in years. The red lines show the limit of capability of Gaia to measure the orbital parameters of a system, with the hatched area being not accessible by Gaia.
In the text

	Fig. B.1 V2 S/N as a function of the number of frames added coherently, for HD 210809, for the different baselines.
In the text

	Fig. B.2 CP error as a function of the number of frames added coherently, for HD 210809, for the triplets of telescopes.
In the text