AN UPDATED LOOK AT BINARY CHARACTERISTICS OF MASSIVE STARS IN THE CYGNUS OB2 ASSOCIATION

While our current approach explores a larger parameter space than KF07, incorporates eccentric orbits, and utilizes the information of 22 close massive binary systems in addition to the radial velocity information of 12 years of spectroscopic observations, there are still several limitations stemming from our modeling approach, our data set, and our current knowledge of Cyg OB2.

Limitations owing to our modeling approach. As in KF07, we still consider all velocity variation as stemming from orbital dynamics rather than atmospheric phenomena such as pulsations. Clark et al. (2010) show that pulsations in evolved early-type stars can significantly affect radial velocity measurements, meaning that this effect may add scatter to radial velocity measurements in roughly one-quarter of our sample. We also still treat higher order systems as binaries. At present, we know of three triple or higher order systems, Cyg OB2 No.5, Cyg OB2 No.8, and MT429 (see references in Table 6 of Paper IV). In addition, we limit characterization of each orbital parameter distribution to a single power law, an assumption that may not hold.

Limitations owing to the data set. As mentioned, the Cyg OB2 radial velocity survey probes binary systems with periods ranging from one day to a few thousand days, meaning we are not sensitive to extremely close (and extremely rare) systems with P < 1 day. Therefore, the MC code only draws from period distributions with periods ranging from 1 to 1000 days. Additionally, the vast majority of spectroscopic observations were obtained June through November when Cygnus is visible to northern latitudes. This, coupled with typical velocity uncertainties of 9.7 km s⁻¹, reduces our sensitivity to long-period, high eccentricity systems (i.e., systems that may pass through periastron and exhibit their largest velocity variation during an unobservable period), systems with extreme mass ratios, and systems with small inclination angles. Even though we have detected and computed solutions for systems with semi-amplitudes as small as 9 km s⁻¹ (Paper VI), we attempt to circumvent these limitations by designating systems in the code with semi-amplitudes smaller than 15 km s⁻¹ as "single."

Another limitation comes in the form of not being able to detect unresolved SB2s. For systems exhibiting highly blended line profiles (even at quadrature), such as long-period binaries with mass and luminosity ratios near unity, radial velocity variations of the blended line profile will be small and possibly undetectable at the resolution of the spectra in this survey, thereby being designated as "single" by the code. Sana et al. (2009) discuss this effect in their analysis of the massive star binary fraction of NGC 6611 and conclude that the effect is minimal (affecting their modeling of the binary fraction by only a few percent). Additionally, there is also a negligible bias toward detecting short-period, high-mass-ratio systems that show up as SB2s as it only takes one observation at quadrature to determine that the system warrants follow-up. These systems, however, have large velocity separations (i.e., large semi-amplitudes) and will get flagged as a "binary" by the code with only a few observations.

Limitations owing to incomplete knowledge of Cyg OB2. A few studies have shown that Cyg OB2 may be contaminated by an older population of stars, meaning that evolutionary effects may be present (e.g., Wright et al. 2010; Drew et al. 2008; Comerón et al. 2008). However, Comerón et al. (2008) point out that the field star contamination is likely composed of F to late B stars at a radius of 1°. This exceeds the radial extent of the region examined in this study and concerns stars of a later type than are included here. On the other hand, some doubt remains in the membership of the 114 OB stars included in this work as recent findings by Stroud et al. (2010) and Linder et al. (2009) provide conflicting distance estimates to Cyg OB2 based on their complete orbital solutions calculated for the eclipsing SB2s, Cyg OB2 No.5 (∼1 kpc), and B17 (1.5–1.8 kpc), respectively. We maintain that the contamination is minor, if it exists, owing to the localized region we examine, but even the removal of a few periods or mass ratios could influence the distributions we observe at present. Similar to this is the assumption that the members of Cyg OB2 have a common age. It has been suggested that even within the core region of Cyg OB2, there has been more than one epoch of star formation (Wright et al. 2010). There are no clear predictions for how a population of close binaries should evolve over the 1–5 Myr lifetime of most massive stars in young clusters, so we have no way to model such effects or even know if they are important.

