The VLTFLAMES Tarantula Survey ^{⋆,}^{⋆⋆,}^{⋆⋆⋆}
VIII. Multiplicity properties of the Otype star population
^{1}
Astronomical Institute Anton Pannekoek, Amsterdam University,
Science Park 904,
1098 XH, Amsterdam, The
Netherlands
email:
h.sana@uva.nl
^{2}
Utrecht University, Princetonplein 5, 3584CC
Utrecht, The
Netherlands
^{3}
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD
21218,
USA
^{4}
Johns Hopkins University, 11100 Johns Hopkins Rd Laurel, Baltimore, MD, USA
^{5}
Astrophysics Research Centre, School of Mathematics and Physics,
Queens University Belfast, Belfast
BT7 1NN,
UK
^{6}
UK Astronomy Technology Centre, Royal Observatory Edinburgh,
Blackford Hill,
Edinburgh, EH9 3HJ, UK
^{7}
Scottish Universities Physics Alliance (SUPA), Institute for
Astronomy, University of Edinburgh, Royal Observatory Edinburgh,
Blackford Hill,
Edinburgh, EH9 3HJ, UK
^{8}
Instituto de Astrofísica de AndalucíaCSIC,
Glorieta de la Astronomía
s/n, 18008
Granada,
Spain
^{9}
Department of Physics and Astronomy, The Open
University, Walton
Hall, Milton
Keynes, MK7 6AA,
UK
^{10} Dept. of Physics & Astronomy, Hounsfield Road,
University of Sheffield, S3 7RH, UK
^{11}
Instituto de Astrofísica de Canarias, C/ Vía Láctea s/n, 38200 La Laguna,
Tenerife,
Spain
^{12}
Departamento de Astrofísica, Universidad de La
Laguna, Avda. Astrofísico Francisco
Sánchez s/n, 38071La Laguna, Tenerife, Spain
^{13}
Institute of Astronomy, University of Cambridge,
Madingley Road,
Cambridge, CB3 0HA, UK
^{14}
ArgelanderInstitut für Astronomie, Universität Bonn,
Auf dem Hügel 71,
53121
Bonn,
Germany
^{15}
ESA/STScI, 3700 San Martin Drive, Baltimore, MD
21218,
USA
^{16}
Armagh Observatory, College Hill, Armagh, BT61 9DG, Northern Ireland,
UK
Received:
17
May
2012
Accepted:
17
September
2012
Context. The Tarantula Nebula in the Large Magellanic Cloud is our closest view of a starburst region and is the ideal environment to investigate important questions regarding the formation, evolution and final fate of the most massive stars.
Aims. We analyze the multiplicity properties of the massive Otype star population observed through multiepoch spectroscopy in the framework of the VLTFLAMES Tarantula Survey. With 360 Otype stars, this is the largest homogeneous sample of massive stars analyzed to date.
Methods. We use multiepoch spectroscopy and variability analysis to identify spectroscopic binaries. We also use a MonteCarlo method to correct for observational biases. By modeling simultaneously the observed binary fraction, the distributions of the amplitudes of the radial velocity variations and the distribution of the time scales of these variations, we constrain the intrinsic current binary fraction and period and massratio distributions.
Results. We observe a spectroscopic binary fraction of 0.35 ± 0.03, which corresponds to the fraction of objects displaying statistically significant radial velocity variations with an amplitude of at least 20 km s^{1}. We compute the intrinsic binary fraction to be 0.51 ± 0.04. We adopt powerlaws to describe the intrinsic period and massratio distributions: f(log _{10}P/d) ~ (log _{10}P/d)^{π} (with log _{10}P/d in the range 0.15−3.5) and f(q) ~ q^{κ} with 0.1 ≤ q = M_{2}/M_{1} ≤ 1.0. The powerlaw indexes that best reproduce the observed quantities are π = −0.45 ± 0.30 and κ = −1.0 ± 0.4. The period distribution that we obtain thus favours shorter period systems compared to an Öpik law (π = 0). The mass ratio distribution is slightly skewed towards low mass ratio systems but remains incompatible with a random sampling of a classical mass function (κ = −2.35). The binary fraction seems mostly uniform across the field of view and independent of the spectral types and luminosity classes. The binary fraction in the outer region of the field of view (r > 7.8′, i.e. ≈117 pc) and among the O9.7 I/II objects are however significantly lower than expected from statistical fluctuations. The observed and intrinsic binary fractions are also lower for the faintest objects in our sample (K_{s} > 15.5 mag), which results from observational effects and the fact that our O star sample is not magnitudelimited but is defined by a spectraltype cutoff. We also conclude that magnitudelimited investigations are biased towards larger binary fractions.
Conclusions. Using the multiplicity properties of the O stars in the Tarantula region and simple evolutionary considerations, we estimate that over 50% of the current O star population will exchange mass with its companion within a binary system. This shows that binary interaction is greatly affecting the evolution and fate of massive stars, and must be taken into account to correctly interpret unresolved populations of massive stars.
Key words: stars: earlytype / stars: massive / binaries: spectroscopic / open clusters and associations: individual: 30 Dor / binaries: close / Magellanic Clouds
Based on observations collected at the European Southern Observatory under program ID 182.D0222.
Full Tables 1–3 are only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr(130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/550/A107
Appendices are available in electronic form at http://www.aanda.org
© ESO, 2013
1. Introduction
Massive binaries are spectacular systems that may have been among the first objects that formed in the early universe (Turk et al. 2009; Stacy et al. 2012). Because the binary fraction among massive stars is high (Mason et al. 2009) and as a large fraction of them are part of short period systems (Sana & Evans 2011), an accurate knowledge of their binary properties is crucial to understand the role of massive stars as a population (Langer et al. 2008; de Mink et al. 2009). On the one hand, the fraction of binaries and their orbital parameter distributions are the tracers of massive star formation and of the early dynamical evolution of their host star cluster. On the other hand, they determine the frequency of evolved massive binary systems, including highmass Xray and double compact binaries (Sadowski et al. 2008), Type Ib/c supernovae (Yoon et al. 2010), short and longduration gammaray bursts (Podsiadlowski et al. 2010).
Among the various multiplicity properties, the most investigated one is the observed fraction of spectroscopic binaries, that sets a lower limit on the true binary fraction. Based on an extensive compilation of the literature, Mason et al. (2009) showed that 51% of the Galactic Otype stars investigated up to 2008 by multiepoch spectroscopy are spectroscopic binaries. This fraction is 56% for objects in clusters and OB associations. The spectroscopic survey of Galactic O and WN stars (Barbá et al. 2010) obtains a similar fraction: 60% of the 240 southern stars investigated indeed display significant radial velocity (RV) variations with an amplitude larger than 10 km s^{1}. Finally, studies of individual young open clusters, such as IC 1805 (De Becker et al. 2006; Hillwig et al. 2006), IC 1848 (Hillwig et al. 2006), IC 2944 (Sana et al. 2011), NGC 2244 (Mahy et al. 2009), NGC 6231 (Sana et al. 2008), NGC 6611 (Sana et al. 2009) and Tr 16 (Rauw et al. 2009) reveal observed binary fractions ranging between 30 and 60%. These cluster to cluster variations are so far compatible with random fluctuations expected from the size of the samples (Sana & Evans 2011).
The intrinsic period and massratio distributions of the Otype binaries have remained poorly constrained until recently. In a pioneering work, Kobulnicky & Fryer (2007) attempted to constrain these distributions using a sample of 900 RV observations of 32 O and 88 Btype stars in Cyg OB2 and a MonteCarlo method to correct for observational biases. Their solution was unfortunately degenerate leading the authors to fix the period distribution to an Öpik law (i.e., a flat distribution in the logarithm of the separation or, equivalently, of the period). They find a very high intrinsic binary fraction of 80 to 100%, although the range of separations considered implies an extrapolation over one to two orders of magnitude outside the sensitivity domain of their observations (see also Sect. 6.2). Based on a smaller sample of 13 O+OB and 37 B+B eclipsing binaries (Harries et al. 2003; Hilditch et al. 2005), Pinsonneault & Stanek (2006) suggested the presence of a population of equalmass (“twin”) binaries albeit Lucy (2006) argued that the considered sample was too small to draw statistically significant conclusions on the presence of a peak close to unity in the massratio distribution.
Sana et al. (2012) used a MonteCarlo method to retrieve the intrinsic multiplicity properties of Otype stars in nearby Galactic open clusters. Based on a sample of 71 Otype objects, they showed that the period distribution does not follow the widely used Öpik law but is overabundant in short period systems and that the massratio distribution is flat in the range 0.1 < M_{2}/M_{1} < 1, ruling out both the presence of a twin population of equal mass binaries and random pairing from a classical initial mass function as a process to associate components in a binary. These results have strong implications for the evolution of massive stars as over twothirds of the stars born as O star will, during their lifetime, interact with a nearby companion in a binary system.
The clusters considered by Sana et al. (2012) have however relatively low masses, in the range 1000−5000 M_{⊙}. In view of the much more energetic environment provided by 10^{5}−10^{6} M_{⊙} clusters, containing up to thousands of massive stars, it is currently unknown whether the binary fraction and the parameter distributions would remain unaffected. Addressing this question has implications in the context of superstarclusters and cluster complexes observed well beyond the Local Group (e.g. Gieles et al. 2010).
The Tarantula Nebula (or 30 Doradus) in the LMC is one of the largest concentrations of massive stars in the Local Group (Walborn & Blades 1997) with a total mass of ~5.5 × 10^{4}M_{⊙} (Andersen et al. 2009). Combined with its youth (the central cluster R136 is approximately 1 to 2 Myr; de Koter et al. 1998) and the presence of several subpopulations, indicative of violent star forming activity in the past tens of millions of years, this complex system is our closest view of a starburstlike region in the Local Universe and an ideal laboratory to investigate the characteristics of massive stars.
By targeting ~800 massive stars the VLTFLAMES Tarantula Survey (VFTS; Evans et al. 2011b, Paper I) provides a nearcomplete census of the 30 Doradus region. The science drivers of this program are presented in e.g. Evans et al. (2011a) and de Koter et al. (2011). The study of the multiplicity properties of this region is an important component of VFTS and to this end a multiepoch observing strategy has been implemented to identify binaries with periods ≲10^{3} days (see Sect. 4.3), i.e. the systems in which the components are expected to interact through Rochelobe overflow. Early serendipitous findings (Taylor et al. 2011; Dufton et al. 2011) clearly emphasize the importance of binaries and the large potential these systems have for our understanding of binary properties and evolution.
In this paper, we use radial velocities (RVs, see Sect. 2) and a variability analysis to investigate the multiplicity properties of the 360 Otype stars in the VFTS sample. We first establish the observed spectroscopic binary fraction in our sample (Sect. 3). Using a MonteCarlo method, we correct for the observational biases and we constrain the intrinsic binary fraction and period and massratio distributions (Sect. 4). Possible correlations of these quantities with spatial distribution, magnitude and spectral type are investigated in Sect. 5. Our findings are discussed in Sect. 6 and summarised in Sect. 7.
2. Observations, sample and RV measurements
2.1. Data overview
The VFTS data and the data reduction process are extensively described in Paper I. Here, we provide a brief description of the VFTS campaign and observational settings related to the present study. Using nine different plate configurations, VFTS has collected multiepoch fibre spectroscopy of about 800 earlytype stars in the Tarantula region, covering a 25′ diameter field of view centered on R136. Each observation consisted of two 30min exposures taken back to back. In this paper, we analyze the Otype star population observed with the Medusa fibres and the low resolution settings LR02 and LR03 of the Giraffe multiobject spectrograph. These data provide continuous coverage of the 3960−5070 Å region at a spectral resolution of about 0.6 Å. The higher resolution data obtained in the Hα region are not used in the present analysis.
Most of the time sampling is provided by the LR02 observations (3960−4565 Å, λ/λΔ ≈ 7000), with five pointings covering time scales of days to months and one additional pointing taken about 300 days after the first observations. The LR03 setup (4500−5070 Å, λ/Δλ = 8500) was observed three times with no specific time constraint. As a result some of the configurations were observed consecutively to LR02 data while others are spread out and provide additional time sampling. For operational reasons, one LR02 plate configuration – field C – was reobserved a seventh time in Oct. 2010, extending the temporal baseline to almost 700 days for 48 O stars in our sample. The signaltonoise ratio (S/N) is however lower than average and this last pointing provides little additional information. We have performed the complete analysis presented in Sects. 3 and 4 with and without these supplementary data and found no significant differences in the results. In the rest of the paper, we ignore that epoch in order to preserve the homogeneity of the temporal sampling for all stars in our data set.
Fig. 1 Top panel: number of RV epochs obtained for each object. One Argus object, VFTS 1022, has been observed at 13 epochs and falls outside the plotted range. Middle panel: number of lines available for absolute RV measurements. Bottom panel: time elapsed between the first and last suitable epochs for RV measurement. The solid line corresponds to objects observed with Medusa while the dashed lines indicate objects observed with the Argus IFU. 

Open with DEXTER 
Fig. 2 Distribution of the Otype stars in the field of view. The population is dominated by the central cluster, NGC 2070, that contains 57% of our sample. NGC 2060, at 6.7′ to the southwest, is slightly older and accounts for 22% of the sample, while the remaining 21% is spread throughout the field of view. The dashed circles, with a 2.4′radius each, show the adopted regions for NGC 2070 and NGC 2060 (see Table 6). 

Open with DEXTER 
The combined effect of chromatic atmospheric diffraction and of seeing may strongly impact the fraction of the starlight that is injected into the Medusa fibres. As a consequence, the total response of the observing system may present significantly different shapes from one exposure to another. Each 30min spectrum is therefore normalised individually using the procedure described in Appendix A that provides a high degree of automatisation and guarantees an homogeneous handling of the data.
Contiguous observations were then combined to improve the S/N, yielding a typical number of six epochs per targets and providing a minimum timebase of about 300 days for most of the stars. Statistics of the number of epochs, diagnostic lines for RV measurements and timebase are provided in Fig. 1.
2.2. Otype stars in VFTS
The sample analyzed in this work includes all Otype stars observed during the VFTS campaign and ignore the Btype stars, the WolfRayet stars and the transitional objects to the WolfRayet domain (the socalled “slashstars”). The multiplicity properties of the B stars are investigated in a separate paper within the VFTS series (Dunstall et al., in prep.).
The Otype stars were classified by comparison, at a resolving power of 4000, with spectral standard stars. The detailed spectral classification is presented in Walborn et al. (in prep.). In total, 332 Otype stars observed with Medusa form the bulk of our sample.
The densest regions of 30 Dor, in and around the core cluster R136, were further observed with the Argus integralfield unit (IFU) and the LR02 Giraffe setup. The RVs of these objects were measured using a similar approach to the one described in Appendix B and the same set of lines and rest wavelengths (see Table B.1 and HénaultBrunet et al. 2012). The Argus RV measurements of 28 additional O stars observed with the Argus IFU have been incorporated in the present analysis. In total, the sample that we analyze thus contains 360 Otype stars and is the largest homogeneous sample of O stars investigated for multiplicity to date.
Figure 2 shows the spatial distribution of the O stars in the 30 Dor field of view. About half of our sample (184 objects) is concentrated in a 2.4′radius around R136, the young and massive cluster at the core of NGC 2070. Seventythree targets are spatially associated to NGC 2060, a slightly older cluster at 6.7′ southwest of R136 and the remaining 103 objects are spread across the field of view.
2.3. Completeness
The list of potential targets to be observed by the VFTS was selected based on a simple magnitude cutoff at V = 17 (see Paper I). The allocation of the MEDUSA fibres to VFTS sources was subject to further constraints related to the physical size of the magnetic buttons supporting the fibres, precluding observation of all sources in the most crowded parts of the fieldofview. The VFTS completeness is about 70% of the total Ostar population in the field of view and, importantly, does not show a bias towards the brighter end of the magnitude distribution. Beside the physical constraint on the fibre allocation, the VFTS sample is complete for the Otype stars, except for the late Otype objects that suffer from a significantly larger extinction – at least one magnitude – than typical in the field of view.
2.4. RV measurements
The radial velocities of our targets are obtained by measuring the Doppler displacement of spectral lines. Because the RVs are a key ingredient of this study, details of the method are provided in Appendix B, including alternative approaches that have been considered. Our method of choice, Gaussian fitting of selected sets of He lines, is prompted by four criteria: (i) robustness, (ii) internal consistency across the full range of O spectral subtypes, and the ability to provide (iii) absolute measurements and (iv) error bars. An extension to the standard linebyline Gaussian fitting is obtained by simultaneously fitting all the data available for a given object (i.e., all the epochs and spectral lines). The different lines fitted at a given epoch are required to provide the same Doppler shift while the shape of the individual lines is kept constant through all the epochs. This approach provides two significant improvements compared to linebyline fitting: a strong reduction of the number of parameters of the fit and a significant gain in robustness in the RVs measured from data with poor S/N.
Special attention has been given to identify a set of lines that provide consistent RV measurements across the range of O spectral subtypes. Lines displaying systematic deviations as a result of, e.g., line blending or wind effects have been discarded (see Appendix B.2 for a detailed discussion). The best consistency is observed between RVs obtained from the He iiλλ4200, 4541, He iλλ4387, 4713 and 4922 lines. Whenever possible, we have thus obtained absolute RV measurements based on these five lines or a subset of these lines depending on the spectral type, wavelength coverage and degree of nebular contamination. The spectral lines used for each object are reported in Table 1 while Table 2 provides the journal of the observations and, for each epoch, the measured RV and 1σ uncertainty.
The typical number of spectral lines used in measuring the RVs is between two and five (Fig. 1). For VFTS 051, 400, 444, 453, 531, 565 and 587, none of these lines were of sufficient quality. RVs for these targets were measured using the He i+iiλ4026 and/or He iiλ4686 lines. Because these are late Otype stars the measured RVs are probably only weakly affected by temperature and/or wind effects but caution should be applied in using the measured RVs in absolute terms. A number of objects displayed doublelined profiles (see Sect. 3.2.1). Measuring RVs for both components sometimes required the use of lines that are not suited to provide absolute RVs. Because SB2s are not suited for absolute RVs unless a full orbit determination can be achieved, we relax our constraints for these objects and include other He lines as indicated in Table 1. For the SB2 cases, Table 2 lists the RV of the earliest component of the pair.
3. Multiplicity analysis
In this section we investigate the multiplicity of the stars in our sample using various approaches that combine variability analysis (Sect. 3.1) and the search for composite sources (Sect. 3.2) revealed either by multiple signatures in their spectrum or by Hubble Space Telescope (HST) imaging. Section 3.3 establishes the observed spectroscopic binary fraction and Sect. 3.4 compares our results with earlier studies. The results of the multiplicity analysis are summarised in Table 3.
Journal of the observations, RV measurements and their 1σ errorbars for the stars in our sample.
Overview of the results of the multiplicity analysis.
3.1. Variable sources
3.1.1. RV variability
To identify objects displaying significant RV variations, we use a statistical test. The null hypothesis of constant RV is rejected if, for a given star, any two RV measurements deviate significantly from one another, i.e.: (1)where v_{i} and σ_{i} are the RVs and their 1σ errors for a given object at epoch i. The confidence threshold of 4.0 provides a probability of about one false variability detection in every 1000 objects, given the sampling and accuracy of our measurements. This corresponds to an expected impact on the measured fraction of RV variable objects of 10^{3}. We find that 165 sources (46%) display significant RV variability.
Keeping all other things equal, Eq. (1) provides a more sensitive criterion to detect variability than a comparison of the v_{i}’s with the average RV of the object. The results of Eq. (1) are also equivalent to those of a variability test based on the goodness of fit of a constant RV model (see e.g. HénaultBrunet et al. 2012) for over 96% of our sample.
3.1.2. Line profile variability
We also search for line profile variability (LPV) using a time variance spectrum (TVS) analysis (Fullerton et al. 1996). The TVS analysis is sensitive to variations in the amplitude, shape and position of the line profile, as well as to inconsistency in the nebular correction between the various epochs and to technical defects such as cosmic rays and continuum localisation. The method also allows us to verify the consistency of the normalisation procedure described in Appendix A. We adapt the original TVS approach of Fullerton et al. (1996) to take into account the known error bars at each pixel, rather than to rely on the overall S/N of the spectrum. This is needed because the S/N is variable as a function of wavelength and depends on the chosen setup. We thus define: (2)where F_{i} and σ_{i} are the observed normalised flux and its 1σ error as a function of wavelength λ at epoch i. ⟨ F ⟩ is the weighted average flux computed over all available epochs. The statistical significance threshold is adopted to be α = 0.01, i.e.: (3)where P_{χ2}(α, N − 1) gives the cutoff value in a χ^{2} distribution with N − 1 degrees of freedom such that the probability that a random variable exceeds P_{χ2}(α, N − 1) is equal to α. In our data, most of the variability detected through the TVS analysis is produced by nonoptimum nebular correction. It results from the combined effects of the limited spectral resolution, the spatially variable level of nebular emission and the fact that FLAMES does not provide a nebular spectrum close to the target.
In addition to residuals of the nebular correction, the TVS detects significant variability in 96 cases. By far, most cases (92) are also associated with significant RV variations so that the TVS analysis does not significantly add to the number of detected variable stars.
3.2. Composite sources
3.2.1. Doublelined binaries
The spectra obtained for each target were visually inspected to search for the signature of doublelined and asymmetric profiles. Despite the limited S/N of the individual epoch spectra, doublelined profiles were clearly identified in 50 cases. Higher multiplicity profiles were not observed. Asymmetric or variable profiles are observed for 13 objects. Candidate SB2s and objects for which asymmetric profiles are suspected represent another 21 cases. For pronounced SB2s, RV measurements were attempted using two Gaussians per profiles. Reliable fits could be achieved in 43 of the 50 cases. In two cases the stronger component – likely the primary, more massive star – shows no RV variation. These objects are thus added to the list of binary candidates although one cannot confirm the physical link between the two component of the pairs.
3.2.2. HST observations
We investigate the targets for visual companions using HST data taken with WFC3/UVIS as the primary camera with ACS/WFC in parallel, both using the F775W filter. These data constitute the first epoch of a proper motion program in the region (GO 12499, PI: Lennon, see de Mink et al. 2012; Sabbi et al. 2012). It provides a ~16′ × 12′ mosaic centered on R136 covering over 90% of the O stars in the Tarantula Survey. We find that 19 O stars observed with Medusa and six with ARGUS have one or several companions that are close and bright enough, for the FLAMES spectra to be a composite of multiple sources. These and an additional five multiple candidates are listed in Table 3.
The HST point spread function (PSF), of ~0.1′, corresponds to a physical separation of about 5 × 10^{3} AU at the distance of 30 Dor. HST thus probes a very different, nonoverlapping, separation range compared to the spectroscopic binaries discussed in Sect. 3.1. Five of the visual pairs identified by HST also show significant RV variations, indicating that at least one of the components of the visual system is also a spectroscopic binary
Fig. 3 Variations of the fraction of systems below and above the critical RV variation amplitude C separating the spectroscopic binaries (solid line) from lowamplitude RV variables (dotted line). The vertical dasheddotted line indicates the adopted threshold of C = 20 km s^{1}. 

Open with DEXTER 
3.3. Observed binary fraction
Of all the diagnostics discussed above, RV variability is by far the most robust test to identify variable sources in our data. Complementary tests (Sects. 3.1.2 to 3.2.2) do not significantly change the number of detected systems. Furthermore the RV variability tests defined by Eq. (1) do not rely on subjective appreciations such as e.g., the visual inspection of the line profile morphology. This is a necessary condition to model the observational biases by means of MonteCarlo simulations (Sect. 4.1). In the following we will thus focus on spectroscopic systems that reveal themselves through their RV variability.
The test provided by Eq. (1) only identifies significant RV variations. To distinguish cases where the observed variability is genuinely caused by orbital motion in a binary system from potential intrinsic variability due to e.g. photospheric or wind effects, we impose a minimum amplitude threshold C on any significant deviations computed from any individual pair of epochs. I.e., an object is considered a spectroscopic binary if at least one pair of RV measurements satisfies simultaneously (4)Figure 3 shows the influence of C on the identified binary fraction. It reveals a kink at about 15 km s^{1} that may mark the transition between a regime were photospheric variability has a significant contribution (ΔRV < 15 km s^{1}) and the regime of genuine spectroscopic binaries. A crude extrapolation of the contribution of the binaries into the lowamplitude ΔRV regime suggests the presence of about 5% of genuine binaries with ΔRV < 15−20 km s^{1}. Similarly, we estimate that photospheric variabilities affect as little as 6 to 7% of our sample.
In Sect. 4 we will correct for the observational biases and consider the fraction of long period systems that potentially fall below the ΔRV cutoff. As a consequence, the exact value of the adopted cutoff is of little importance but has to be large enough to avoid a significant number of false detections due to measurement errors and/or photospheric variability.
Because photospheric variations in supergiants can mimic variations with amplitudes of up to 20 km s^{1} (e.g. Ritchie et al. 2009), we conservatively adopt C = 20 km s^{1}. Eleven per cent of the sample (39 stars) are thus identified as variable but feature an amplitude of the significant RV variations that remains below 20 km s^{1}. These objects might be long period binaries or stars with photospheric variability and, as discussed above, we estimate an equal contribution of both categories.
Given Eq. (4) with C = 20 km s^{1}, the observed spectroscopic binary fraction identified through RV variation of the component with the strongest line – likely the primary – is thus 35.0 ± 2.5%, where the 1σ error bar gives the statistical error that results from the size of the sample (Sana et al. 2009).
3.4. Comparison with earlier results
Bosch et al. (2009, BTT09 reported on a similar study of 54 O and early B stars in NGC 2070. They obtained R ~ 3700 sevenepoch spectroscopy covering a wavelength range equivalent to the coverage of our LR02+LR03 settings. Using the externaltointernal velocity dispersion ratio σ_{E}/σ_{I} and a threshold of 3 they identified 17 RV variables and nine additional SB2 binaries and binary candidates for a maximum detected binary fraction of 0.48. The results of BTT09 thus suggest a larger binary fraction than the one we observed in VFTS. To understand this difference, we perform a detailed comparison of the VFTS and BTT09 approaches.
We first apply our variability and binary criteria on BTT09’s RV measurements. We detect 24 RV variable objects among the 48 objects listed in their Table 2. The difference comes from the fact that the ratio σ_{E}/σ_{I} is artificially reduced if one measure has a significantly worse accuracy that the others. All of the stars identified as RV variables have peaktopeak amplitudes in excess of 20 km s^{1} and thus qualify as spectroscopic binaries according to the criteria of Eq. (4). Applying our method on the BTT09 RV data set, we measure a spectroscopic binary fraction of 50.0% while BTT09 obtained 35.4% by using the σ_{E}/σ_{I} ratio. Two of the seven binaries not identified as RV variable by BTT09 are still accounted for in their total binary count because of asymmetry in their line profiles. Even though we require a stricter (4σ) threshold, we conclude that our binary criteria are slightly more sensitive than the one used by BTT09 based on σ_{E}/σ_{I} > 3.
As a second step in comparing our results with those of BTT09 we compute the binary fraction from the overlapping sample between BTT09 and the present work, which are 34 of the 54 objects studied by BTT09. Of these 34 objects, BTT09 reported RVs for 30 of them among which they identified 12 binaries (40%) through RV variations. Among the same set of 30 stars, VFTS identify 13 binaries (43%) and 6 (20%) lowamplitude variable objects, with ΔRV ranging from 7 to 19 km s^{1}.
On the overlapping sample, the two studies agree well, with VFTS being slightly more sensitive due to the improved RV accuracy and statistical criteria. The significantly lower binary fraction observed in the whole VFTS sample relative to the work of BTT09 is attributed to the focus of BTT09 on the brightest stars. We will show in Sect. 5.2 that fainter Otype objects have an intrinsically lower binary fraction.
Fig. 4 Histogram of the peaktopeak amplitude RV variations of each object. Only significant variations according to Eq. (4) are considered. Therefore the first bin is empty. 

Open with DEXTER 
4. Intrinsic multiplicity properties
The detected binary population provides a minimum estimate of the true binary population in the Tarantula region. The number of undetected binaries depends on the sensitivity of our campaign in terms of RV accuracy and time sampling and on the intrinsic distribution of the orbital parameters. The problem of correcting for the observational biases is thus illposed and we employ here a MonteCarlo approach to estimate the impact of the observational biases and to provide constraints on the intrinsic binary fraction and distributions of orbital parameters.
4.1. MonteCarlo method
In our MonteCarlo approach, we simulate massive star populations with intrinsic multiplicity properties randomly drawn from adopted parent distributions. Accounting for observational biases due to sampling and measurement uncertainties we search for sets of distributions that reproduce the properties of the observations in three aspects: the observed binary fraction, the peaktopeak amplitude ΔRV of the RV variations (Fig. 4) and the minimum time scale ΔHJD for significant RV variation to be observed (Fig. 5).
The method is detailed in Appendix C. In short, simulated and observed distributions are compared by means of KolmogorovSmirnov (KS) tests. The detected binary fraction predicted by the simulations is compared to the observed fraction using a Binomial distribution (Sana et al. 2012). A global merit function (Ξ′) is constructed by multiplying the two KS probabilities obtained for the ΔRV and ΔHJD distributions and the Binomial probability: (5)where is the simulated fraction of detected binaries, N_{bin} is the detected number of binaries in the observations and N is the population size.
Following Kobulnicky & Fryer (2007) and Sana et al. (2012), we use powerlaws to describe the intrinsic distribution functions of orbital parameters: f(log _{10}P/d) ~ (log _{10}P)^{π}, f(q) ~ q^{κ} and f(e) ~ e^{η}. The exponents π, κ and η are then varied to explore different shapes for the parameter distributions. Because the global merit function is insensitive to the eccentricity distribution (see Fig. C.1), we cannot constrain η and we adopt the eccentricity distribution measured using the O star populations of Galactic young open clusters (Sana et al. 2012), i.e. f(e) ~ e^{0.5}. We thus explore the behaviour of Ξ′ in the threedimensional parameter space defined by the intrinsic binary fraction and distributions of the orbital periods and massratios. The exploration domains and mesh sizes are given in Table 5.
4.2. Results
Figure 6 displays the merit function (Ξ′) projected on planes defined by the various pairs of parameters. It shows that the optimum is well defined within the explored ranges. The best agreement between the observed and simulated distributions is obtained for π = −0.45 ± 0.30, κ = −1.00 ± 0.40 and a intrinsic binary fraction of f_{bin} = 0.51 ± 0.04. The quoted errors have been estimated as the 1σ confidence interval on the retrieved parameters of 50 synthetic data sets randomly drawn from parent distributions close to our best fit distributions (see Appendix C.3 and the lower part of Table C.1). The comparison between the simulated and observed cumulative distribution functions (CDFs) of the peaktopeak RV amplitudes and of the variability time scales is presented in Fig. 7.
Fig. 5 Cumulative distribution of the minimum time difference ΔHJD for a binary to display significant RV variations according to Eq. (4) and for different RV amplitude thresholds (C) in the RV variation amplitudes. The cumulative distributions are normalised to the total number of binaries. 

Open with DEXTER 
Overview of the results of various variability and multiplicity criteria used in this paper.
Overview of properties of the MonteCarlo grid.
Fig. 6 Projections of the global merit function Ξ′ on the various pairs of planes defined by π, κ and f_{bin}. The absolute maximum is indicated by a cross (+). Contours indicate loci of equalvalues of the merit function, expressed as a fraction of its absolute maximum (see legend). 

Open with DEXTER 
Fig. 7 Comparison between the observed (crosses) and simulated (lines) cumulative distributions of the peaktopeak RV amplitudes (top row) and of the variability time scales (bottom row). In the left (resp. right) hand panels, we vary the exponent π of the period distribution (resp. κ of the massratio distribution) by ±1σ. The upperleft legends indicate the values of π, κ and η considered in each panel. The bottomright legends give the overall VFTS detection probability for the adopted parent distributions. 

Open with DEXTER 
While the overall shapes of the observed and simulated functions are well matched, some features are not fully reproduced by the simulations. The CDF of the peaktopeak RV amplitudes presents an overabundance of systems with ΔRV in the range 300 to 400 km s^{1}. This could be reproduced by increasing the fraction of short period systems, i.e. a more negative π exponent, but the CDF of the variability time scale would then not be properly reproduced. Alternatively, a flatter distribution of massratios would help to improve the CDF(ΔRV) fit in the 300 to 400 km s^{1} range but would overestimate it in the lower amplitude range (ΔRV < 100 km s^{1}).
We also attempt to reproduce the overabundance of systems with ΔRV in the range 300 to 400 km s^{1} by including specific populations of short period systems, of equal mass systems, of highmass primaries but none of these, nor any combinations of them, were able to reproduce the observed localised bump.
The bump in the CDF(ΔRV) curve comprises less than 15 systems and can also result from statistical fluctuations. Indeed, the individual P_{KS} probabilities are already large: P_{KS}(ΔRV) ≈ 0.4 and P_{KS}(ΔHJD) ≈ 0.8. We thus consider the present fit to provide an adequate representation of the data.
The value of the period distribution index is remarkably close to the value π_{GOC} = −0.55 ± 0.22 found in the O star population of Galactic open clusters (GOCs, Sana et al. 2012). The intrinsic binary fraction is 51% in VFTS compared to 69 ± 9% in the GOCs. Finally, the massratio distribution obtained in the present study favours lower mass companions while the mass ratios in the Sana et al. study were equally distributed (κ_{GOC} = −0.1 ± 0.6). Both studies do however agree within 2σ. The implications of these results are discussed in Sect. 6.
4.3. Sensitivity domain of the VFTS campaign
We use the method of Sana et al. (2009) to estimate, given the distributions of the orbital parameters that we find, the sensitivity of our campaign with respect to the orbital period, the mass ratio, the eccentricity and the primary mass.
Compared to the Sana et al. (2009) approach, we now include the measurement uncertainties obtained at each epoch in the simulations and we implement the binary detection criteria given by Eqs. (1) and (4) to consistently follow the binary detection strategy laid out in this paper.
The detection probability curves in Fig. 8 are projections from the multidimensional space onto each of the four dimensions investigated. The binary detection probability is a strong function of orbital period. It moderately depends on the mass ratio except for systems with a large mass difference. The detection probability is also weakly affected by the primary mass as well as by the eccentricity. Because we have not taken into consideration the fact that high eccentricity systems may preferentially have a long period, we underestimate the detection rate at very high eccentricities (e > 0.7).
Fig. 8 Overview of the binary detection probability achieved by the VFTS campaign as a function of the orbital period, the mass ratio, the eccentricity and the primary mass. (See Sect. 4.3 for more details.) 

Open with DEXTER 
5. Binary fraction
In this section, we discuss how the binary fraction correlates with the spatial distribution, brightness, spectraltype and luminosity class of the stars in our sample. Supplementary figures supporting this discussion are provided in Appendix D.
Observed fraction of RV variables and binaries among different O star populations in 30 Dor.
5.1. Spatial variations
We search for variations in the observed binary fraction throughout the field of view using various techniques. First we divide the stars in three samples based on their spatial location: stars within a 2.4′ radius from the centre of NGC 2060 and NGC 2070, and stars outside these two clusters. NGC 2060 tends to display a slightly larger fraction of lowamplitude RV variable stars (Table 6). This suggests a larger number of objects affected by photospheric variability, likely supergiants, reflecting the older age of the cluster, hence the larger fraction of supergiants (see also Fig. 13). When the sample is limited to stars with large amplitude RV variations, the genuine binaries, there are no obvious spatial differences between the various populations.
We also look for spatial variations in the field of view using a 2D smoothing. We do not identify significant variations save for two trends. There are slightly fewer binaries towards the southern and, to a lesser extent, northern parts of the field of view, as well as some enhancement in between NGC 2070 and NGC 2060. These trends are significant at the 95% level (≈1.7σ), but do not reach the 2σ confidence level.
Finally, we also build the radial cumulative distribution of the binary fraction, starting both from the inside and from the outside (Fig. 9). The centre of the field of view is taken at the core of R136, as defined in Table 6. Overall, there is no significant variability but for the edges of the field. The binary fraction for the 22 outermost objects is 0.14 ± 0.07 which is just at 3σ from the average binary fraction in the whole sample. A visual impression of the lack of binaries in the outer regions can be obtained by comparing the locations of the single and binary stars as shown in Fig. 10.
Fig. 9 Cumulative radial distribution of the observed binary fraction as a function of the distance to the centre of the field. The black solid line shows the distribution computed insideout, and the grey/red line shows the distribution computed outsidein. The dashed curves indicate the ± 1σ confidence interval on both curves. 

Open with DEXTER 
Fig. 10 From left to right, distributions of the RV constant, likely single stars, the binary candidates (lowamplitude RV variables) and the binary stars in the 30 Dor field of view. 

Open with DEXTER 
5.2. Brightness dependence
Objects in our sample span five magnitudes in the V band. We use the Vband photometry published in Paper I to investigate any variations of the binary fraction with magnitude. We also employ K_{s}band photometry from the VISTA Magellanic Clouds (VMC) Survey (Cioni et al. 2011) to the same end. Specifically, we use PSFfitting photometry of the 30 Dor VMC Survey tile from Rubele et al. (2012), with crossmatches to the VFTS sources listed by Zaggia et al. (in prep.). K_{s}band magnitudes from the PSFfitting analysis are unavailable for 15 VFTS stars, either due to crowding, blending or saturation issues. VMC photometry is also missing for the Argus targets and these are henceforth excluded from the present analysis. Where possible, the VMC Survey results are supplemented with 2MASS photometry using the crossmatches given in Paper I.
The analysis of both bands provides similar results. Because of the variable extinction in the region (Maíz Apellániz et al., in prep.) we focus here on the K_{s} data only. Figure 11 displays the magnitude distributions of the VFTS O stars and the observed binary fraction in various magnitude bins. The latter seems relatively uniform for objects brighter than K_{s} = 15.5 mag. For dimmer objects, the observed binary fraction starts to decrease and reaches 16% at the magnitude cutoff of our survey: K_{s} ≈ 16.5 mag, corresponding to V = 17 mag.
The apparent lower binary fraction for fainter stars can be explained by the combination of two observational effects:

(i)
spectra from fainter stars have a lower S/N, thus RV measurements are on average less accurate and small RV variations are more difficult to detect;

(ii)
binaries are typically brighter than single stars of the same spectral type by up to 0.75 mag, depending on the brightness of the companion. Binary systems thus tend to avoid the faintest bins, which then appear to present a lower binary fraction.
Fig. 11 Top panel: number of objects (dashed line) and of binaries (solid line) per K_{s}band magnitude bin. Bottom panel: observed binary fraction in corresponding magnitude bins (solid line) and intrinsic binary fraction (dashed line). Error bars are plotted whenever a bin contains more than one star. Figure D.2 provides similar figures showing the dependence with luminosity class. 

Open with DEXTER 
Number of objects, observed binary fraction detection probability and intrinsic binary fraction as a function of the K_{s} magnitude.
To investigate the effect of (ii) we simulate the magnitude distribution of an Otype population madeup of 50% single stars and 50% binaries. We find that the intrinsic binary fraction in the brightest bins was typically close to 60% but drops to 20−25% in the last bin, in excellent agreement with the intrinsic binary fraction measured in Table 7. Adding a slightly variable extinction easily allows us to extend the effect over the last two bins. If early Btype stars with a similar binary fraction are included in the simulations the magnitude bins at K_{s} ~ 15.5−16.5 do not show such a drop.
We thus conclude that the variation of the intrinsic binary fraction with magnitude is simply a consequence of the fact that our sample is defined by a cutoff in spectral type, i.e. the Otype subsample of the VFTS is not magnitude limited, save for the few stars with severe extinction. We also conclude that O star binary fractions obtained from magnitudelimited samples can be significantly overestimated, as such a selection criterion picks up extra binaries near the magnitude limit.
This depletion of binaries for fainter objects is actually responsible for the apparent difference between our overall results and the observed binary fraction obtained by Bosch et al. Our observed binary fraction among objects with K_{s} < 15.5 is 41.4%. Based on the results of Table 7, we estimate an intrinsic binary fraction of 61.4% in that brightness range.
Fig. 12 Top panel: number of objects (dashed line) and of binaries (solid line) per O spectral subtype. Bottom panel: observed binary fraction as a function of the spectral subtypes. Error bars are plotted whenever a bin contains more than one star. Figure D.3 provides similar figures for the different luminosity classes. 

Open with DEXTER 
Fig. 13 Spatial distribution and multiplicity properties of the supergiants in the field of view. 

Open with DEXTER 
5.3. Spectral subtype and luminosity class dependence
To investigate possible variations of the detected binary fraction with the spectral subtype and luminosity class, we ignore all objects which spectra have a too poor S/N or are too heavily contaminated by nebular lines to yield decent spectral classification. We also ignore the Argus sources that have a more limited spectral coverage and thus a less accurate spectral classification, and the seven Onfp stars, that present peculiar spectra as noted by the qualifier “p” appended to their spectral type. We are left with 295 out of 328 Medusa objects (i.e., 90% of the Medusa sample). These objects are distributed among the luminosity classes as follows: 30 III objects, 77 III objects and 186 IVV objects. Figure 12 shows the observed binary fraction as a function of spectral subtype and Fig. D.3 splits Fig. 12 according to luminosity class. Table 8 indicates the overall lowamplitude RV variable and binary fraction in each category. Most of the fluctuations are compatible with the sample sizes so that the binary fraction seems homogeneous across the spectral subtypes in the various samples. With two binaries and one lowamplitude RV variable object in a sample of seven, the binary fraction of the Onfp stars is also compatible with that of the main sample.
Observed fractions of lowamplitude RV variables and binaries as a function of the luminosity class.
The O6 stars display a significantly lower binary fraction than the average of the population. Inspection of Fig. D.3 reveals that the drop in the binary fraction at O6 is due to the dwarfs, with only one out of 13 O6 V stars being detected as a binary, thus a binary fraction of 0.08 ± 0.07.
The binary fraction among supergiants (23%) is smaller than among dwarfs and giants (~40%, Table 8). Interestingly, these numbers are roughly in line with expectations from binary evolution. The larger radii of I and II stars may imply that the shortest period binaries (P < 2−3 d) have already interacted. Because postinteraction systems are most likely to appear as single stars in RV studies (de Mink et al., in prep.), the fraction of detected binaries is thus expected to be lower.
With only one binary out of 12 objects (corresponding to an observed binary fraction of 0.08 ± 0.08) but a large fraction (42%) of objects showing lowamplitude RV variability, the O9.7 I and II stars are actually responsible for most of the differences between the binary fraction among supergiants and among dwarfs and giants. Figure 13 compares the spatial location of the O9.7 III objects with the other class III objects in the field. It reveals that the O9.7 supergiants preferentially belong to NGC 2060 while the hotter supergiants are distributed according to the ratio of the populations of NGC 2060 and NGC 2070. The data associated with the O9.7 supergiants are of a high quality, therefore exclude the possibility that the sample suffers from a higher detection threshold than the other supergiants. It may thus be interesting to speculate that the O9.7 III population in NGC 2060 contains a large fraction of apparently single postinteraction stars.
5.4. Summary
The binary fraction displays a high degree of homogeneity across the different populations. Only three subcategories deviate significantly from the mean:

i.
the overall binary fraction is lower in the outer region (r > 7.8′) of the field of view;

ii.
the binary fraction is lower for the fainter stars in the sample (K_{s} > 15.5), an apparent effect that results from observational biases and the definition of the sample in terms of spectral type;

iii.
the O9.7 III stars, mostly located in NGC 2060, are predominantly single, though a relatively large fraction show lowamplitude RV variability.
Fig. 14 Massratio distributions obtained using different pairing mechanisms (solid and dashed lines) compared to those obtained with a powerlaw density distribution with different κ values (dotted lines). 

Open with DEXTER 
6. Discussion
6.1. Nature of the pairing process
We investigate two different pairing processes to associate the components of the binary systems and we compare the results of the massratio distributions of simulated populations with the one derived from the observations.

i.
Random pairing: primary and secondary masses are independently drawn, between 0.05 and 80 M_{⊙}, from a mass function with a slope of −2.35. Stars are randomly associated in pairs;

ii.
Truncated secondary mass function: the primary masses are drawn as above. The secondary masses are drawn from a mass function truncated at the primary mass.
Figure 14 illustrates the massratio distributions obtained. None of the tested pairing processes produces a distribution function that resembles f(q) ∝ q^{κ} with − 1 < κ < 0. It also shows that, as expected, random pairing results in κ = −2.35 and that the truncated secondary mass function comes close to κ = 1.0.
6.2. On the universality of the multiplicity properties
Our results can be discussed in the broader context laid out by recent studies of two different O star samples in the Milky Way. First, Kobulnicky & Fryer (2007) attempted to model the multiplicity properties in Cygnus OB2 using the results of a multi epoch RV campaign targeting 120 OB stars (900 epochs in total). Because they only adjusted the χ^{2} distribution, the authors obtained a degenerate solution in the exponent of the period and massratio distributions and they chose to fix π at 0.0. The separation regime considered in their modelling extends at least up to 100 AU and, for some of their simulations, to 10 000 AU, i.e., quite beyond the sensitivity range of their spectroscopic observations. In the considered parameter ranges, Kobulnicky & Fryer obtained a large binary fraction, at least 80%, and κ values ranging from − 0.6 to 1.0. Because of the different approaches and the different parameter ranges considered, direct comparison with the current results is hazardous.
Second, Sana et al. (2012) studied the multiplicity properties of the O star populations in nearby young open clusters. They performed a similar MonteCarlo analysis as done here. Most of the binaries identified in their sample have constraints on their orbital parameters. This allowed them to test directly the period and massratio distributions. They obtained an intrinsic binary fraction of 0.69 ± 0.09 and values of π = −0.55 ± 0.22 and κ = −0.1 ± 0.6 for the powerlaw exponents of the period and massratio distributions.
The period distributions obtained in both studies are in excellent agreement. The massratio distribution in Sana et al. is almost flat while that of the VFTS sample is more abundant in lower massratios. Both results are not inconsistent within 2σ and, as discussed in Appendix C.2, the present approach is only weakly sensitive to κ. The main difference resides in the intrinsic binary fraction which seems significantly lower in the Tarantula region. This aspect will be discussed in Sect. 6.3.
6.3. Evolutionary implications
Assuming that today’s distributions are representative of the distributions at birth and using a series of assumptions on critical period and massratio values that delimit various binary interaction scenarios, the integration of the distribution functions of Sect. 4.1 provides an estimate of the evolutionary fate of the 30 Dor O star population. For comparison purposes, we adopt the same set of assumptions as in Sana et al. (2012), namely: the maximum period at which significant binary interaction takes place is 1500 d. Binaries with periods less than 6 days will initiate masstransfer while the primary is still on the main sequence (Case A mass transfer). Among these systems, binaries with periods less than two days will merge. Between 2 and 6 days, binaries with massratios less than 0.65 will also merge. Similarly a small fraction (5%) of case BC mass transfer (i.e., mass transfer occurring after the main sequence) will lead to coalescence. In the remaining interaction cases, the primary will be stripped from most/all of its envelope while the secondary will gain mass and angular momentum and will be spun up to critical rotational velocities.
Under these assumptions and given an intrinsic binary fraction of 0.51 in the range 0.15 < log _{10}P/d < 3.5, we estimate that 53% of all the stars born as Otype star belong to a binary system with a period less than 1500 d. The evolution of these stars is strongly affected by binary interaction: 18% of the O stars will merge with a companion, 27% will be stripped from their envelope and 8% are expected to be spun up. About one third of the binary interaction thus results in coalescence of the two companions.
Compared to the fractions of stars in the various binary evolution channels derived in the Milky Way (Sana et al. 2012), we obtain from the VFTS population an equivalent merger rate but lower rates of stripping and spinup. In the above, we have implicitly assumed that today’s distributions are representative of the distributions at birth. In an environment such as 30 Dor, this assumption may not hold. The presence of different populations in the field, some already quite old, and the fact that at least 5% of the O star population is running away, either because of dynamical interaction or supernova kick, suggest that a fraction of the current O star population has already undergone significant dynamical and/or evolutionary interaction.
In particular, today’s binary fraction might be dragged towards lower values by the presence of runaways and postinteraction objects. As mentioned already the latter are indeed expected to be predominantly single (de Mink et al., in prep.). Runaways and postinteraction objects would need to represent about 20% of the current 30 Dor population to reconcile the intrinsic binary fraction in 30 Dor (51%) with that of the Galactic open cluster population (69%, Sana et al. 2012). Establishing the star formation and dynamical history of the Tarantula region would constitute a major step in deciding whether such a high fraction is plausible.
It remains that binary interaction effects have a critical influence on the evolutionary path of more than half the stars born as Otype stars. These effects need to be considered to better interpret massive star populations seen in integrated light and to accurately explore the frequency of progenitors of compact objects, highmass Xray binaries, hydrogenpoor supernovae and longduration gammaray bursts.
7. Conclusions
In this paper, we investigate the multiplicity properties of a sample of 360 Otype stars observed by multiepoch spectroscopy in the framework of the VLTFlames Tarantula Survey. Fortysix per cent of our sample present significant RV variations and, for 35%, the detected variations have an amplitude larger than 20 km s^{1} and are very likely spectroscopic binaries.
The observed binary fraction presents a high degree of homogeneity across the field of view and among the various spectraltypes and luminosity classes. Three subgroups however show a significantly lower binary fraction:

i.
the binary fraction is lower in the outer region (r > 7.8′; ≈117 pc) of the field of view. This outer population may be dominated by single stars ejected from the core of the region;

ii.
the binary fraction is lower for the fainter stars in the sample (K_{s} > 15.5), which results from the sample definition and traces the brightness separation between the O and B stars ;

iii.
the O9.7 III stars are predominantly single and belong to NGC 2060, the older Ostar cluster in the field. Though speculative at the moment, a large fraction of the O9.7 III stars may be postinteraction stars.
We use a MonteCarlo method to correct for the observational biases and to constrain the intrinsic multiplicity properties of the Otype star population. We obtain an intrinsic binary fraction of f_{bin} = 0.51 ± 0.04. The most likely period distribution, f(log _{10}P/d) ~ (log _{10}P/d)^{0.45} in the interval log _{10}P/d ∈ [0.15,3.5], favours shorter period systems compared to a flat distribution in log _{10}P/d. Similarly the massratio distribution, f(q) ~ q^{1.0} in the interval q ∈ [0.1,1.0], favours lowermass companions but we note that our method only provides weak constraints on the q distribution.
The period distribution obtained here is strikingly similar to the one obtained for Galactic Otype binaries, possibly hinting at a universal property among massive binaries. The intrinsic fraction of binaries with periods less that 10^{3.5} d (0.51 ± 0.04) seems lower than obtained for the Galactic Otype binaries (0.69 ± 0.09) albeit the two results agree within 2σ. This could imply an environmental effect on the outcome of star formation. It could also reflect the fact that the binary properties in the Tarantula region have already been significantly affected by dynamical and/or stellar evolution that would predominantly lead to either a merger event or disrupt the binary. Quantitatively understanding if and how binary evolution and dynamical interaction have affected the multiplicity properties in the Tarantula region would help to decide whether the observed differences are nature or nurture. A key ingredient in this process will be a better understanding of the star formation history of the region.
Our results emphasize again that multiplicity is one of the main characteristics of massive stars and that it needs to be taken into account to properly predict the evolution of entire populations of massive stars. Dynamical interactions and evolutionary effects that could have possibly affected our sample would be expected to, predominantly, decrease the observed number of binaries; our derived binary fraction can therefore be seen as a lower
limit to the binary fraction at birth. This implies that the most frequent final product of highmass star formation is an O + OB binary and that massive star formation theories not only have to explain the formation of high mass stars but also have to reproduce their multiplicity properties.
Acknowledgments
S.d.M. acknowledges NASA Hubble Fellowship grant HSTHF51270.01A awarded by STScI, operated by AURA for NASA, contract NAS 526555.S.E. STScI is operated by AURA, Inc., under NASA contract NAS526555. J.M.A. acknowledges support from the Spanish Government Ministerio de Educación y Ciencia through grants AYA201015081 and AYA201017631 and the Consejería de Educación of the Junta de Andalucía through grant P08TIC4075. AH ackowledges support by the Spanish MINECO under grants AYA201021697C0504 and ConsoliderIngenio 2010 CSD200600070, and by the Canary Islands Government under grant PID2010119. M.G. acknowledges financial support from the Royal Society and VHB acknowledges support from the Scottish Universities Physics Alliance (SUPA) and from the Natural Science and Engineering Research Council of Canada (NSERC). Finally, H.S. acknowledges support from the SARA Computing and Networking Services.
References
 Andersen, M., Zinnecker, H., Moneti, A., et al. 2009, ApJ, 707, 1347 [NASA ADS] [CrossRef] [Google Scholar]
 Barbá, R. H., Gamen, R., Arias, J. I., et al. 2010, in Rev. Mex. Astron. Astrofis. Conf. Ser., 38, 30 [Google Scholar]
 Bosch, G., Terlevich, E., & Terlevich, R. 2009, AJ, 137, 3437 [NASA ADS] [CrossRef] [Google Scholar]
 Cioni, M.R. L., Clementini, G., Girardi, L., et al. 2011, A&A, 527, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Conti, P. S., Leep, E. M., & Lorre, J. J. 1977, ApJ, 214, 759 [NASA ADS] [CrossRef] [Google Scholar]
 De Becker, M., Rauw, G., Manfroid, J., & Eenens, P. 2006, A&A, 456, 1121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de Koter, A., Heap, S. R., & Hubeny, I. 1998, ApJ, 509, 879 [NASA ADS] [CrossRef] [Google Scholar]
 de Koter, A., Sana, H., Evans, C. J., et al. 2011, J. Phys. Conf. Ser., 328, 012022 [NASA ADS] [CrossRef] [Google Scholar]
 de Mink, S. E., Cantiello, M., Langer, N., et al. 2009, A&A, 497, 243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de Mink, S. E., Lennon, D. J., Sabbi, E., et al. 2012, in Am. Astron. Soc. Meet. Abstr., 219, 151.13 [Google Scholar]
 Dufton, P. L., Dunstall, P. R., Evans, C. J., et al. 2011, ApJ, 743, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Evans, C., Taylor, W., Sana, H., et al. 2011a, The Messenger, 145, 33 [NASA ADS] [Google Scholar]
 Evans, C. J., Taylor, W. D., HénaultBrunet, V., et al. 2011b, A&A, 530, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fullerton, A. W., Gies, D. R., & Bolton, C. T. 1996, ApJS, 103, 475 [NASA ADS] [CrossRef] [Google Scholar]
 Gieles, M., Sana, H., & PortegiesZwart, S. F. 2010, MNRAS, 402, 1750 [NASA ADS] [CrossRef] [Google Scholar]
 Harries, T. J., Hilditch, R. W., & Howarth, I. D. 2003, MNRAS, 339, 157 [NASA ADS] [CrossRef] [Google Scholar]
 HénaultBrunet, V., Evans, C. J., Sana, H., et al. 2012, A&A, 546, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hilditch, R. W., Howarth, I. D., & Harries, T. J. 2005, MNRAS, 357, 304 [NASA ADS] [CrossRef] [Google Scholar]
 Hillwig, T. C., Gies, D. R., Bagnuolo, Jr., W. G., et al. 2006, ApJ, 639, 1069 [NASA ADS] [CrossRef] [Google Scholar]
 Kobulnicky, H. A., & Fryer, C. L. 2007, ApJ, 670, 747 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Langer, N., Cantiello, M., Yoon, S.C., et al. 2008, in Massive Stars as Cosmic Engines, eds. F. Bresolin, P. A. Crowther, & J. Puls, IAU Symp., 250, 167 [Google Scholar]
 Lucy, L. B. 2006, A&A, 457, 629 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mahy, L., Nazé, Y., Rauw, G., et al. 2009, A&A, 502, 937 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mason, B. D., Hartkopf, W. I., Gies, D. R., Henry, T. J., & Helsel, J. W. 2009, AJ, 137, 3358 [NASA ADS] [CrossRef] [Google Scholar]
 Pinsonneault, M. H., & Stanek, K. Z. 2006, ApJ, 639, L67 [NASA ADS] [CrossRef] [Google Scholar]
 Podsiadlowski, P., Ivanova, N., Justham, S., & Rappaport, S. 2010, MNRAS, 406, 840 [NASA ADS] [Google Scholar]
 Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rauw, G., Nazé, Y., Fernández Lajús, E., et al. 2009, MNRAS, 398, 1582 [NASA ADS] [CrossRef] [Google Scholar]
 Ritchie, B. W., Clark, J. S., Negueruela, I., & Crowther, P. A. 2009, A&A, 507, 1585 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rubele, S., Kerber, L., Girardi, L., et al. 2012, A&A, 537, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sabbi, E., Lennon, D. J., Gieles, M., et al. 2012, ApJ, 754, L37 [NASA ADS] [CrossRef] [Google Scholar]
 Sadowski, A., Belczynski, K., Bulik, T., et al. 2008, ApJ, 676, 1162 [NASA ADS] [CrossRef] [Google Scholar]
 Sana, H., & Evans, C. J. 2011, in IAU Symp. 272, eds. C. Neiner, G. Wade, G. Meynet, & G. Peters, 474 [Google Scholar]
 Sana, H., Gosset, E., Nazé, Y., Rauw, G., & Linder, N. 2008, MNRAS, 386, 447 [NASA ADS] [CrossRef] [Google Scholar]
 Sana, H., Gosset, E., & Evans, C. J. 2009, MNRAS, 400, 1479 [NASA ADS] [CrossRef] [Google Scholar]
 Sana, H., James, G., & Gosset, E. 2011, MNRAS, 416, 817 [NASA ADS] [CrossRef] [Google Scholar]
 Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Stacy, A., Greif, T. H., & Bromm, V. 2012, MNRAS, 422, 290 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, W. D., Evans, C. J., Sana, H., et al. 2011, A&A, 530, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Turk, M. J., Abel, T., & O’Shea, B. 2009, Science, 325, 601 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Walborn, N. R., & Blades, J. C. 1997, ApJS, 112, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Yoon, S.C., Woosley, S. E., & Langer, N. 2010, ApJ, 725, 940 [NASA ADS] [CrossRef] [Google Scholar]
Online material
Appendix A: Normalisation
In this section, we describe the semiautomatic procedure used to normalise the spectra discussed in this paper. The adopted algorithm uses all continuum points after automatically rejecting stellar and interstellar lines (both in emission and absorption) and cosmetic defects (bad columns, cosmic rays).
The global trend in the spectrum is first removed using a onedegree polynomial fitted to the continuum, resulting in a spectrum scaled around unity. A polynomial of degree 4 to 8 is then fitted to the scaled continuum. The error bars on the flux in each pixel are in principle provided by the ESO pipeline. By a comparison of the error spectra with the empirical S/N in continuum regions, we find that the pipeline actually overestimates the error spectra by a factor 1.2, which we correct for before using the pipelineprovided errors to weight the fit. The normalisation procedure proceeds as follows:

i.
the spectrum is smoothed using a 21pixel wide median filter;

ii.
the spectral lines are removed by comparing the results of median and max/min filters applied on the initial spectrum, accounting for the spectrum S/N;

iii.
a first guesssolution of the continuum position is obtained using the smoothed spectrum (ii);

iv.
a secondguess solution is obtained using the smoothed spectrum (ii) after σclipping around (iii);

v.
the final solution is obtained using the original spectrum (i.e., no smoothing) after σclipping around (iv). Manual exclusion of specific regions is allow at this stage to incorporate a priori knowledge on unsuitable continuum regions (e.g., broad diffuse interstellar bands, broad emission regions);

vi.
the spectrum (v) is scaled by the median value computed over its continuum regions. This allows to correct for slight global under/over estimates of the average continuum level which occurs when there are many absorption/emission lines as the line wings may induce a small but systematic bias. This empirical correction is of the order of a fraction of a per cent;

vii.
the normalisation function is applied to the error spectrum.
Appendix B: Radial velocity measurements
Appendix B.1: Methodology
The choice of the RV measurement method is guided by the need to provide consistent, unbiased and homogeneous absolute RVs with representative error bars over the full range of O spectral subtypes and given the specificity of the data set. In our case, only a few lines are suitable for RV measurements. The number of epochs, the number of lines and their quality change from one star to another as a function of the star’s spectral type and brightness and the degree of nebular contamination. The S/N also varies significantly from one star to another, from one epoch to another and over different spectral regions.
Three RV techniques were considered: line moments, crosscorrelation and Gaussian fitting. Tests were performed both on synthetic data and on a small set of objects taken from our survey. Synthetic spectra for a range of temperatures, gravities, projected rotational velocities (vsini) and S/N representative of the VFTS O star data were generated using fastwind (Puls et al. 2005) and degraded to the resolution of our survey. We briefly discuss the pros and cons of the three approaches below.
Line moments turn out to be very sensitive to residuals of the nebular correction and the accuracy of the method is not explored further. Auto crosscorrelation (Dunstall et al., in prep.) provides accurate RVs and is rather robust against residuals of the nebular correction. However, for stars with less than four lines suitable for RV measurement, crosscorrelation has intrinsic difficulties in providing accurate error bars. Gaussian fitting provides slightly less accurate RVs than crosscorrelation, but with more reliable error bars for stars with only few lines of sufficient quality. However, the performance of Gaussian fitting degrades strongly with the quality of the lines (i.e., low S/N or high vsini).
To improve the performance of the Gaussian fitting, we developed a tool to allow us to fit all the available lines at all epochs (and both wavelength setups) simultaneously. We assume that the line profiles are constant and well described by Gaussians. Each considered line is required to have the same fullwidth at halfmaximum (FWHMs) and amplitude (A) at all epochs. All the lines of a given epoch are further required to display the same RV shift. If L is the number of considered lines and N the number of epochs, the number of parameters in the fit is thus : N + 2 × L. Typical values for N are 18 and 6, depending on whether consecutive exposures are taken individually or summed up. Values for L vary between 2 and 5. This yields a maximum number of parameters of the order of 28 in the present approach, to be compared with the 3 × N × L = 270 parameters of the standard approach that adjusts each line separately. In essence the proposed method allows the line profile parameters (FWHMs and amplitudes) to be bootstrapped by the best quality spectra, improving the RV measurements of the lowest S/N data.
For simplicity and robustness in the SB2 cases (see below), we have chosen to fit the full line profile using Gaussians. Intrinsic He profiles are however not always well represented by Gaussian shapes. For slow rotators (v_{rot}sini ≤ 80 km s^{1}), the shape of the line is dominated by the instrumental profiles that is well approximated by a Gaussian. The profile of moderately rotating stars (v_{rot}sini ≤ 300 km s^{1}) is also well matched by a Gaussian, but with a small difference in the line wings. For fast rotators, rotationally broadened line profiles deviate significantly from the Gaussian shape. However, because we only attempt to measure the position of the centre of the lines and because the fitted lines are symmetric to within firstorder, no significant bias is expected from this approximation.
The fitting method is based on leastsquare estimates where the merit function accounts for the error spectrum to provide individual error bars for each pixel. The optimisation relies on a modified LevensberqMarquardt algorithm that requires a first estimate of the parameters. For constant RV stars and singlelined (SB1) binaries, the method is found to be robust against the initial guesses and we have taken v = 270 km s^{1}, FWHM = 3.0 Å and A = 0.1 to initialize the fit. The allowed ranges for the various parameters are as follow: v ∈ [−300, + 700] km s^{1}, A ∈ [0.0,0.6] and FWHM ∈ [0,10] Å. Only absorption lines were considered. For the doubledlined spectroscopic binaries (SB2) in our sample, we follow the same approach with two Gaussian components per line profile. For some of the objects the initial guesses needed to be iterated once or twice to improve the quality of the SB2 fit.
Appendix B.2: Line selection
In this section, we discuss the choice of spectral lines suitable for relative and absolute RV measurements. The spectra of Otype stars are dominated by H i, He i and He ii lines. Hydrogen lines are generally poorer RV indicators than (most of) the He lines because they are broader and more sensitive to wind effects. The hydrogen lines of most of the VFTS O stars are also strongly affected by nebular emission. Metal lines (C, N, O, Si, Mg) are also present in O star spectra but are typically weaker than He lines. The presence of a given metal line is often limited to a small range of spectral subtypes and does not allow for an homogeneous approach across the O star domain. Our analysis is thus focused on the strong He i and He ii lines from 4000 to 5000 Å, i.e. He iλλ4387, 4471, 4713, 4922, He iiλλ4200, 4541, 4686 and the He i+ii blend at λ4026. In the rest of this section, we discuss the respective merits of each of these lines. Two aspects are important to consider:

the internal consistency, i.e. RV measurements from various lines of a given star should, coherently, reflect the systemic velocity of the star;

the homogeneity with respect to the full sample, i.e. the selected set of lines should provide RV measurements that can be compared to RVs from other stars even if only a subsample of the lines are available. This is important as the spectrum of early O stars is dominated by He ii lines and those of late O stars by He i lines.
Of all helium lines, stellar winds mostly affect He iiλ4686, though not all spectral subtypes or luminosity class are equally affected by the filling in of the photospheric absorption line by wind emission. We find that RVs obtained from He iiλ4686 are reliable for late and mid mainsequence subtypes. They are biased towards negative values by up to 25 km s^{1} in early mainsequence stars. RVs from giants and supergiants are blueshifted even for lateO stars and the line cannot be used as a probe for absolute RVs.
We also compare RV measurements obtained individually from different lines in a representative subsample of the O stars, covering the complete O spectral subrange. Only stars displaying constant RVs are considered for this test. Because the line is observed both in the LR02 and LR03 setups and remains of a reasonable strength even in the latest O subtypes, we use the He iiλ4541 line as a reference for the comparison. Figure B.1 shows the difference between RVs obtained for a given line and those obtained from He iiλ4541 plotted against the spectral subtype. Temperature effects on He i+iiλ4026 are clearly illustrated. Wind effects turn out to have a limited impact on He iiλ4686 measurements mostly because our sample is dominated by mainsequence stars and because many single supergiants present lowamplitude RV variations and are thus excluded from the comparison.
RVs from the He iλ4471 line show large discrepancies, which we attribute to the combination of the following factors: He iλ4471 is a triplet transition, with a forbidden component. It is quasimetastable, similar to but less extreme than He iλλ3889, 5876, 10830 and is highly susceptible to wind effects. For late spectral subtypes it is blended with O ii in its blue wing and, among the measured He i lines, it suffers the most from the nebular contamination.
For late spectral types, He iiλ4200 is blended with N iiiλ4200.07 (i.e., about 17 km/s redwards from He iλ4200), which can explain some of the deviations seen in the corresponding panel in Fig. B.1. For late O stars, the uncertainties associated with He iiλ4200 RVs are large anyway and the line has a very limited weight in the final RV measurements obtained from simultaneously fitting all available lines (Fig. B.2).
Fig. B.1 Difference between the RVs measured from individual He lines and those measured from He iiλ4541 as a function of the spectral subtype. The panel at the bottom right corner compares RVs obtained using all available He lines (except He i+iiλ4026 and He iiλ4686) with the He iiλ4541 RVs. The dashed line indicates the average of all comparisons weighted by their uncertainties. Only stars with constant RVs and data points with σ_{ΔRV} < 20 km s^{1} have been included. 

Open with DEXTER 
Fig. B.2 Comparison of the RVs measured using simultaneously all suitable He i and He ii lines (left and right panel, respectively) with the RVs measured using all suitable He lines. The He i+iiλ4026 and He iiλ4686 lines have not been used in the measurements. 

Open with DEXTER 
The best consistency is observed between the lines He iiλλ4200, 4541, He iλλ4387 and 4713. He iλ4922 displays good consistency over the range of O spectral subtypes but provides RVs that are on average shifted by 3 km s^{1} compared to the other lines. We mitigate this by modifying the effective rest wavelength of the transition.
As a conclusion, up to five lines can be used for absolute RV measurements (He iiλ4200, 4541, He iλλ4387, 4713, 4922) depending on the spectral types, S/N and nebular contamination. He i+iiλ4026, and, for relatively weak wind stars, He iiλ4686 can be used for relative measurements.
Rest wavelengths (λ_{0}) used to computed the RVs.
Appendix C: A MonteCarlo method to constrain the intrinsic orbital parameter distributions
Fig. C.1 Comparison between the observed (crosses) and simulated (lines) cumulative distributions of the peaktopeak RV amplitudes (top row), variability time scales (middle row) and χ^{2} (bottom row). Lefthand panels vary the exponent π of the period distribution. Panels in the central column vary κ, the exponent of the massratio distribution and the righthand panels vary η, the exponent of the eccentricity distribution. π, κ and η values are indicated in the upperleft legend of each panel. The bottomright legends indicate the overall VFTS detection probability for the adopted parent distributions. 

Open with DEXTER 
In this section, we provide the details of the MonteCarlo method used to constrain the intrinsic orbital parameter distributions. The general philosophy resembles that of Kobulnicky & Fryer (2007), in the sense that it relies on a MonteCarlo modeling of the orbital velocities. On the one hand, we follow Kobulnicky & Fryer in using KS tests to compare simulated and observed distributions and the diagnostic distributions based on ΔRV and χ^{2} (see Sect. C.3) are similar to, respectively, the V_{h} and V_{rms} used by Kobulnicky & Fryer. On the other hand, our method fundamentally differs in various aspects. First we use the specific measurement errors provided by the RV analysis throughout the complete MonteCarlo analysis. Second, we do not apply external constraints such as fixing the period distribution or using the fraction of Type Ib/c supernovae corresponding to each model. We thus attempt to simultaneously constrain the intrinsic period and massratio distributions and the intrinsic binary fraction.
The core of the method is identical to that described in Sana et al. (2012) but differs in the adopted merit function. Sana et al. compared the observed orbital parameter distributions while we will use distributions constructed from the RV database. Indeed most of the detected binaries in 30 Dor have too few RV measurements to compute a meaningful orbital solution, precluding the direct approach used by Sana et al. (2012).
Appendix C.1: Method overview
As discussed in the main text, we use power laws to describe the intrinsic distributions of orbital parameters: f(log _{10}P/d) ~ (log _{10}P/d)^{π} , f(q) ~ q^{κ} and f(e) ~ q^{η}, the exponents of which are left as free parameters. For each combination of π, κ and f_{bin} in our grid, we draw populations of N = 360 primary stars using a Kroupa mass function (Kroupa 2001) between 15 and 80 M_{⊙}, covering thus the expected mass range of Otype stars. Each star is assigned an observational sequence (observing epochs and RV accuracy at each epoch) from our VFTS sample. A fraction of these stars are randomly assigned to be binaries following a binomial with a success probability f_{bin}. The orbital parameters of the binaries are randomly drawn as follows:

the period and massratio are taken from the powerlaw distributions given π and κ;

–
the eccentricity is taken from a powerlaw distribution with an index given by η = −0.5 (see discussion in Sect. 4.1);

the inclination and periastron argument are taken assuming random orientation of the orbit in space;

the time of periastron passage is assumed to be uncorrelated with the observational sampling and is thus taken from a uniform distribution.
The binary detection criteria of Eq. (4) are then applied and the simulated detected binaries are flagged. Simulated observed distributions (see Sect. C.3) are built based on the simulated RV measurements of the detected systems and the simulated measured binary fraction is computed.
For a given combination of π, κ and f_{bin} the process is repeated 100 times to build the parent statistics. The simulated parent distributions are then compared to the observed ones by means of a onesided KS test. We also compute the binomial probability to detect N_{bin} binaries among a population of N stars given the success probability predicted by the simulations, i.e. . Adopting N_{bin} to be the actual number of detected Otype binaries in VFTS allows us to estimate the ability of the assumed intrinsic parameters to reproduce the observed binary fraction while taking into account the finite size of the sample.
Appendix C.2: Diagnostic distributions
In this section, we specifically discuss the empirical distributions used as diagnostics. The choice of the merit function will be discussed in the next section. We investigate four diagnostic quantities:

i.
the distribution of, for each detected binary, the amplitude of the RV signal (ΔRV);

ii.
the distribution of, for each detected binary, the χ^{2} of a constant RV model;

iii.
the distribution of, for each detected binary, the minimum time scale for significant RV variations to be seen (ΔHJD);

iv.
the detected binary fraction.
Appendix C.3: Suitable merit functions
At each point of the grid defined by the the intrinsic values π, κ and f_{bin}, the simulated distributions (i) to (iii) are individually compared to the observed ones using KStests while the simulated and observed binary fraction are compared using a binomial statistics. Because we aim at matching all the diagnostic distributions simultaneously, we also combine the individual probabilities in a global merit function.
We run several test to assess different merit functions (i.e., different ways of combining the individual probabilities) as well as to validate our method. We generate synthetic data from known input distributions, i.e. known π, κ and f_{bin}, and we apply our code using different merit functions. Input and output distributions are then compared to check the suitability of the various merit functions as well as the general accuracy of the method. The process has been repeated using 50 synthetic data sets in each case. The synthetic data have been generated such that they share the same observational properties (sampling, measurement accuracy) as the VFTS Medusa data. Table C.1 provides an overview of the results of the tests. Because of the computational load of such tests, we use a slightly coarser grid with steps of 0.02 in f_{bin}. Explored ranges in f_{bin}, π and κ and steps in the power law indexes π and κ are identical to those quoted in Table 5. In addition to appraising the merit function, we also test the need to apply a minimal threshold C for the detected significant RV signal.
Overview of the test results using different merit functions (Col. 2).
The first striking result is that most merit functions can recover the input binary fraction within a few percents albeit with various precision. Some of the best results are obtained when including the B(N_{obs},N,f_{sim}) term in the merit function.
The index π of the period distribution is poorly constrained whenever one uses no detection threshold (i.e., C = 0 km s^{1}). This indicates that the detection threshold is not only useful to reject false detections due to intrinsic profile variability but is also a critical ingredient of the method. Interestingly, π is reasonably well constrained by the ΔHJD distribution alone, i.e. merit function (l), but a better overall result is obtained when used in conjunction with other diagnostic quantities.
The ΔRV and χ^{2} distributions are useful to refine the κ value but the uncertainties remain large. In essence, the ΔRV
and χ^{2} distributions carry similar information. Because merit functions using ΔRV tend to behave slightly better than those using χ^{2}, and because the χ^{2} values are by definition much more sensitive to the exact knowledge of the measurement errors, we select our final merit function as given by the product of the P_{KS} probabilities computed from the ΔHJD and ΔRV distributions and of the Binomial probability (Eq. (5)), hence merit function (u) in Table C.1.
Overall, the retrieved values of π and κ tend to be smaller than their input values but the correct indexes concur with the 1σ confidence intervals quoted in Table C.1. These confidence intervals can further be used to estimate the formal uncertainty of the method. Given that the best representation of the data are obtained with f_{bin} = 0.51, π = −0.45 and κ = −1.00 (see Sect. 4.2), the bottom part of Table C.1 applies and we adopt σ_{fbin} = 0.04, σ_{π} = 0.3 and σ_{κ} = 0.4.
Appendix D: Supplementary figures
Fig. D.1 K_{s}band magnitudes as a function of the spectral subtypes for luminosity classes VIV (left panel), III (middle panel) and III (right panel). Different symbols identify the single stars, the lowamplitude RV variables and the binaries. Symbols for single and binaries of a given spectraltypes have been shifted, compared to one another, by a small amount amount along the xaxis to improve the clarity of the plots. 

Open with DEXTER 
Fig. D.2 Same as Fig. 11 for the luminosity classes VIV (left panels), III (middle panels) and III (right panels). 

Open with DEXTER 
Fig. D.3 Same as Fig. 12 for the luminosity classes VIV (left panels), III (middle panels) and III (right panels). 

Open with DEXTER 
All Tables
Journal of the observations, RV measurements and their 1σ errorbars for the stars in our sample.
Overview of the results of various variability and multiplicity criteria used in this paper.
Observed fraction of RV variables and binaries among different O star populations in 30 Dor.
Number of objects, observed binary fraction detection probability and intrinsic binary fraction as a function of the K_{s} magnitude.
Observed fractions of lowamplitude RV variables and binaries as a function of the luminosity class.
All Figures
Fig. 1 Top panel: number of RV epochs obtained for each object. One Argus object, VFTS 1022, has been observed at 13 epochs and falls outside the plotted range. Middle panel: number of lines available for absolute RV measurements. Bottom panel: time elapsed between the first and last suitable epochs for RV measurement. The solid line corresponds to objects observed with Medusa while the dashed lines indicate objects observed with the Argus IFU. 

Open with DEXTER  
In the text 
Fig. 2 Distribution of the Otype stars in the field of view. The population is dominated by the central cluster, NGC 2070, that contains 57% of our sample. NGC 2060, at 6.7′ to the southwest, is slightly older and accounts for 22% of the sample, while the remaining 21% is spread throughout the field of view. The dashed circles, with a 2.4′radius each, show the adopted regions for NGC 2070 and NGC 2060 (see Table 6). 

Open with DEXTER  
In the text 
Fig. 3 Variations of the fraction of systems below and above the critical RV variation amplitude C separating the spectroscopic binaries (solid line) from lowamplitude RV variables (dotted line). The vertical dasheddotted line indicates the adopted threshold of C = 20 km s^{1}. 

Open with DEXTER  
In the text 
Fig. 4 Histogram of the peaktopeak amplitude RV variations of each object. Only significant variations according to Eq. (4) are considered. Therefore the first bin is empty. 

Open with DEXTER  
In the text 
Fig. 5 Cumulative distribution of the minimum time difference ΔHJD for a binary to display significant RV variations according to Eq. (4) and for different RV amplitude thresholds (C) in the RV variation amplitudes. The cumulative distributions are normalised to the total number of binaries. 

Open with DEXTER  
In the text 
Fig. 6 Projections of the global merit function Ξ′ on the various pairs of planes defined by π, κ and f_{bin}. The absolute maximum is indicated by a cross (+). Contours indicate loci of equalvalues of the merit function, expressed as a fraction of its absolute maximum (see legend). 

Open with DEXTER  
In the text 
Fig. 7 Comparison between the observed (crosses) and simulated (lines) cumulative distributions of the peaktopeak RV amplitudes (top row) and of the variability time scales (bottom row). In the left (resp. right) hand panels, we vary the exponent π of the period distribution (resp. κ of the massratio distribution) by ±1σ. The upperleft legends indicate the values of π, κ and η considered in each panel. The bottomright legends give the overall VFTS detection probability for the adopted parent distributions. 

Open with DEXTER  
In the text 
Fig. 8 Overview of the binary detection probability achieved by the VFTS campaign as a function of the orbital period, the mass ratio, the eccentricity and the primary mass. (See Sect. 4.3 for more details.) 

Open with DEXTER  
In the text 
Fig. 9 Cumulative radial distribution of the observed binary fraction as a function of the distance to the centre of the field. The black solid line shows the distribution computed insideout, and the grey/red line shows the distribution computed outsidein. The dashed curves indicate the ± 1σ confidence interval on both curves. 

Open with DEXTER  
In the text 
Fig. 10 From left to right, distributions of the RV constant, likely single stars, the binary candidates (lowamplitude RV variables) and the binary stars in the 30 Dor field of view. 

Open with DEXTER  
In the text 
Fig. 11 Top panel: number of objects (dashed line) and of binaries (solid line) per K_{s}band magnitude bin. Bottom panel: observed binary fraction in corresponding magnitude bins (solid line) and intrinsic binary fraction (dashed line). Error bars are plotted whenever a bin contains more than one star. Figure D.2 provides similar figures showing the dependence with luminosity class. 

Open with DEXTER  
In the text 
Fig. 12 Top panel: number of objects (dashed line) and of binaries (solid line) per O spectral subtype. Bottom panel: observed binary fraction as a function of the spectral subtypes. Error bars are plotted whenever a bin contains more than one star. Figure D.3 provides similar figures for the different luminosity classes. 

Open with DEXTER  
In the text 
Fig. 13 Spatial distribution and multiplicity properties of the supergiants in the field of view. 

Open with DEXTER  
In the text 
Fig. 14 Massratio distributions obtained using different pairing mechanisms (solid and dashed lines) compared to those obtained with a powerlaw density distribution with different κ values (dotted lines). 

Open with DEXTER  
In the text 
Fig. B.1 Difference between the RVs measured from individual He lines and those measured from He iiλ4541 as a function of the spectral subtype. The panel at the bottom right corner compares RVs obtained using all available He lines (except He i+iiλ4026 and He iiλ4686) with the He iiλ4541 RVs. The dashed line indicates the average of all comparisons weighted by their uncertainties. Only stars with constant RVs and data points with σ_{ΔRV} < 20 km s^{1} have been included. 

Open with DEXTER  
In the text 
Fig. B.2 Comparison of the RVs measured using simultaneously all suitable He i and He ii lines (left and right panel, respectively) with the RVs measured using all suitable He lines. The He i+iiλ4026 and He iiλ4686 lines have not been used in the measurements. 

Open with DEXTER  
In the text 
Fig. C.1 Comparison between the observed (crosses) and simulated (lines) cumulative distributions of the peaktopeak RV amplitudes (top row), variability time scales (middle row) and χ^{2} (bottom row). Lefthand panels vary the exponent π of the period distribution. Panels in the central column vary κ, the exponent of the massratio distribution and the righthand panels vary η, the exponent of the eccentricity distribution. π, κ and η values are indicated in the upperleft legend of each panel. The bottomright legends indicate the overall VFTS detection probability for the adopted parent distributions. 

Open with DEXTER  
In the text 
Fig. D.1 K_{s}band magnitudes as a function of the spectral subtypes for luminosity classes VIV (left panel), III (middle panel) and III (right panel). Different symbols identify the single stars, the lowamplitude RV variables and the binaries. Symbols for single and binaries of a given spectraltypes have been shifted, compared to one another, by a small amount amount along the xaxis to improve the clarity of the plots. 

Open with DEXTER  
In the text 
Fig. D.2 Same as Fig. 11 for the luminosity classes VIV (left panels), III (middle panels) and III (right panels). 

Open with DEXTER  
In the text 
Fig. D.3 Same as Fig. 12 for the luminosity classes VIV (left panels), III (middle panels) and III (right panels). 

Open with DEXTER  
In the text 