© ESO, 2017

1. Introduction

Massive stars are key components of galaxies and, despite their importance in modern astrophysics, they remain incompletely understood. Direct measurements of some critical parameters remain scarce (Gies 2003, 2012). In particular, masse estimates with an accuracy better than 5% are still needed to challenge stellar evolution models (Schneider et al. 2014).

So far, double-lined spectroscopic eclipsing binaries (SB2E) offer the most reliable constraints on the absolute masses. In the best cases, the formal precision on the masses can be as low as 1%. Unfortunately, because of the decreasing alignment probability with increasing separation, the few known SB2E systems only give access to tight, short period binaries (P< 10 d, typically). The complex influence of tidal forces on the internal structure is not negligible. Additionally, it is sometimes difficult to strictly exclude a prior mass-transfer event that could modify the amount of mass in each individual star. This is especially true because the evolutionary status of massive stars is not always easy to determine. These phenomena add significant degrees of uncertainty to the difficult task of calibrating evolutionary tracks from short period binaries (Schneider et al. 2014; Almeida et al. 2015).

In this context, being able to measure the true masses of longer-period, non-eclipsing systems could bring two fundamental advantages: (1) they alleviate the limitations imposed by the rareness of eclipsing systems; and (2) they avoid the uncertainties related to the effects of tides and mass-transfer given the significantly larger physical separations. Long-period binary systems can be characterised through radial-velocity up to periods of about ten years (Sana & Evans 2011). However the inclination of their orbital plane with respect to the line of sight needs to be obtained by other means. Optical long baseline interferometry (OLBI) has long been recognised as a promising solution (Boyajian et al. 2007; Kraus et al. 2009), although hampered by sensitivity and efficiency issues.

In recent years, the OLBI techniques underwent major improvements in performances, thanks to infrastructure and instrumental evolutions. From a handful of O-type stars (e.g. γ²Vel, θ¹ Ori C, η Car), the accessible sample increased to more than one hundred (Sana et al. 2014). Many massive binaries with periods in the range 100 to 5000 days at 1 to 2 kpc (a typical distance for nearby O stars) are now accessible by both radial-velocity studies and by the Very Large Telescope Interferometer (VLTI). A quantitative illustration of the overlap and complementarity between the techniques can be found in Fig. 1 from Sana & Le Bouquin (2010). This new situation has triggered a flurry of results in the field of massive stars (Millour et al. 2009; Sana et al. 2011, 2013; Sanchez-Bermudez et al. 2013; De Becker et al. 2012; Sanchez-Bermudez et al. 2014).

In this context, we initiated a follow-up of long-period massive binaries that could be spatially resolved by optical interferometry. Most of our targets were selected from the catalogue of Sana et al. (2012) with proven or suspected spectroscopic double-lined signature (SB2). By combining the SB2 radial velocity amplitudes with the size and inclination of the relative orbital motion on the sky, it is possible to determine the individual masses of each component and the distances to the systems (though see difficulties below). Our final aim is to confront evolutionary models and existing theoretical calibration laws such as that in Martins et al. (2005). Additionally, a complete knowledge of the orbital elements is of particular interest for specific objects, such as non-thermal radio emitters and hierarchical triples. In the long run, obtaining a good sample of the latter systems may indeed help understanding why over one third of massive systems are actually formed by three or more stars (Sana et al. 2014), and to constrain star formation theories.

This paper presents the observational results of the interferometric follow-up, and discusses the consequences with respect to the final goal of calibrating models. The second section details the sample, observations and analyses. The third section presents the results for each individual object. The fourth section discusses the potential of these results to provide reliable masses, with a special focus on perspectives offered by Gaia. The paper ends with brief conclusions.

2. Observations and analyses

2.1. The sample

The sample of stars is presented in Table 1. We focus on long-period massive binaries with proved or suspected SB2 signature. The sample is restricted to southern targets with magnitude H< 8 to match the location and sensitivity limit of the PIONIER instrument with the 1.8 m Auxiliary Telescopes. We add three non-thermal radio emitters that were recently resolved as long-period binaries.

The initial sample also included HD 152234, an O9.7I with a known spectroscopic period of ≈ 125 days (Sana et al. 2012). However the spectroscopic companion was only marginally resolved in our first interferometric observations (Sana et al. 2014). The obtained accuracy with the current infrastructure would not allow us to reach our science goal for this system, and so we stopped interferometric monitoring of this object.

The individual masses listed in Table 1 are computed from the spectral types using the calibration from Martins et al. (2005). Caution should be applied when using these masses, since their uncertainty is as high as 35 to 50%. The reference for the spectral types (and/or directly the masses) are given in Col. 4. The supposed association and corresponding distances are reported in Table 1, with references in Col. 7.