Binary surveys that examine the orbital period distribution probe widely varying ranges in period/separation and stellar mass, and they find widely varying characterizations of its functional form. For solar-type stars with periods from 1 to 10,000+ days, Duquennoy & Mayor (1991) find a log-normal distribution. Among massive stars, most studies find a uniform distribution in log space (β = 0), sometimes known as an Öpik (1924) distribution. For example, Garmany et al. (1980) find β ∼ 0 for a large number of Galactic O-star binaries in the northern hemisphere. Kouwenhoven et al. (2007) also find an Öpik's distribution for a survey of visual and spectroscopic binaries in Sco OB2 composed of A- and B-type primaries (for 0.7 days ⩽P ⩽ 3 × 10⁸ days). Zinnecker & Yorke (2007) speculate that rather than a simple power-law distribution for O-star binaries, there may be a concentration of systems having periods between two and five days. The O-star binaries of NGC 6231 show this trend (Sana et al. 2008). However, Sana & Evans (2011) conclude that a broken Öpik distribution provides a better description of the period distribution for spectroscopic O-binaries in nearby clusters, where the broken distribution could better explain the higher concentration of systems with periods less than 10 days.

The observed orbital period, eccentricity, and mass ratio distributions for the 22 known close massive binaries in Cyg OB2 are shown in Figure 1. The orbital period and mass ratio distributions from Garmany et al. (1980) are also shown and represented by the dotted lines in the upper and lower panels, respectively. The Cyg OB2 and Garmany et al. (1980) period distributions show a concentration between three and five days (with a peak around four to five days for Cyg OB2). A possible explanation for such a concentration might be observational bias owing to the average observing run lasting a few days to a week. Three systems in Cyg OB2 with periods between three and five days originate from other studies (MT421, MT429, and B17) where this may be the case. However, the remaining five systems in this period range originate from our study which has a time coverage of up to 12 years and individual observing run durations of one to three weeks. This allows for a sensitivity to orbital periods up to 25 days without appreciable bias. This abundance of three to five day systems is not limited to the SB2s having the largest mass ratios either. Three of the nine SB1s in Cyg OB2 also have periods that are in the four to five day range. This is illustrated in Figure 2, which shows eccentricity versus orbital period (upper panel) and mass ratio (lower panel) versus orbital period. The SB1s and SB2s are the open and filled circles respectively, and the asterisks represent the orbital periods and eccentricities for the massive star binaries (primaries of type B3 and earlier) in the Ninth Spectroscopic Binary Catalog (Pourbaix et al. 2004). One can easily see the concentration of points around three to five days and a peculiar shortage of systems with periods from 7 to 14 days in the Cyg OB2 data. An initial hypothesis for this shortage of points might be that completeness plays a role. However, with multiple studies with varying observing cadences contributing to the catalog of binaries in Cyg OB2, we maintain that a more likely observation would be to see a continually diminishing number of binary detections with increasing period as completeness declines. Given the existence of this excess in periods around three to five days, we caution that a simple power law may be insufficient to describe the distribution of orbital periods for early-type stars in this cluster.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 3, we show the results from the MC analysis of the 22 close OB binaries in the "complete" sample. The four panels display the 1σ (68.3%), 2σ (95.4%), and 3σ (99.7%) contours relative to the maximum for the average probability of a generation as a function of two parameters for fixed best-fit values of the remaining two parameters. For example, the upper left panel shows the average probability as a function of α and γ for fixed best-fit values of β (0.2) and B.F. (45%). The upper right panel similarly shows the average probability as a function of β and γ for fixed best-fit values of α (−0.1) and B.F. (45%). The cross in each panel marks the location of the best fit for each parameter (i.e., the maximum of the average probability). We find that the period distribution is confined at the 1σ level to −0.3 < β < 0.5, with a best fit of β = 0.2 ± 0.4. This points to a nearly log-flat distribution, consistent with Öpik's law (Öpik 1924), consistent with the previous O-star findings of Garmany et al. (1980), Kouwenhoven et al. (2007), Sana & Evans (2011), and KF07. Figure 4 shows the cumulative period distribution of the Cyg OB2 data (diamonds) and the best-fit cumulative period distribution from the MC code (solid curve). The hump in the distribution at log(P) ∼ 0.65 or P ∼ 4.5 days greatly influences the fit shown, and it is clear that a power law of β = 0 might be a better fit if this hump were not present or less pronounced. Indeed, no single power law provides a good fit to the cumulative period distribution. We address the peculiar nature of the period distribution in Section 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows the survey completeness as a function of orbital period for the best-fitting power-law indices. We use the MC code to estimate the survey completeness by tallying the fraction of synthetic systems deemed detectable in our survey given the time sampling and velocity semi-amplitude limit of K_MIN  =  15 km s⁻¹. The survey is 88% complete at P  =  10 days and ≳ 83% complete at P  =  30 days. Our completeness analysis implies that, given 14 detected binaries with periods P  <  10 days, there are an additional one to two undetected binary systems in this range. Statistically, the MC code shows that all of these would be undetectable by the survey because they have some combination of low inclination, small mass ratio, large eccentricity, and larger radial velocity uncertainties. If, for the sake of argument, we ignore the excess of short-period systems and assume that the β  ∼  0 power law holds, then the ∼15 binaries in the 1–10 day period range would imply an equivalent number between 10 and 100 days. There are eight known systems in this range, leaving an additional seven putative binaries as undetected. The mean completeness over this range is ∼82%, suggesting that ∼12 of the 15 would be detectable by our survey. This implies that the ongoing Cyg OB2 survey should uncover as many as ∼5 additional binaries with periods between 10 and 100 days, given sufficient long-term data. For periods longer than 100 days, where there are presently only two systems with orbital solutions, Schulte 9 (Nazé et al. 2010) and MT070 (Paper VI), detection becomes more difficult. In the case of Schulte 9, with a period of ∼2.355 years and an eccentricity of nearly 0.7, even an observing campaign such as this one experienced difficulty in uncovering this massive "twin" system. All of the foregoing discussion, however, assumes an Öpik period distribution, and we have already shown in Figure 4 that no single power law adequately reproduces the data at the shortest periods, P < 14 days.

