A relation between the radial velocity dispersion of young clusters and their age Evidence for hardening as the formation scenario of massive close binaries

The majority of massive stars ( > 8 M (cid:12) ) in OB associations are found in close binary systems. Nonetheless, the formation mechanism of these close massive binaries is not understood yet. Using literature data, we measured the radial-velocity dispersion ( σ 1D ) as a proxy for the close binary fraction in ten OB associations in the Galaxy and the Large Magellanic Cloud, spanning an age range from 1 to 6Myr. We ﬁnd a positive trend of this dispersion with the cluster’s age, which is consistent with binary hardening. Assuming a universal binary fraction of f bin = 0.7, we converted the σ 1D behavior to an evolution of the minimum orbital period P cuto ﬀ from ∼ 9.5years at 1Myr to ∼ 1.4days for the oldest clusters in our sample at ∼ 6Myr. Our results suggest that binaries are formed at larger separations, and they harden in around 1 to 2Myr to produce the period distribution observed in few million year-old OB binaries. Such an inward migration may either be driven by an interaction with a remnant accretion disk or with other young stellar objects present in the system. Our ﬁndings constitute the ﬁrst empirical evidence in favor of migration as a scenario for the formation of massive close binaries.

The first effort to characterize the binarity properties of a sample of O stars in compact H ii regions was performed by Apai et al. (2007). They did a multi-epoch (two to three epochs) radial velocity (RV) study of a sample of 16 embedded O stars in seven massive star-forming regions. They identified two close binary stars based on their RV variations (∼90 km s −1 ) and measured an RV dispersion (σ 1D ) of 35 km s −1 for the whole sample and 25 km s −1 when excluding the two close binaries. After pioneering studies to spectroscopically characterize single massive young stellar objects (mYSOs) such as those carried out by Bik et al. (2006Bik et al. ( , 2012, Ochsendorf et al. (2011), Ellerbroek et al. (2013), and Ramírez-Tannus et al. (2017) performed a single-epoch VLT/X-shooter spectroscopic study of a sample of eleven candidate mYSOs in the very young giant H ii region M 17 ( 1 Myr). The stars range in mass from 6−25 M which is the mass range that dominates the samples from which multiplicity characteristics of 2−4 Myr old main-sequence OB stars are derived Kobulnicky et al. 2014). The measured radial-velocity dispersion of these mYSOs is σ 1D = 5.6 ± 0.2 km s −1 . In low density clusters, such as M 17, σ 1D of a single epoch is strongly dominated by the orbital properties of the binary population. For example, if a given cluster has several close binaries of similar masses, one would expect the individual radial velocities of the stars to differ significantly from each other and, therefore, for σ 1D to be large. For 2−4 Myr clusters, with binary fractions >0.5 and minimum periods of ∼1.4 days, a dispersion of 30 to 50 km s −1 is typical (e.g., Kouwenhoven et al. 2007;Sana et al. 2008Sana et al. , 2012Sota et al. 2014;Kobulnicky et al. 2014). The latter is in stark contrast with our observation of M 17, suggesting a lack of close massive binaries in this region.
In Sana et al. (2017), two scenarios are explored that may explain the low σ 1D observed in M 17: a small binary fraction A&A 645, L10 (2021) f bin and/or a lack of short-period binaries. They conclude that the observed dispersion can be explained either if f bin = 0.12 +0.16 −0.09 or if the minimum orbital period P cutoff > 9 months. Parent populations with f bin > 0.42 or P cutoff < 47 days can be rejected at the 95% significance level. Since it is unlikely that the binary fraction for M 17 would be so far below that of other clusters, this very interesting result suggests that massive binaries form in wide orbits that migrate inward over the course of a few million years. In this Letter, we refer to the generic mechanism of shrinking binary periods as the migration scenario. One strong test for this scenario is to compare the velocity dispersion observed in clusters spanning a range of ages. If the binary orbits harden with time, one would expect σ 1D to increase as the cluster age increases.
In Ramírez-Tannus et al. (2020), VLT/KMOS spectra of around 200 stars in three very young clusters (M 8, NGC 6357, and G333.6−0.2) were obtained. Introducing an automatic method to classify the spectra, the effective temperatures, and luminosities of the observed stars were characterized in order to place them in the Hertzsprung-Russell diagram (HRD). The age and mass range of the observed populations was constrained by comparison to MESA evolutionary tracks obtained from the MIST project (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015Dotter 2016;Choi et al. 2016). The main sequence stars in M 8 have masses between ∼5 and ∼70 M and the age of this cluster is between 1 and 3 Myr. In G333.6−0.2, the main sequence population ranges in mass between ∼5 and ∼35 M and the estimated age of this region is <3 Myr. The main sequence stars in NGC 6357 have masses between ∼10 and ∼100 M and their ages range from 0.5−3 Myr.
The goal of this Letter is to provide a first test of the migration scenario for the formation of massive close binaries. We aim to study a possible age evolution of σ 1D , and P cutoff , to constrain a timescale for binary hardening assuming a universal binary fraction ( f bin = 0.7; Sana et al. 2012). We base our analysis on clusters younger than 6 Myr to ensure that neither secular evolution nor the effect of binary interactions (Wellstein & Langer 1999;de Mink et al. 2007) affect our results significantly. In Sect. 2 we measure the radial velocities of the high-mass stars in M 8 and NGC 6357 and calculate their σ 1D . Next, we compare our findings with those presented by Sana et al. (2012) for Galactic clusters of 2−4 Myr, with those from Zeidler et al. (2018) for Westerlund 2 (Wd2), with those from Hénault-Brunet et al. (2012) for R136 in the Large Magellanic Cloud, and with those from Ramírez-Tannus et al. (2017) for the very young massivestar forming region M 17. This reveals a temporal behavior of σ 1D (Sect. 3) that is converted into an evolution of the minimum binary period (Sect. 4), as binary motion is dominating the velocity dispersion of young massive clusters. In Sect. 5 we discuss and conclude this work.

Observations
The sample studied in this Letter consists of the OB stars in M 8 and NGC 6357. The data acquisition and reduction are described in detail in Ramírez-Tannus et al. (2020). In short, we obtained around 200 H and K-band intermediate resolution spectra (with spectral resolution power, λ/∆λ, between 6700 and 8500, i.e., 30 < ∆v < 40 km s −1 ) of stars in the abovementioned giant H ii regions with VLT/KMOS (Sharples et al. 2013). The final samples of massive stars consist of 16 stars in M 8, 22 in NGC 6357, and four in G333.6−0.2. We discarded G333.6−0.2 from our analysis because there are not enough stars with RV measurements to calculate σ 1D . A description of the age and mass range determination can be found in Ramírez-Tannus et al. (2020), and Appendix A presents a detailed discussion about the accuracy of the age determination.
The radial velocity (RV) of the intermediate to high-mass stars was obtained by measuring the Doppler shifts of a suitable set of photospheric lines. Tables B.1 and B.2 list the RV obtained for each star together with its error and the spectral lines used in our analysis.
The RV-fitting approach is similar to the one adopted by Sana et al. (2017). First, for the profile fitting, we adopted Gaussian profiles. Second, we clipped the core of diagnostic lines that were still contaminated by residuals of the nebular emission. Third, we simultaneously fit all spectral lines available, thereby assuming that the Doppler shift is the same for all lines (see Sect. 2 and Appendix B of Sana et al. 2013). Figure C.1 shows the radial-velocity distribution for the two regions. We calculated the errors of the histogram bins by randomly drawing RV values from a Gaussian centered at each measured RV and with a sigma corresponding to the measurement error; we repeated that process 10 5 times. The value shown for each bin is the mean of all the RVs in that bin's measurements and the error bar corresponds to the standard deviation. We obtained σ 1D by calculating the weighted standard deviation of the measured RVs. The weighted mean and standard deviation are listed in the top-left corner of each histogram in Fig. C.1. The measured σ 1D for M 8 and NGC 6357 are 32.7 ± 2.6 and 26.9 ± 1.3 km s −1 , respectively.

Velocity dispersion versus cluster age
Based on single-epoch radial-velocity measurements of young massive stars in M 17, Sana et al. (2017) conclude that this young (∼1 Myr) star-forming region hosts only a few close binary systems. This is in contrast to the observation that most massive stars in somewhat older clusters are in close binaries. Under the hypothesis that massive stars form in wide binary systems that harden their orbits within the first million years of evolution, one would expect that during those initial couple of million years the radial-velocity dispersion, σ 1D , increases with time.
We plotted σ 1D of the clusters studied in this Letter and compare it with the results by Sana et al. (2012, NGC 6231, IC 2944, IC 1805, IC 1848, Zeidler et al. (2018, Wd2), Hénault-Brunet et al. (2012, R136), and Sana et al. (2017, M 17). R136 and Wd2 are very relevant for our study given their relatively young age (1−2 Myr) which is in between the age measured for M 17 ) and that of the somewhat older clusters . In order to compare the multi-epoch RV data provided by Sana et al. (2012) with the single epoch data of M 17, M 8, NGC 6357, and Wd2, we drew the RV measured for each star in a given cluster in a random epoch and we computed the RV dispersion. We repeated this procedure 10 5 times and then calculated the most probable σ 1D and its standard deviation. The σ 1D obtained for each cluster is listed in the third column of Table 1. The second and seventh columns show the age of the clusters and the respective references.
In Fig. 1 we show σ 1D versus age of the clusters. We performed an orthogonal distance regression (ODR ;Churchwell 1990) to the data and find a positive correlation between the age of the clusters and σ 1D . The solid black line represents the best fit to the data and the gray area represents the 1-σ error on the fit. The Pearson coefficient for the observed relation is 0.7, which indicates a strong positive correlation. Nevertheless, this coefficient does not take into account the errors in the parameters. To test the validity of our results we performed two Monte Carlo tests whose results are shown in Fig. 2. The left panel shows the distribution of Pearson coefficients obtained from drawing random points centered on our data (age, σ 1D ) with a standard deviation equal to our error bars. The right panel shows the L10, page 2 of 7 M. C. Ramírez-Tannus et al.: A relation between the radial velocity dispersion of young clusters and their age  probability (2%) that a random distribution causes the observed coefficient. Even though the scatter is substantial, we conclude that for the present data there is a positive correlation between σ 1D and the age of the clusters. Although the number of clusters considered here remains limited and there are sizable uncertainties as to the individual age determinations, our results indicate an increase in the fraction of close binary systems as the clusters get older. M 17 seems to be a unique cluster given its very young age and extremely low σ 1D . Trumpler 14 in this sense seems comparable, being ∼1 Myr old and having no evidence for close binaries amongst its O stars (Sana & Evans 2011). We do not include this cluster in the analysis because there are not enough stars with RV measurements to calculate σ 1D . Pearson r 2% Fig. 2. Distribution of Pearson coefficients after randomly drawing 10 5 samples from the data shown in Fig. 1. Left: 2D Gaussians centered at our data points and with σ equal to our error bars. The green, red, and blue areas show the 68, 95, and 99% confidence intervals, respectively. The black line shows the observed coefficient. Right: 2D Gaussians centered at random locations of the parameter space (age between 0 and 7 Myr and σ 1D between 0 and 80 km s −1 ) and with σ equal to our error bars. The blue shaded area represents the probability (2%) that the observed coefficient is caused by a random distribution.

Physical parameters
Our aim is to characterize the multiplicity properties of observed clusters based on the observed σ 1D . This section presents the results of Monte Carlo population syntheses computed with different underlying multiplicity properties.

Binary fraction and minimum period
We focus on the effect that the binary fraction, f bin , and minimum orbital period, P cutoff , have on the observed σ 1D . The methodology is similar to that used in Sana et al. (2017). A parent population of stars is generated with a certain f bin and P cutoff . The binary star systems in this population are described by their orbital properties. These are the primary mass, M 1 , mass ratio, q, period, P, and eccentricity, e. For the multiplicity properties that we do not vary, we adopt the values from Sana et al. (2012) for Galactic young clusters as the basic properties for our population. For the primary star, we adopt a Kroupa mass distribution (Kroupa 2001). The mass ratio distribution is uniform with 0.1 < q < 1. The probability density function of the period is described as pdf(P) ∝ (log P) −0.5 , with the period in days and log P cutoff < log P < 3.5. The eccentricity distribution depends on the period of the binary system. For P < 4 days, we assume circular orbits; for periods between 4 and 6 days, the eccentricities are sampled from pdf(e) ∝ e −0.5 , with 0 ≤ e < 0.5; for periods longer than 6 days the same distribution is used, but with 0 ≤ e < 0.9. In order to calculate the radial velocity of the binary star systems, they are given a randomly generated inclination, i, longitude of periastron, ω, and eccentric anomaly, E. The latter is determined by generating a random mean anomaly and numerically solving Kepler's equation to find the corresponding eccentric anomaly using Brent's method (Brent 1973) 1 . This gives all of the information required to calculate the binary component of the radial velocity of the primary star through with K 1 being the amplitude of the orbital velocity of the primary star. The binary component of the radial velocity is added to the velocity due to cluster dynamics, σ dyn . This results in a population of stars with radial velocities based on either cluster dynamics only (in the case of f bin = 0) or both cluster dynamics and binary orbits. The secondary stars are assumed to be undetected. We adopt a cluster velocity dispersion for all clusters of 2 km s −1 , which corresponds to the typical value found by Kuhn et al. (2019) who measured σ dyn for several massive clusters in the Milky Way. Adopting σ dyn = 1 or σ dyn = 5 km s −1 has no significant effect on our results. Sana et al. (2012) also conclude that the dynamical dispersion of young clusters is negligible with respect to the dispersion due to binary motions. Each Monte Carlo run consists of a generated parent population of 10 5 stars, which is sampled 10 5 times to simulate the observation of a cluster. The mass distribution of the stars in the parent distribution are matched to the mass distribution in the observed cluster. The samples contain the same number of stars and measurement accuracy as the observed sample. We then constructed density distributions of σ 1D for each parent population. Each grid point in Fig. 3 shows the median σ 1D obtained for a simulated cluster of 20 stars sampled from a parent population with masses between 15 and 60 M , where we varied f bin , and P cutoff . The color bar corresponds to σ 1D in km s −1 . The dashed lines show the σ 1D contours for 5.5, 15, 25.3, 30.9, 50.3, and 66 km s −1 , which correspond to the σ 1D measured for our clusters. As the number of stars and the mass ranges are different for each cluster, the contours do not represent the exact way in which we measured P cutoff but are meant to show examples of the trends that a certain σ 1D follows in this diagram. A comparison of our results with Sana et al. (2017) is shown in Appendix D.

Timescale of binary hardening
Assuming that the binary fraction of massive stars is consistent with that observed in OB stars in young open clusters ( f bin = 0.7; Sana et al. 2012), but that binary stars are born in wider orbits, we can estimate P cutoff which best represents the observed σ 1D for each cluster. For each cluster, we kept the distribution of orbital properties fixed and we adjusted the sample size and Wd2 P cutoff (t)=10 (5 ± 1) exp( t/(0.19 +0.06 0.04 )) + 1.40 Fig. 4. P cutoff versus age of the clusters. The error bars represent the P cutoff corresponding to the distributions that represent each σ 1D within their 68% confidence range. The blue line and shaded region show the fit to the data and its 1-σ error bars. mass ranges in accordance to those of the observed samples (see Table 1). This allowed us to make a similar plot as in Fig. 1 but in terms of P cutoff . Figure 4 shows the estimated P cutoff for each cluster as a function of cluster age. We determined that, assuming f bin = 0.7, the most likely P cutoff to explain the observed σ 1D of M 17 is around 3500 days and that a P cutoff of 665 days and 126 days can be rejected at the 68% and 95% confidence levels, respectively. For the clusters with the largest σ 1D in our sample, IC 1805, IC 1848, and NGC 6231, the cutoff period that best explains the observed σ 1D is 1.4 days, which is the P cutoff adopted by Sana et al. (2012).
In order to get a first estimate of the binary hardening timescale, we fit a function of the form P(t) = P 0 e −t/t 0 + c to the data in Fig. 4, where P 0 is the minimum period at the moment of binary formation and t 0 corresponds to the e-folding time. As P 0 is very uncertain, we assumed a typical value of 10 5 days, corresponding to an initial separation of ∼100 AU for a pair of 10 M stars. In order to find the range of t 0 , we varied P 0 from 10 4 to 10 6 days. The resulting fit and corresponding parameter ranges are shown with the blue line and shaded area in Fig. 4. We obtain a typical e-folding time of t 0 = 0.19 +0.06 −0.04 Myr. To harden an orbit from 3500 to 1.4 days, ∼8 e-foldings are needed, this implies that a typical binary hardening timescale is ∼1.6 Myr. Assuming a binary system with two 10 M stars, this would imply a mean hardening rate of ∼7.5 AU yr −1 . For M 17, we discarded periods lower than 665 d with a 68% confidence level. Following the same argument as before, the typical timescale to harden an orbit from 665 to 1.4 days would be ∼1 Myr equivalent to a mean hardening rate of ∼4 AU yr −1 .

Discussion and conclusions
In this Letter we present observational evidence that the 1D velocity dispersion of massive stars in young clusters (σ 1D ) increases as they get older (Fig. 1). Additionally, we performed Monte Carlo simulations which allowed us to convert the measured values of σ 1D to physical parameters. Assuming that stars are born with the binary fraction representative of OB stars in 2−4 Myr clusters ( f bin = 0.7), we calculated the cutoff period that would correspond to each of the observed populations. From Fig. 4 we can conclude that the orbits would harden L10, page 4 of 7 in 1−1.5 Myr. It is worth pointing out that the fit relies heavily on the M 17 point. If the binary fraction of this cluster were lower than 0.7, the timescales derived would be larger. Nevertheless, the qualitative results of this Letter remain the same. Sana et al. (2017) propose that the inward migration process could be driven either by the interaction with the remnants of an accretion disk or with other young stellar objects. The hardening of the orbits would then stop either when the disk is destroyed, or when the third body is pushed far out or ejected from the orbit. For a typical system of two 10 M stars, a total angular momentum of ∼6.5 × 10 47 kg m 2 s −1 (∼3.5 × 10 47 kg m 2 s −1 ) needs to be transferred for the binary to harden from 3500 (665) to 1.4 days. Rose et al. (2019) explored the possibility of the orbit hardening via the Eccentric Kozai-Lidov (EKL) mechanism, where a third companion perturbs the orbit of a binary system. They find that, beginning with a cutoff period of 9 months the EKL mechanism is insufficient to reproduce the population of short period binaries observed by Sana et al. (2012). They suggest that type II migration (Lin & Papaloizou 1986) might explain the tightening of binaries in such a time period. Moe & Kratter (2018) find that the main mechanism to harden binaries should be the dynamical disruption of coplanar triples that initially fragmented in the disk in combination with energy dissipation within this disk. Other mechanisms include the combination of EKL oscillations with tidal friction both during the pre-main sequence and the main sequence (Bate et al. 2002;Bate 2009). Given the densities of the clusters studied, stellar encounters are uncommon and, therefore, mechanisms such as binary-binary or single-binary interactions should not play a significant role in the formation of close binaries (Fujii & Portegies Zwart 2011).
Our findings affect predictions of binary population synthesis models that follow the evolution of an ensemble of binary systems subject to, among others, a distribution function for the zero-age orbital periods. Such models (e.g., Sana et al. 2012;Schneider et al. 2015) predict that a small fraction of systems interact within a timescale of 3 Myr, almost invariably resulting in a merging of the two components thus creating a blue straggler. Specifically, Schneider et al. (2015) report that of the population within two magnitudes in brightness of the main-sequence turnoff, only 1%, 2%, and 10% is a blue straggler after 1, 2, or 3 Myr. Therefore, unless one is specifically interested in early blue straggler formation, the temporal evolution of orbital properties in the first few million years reported here should not affect predictions made by binary population synthesis models that rely on initial conditions defined at the zero-age main sequence only. Except for the caveat mentioned, adopting initial conditions for the massive star population as reported on by Sana et al. (2012), for example, remains a justified approach.