2.2. Observations and data reduction

Most interferometric data were obtained with the PIONIER¹ combiner (Le Bouquin et al. 2011) and the four 1.8 m Auxiliary Telescopes of the VLTI (Haguenauer et al. 2010). To calibrate the fringe visibilities and closure phases, we interleaved observations of KIII stars between those of the O-type stars. These reference stars were found with the SearchCal² software (Bonneau et al. 2011). These K giants have apparent diameters smaller than 0.5 mas therefore providing calibration with 5% accuracy even on the longest baselines of VLTI. Data were reduced and calibrated with the pndrs³ package (Le Bouquin et al. 2011). Each observing block (OB) provides five consecutive files, each containing six visibilities and four closure phases dispersed over three spectral channels across the H band. The dataset is complemented by about ten observations performed with the AMBER combiner and three 8 m Unit Telescopes, mostly before 2011. These data were reduced and calibrated with the amdlib-3 package (Tatulli et al. 2007; Chelli et al. 2009). Finally, two 2016 observations were obtained during the science verification of the GRAVITY instrument (Eisenhauer et al. 2011). They were reduced and calibrated by the standard pipeline (Lapeyrere et al. 2014). The description of the GRAVITY instrument and its first results will be presented in detail in a forthcoming paper from the GRAVITY Collaboration (corresponding author: F. Eisenhauer, e-mail: eisenhau@mpe.mpg.de).

2.3. Resolved astrometry fitting

The interferometric observations were adjusted with a simple binary model composed of two unresolved stars. This is a valid assumption because the expected diameters of the individual components (< 0.2 mas) are unresolved, even with the longest VLTI baselines. The free parameters are the two coordinates of the apparent separation vector (East and North, in milliarcsecond) and the flux ratio f_H between the secondary and the primary. The latter is considered constant over the H band and over the different epochs. Given the long orbital periods of these binaries (P ≈ 0.5 − 15 yr), we neglected the orbital motion in the course of a single night.

We here report ≈130 observations, a tenth which have already been presented in Sana et al. (2014). The journal of the observations and the astrometric measurements are available at the CDS. It includes the best-fit binary position and its associated uncertainty ellipse for each modified Julian Date (MJD). The total execution time for a single resolved astrometric observation is about 30 min, including calibrations.

2.4. Orbit determination

We use a Levenberg-Marquardt method to find the best orbital solutions. A more complex search is not necessary because the temporal sampling is dense enough for the period to be determined unambiguously. We adjust simultaneously the resolved astrometric positions and the radial velocities, when available. This is similar to what is presented in Le Bouquin et al. (2013). There are ten independent parameters to fit, namely the orbital period P, the time of periastron passage T, the size of the apparent orbit a, the orbital eccentricity e, the inclination i, the longitude of ascending node Ω, the argument of periastron of the secondary ω, the semi-amplitudes K_a and K_b of the primary and secondary radial velocity curves, and the systemic heliocentric offset velocity g.

The ascending node is the node where the motion is directed away from the observer, that is when the passing star has a positive radial velocity. The argument of periastron is the angle between the ascending node and the periastron of the secondary, measured in the plane of the true orbit and in the direction of the motion of the secondary. The longitude of ascending node, measured in the plane of the sky, is counted positively from North to East (Ω = 90 deg). An inclination larger than 90 deg indicates an astrometric orbit covered clockwise.

We suppose the primary component for radial velocity (most massive) corresponds to the astrometric primary (brightest in H band). When only SB1 radial velocities are available, we assume that they represent the radial velocities of the primary. When no radial velocities are available, we enforce 0 < Ω < 180 deg (using Ω + 180 deg and ω + 180 deg gives the same fit to the astrometric orbit alone).

2.5. Uncertainties

The uncertainties are propagated numerically. A series of adjustments is performed while the data are corrupted by a random additive noise corresponding to their uncertainties. The quoted uncertainty on the orbital parameter is the root mean square of these various realisations. This method also allows us to explore the possible correlation between the parameters. A systematic uncertainty of 1.5% should be added on the size of the apparent orbit a. It comes from the calibration of the effective wavelength of the PIONIER instrument that was used to obtain almost all resolved astrometric points. This translates into a “plate scale error” through the interferometric equation.

2.6. Physical quantities