Zoom In Zoom Out Reset image size

Observationally, binaries with solar-type primaries and periods longer than a year have a broad distribution of orbital eccentricities with a median value of 0.55 (Duquennoy & Mayor 1991). For periods shorter than this, these systems have a higher concentration near zero eccentricity stemming from circularization of the orbits. This is also shown to be true for O-star binaries with periods as long as 10,000 days (Sana & Evans 2011). Although the method by which binary orbits are circularized and the critical period below which all are circular is unclear from existing data, a standard eccentricity versus period scatter diagram, like Figure 2, illustrates that binaries with short periods do tend toward nearly circular orbits. The figure shows a decline of mean eccentricity for periods of less than 10 days. Giuricin et al. (1984) found that nearly all (92%) early-type binaries with periods shorter than two days had an eccentricity less than 0.05. They also found that the eccentricities of nearly all of the early-type binaries in their survey could be explained by the Zahn (1977) tidal theory, which describes the circularization time for an equal mass binary with an initial eccentricity of e < 0.3 as being related to the main sequence lifetime of that binary. Khaliullin & Khaliullina (2010) provide a modified version of the Zahn (1977) theory and address the evolution of the orbital parameters over the circularization time of the orbit. The circularization timescales are, on average, a third of those predicted by Zahn (1977) and in better agreement than the predictions of Tassoul (1988) with the massive star binary statistics of Giuricin et al. (1984). The predicted circularization time for an early-type (pair of B0V stars), nearly equal-mass binary at a separation comparable to the Cyg OB2 early B-type binaries (e.g., asin i  ∼ 0.16 AU for MT720) is close to the main-sequence lifetime of an early B star. However, a study by Witte & Savonije (2001) found that there are a significant number of early-type binaries with periods up to 10 days with an eccentricity near zero that cannot be explained with the Zahn (1977) theory and are unlikely to have been primordially circular. Krumholz & Thompson (2007) have proposed a method by which orbits become significantly circularized during the star formation process, thereby diminishing the conflict between circularization time and the theoretical life span of massive stars. However, the theory assumes the binary is formed concurrently with the formation of the components, and that is a highly debated topic in massive star formation theory (Zinnecker & Yorke 2007).