Once the orbital parameters are known, additional physical quantities can be derived. We define two strategies: one fully independent of the radial velocities and one fully independent of the adopted distance. In the first approach, we impose the distance d from external constraints such as membership to a known cluster or association (see Table 1). It allows to compute the physical size of the orbit a_p = da. The total mass M_t of the system is then derived with Kepler’s third law: ${\mathit{M}}_{\mathrm{t}}\mathrm{=}\frac{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}{\mathit{a}}_{\mathrm{p}}^{\mathrm{3}}}{\mathit{G}\hspace{0.17em}{\mathit{P}}^{\mathrm{2}}}\mathit{,}$(1)where G is the gravitational constant. In the second approach, we use the radial velocity amplitudes, when available, to determine the physical size of the orbit, ${\mathit{a}}_{\mathrm{p}}\mathrm{=}\frac{\mathit{P}\hspace{0.17em}\mathrm{(}{\mathit{K}}_{\mathit{a}}\mathrm{+}{\mathit{K}}_{\mathit{b}}\mathrm{)}}{\mathrm{2}\mathit{\pi }}\hspace{0.17em}\frac{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}}{\sin \mathrm{\left(}\mathit{i}\mathrm{\right)}}$(2)and thus the total mass from Eq. (1). The distance of the system is d = a_p/a and the individual masses are given by: ${\mathit{M}}_{\mathit{a}}\mathrm{=}\frac{{\mathit{K}}_{\mathit{b}}}{{\mathit{K}}_{\mathit{a}}\mathrm{+}{\mathit{K}}_{\mathit{b}}}\hspace{0.17em}{\mathit{M}}_{\mathrm{t}}\mathrm{;}{\mathit{M}}_{\mathit{b}}\mathrm{=}\frac{{\mathit{K}}_{\mathit{a}}}{{\mathit{K}}_{\mathit{a}}\mathrm{+}{\mathit{K}}_{\mathit{b}}}\hspace{0.17em}{\mathit{M}}_{\mathrm{t}}\mathit{.}$(3)

3. Results

The results for each binary of the sample are displayed in Appendices A and B. Here we discuss the specificities of each system.

3.1. HD 54662

This massive binary is a known runaway star producing a bow shock as its winds interact with the interstellar medium (Peri et al. 2012). It was first recognised as a promising target for mass determination by Boyajian et al. (2007). The three-dimensional orbit is now fully constrained by the resolved astrometric observations. However the orbital period of almost six years is completely different from the preliminary estimation from radial velocities (557 d). It illustrates that resolved astrometry is less prone to phase ambiguities. We have attributed the historical SB1 radial velocities to the primary as they closely match its SB2 amplitude. We did so to confirm that the historical spectroscopic binary corresponds to the interferometric binary.

The total mass is in agreement with expectations from spectral types. However the individual masses are still uncertain because of the lack of accurate SB2 radial velocity measurements. The preliminary estimates for K_a and K_b from Boyajian et al. (2007) also lead to unrealistic masses, and should then be considered with caution.

Finally, the eccentricity is rather small considering the period of six years.

3.2. HD 93250

This target is a non-thermal radio emitter that was first spatially resolved as a binary by Sana et al. (2011). Here we present additional observations that fully constrain its astrometric orbit. The projected semi-major axis of 1.2 mas is at the limit of the spatial resolution of the VLTI and some measurements have large relative errors. Accordingly, it has the shortest period of our sample (P = 194.3 d). The orbital period is still long enough to avoid tidal interaction and circularisation. The total mass of the system points toward a pair of similar high mass stars (O4+O4), as do the measured flux ratio near unity. The system is also known for the presence of non-thermal radio emission (colliding winds).

There is no radial velocity variation reported for this star (Rauw et al. 2009), but given its short period it was worth obtaining a full orbit to guess the expected signal. However, because of the small inclination, the expected radial velocity variations due to the orbital motions are about 20 km s^-1. It is probably impossible to disentangle this from the intrinsic line width, as first suggested in Sana et al. (2011). It illustrates well the detection biases affecting binaries with (quasi) identical components.

While HD 93250 is inadequate for direct mass determination from resolved astrometry and radial velocities, constraints on the orbital parameters are valuable to properly interpret the non-thermal radio emission. A significant eccentricity was also required to explain the variations observed in the X-ray flux and is confirmed by this study (e = 0.22 ± 0.01).

3.3. HD 150136

This is a known hierarchical triple system and one of the most massive stellar systems known (Mahy et al. 2012). The first spatially resolved observations of the outer pair were discussed by Sana et al. (2013) and Sanchez-Bermudez et al. (2013). Here we present additional observations and an improved orbit. The new orbit has a slightly longer period (3065 d vs. 3008 d), slightly lower eccentricity (0.68 vs. 0.73) and slightly different absolute masses. All stellar and orbital parameters are within 1σ of the previously published values. The accuracies are, however, improved. In Table B.3, component a corresponds to the inner 2.7 d binary while component b is the O6.5-7 V outer companion.