In Figure 6, we show the MC code best fit to the Cyg OB2 eccentricity distribution. As in Figure 4, diamonds show the normalized cumulative distribution (left Y-axis) and the solid curve shows the normalized cumulative distribution for γ = −0.6. Figure 3 (upper left and right panels) shows that the index γ is confined, at the 1σ level, to −1.0 < γ < −0.4 with a best fit of γ = −0.6 ± 0.3. The steep negative power law is required to reproduce the abundance of systems with low eccentricity. A K-S test shows that the γ = −0.6 power law is compatible (i.e., the observed sample and the synthetic sample are consistent with being drawn from the same parent distribution) with the data at the ∼90% level. The absence of eccentricities with e > 0.55 is similar to what is observed for short-period binaries in other OB clusters (e.g., NGC 6231; Sana et al. 2008). Highly eccentric systems are apparently rare in this period regime, either because circularization mechanisms are efficient or possibly because they are more easily disrupted during encounters with other cluster stars. Finally, we note that the eccentricity distribution within Cyg OB2 is also compatible with that of OB binaries in the Ninth Spectroscopic Binary Catalog at the ∼50% level.

Zoom In Zoom Out Reset image size

We show an estimate of the survey completeness as a function of eccentricity for two different period ranges in Figure 7. The solid and dot-dashed curves represent the completeness for P  <  1000 days and P  <  30 days, respectively. The plot indicates that the survey is ∼75% complete at e ∼ 0.65 for binaries with P < 1000 days and ∼88% complete at e ∼ 0.75 for binaries with P < 30 days. For P < 30 days, the survey is likely to have detected the vast majority of systems, even those with large eccentricities, if they exist.

Zoom In Zoom Out Reset image size

The nature of the mass ratio distribution is one of the more highly debated and unresolved questions pertaining to massive binaries and their formation. The distribution is difficult to characterize because of the unavoidable bias of observational studies, including this one, against low-q systems. The general consensus is that short-period binaries have a different mass ratio distribution compared to the long-period binaries. Goldberg et al. (1991) and Duquennoy & Mayor (1991) show this with lower-mass binaries, but they also point out that the period range over which this change occurs is still ill-defined. Most studies on the mass ratio distribution are for mid- to late-type stars. A few examples include Duquennoy & Mayor (1991), who find α ⩽ −2 (for q > 0.1 and periods 1–10,000+ days) for a large survey of solar-type spectroscopic binaries, Kouwenhoven et al. (2007), who find α = −0.4 to −0.5 for a large survey of spectroscopic and visual A- and B-type primaries in Sco OB2, Fisher et al. (2005), who find α ∼ 0 for all the SB1s and α > 0 for all the SB2s in the solar neighborhood (for periods 1–500 days), and Hogeveen (1992), who find α ≃ 0 for mass ratios of q < 0.3 and α ⩽ −2 for mass ratios of q ≳ 0.3 in binaries of B-type and later (for periods <1 day to 10,000+ days). Among the massive stars, studies such as Pinsonneault & Stanek (2006) conclude from their study of 21 eclipsing massive binaries (P < 5 days) in the Small Magellanic Cloud that the secondary mass function for massive binaries does not follow a Salpeter (1955) slope. Instead, the authors find a flat distribution of mass ratios and predict that roughly 45% of all massive systems are part of a "twin" system with a mass ratio of q > 0.95. Sana & Evans (2011) have refuted this proposition, but also find a preference for more massive companions among O-stars in clusters.

The distribution of observed mass ratios in Cyg OB2, shown in Figure 1, appears relatively flat, except for a slight preference for mass ratios between 0.7 and 0.9. This may be an artifact of using a standard histogram for small numbers. Not surprisingly, the upper end of the distribution is populated exclusively by SB2s, while the lower end of the distribution is almost entirely populated by SB1s. For five of these single-lined systems, we calculated the expected orbital inclination, 〈i〉, between i_MIN and 90°, using the inclination lower limits estimated from Paper II, Paper III, and Paper VI. We calculate the most probable mass ratios to be 0.27, 0.30, 0.03, 0.22, and 0.11 for MT059, MT145, MT174, MT258, and MT267, respectively. These are close to the lower-limit values for q because 〈i〉 is not far from the upper limit on i. There is a possibility, of course, that for one or more of these systems, the orbital inclination could be less than 〈i〉, and therefore, the system(s) would have a more massive companion. In such a case, our present estimates may incorrectly inflate the number of systems with mass ratios between 0.05 and 0.30, and the true distribution would not be as flat as currently shown in Figure 1.

From the MC code, we find a best fit of α = 0.1 ± 0.5 for the mass ratio distribution. Figure 8 shows the normalized cumulative distribution of mass ratios (diamonds), along with the corresponding simulated cumulative distribution for α = 0.1 (thick curve). The 1σ contours (top left and bottom left panels in Figure 3) show a range of −0.3 < α < 0.7. Such a shallow positive power law indicates that mass ratios in the cluster are quasi-uniformly distributed between q = 0.005 and q = 1.0 with a preference for larger mass ratios. This type of distribution implies that the secondaries are not drawn from a Salpeter (1955) stellar initial mass function, in contrast to lower-mass binary populations which exhibit an abundance of low-mass ratios. Furthermore, the histogram of mass ratios in Figure 8 is also inconsistent with the "twin-heavy" population suggested by Pinsonneault & Stanek (2006), whereby 45% of systems are predicated to have q > 0.95.

Zoom In Zoom Out Reset image size

Figure 9 shows the completeness of the survey as a function of mass ratio. For systems with P < 30 days (dot-dashed curve), the survey is at least 90% complete down to q ∼ 0.25, but there is a sharp drop in the fraction of systems detected for q ≲ 0.25. The average completeness is 10%–20% lower, overall, for longer period systems with P  <  1000 days (solid curve). Our measurement of MT174 (q  =  0.03; K₁  =  9 km s⁻¹) below our nominal semi-amplitude detection threshold of 15 km s⁻¹ provides evidence that such low-amplitude, low-q systems exist and are detectable in our survey. Figure 9 shows that incompleteness is an issue only at the lowest mass ratios. Any attempt to correct for incompleteness at q < 0.25 would add low-q systems to the histogram and serve to flatten the distribution.

Zoom In Zoom Out Reset image size

Table 2 shows a selection of the larger surveys that investigated the O-star binary fraction within individual clusters/associations. Column 1 contains the common cluster identification, Column 2 is the number of O stars included in the survey, Column 3 is the number of early-B stars included in the survey, Column 4 is the lower limit on the O-star binary fraction, Column 5 is the lower limit on the massive star binary fraction (O and early-B primaries), and Column 6 references the original works. The observed O-star binary fractions (designated b.f. rather than B.F. here to distinguish from the intrinsic binary fraction) within the clusters listed range widely from 0% to 67% (Column 4). More inclusively, the massive star b.f. (with primaries B3 and earlier; Column 5) varies between 4% (NGC 330 in the SMC) and 53% (NGC 6231). In Paper I, we find an initial lower limit on the binary fraction of 30% based on the χ² PDF for radial velocities from the first six years of observations of Cyg OB2, which would place it near the middle of the binary fraction range among massive stars. KF07 find a more probable binary fraction of B.F.  >  70% based on their MC analysis of the raw radial velocities, which would place it closer to that of NGC 6231. However, based purely on the binary detections of Papers II–IV and VI, the hard lower limit on the binary fraction among massive stars in Cyg OB2 is 21% (i.e., 24/114).