The three-dimensional orbit of the 8.4 year outer system is now uniquely constrained thanks to our long term follow-up. Still, it would benefit from few additional spatially resolved observations in the coming years to further reduce the uncertainties on the absolute masses.

Regarding the inner, short period binary: accounting for the minimal masses from Mahy et al. (2012) and our new total mass of 88 M_⊙, we obtain an inclination of i = 53° and thus individual masses of 54 and 34 M_⊙. Again, these values are at the lower bound of the 1σ confidence interval from the preliminary study of Sana et al. (2013).

The distance estimate from the combined fit favours the new distance to the NGC 6193 cluster (1.15 kpc) from Kharchenko et al. (2005) compared to the older one (1.32 kpc) from Herbst & Havlen (1977).

3.4. HD 152233

The well constrained three-dimensional orbit reveals an almost edge-on eccentric binary. The measured properties correspond with those of the long period SB2 binary discussed by Sana et al. (2008). The total mass of the system (≈ 45 M_⊙) is significantly lower than the expectation from the spectral types (38 + 30 M_⊙, Table 1). Together with the small flux ratio, it points towards a lower mass secondary, most likely of a later spectral type as we discuss further below.

The individual masses (≈ 2M_⊙) and distance (700 pc) derived from the combined astrometric and RV fit are unrealistic. They cannot be reconciled with massive stars, even considering large uncertainties. A visual inspection of the UVES spectra shows that the spectral lines are never fully disentangled by the radial velocity shifts. The radial velocities of the secondary may be systematically underestimated in the spectral analysis, thus lowering the total mass and the distance of the system. If we impose a distance of ≈ 1.5 kpc (thus somehow imposing the total mass to be ≈ 45 M_⊙), the fit of the radial velocities of the secondary yields a systematic residual of about 15 km s^-1, that is of similar amplitude as the uncertainties. The individual masses rise to 31 and 12 M_⊙. These are more realistic values, also compatible with the small flux ratio measured. We conclude that the radial velocities are probably faulty.

Given the high inclination, the closest apparent approach is 0.24 ± 0.05 mas. However there is little chance that this system could be a long-period eclipsing binary because the expected apparent stellar diameters are about 0.17 mas.

3.5. HD 152247

As seen in the pair diagram, the three-dimensional orbit is well constrained. It corresponds to the long period, highly eccentric spectroscopic binary discussed in Sana et al. (2012).

The total mass is compatible with the masses expected from spectral types. The measured flux ratio agrees with the expectation for late O-type stars with luminosity classes III and V. A small amount of accurate SB2 data should allow to obtain individual mass estimates independently from the distance. Attributing the reported SB1 radial velocities to the primary, we expect SB2 velocity semi-amplitudes of ≈ 61 km s^-1 for the secondary (fainter and less massive). Thus HD 152247 is confirmed as a good candidate for precise mass determination.

3.6. HD 152314

The three-dimensional orbit is well constrained. It corresponds to the long period spectroscopic binary discussed in Sana et al. (2012).

The total mass and the measured flux ratio are roughly compatible with the masses expected from spectral types. Attributing the reported SB1 radial velocities to the primary, we expect SB2 velocity semi-amplitudes of ≈ 30 km s^-1 for the B-type secondary.

3.7. HD 164794 (9 Sgr)

The three-dimensional orbit is well constrained. Our combined fit confirms that the resolved pair corresponds to the long-period SB2 presented in Rauw et al. (2012, 2016). We had to introduce a constant RV shift of 10.0 km s^-1 between the two components in order to achieve a correct fit of the SB2 radial velocities. This is not unusual for very massive stars and is thought to be an effect of the stellar winds. The combined fit is convincing but this result is the most puzzling of our sample.

First, placing the system at a distance of d ≈ 1.25 kpc makes the total mass compatible with the expected spectroscopic one (≈ 80 M_⊙). Thus, in the long controversy about the distance to NGC 6530, and assuming the total mass from spectroscopy is correct, our observations confirm the distance from Prisinzano et al. (2005) and Kharchenko et al. (2005). This is significantly lower than the 1.78 ± 0.08 kpc assumed by Rauw et al. (2016) and based on Sung et al. (2000).