We show the observed χ² PDF for the radial velocities of Cyg OB2 stars in Figure 10. There are 53 having P(χ², ν)  < 5% and 43 with P(χ², ν)  < 1%. These systems encompass the 24 known binary systems and additional candidates for radial velocity variables. The number of variables has increased from Paper I, despite removing some low-S/N spectra that may have caused extraneous velocities. This is understandable given the additional data acquired on the sample, especially on long-period binaries like MT070 (P = 7.58 years). Figure 11 shows the MC code best fit to the cumulative χ² PDF for radial velocities in Cyg OB2 using the "unbiased" sample of 110 stars. The solid curve represents the cumulative χ² PDF produced by the MC code with best-fit indices of α = 0.1 (q_MIN = 0.005), β = 0.4 (P_MAX = 1000 days), and γ = −0.6 (e_MIN = 0.0001, e_MAX = 0.9). The open diamonds represent the normalized cumulative χ² PDF for the unbiased sample of OB primaries. The best MC code fit results in a binary fraction of B.F. = 45% ± 8% if all binary systems are assumed to have P < 1000 days. Simulated binary populations with binary fractions of 40%–50% in increments of 1% and with fixed α, β, and γ, were individually examined for better agreement with the data. These additional simulations show that a binary fraction of 44% yields a slightly better fit to the data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 12 shows the results of the MC analysis for the 18 OB close binaries comprising the "unbiased" sample. Similar to Figure 3, the four panels provide the 1σ, 2σ, and 3σ probability contours for α, β, γ, and binary fraction. The three panels depicting probabilities for power-law indices (upper left, upper right, lower left) are nearly identical to the results from the MC analysis of the "complete" sample shown in Figure 3. The best-fit β is slightly larger (0.4 instead of 0.2) owing to the exclusion of several shorter-period binaries from the unbiased sample. The best-fit α and γ are the same, just with larger uncertainties. The contours in the lower right panel show that the binary fraction is confined to 37% <B.F. < 49% at the 1σ level and 30% <B.F. < 60% at the 3σ level for P_MAX = 1000 days. This result is consistent with the range of massive binary fractions for Galactic clusters listed in Table 2. In order to compare our current results to the MC analysis of KF07, who explored a larger range of periods (separations of 0.02–10,000 AU corresponding to P  ∼ 1–10⁷ days for two 20 M_☉ stars), we also performed an MC run with P_MAX = 10⁷ days. In both analyses, 10,000 AU is considered a reasonable upper limit for bound systems. In this case, we find that the binary fraction must be ∼100% in order to replicate our observed χ² PDF. This is in rough agreement with the conclusions of KF07, but based on extrapolating the power-law distributions of mass ratios and periods determined from close binaries over many orders of magnitude to characterize loosely bound systems.

Zoom In Zoom Out Reset image size

Is an extrapolation of the power-law index β to periods as long as 10⁷ days warranted when our radial velocity survey is sensitive only to systems with periods less than about 10 years? Caballero-Nieves et al. (2011) conducted an adaptive optics and HST Fine Guidance Sensor survey of 75 Cyg OB2 systems and detect highly probable companions around 31 OB primaries at separations >200 AU (P > 630 years). Of these, nine (9/75 = 12%) have probable periods in the range 10³–10⁴ years (S. Caballero-Nieves 2012, private communication).⁴ If we adopt P_MAX = 10⁴ years (∼365, 250 days), then the MC code provides a best-fit binary fraction of ∼90%, and this extrapolation of the β ≃ 0.2 power law predicts approximately 15 binaries with orbital periods in each 1 dex interval (i.e., 15 with periods between 10 years and 100 years, 15 between 100 years and 1000 years, etc.) given the 110 stars in our unbiased sample. This translates into about 14% of the binary systems falling within each logarithmic period interval. Hence, the predictions of an extrapolated power law appear consistent with the emerging high angular resolution observations at larger separations. We therefore suggest that the binary fraction among massive stars may be as high as 90% ± 10% if periods as long as 10⁴ years are allowed.

In addition to exploring the effect of varying P_MAX, we also examined the effect of varying q_MIN and the semi-amplitude detection threshold, K_MIN, in the MC code. As K_MIN is changed from 15 km s⁻¹ to 20 km s⁻¹ or larger, a greater fraction of synthetic systems are deemed undetectable by the code, and the detected systems have smaller separations and/or larger companions. Therefore, we expect some combination of increase in α (larger mass ratios), decrease in β (shorter periods), and increase in binary fraction. We find insignificant change (i.e., α, β, and γ are the same) when K_MIN = 20 km s⁻¹. A significant change in the indices is not observed until K_MIN is increased to much larger values of ∼50 km s⁻¹. We are also confident our detection threshold is much closer to 15 km s⁻¹, as evidenced by the measurement of MT070 (H. A. Kobulnicky et al. 2012, in preparation). We conclude that the minimum semi-amplitude for binary detection is appropriately set at 15 km s⁻¹.

In the case of increasing q_MIN, the lowest mass ratio simulated in the code, we primarily expect to see a decrease in α and a decrease in the binary fraction. Binaries will have more massive companions and semi-amplitudes will be larger, meaning that a smaller fraction of binary systems will be deemed undetectable in the code. We observed only a small decrease in binary fraction in runs with q_MIN = 0.1 (a few percent and still consistent with our initial results). We also find only a small increase in α and no change in β or γ. However, increases in α become larger when q_MIN > 0.1. For example, α changes from 0.1 to −0.4 for the "complete" sample when q_MIN = 0.2; however, such a large q_MIN is probably unphysical given the existence of MT070 (q ⩾ 0.03) and MT267 (q ⩾ 0.11). These results, for both tests, indicate that our MC results are not sensitive to small changes in the adopted values of the mass ratio lower limit or semi-amplitude threshold.

The presence of an excess of binary systems with periods around three to five days (with a peak at four to five days) and a shortage of systems with periods around 7–14 days is not only observed in Cyg OB2. Sana et al. (2008) observe a similar occurrence in NGC 6231, and the excess and shortage of systems occur at roughly the same ranges in period. Systems in Cyg OB2 that lie within the excess encompass a range of spectral types, mass ratios, and binary types (i.e., eclipsing, SB1 and SB2) and are scattered randomly across the cluster, having no preferred location either in dense or sparse regions. We do not believe that this peculiar distribution can be attributed to observational biases given the extended time coverage of this survey and the different observing cadences of the observing runs that have contributed to the Cyg OB2 binary catalog. While the existence of an excess of systems with periods of around three to five days seems well established in the literature and in our present data, the possibility of a shortage in the 7–14 day range has not been previously noted. A comparison between the Cyg OB2 cumulative period distribution with the data from Garmany et al. (1980) reveals a flattening in this 7–14 day range (though we should be cautious given the use of a standard histogram and small numbers). Their survey (of 28 systems) contains three in this 7–14 day range while ours (22 systems) contains one.

While we find that no single power law can adequately describe the log-period distribution in Cyg OB2, it is unlikely, given any power-law distribution, to find so few systems between 7 and 14 days. We examined the probability that only a single system would be observed in the 7–14 day period range given a random drawing of 22 systems with 1.46 days < P < 25.2 days (for β = −1 to 1). Random samples were drawn 10,000 times and the final percentages were adjusted for a 90% completeness level (the mean level as measured in Figure 5 for the 7–14 day range). The likelihood of detecting a single system for β = 0 is ∼2%. However, if the power law describing the log-period distribution is assumed to be closer to −1, then the probability rises quickly to ∼13% (with ∼5% for β = −0.5, and <1% for β = 0.5 and 1). This tells us that if the power law(s) that describe the log-period distribution have β ⩾ −0.5, then the shortage is a real phenomenon at the 2σ level.

In an attempt to replicate the Cyg OB2 period distribution with a physically motivated prescription, we made a version of the MC code to take 70% of the systems assigned a period between 7 and 14 days and reassign them assuming the same power law describing the rest of the distribution to periods between 3 and 5 days. The 70% reassignment fraction reflects a 90% completeness level and detection of one out of a probable three or four systems in this range. This approach approximates what would happen if binary systems with periods 7 and 14 days migrated to shorter periods. The best fit using this approach is shown in Figure 13 and comes from reassigning systems with a period between 7 and 13.5 days to a period between 4 and 5 days (corresponding to the peak of the excess in the Cyg OB2 distribution) with only a slight increase in β (0.5 and within the uncertainty of our results). α, γ, and the binary fraction are the same as in Section 4. The normalized cumulative distribution for the Cyg OB2 data is represented by the diamonds and the solid curve represents the normalized cumulative period distribution for the best fit using this version of the MC code. Though this fit requires a slightly more positive slope for the initial period distribution, the overall fit to the data is visually better than in Figure 4, and the multiple changes in the slope are also replicated well. A two-sided K-S test yields a compatibility of greater than 99% (compared to ∼60% with the simple power-law approach).

Zoom In Zoom Out Reset image size

A plausible explanation for such period migration could be the mechanism proposed by Krumholz & Thompson (2007) whereby the circularization of the orbit occurs primarily during the formation process instead. In Figure 14, we show the evolution of orbital period/semi-major axis and mass ratio for a binary system with a total mass of 30 M_☉ as the primary loses mass to the secondary. The plotted curves show tracks of constant angular momentum as systems evolve toward shorter periods with mass ratios approaching unity, assuming that mass and angular momentum are conserved,