Although it makes no assumption on the distance, the combined fit between radial velocities and the resolved astrometry does not solve the issue. The deduced individual masses and the distance are even smaller. We note that, fitting the same RV data, Rauw et al. (2016) obtain similar radial velocity amplitudes. They conclude that the masses were compatible with expectation, but assuming an inclination of 45 ± 1 deg. Looking at the resolved astrometry in Fig. A.7, it seems hard to significantly lower the inclination of our best-fit of 86 deg. In fact, such a low inclination is definitely ruled out by the VLTI observations. Thus we are left with a clear discrepancy between the masses expected from spectral types, the best-fit semi-amplitudes of the radial velocity curves and the distance to the system.

To further explore, we ran a combined fit of astrometry and radial velocity, but imposing the distance to be 1.25 kpc (thus somehow imposing the total mass to be ≈ 75 M_⊙). The agreement to the radial velocities is obviously significantly degraded, but the systematic residuals are never larger than 20 km s^-1. Altogether, we believe this is our best proxy for this system (M_a ≈ 40 M_⊙; M_b ≈ 33 M_⊙; d ≈ 1.25 kpc; i ≈ 86 deg). A careful re-analysis of the radial velocity amplitudes taking into account the possible systematic biases is required to understand the disagreement among the observations.

Finally, even if the orbit is almost edge-on, the system is unlikely to be eclipsing. The closest apparent approach is 0.39 ± 0.1 mas, for estimated stellar diameters of 0.14 mas.

3.8. HD 167971

This is a known hierarchical triple system in which we observe the orbital motion of the third component around the short period eclipsing pair. This outer pair is also the longest-period system in our sample. It was first resolved by De Becker et al. (2012). Here we provide additional interferometric observations.

Any period between 15 and 25 yr could easily match the astrometric observations alone. Our temporal coverage is still not sufficient to obtain a definite estimate of the period, and thus of any of the orbital elements. This should be easily tackled with additional interferometric points in the coming years.

Interestingly, the close eclipsing pair is intensively studied and modelled by Ibanoglu et al. (2013). By measuring the times of minimum light, they infer a period of P = 21.7 ± 0.7 yr for the outer pair. This value is compatible with the long term radio variability, interpreted as the signature of a colliding wind region between the winds of the inner binary and that of the third component (Blomme et al. 2007). We thus add these constraints (P = 21.7 ± 0.7 yr, e = 0.53 ± 0.05) in our fit of the astrometric orbit. The total mass of the system inferred from the orbit is in good agreement with the values given in Table 1, even though one has to notice the large uncertainty on this quantity quoted in Fig. A.8.

The inclination of the outer orbit (180 − i = 35 deg) is very different than that of the eclipsing pair (73 deg).

We found an average flux ratio of 0.78, similar to the one first reported in De Becker et al. (2012). Our interferometric measurements alone could not lift the ambiguity in the association of the resolved components to the close eclipsing binary and to the distant third star.

3.9. HD 168137

The three-dimensional apparent orbit corresponds to the long period binary discussed from radial velocities by Sana et al. (2009). The preliminary orbit has an impressive eccentricity of 0.90. We checked various ways to explain this eccentricity with biases in the observations but could not find any plausible explanations. Additional observations in the coming years will firmly establish this peculiar eccentricity, and lift any remaining degeneracies. The next periastron passage is expected in July 2020.

The total mass is unconstrained because of the still poor precision on the size of the orbit (± 50%). Again, a few additional interferometric observations will easily tackle this issue. Moreover the few existing SB2 radial velocities are very difficult to model, except for the well separated spectra. The other measurements might be a mixture of the two components, thus cancelling each other out. This is especially possible because the two stars have similar brightness (f_H = 0.77) and similar line depths and widths (see Fig. 5 in Sana et al. 2009). We have ignored these radial velocities for the fit. The radial velocity elements are consequently degenerated, and the individual masses are unconstrained.

The most interesting result of this system is its huge eccentricity that could potentially make it a long-period eclipsing system. The expected apparent closest approach is 0.15 ± 0.13 mas. This value thus has some probability of being lower than 0.1 mas, while the stellar diameters should be of this size. Altogether this system deserves additional monitoring to better determine its inclination, eccentricity and time of closest approach (so far known within ± 200 days).

3.10. CPD-47 2963

This non-thermal radio emitter was first spatially resolved as a binary in the course of the SMASH+ survey (Sana et al. 2014). Here we present additional observations and the first orbit. The three-dimensional orbit is constrained but would benefit from few additional interferometric observations to lift the remaining degeneracy between Ω and ω and between a and i. The total mass is at the lower bound of the expectation from the spectral type, although still uncertain.