$$\begin{equation} {a \over a_0} = \left({{q_0} \over {q}} \right)^2 \left({{1+q} \over {1+q_0}}\right)^4, \end{equation} \tag{ 1 }$$

(e.g., Equation (8) of Krumholz & Thompson 2007). Here, a₀ and q₀ are the initial semi-major axis and mass ratio while a and q are the final semi-major axis and mass ratio after mass transfer. The two solid curves are tracks for systems with initial mass ratios q = 0.5 and initial periods of 7 and 14 days. These systems may end up having mass ratios near unity, but their final periods are longer than five days. The dotted curves are tracks for systems that begin with mass ratios q = 0.3 and periods of 7 and 14 days. These systems may end up with mass ratios near unity and periods in the range three to five days. This is precisely where we observe an excess of systems with short periods in Cyg OB2, and we also find a relative abundance of systems with mass ratios q = 0.8–0.9—not exactly "twins," but systems where significant mass transfer may have taken place. This figure indicates that, if the Krumholz & Thompson (2007) mechanism is the operative one producing an abundance of short-period, high-q systems, then some fraction of the population of 3–5 day systems should originate as systems with q ≃ 0.3 with 7–14 day periods. While the evidence for "missing" 7–14 day systems in Cyg OB2 is compelling, we cannot claim, on the basis of the present data, to see the corresponding lack of systems with q ≃ 0.3. Our survey is sufficiently incomplete in this mass ratio range that observational biases preclude any firm conclusion.

Zoom In Zoom Out Reset image size

The Cyg OB2 radial velocity survey originated as an observational effort to inform theoretical models predicting the frequency of energetic events like supernova and γ-ray bursts that may have close binaries as progenitors. A key ingredient of these modeling efforts is the fraction of massive stars that might interact through mass transfer or undergo a common envelope phase, especially during the post-main-sequence evolution of either star in the system. Red supergiants (RSGs) are the descendants of massive stars, so we adopt the RSG radii deduced by Levesque et al. (2005) which range between 200 R_☉ and 1500 R_☉. Their Figure 3 shows stellar evolutionary tracks indicating that stars in the mass range 15–30 M_☉ may reach or exceed 1000 R_☉, while stars with initial masses over 10 M_☉ can generally be expected to reach 500 R_☉. We adopt 500 R_☉ as a maximum evolved radius for stars in our sample, meaning that systems with semi-major axes smaller than this value will interact strongly, probably through the formation of a common envelope. Other systems may interact via mass transfer after Roche lobe overflow. Eggleton (1983) provides an analytic description of the Roche lobe radius as a function of mass ratio. His tabulated values show that the Roche lobe radii for binary systems with mass ratios q = 0.2–1.0 fall in the range 0.38a–0.55a, where a is the separation of the components. Consequently, we adopt 1000 R_☉ (4.7 AU) as a conservative threshold separation, below which systems are deemed to interact during the evolution of the system. In our MC code, we flag all systems with periastron distances less than 1000 R_☉ as capable of interaction. Assuming α = 0.1, β = 0.2, γ = −0.6, and P_MAX = 10⁴ years, this predicts that 45% of the binary systems will have sufficiently close companions to interact during the binary's lifetime (37% if we adopt the slightly larger β = 0.5 from Section 5.1).

The luminous contribution from secondary stars, which are unresolved in photometric surveys, even with the HST or adaptive optics, is not negligible. Consequently, measuring accurate distances to star clusters using spectroscopic parallax, main-sequence fitting, or similar standard-candle techniques is systematically in error if close companions are neglected. As part of our MC code, we tallied the mean V-band light contribution from primary and secondary stars by using a zero-age main-sequence (ZAMS) relation between stellar mass and absolute V magnitude (Martins et al. 2005; Drilling & Landolt 2000). For the best-fitting power-law indices described above, allowing for periods as long as 10⁴ years, we find that the secondary stars contribute 16% of the total stellar light at the V band. This is equivalent to a systematic error in the distance modulus of 0.16 mag and a systematic distance error of about 8%, in the sense that distances derived without correcting for luminous secondaries will be too small. We expect the secondary luminosity fractions in evolved populations to vary stochastically about this 16% ZAMS value as the primary and secondary components throughout the population become giants and supergiants. The intricacies of mass exchange between components will further limit any simple description of the time evolution in the nominal ZAMS value.