CPD-47 2963 has been reported as SB1 with a preliminary period of 59 days and a semi-amplitude of 9 km s^-1 (OWN survey pre-publication, Sota et al. 2014). Considering the eccentricity and the period of the interferometric binary (e = 0.66, P = 655 d), a stable hierarchical system seems unlikely according to the criterion quoted by Tokovinin (2004): ${\mathit{P}}_{\mathrm{out}}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{out}}{\mathrm{)}}^{\mathrm{3}}\mathit{>}\mathrm{5}\hspace{0.17em}{\mathit{P}}_{\mathrm{in}}\mathit{.}$(4)It is more plausible that all observations refer to the same pair.

4. Discussion

4.1. Orbit of non-thermal radio emitters

Non-thermal (synchrotron) radio emitters belong to the category of the so-called particle accelerating colliding-wind binaries (PACWBs), a first catalogue of which is published by De Becker & Raucq (2013). These objects are of great interest for understanding non-thermal physics. Particle acceleration arises from the strong hydrodynamic shocks in the colliding-wind region, providing an opportunity to investigate particle acceleration in environments different from, for example, supernova remnants.

PACWBs are intrinsically variable sources on a timescale corresponding to the orbital period, with an observed variable behaviour depending on the orientation of the orbit. The determination of their orbital elements is thus especially important.

Short period (i.e. a few days) systems are not suitable for the investigation of non-thermal processes. First of all, the acceleration of electrons in short period binaries is expected to be severely inhibited due to the highly efficient cooling by inverse Compton scattering of photospheric photons (De Becker 2007). Second, the large dimensions of the radio-sphere (up to several hundreds of R_⊙, see e.g. Mahy et al. 2012) also prevent any potential synchrotron radio emission produced in the colliding-wind region to escape if the two stars are too close to each other. Long-period systems are thus better suited, but a proper knowledge of their orbit is often lacking. Those systems are challenging to follow-up with radial velocities due to large amounts of time needed and to the small amplitudes of the radial velocity variations.

Our sample includes one relatively short-period (HD 93250) and three long-period non-thermal radio emitters. HD 164794 and CPD-47 2963 have their three-dimensional orbits fully mapped. HD 167971 will be completed soon. The updated number of PACWBs with at least partly determined orbits is now 14, compared to the 11 initially mentioned by De Becker & Raucq (2013). The full determination of their orbits is necessary to define adequate strategies for future high angular resolution radio observations using very long baseline interferometry (VLBI). It permits to select adequate orbital phases for observation, and allows for an adequate interpretation of radio images of the synchrotron emission region afterwards. Finally, potential future observations of PACWBs in gamma-rays will require a careful selection of the orbital phases as well.

4.2. Periods, inclinations and eccentricities

Eight binaries of the sample fall within the “interferometric gap” defined by Mason et al. (1998) where the periods are too long to be easily catched in radial velocity snapshot surveys but the apparent separations are too small to be resolved by imaging surveys. However most of them still have detectable radial velocity variations (except HD 93250 which remained hidden because of its small inclination and nearly twin components).

Our target selection is biased toward large inclination because we favoured systems with detected SB2 lines (at least tentatively). Accordingly, the seven binaries with reported radial velocity variations have inclinations less than 15 deg from edge-on. However, given the long periods involved, only one of them has some chance of being eclipsing.

Figure 1 shows the period-eccentricity diagram. Nine of the ten binaries in our sample have a well-determined eccentricity, and this parameter is partially constrained for the remaining one. It is known that the radial velocity technique requires a good sampling near the periastron passage in order to constrain orbits with high eccentricity. Resolved astrometry is more robust to high eccentricity because it measures the integrated orbital motion, and not its derivative.

It is noteworthy that the eccentricities are distributed over a large range. For example, there are two systems with periods near 2000 days: one of them has the highest eccentricity of the sample (0.90) while the other has the lowest eccentricity (0.06). Interestingly, these two objects fall in parts of the diagram that were left empty before. They completely break the correlation between period and eccentricity tentatively observed in previous studies (see open symbols in Fig. 1). This finding suggests that the current lack of long-period systems with low (e< 0.25) or high (e> 0.75) eccentricities results from observational biases. In the former case, the small radial velocity variations are easily hidden by the typical rotational broadening of the stellar lines in massive stars (see the discussion in Sana et al. 2011). In the later case, the radial velocity variations potentially occur over a small fraction of the orbit. The difficulty of catching the periastron passage and the long timescales between two periastron passages make these systems challenging to detect and characterise.

The fact that such a variety of eccentricities co-exist indicates either that massive star formation produces a wide range of systems, or that several binary formation mechanisms are at play.

4.3. Increasing the interferometric sample

Recently, Sana et al. (2014) resolved about 20 additional systems within a similar range of separations. The expected periods should be similar, and a long-term interferometric follow-up has been initiated. Another interesting target is Herschel 36 (Aa-Ab), which has a period of 500 d and was recently spatially resolved by VLTI as well (Sanchez-Bermudez et al. 2014).

It would be interesting to reach systems with slightly shorter periods, because they are more frequent and have larger radial velocity amplitudes. At VLTI, the limitation comes from the length of the interferometric baselines (150 m, 1 mas at H band). The next generation GRAVITY/MATISSE instruments will not change the game, but simply participating in this follow-up initiated with PIONIER. The only way to improve for VLTI is to use shorter wavelengths.

The Center for High Angular Resolution Array (CHARA) has significantly longer baselines (330 m, 0.5 mas at H band, Ten Brummelaar et al. 2008). However, efficient follow-up of binaries are nowadays limited to stars brighter than 5 mag. This should improve in the coming years with the arrival of new detectors and the installation of adaptive optics (Ten Brummelaar et al. 2016).

4.4. Preliminary masses

The distances to the systems are known with a precision of 5 to 15%. For most binaries, it dominates the error budget of Eq. (1). The total masses are thus known with a precision of 15 to 50%. Only HD 167971 and HD 168137 have the precision on their total mass limited by the knowledge of the period and/or size of the orbit. For these two stars, this will be tackled in the coming years with some additional interferometric observations at VLTI.

The preliminary masses estimated from the SB2 radial velocities are puzzling. One system has realistic masses (HD 150136), one has masses obviously overestimated (HD 54662) and two have masses unrealistically low (HD 152233 and HD 164794). We speculate that the measured SB2 radial velocities are affected by systematic biases, phased with the orbital motion. The origin of such biases can be twofold:

• An astrophysical effect makes the line-forming re-gion appear differently depending on the orbital phase.A known example is the Struve-Sahade effect (Linderet al. 2007; Palate& Rauw 2012). However this hypothesisis hard to support because of the very large physical sepa-rations of the binaries in our sample. Another effect couldbe related to the wind-wind collision in eccentric bina-ries. In HD 164794 for instance, Rauwet al. (2012) described that var-ious primary lines might behave differently around periastronpassage. Depending on the relative weighting of the lines in thefinal orbital solution, this could bias the estimation of the radialvelocity amplitude that critically depends on the measurementsat periastron.

• The spectral modelling is corrupted by the entanglement of the lines of each component, which varies according to the orbital phase. This is obviously the case when the true line profiles of each component are incorrectly modelled. This second hypothesis is further supported by the fact that the only system for which the masses are found to be in agreement is also the only SB3. The orbital motion of the inner pair introduces large radial velocity shifts (>200 km s^-1). This is well illustrated in Fig. 1 from Mahy et al. (2012). These enhanced shifts permit them to disentangle the line profile of the third component, and thus to model it adequately in most blended spectra.

4.5. Perspectives for accurate SB2 velocities

The limiting factor for mass determination of these long-period binaries is the accuracy of the SB2 radial velocity amplitudes (Sect. 4.4), and especially their sensitivity to the adopted line profiles. A possible mitigation is to directly fit the entire set of spectra by a combination of stellar line profiles, shifted accordingly to the known orbital parameters. The three remaining free parameters, namely K_a, K_b and g, can be estimated without extracting radial velocities for each phase (note that g may have to be be left free for each individual line). This approach makes optimal use of the information included in the blended spectra. Varying the stellar line profile parameters allows to investigate the final mass accuracies, including systematics.

To reduce this uncertainty, one should input the best possible individual line profiles. Future spectroscopic observations should focus on the predicted phases of maximal velocity, to provide the best possible disentangled spectra. Table 2 summarises the expected maximal velocity shifts. Five binaries have predicted shifts of the order of 100 km s^-1 or larger, thus enabling decent disentangling and model-independent masses. Disentangling will be much more difficult for the five other systems.

4.6. Perspectives with Gaia

At the considered range of distances, Gaia will deliver exquisite parallax accuracy (<2%), especially because the binary effect can be adequately subtracted thanks to the knowledge of the astrometric orbital parameters presented in this study. Combined with the size of the apparent orbits and the periods, we expect the total masses to be measured with an uncertainty of 5%, following Eq. (1). For most systems, the error budget will be equally shared between the interferometric resolved astrometry and the Gaia distance. These measurements of the total masses will be completely independent from radial velocities.

Gaia will also constrain the orbital photocentre displacement μ_G. This quantity can be extracted even if the orbital period is poorly sampled by Gaia because of the knowledge of all orbital parameters. It permits calculating the mass ratio q: $\mathit{q}\mathrm{=}\frac{{\mathit{\mu }}_{\mathit{G}}{\mathit{f}}_{\mathit{G}}\mathrm{+}{\mathit{\mu }}_{\mathit{G}}\mathrm{+}\mathit{a}{\mathit{f}}_{\mathit{G}}}{\mathit{a}\mathrm{-}{\mathit{\mu }}_{\mathit{G}}{\mathit{f}}_{\mathit{G}}\mathrm{-}{\mathit{\mu }}_{\mathit{G}}}$(5)where a is the (known) size of the apparent orbit and f_G the flux ratio in the Gaia G band. Interestingly, f_G is well approximated by the (known) flux ratio in the H band f_H because hot stars have their optical and near-infrared spectra in the Rayleigh-Jeans regime.

In order to test feasibility, we computed the expected photocentre displacement μ_G from the orbits derived in this study. This was done by inverting Eq. (5), guessing the expected mass ratio. We used two approaches, both relying on the calibrations from Martins et al. (2005). We first converted the observed flux ratio f_H into mass ratio. For a pair of main-sequence stars, the calibration is well approximated by $\mathit{q}\mathrm{=}{\mathit{f}}_{\mathit{H}}^{\mathrm{0.7}}$. The exponent is even lowered for a massive, evolved primary with a less-massive main sequence secondary. In a second approach, we simply use the expected mass ratio from the individual spectral types when available (see Table 1). The computed μ_G are reported in Fig. 2. The predicted displacements range from 0.1 mas to 1 mas. This is a factor ten to a hundred larger than the Gaia performances (assuming the orbital photocentre shift is measured with a similar accuracy than the parallaxes).

Accounting for 5% uncertainty on the total mass (shared between the interferometric resolved astrometry and the Gaia distance, see Eq. (1)), and for 1 to 10% uncertainty on the mass-ratio (shared between the interferometric flux ratio and the photocentre shift amplitude from Gaia, see Eq. (5)), we expect to measure the individual masses with an uncertainty of 5 − 15%, independently of radial velocities.

5. Conclusion

In this paper we report ≈130 observations of spatially resolved astrometry on ten massive O-type binaries, using OLBI. This dataset is used to derive the elements of the apparent orbits. Combined with the distance, it provides the total mass of the system. We also compute preliminary individual component masses for the four systems with available SB2 radial velocities. The conclusions are the following:

• Nine over ten binaries have their three-dimensional orbits wellconstrained. One requires additional OLBI observations to liftremaining degeneracies.

• We provide orbital elements for four colliding wind, non-thermal radio emitters. They constitute valuable targets for future high-angular resolution radio imaging.

• The large range of eccentricity suggests that the current lack of long-period systems with low (e < 0.25) or high (e > 0.75) eccentricities results from observational biases.

• For the most-studied systems, we find a clear conflict between the astrometric orbit and the SB2 radial velocity amplitudes. We speculate that the SB2 amplitudes are affected by systematic biases, that could represent an intrinsic limitation for estimating dynamical masses from OLBI+SB2 or Gaia+SB2.

• Our results can be combined with future Gaia parallaxes and orbital photocentre displacements to measure the masses of the individual components with an accuracy of 5 to 15%, independently of the radial velocities.

---

Acknowledgments

PIONIER is funded by the Université Joseph Fourier (now Université Grenoble Alpes, UGA), the Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), the Agence Nationale pour la Recherche (ANR-06-BLAN-0421, ANR-10-BLAN-0505, ANR-10-LABX56, ANR-11-LABX-13), and the Institut National des Sciences de l’Univers (INSU PNP and PNPS). Its integrated optics beam combiner is the result of a collaboration between IPAG and CEA-LETI based on CNES R&T funding. This work is based on observations made with the ESO telescopes. It made use of the Smithsonian/NASA Astrophysics Data System (ADS) and of the Centre de Données astronomiques de Strasbourg (CDS). All calculations and graphics were performed with the freeware Yorick. Some of the observations were acquired in the framework of the Belgium Guaranteed Time of VLTI (VISA project). The authors warmly thank everyone involved in the VLTI project. The authors warmly thank the entire GRAVITY consortium. GRAVITY is developed in a collaboration by the Max Planck Institute for Extraterrestrial Physics, LESIA of Paris Observatory and IPAG of Université Grenoble Alpes/CNRS, Max Planck Institute for Astronomy, the University of Cologne, the Centro Multidisciplinar de Astrofísica Lisbon and Porto, and the European Southern Observatory.

References

Appendix A: Individual results

Appendix B: Additional tables

All Tables

All Figures