ton, Kisseleva-Eggleton & Dearborn 2007;Tokovinin 2008;Griffin 2012;Moe & Di Stefano 2017). Similarly, massive binaries may often reside in triple configurations (Zinnecker & Yorke 2007;Raghavan et al. 2010). While the fraction of massive stars in triples is not well known, the average number of companions per OB primary may be at least three times higher than that of low-mass stars, which have 0.5 companions on average (Preibisch, Weigelt & Zinnecker 2001;Grellmann et al. 2013).

Dynamical stability arguments imply that these triple systems must have a hierarchical configuration: the third body orbits the inner binary on a much wider outer orbit. In this configuration, coherent gravitational perturbations from the outer body influence the long-term evolution of the inner orbit. The orbits can be treated as massive wires that torque each other, where the line-density of each wire is inversely proportional to the orbital velocity. In this orbit-averaged, or secular, approximation, the semimajor axis ratio remains constant. The gravitational potential can be expanded in terms of this ratio a 1 /a 2 , where a 1 (a2) is the semimajor axis of the inner (outer) orbit (Kozai 1962;Lidov 1962). The hierarchical configuration makes this expansion possible by ensuring that a 1 /a 2 is a small parameter. The lowest order, or quadrupole level, of approximation is proportional to (a 1 /a 2 ) 2 . The next level of approximation is called the octupole level (see Naoz 2016, a recent review).

This secular dynamical process is referred to as the Eccentric Kozai-Lidov (EKL) mechanism.

While quasi-secular systems may also exist, we restrict ourselves to hierarchical systems. This focus is motivated by G-type star observations, which indicate that the vast majority of systems are hierarchical (e.g. Tokovinin 2008). Non-secular effects dominate the dynamics of marginally hierarchical systems, where, for example, the EKL time-scale is approximately the outer orbital period. These systems may undergo eccentricity oscillations on an orbital timescale (e.g. Antognini et al. 2014;Antonini, Murray & Mikkola 2014), and these oscillations may have a higher amplitude than the secular case (e.g. Katz & Dong 2012;Bode & Wegg 2014;Grishin, Perets & Fragione 2018). Assuming that the orbital configurations of massive triples parallel those of G-dwarfs, we underestimate the number of binaries that cross the Roche limit (RL) by a few per cent.

The applications of hierarchically configured stellar triples have been explored extensively in the literature. For example, triples have important implications for the evolution of sun-like stars (e.g. Harrington 1969;Mazeh & Shaham 1979;Kiseleva, Eggleton & Mikkola 1998;Fabrycky & Tremaine 2007;Naoz & Fabrycky 2014;Bataille, Libert & Correia 2018;Moe & Kratter 2018). Triple dynamics also represents a key formation mechanism for bluestragglers (e.g. Perets & Fabrycky 2009;Naoz & Fabrycky 2014;Fragione 2019). Furthermore, short period compact object binaries, including black hole, neutron star, and white dwarf binaries, may result from hierarchical triple dynamics (Thompson 2011;Antonini & Perets 2012;Katz & Dong 2012;Naoz et al. 2016;Hoang et al. 2018;Toonen et al. 2018;Fragione & Loeb 2019). The evolution of binaries with a distant massive companion proves essential to understanding mergers within galactic nuclei (Antonini & Perets 2012;Prodan, Antonini & Perets 2015;Stephan et al. 2016Stephan et al. , 2019;;Fragione et al. 2018;Hoang et al. 2018Hoang et al. , 2019;;Fragione 2019). This work explores the implications of triple dynamics for massive star systems, which connect to transient phenomena such as mergers.

We examine the impact of EKL-driven dynamical evolution on the properties of massive binaries embedded in triples. The birth orbital configurations of massive binaries are not well constrained. We therefore explore different possibilities for the period and eccentricity birth distributions with the objective of discerning the birth properties from the final results. We test a broad range of initial distributions, tidal recipes, and efficiencies to identify the manner in which observed period and eccentricity distributions map back on to the initial distributions. While point-mass dynamics coupled with stellar evolution has many interesting consequences (e.g. Pejcha et al. 2013;Shappee & Thompson 2013;Michaely & Perets 2014;Stephan et al. 2016Stephan et al. , 2019;;Fragione & Loeb 2019), we focus on 10-Myr-old systems, before significant stellar evolution has taken place. Our simulations include tidal dissipation and general relativistic precession.

We run 38 sets of Monte Carlo numerical simulations with a variety of initial conditions and tidal recipes. Each simulation has 1000 systems, bringing the total number of realizations to 38 000. Section 2 reviews our methods and initial conditions. In Section 3, we present our results and analyse the statistical distributions of the inner orbital periods and eccentricities. Section 3.4 investigates traces of the initial conditions in the final distributions, traces which persist even after 10 Myr of EKL-driven evolution. Finally, we consider our simulated results in the context of observations and discuss the implications in Section 3.5. Our simulations suggest that the EKL mechanism can re-distribute the orbital periods to match observed distributions, while the eccentricity distribution retains its original shape.

We solve the secular equations of motion for the hierarchical triple to the octupole-level of approximation, as described in Naoz et al. (2013a). Stars with masses m 1 and m 2 compose the inner binary, while the tertiary body with mass m 3 and the inner binary form an outer binary. The inner (outer) orbit has the following parameters: a 1 (a 2 ), e 1 (e 2 ), ω 1 (ω 2 ), and i 1 (i 2 ) for the semimajor axis, eccentricity, argument of periapsis, and inclination with respect to the total angular momentum, respectively. We define i tot as i 1 + i 2 . We include general relativistic (GR) precession following Naoz et al. (2013b), who demonstrate that GR precession can suppress or facilitate eccentricity excitations in different parts of phase space.

We limit our simulations to 10 million years to allow for comparisons with observed massive stars in young stellar clusters (e.g. Sana et al. 2012). We treat the mass and radius of each star as constant. Over 10 million years, a 20 M star, the largest mass allowed in these simulations, will lose about three per cent of its mass, and the semimajor axis of its orbit will expand by 3 per cent. Thus, while the interplay between the point-mass dynamics and stellar evolution has interesting consequences (e.g. Stephan et al. 2016), it has negligible effect on our calculations.

We include tidal dissipation, which acts to circularize and tighten the inner binary. Following Naoz & Fabrycky (2014), we adopt the tidal evolution equations of Eggleton & Kiseleva-Eggleton (2001). These equations implement the equilibrium (E) tide model of Eggleton, Kiseleva & Hut (1998). Eggleton & Kiseleva-Eggleton (2001) relate the viscous time-scale t v of a star to the tidal dissipation time-scale t F using the Love parameter k L, j , where j = 1, 2 for masses m 1 and m 2 , respectively. The Love parameter quantifies the quadrupolar deformability of a star. For the primary star with mass m 1 ,

(1)

We use k L, j = 0.014 for an n = 3 polytrope (Eggleton & Kiseleva-Eggleton 2001).

The spin of the stars influences the tidal evolution of the binary. As we discuss below, our tides are not very efficient, and thus the spin period may only be important to the lower portion of the mass distribution (∼1 M ). To accommodate this contribution, we choose the primary and secondary stars to have spin periods of 4 d. This 4 d estimate is consistent with magnetic braking evolution up to the Sun's current spin rate (Dobbs-Dixon, Lin & Mardling 2004). We note that the expected period of massive stars ranges between ∼2 d (e.g. Bouvier 2013) to ∼0.5 d (e.g. Hurley, Pols & Tout 2013;Stephan et al. 2019). We verify that a choice of the smaller 1 d spin gives results (not shown to avoid clutter) consistent with those depicted in this paper, which use a 4 d spin period.

We also implement a prescription for radiatively damped dynamical (RDD) tides (Zahn 1977) for stars of mass greater than 1.5 M . Following Hurley, Tout & Pols (2002) and Zahn (1977), we express

.

(2)

We note that in a comparison with numerical calculations, Chernov (2017) suggests that this tidal prescription underpredicts the tidal efficiency in short-period binaries. Claret & Cunha (1997) find that Zahn's formalism fails to account for some observed circularized systems. The treatment of stellar tides has invoked much debate (Langer 2009). Noting this lack of accord, Yoon, Woosley & Langer (2010) include a parameter for the rate of tidal synchronization to test for its effects. Similarly, we opt to use the equilibrium tide model for all stars in select simulations, allowing us to treat the viscous time-scale as a free parameter and tune the tidal efficiency. We label our simulations according to tidal prescription and efficiency as detailed in Table 1. Both Hurley et al. (2002) and Zahn (1977) note the inefficient nature of RDD tides compared to E tides. We use E tides with t v = 0.005, an unrealistically small value, in a few simulations for comparative purposes.

To explore the influence of initial parameter distributions on the results, we use several different probability distributions to draw initial conditions for the Monte Carlo simulations. We label simulations according to their birth distributions.

The simulations fix the mass of the primary m 1 at 10 M . We then draw m 2 and m 3 from uniform distributions of the mass ratios q 1 = m 2 /m 1 and q 2 = m 3 /(m 1 + m 2 ), where 0.1 < q < 1. We discuss simulations that use a Kroupa mass function in Appendix A. A power law of the form R ∝ M α derives the stellar radii from the masses. For most simulations, we use the ZAMS mass-radius relation R = 1.01M 0.57 (Demircan & Kahraman 1991). 1 We draw the inclinations cos i 1 (cos i 2 ) of the inner (outer) orbit from a uniform distribution. In all simulations, we include a spin-orbit misalignment for the primary (secondary) star ψ 1 (ψ 2 ) from a uniform distribution. These choices stem from an agnostic view of the formation mechanism of these massive binaries. Like Naoz & Fabrycky (2014), we adopt an isotropic distribution because it allows for chaotic processes during formation.

Observations of young (∼ few Myrs old) massive clusters suggest power-law distributions for period and eccentricity (with different indices). The most conservative choice of initial conditions assumes that the final distributions remain unchanged from the birth distributions. The most agnostic choice assumes uniform initial distributions for period and eccentricity. Overall, we test six birth distribution combinations of eccentricity and period. Below, the abbreviation for each initial condition is given in parentheses. We use these abbreviations in the simulation labels.

We use three possible initial distributions for eccentricity. For populations of massive stars, observations indicate that the probability density distribution is either f(e) ∝ e -0.5 ('e05') (e.g. Sana & Evans 2011;Sana et al. 2012) or uniform ('e 1 ') (e.g. Almeida et al. 2017). We supplement these simulations with ones that assume a thermal eccentricity distribution ('e 2 ').

For the period distribution, we adopt a probability density of the form f(log P/days) ∝ (log P/days) -0.5 ('PS') (e.g. Sana et al. 2012). Since this observed probability density does not necessarily represent the birth distribution, we also use a uniform distribution in logarithmic space ('PU') for certain simulations. Additionally, Sana et al. (2017) present observations of 11 young (< 1 Myr) massive binaries in the open cluster M17 with a low radial-velocity dispersion. They offer multiple explanations for the peculiar velocity dispersion. In one of these explanations, the massive binaries form at larger periods, with some minimum birth separation, and later tighten and circularize. Specifically, Sana et al. (2017) suggest that a lower period cut-off of 9 months is consistent with the velocity dispersion in M17 within 1σ . To explore this intriguing possibility, we include a lower period cut-off of 9 months in multiple ('gap') simulations.

To ensure that the secular approximation is valid, we apply three criteria to the initial conditions (Naoz & Fabrycky 2014). The Mardling & Aarseth (2001) criterion ensures the long-term stability of the system, where i tot denotes the mutual inclination between the inner and outer orbits in radians:

(1e 2 ) 6/5 (1 -0.3i tot ) .

(3)

We note that this stability criterion does not consider the system's lifetime (see e.g. Mylläri et al. 2018), and therefore underpredicts the number of systems that may undergo EKL oscillations during the 10-Myr integration time.

The second criterion compares the amplitudes of the octupole and quadrupole terms to verify that the perturber is weak, i.e. that higher order terms have a negligible effect. This criterion stipulates that , the ratio of the octupole and quadrupole amplitudes, remains small:

The last criterion requires that the inner orbit falls outside the RL to avoid an early merger, before secular effects occur. We use the RL as defined by Eggleton (1983):

where j = 1, 2 and R i denotes stellar radius. All of our systems begin with a 1 (1

The stability criteria affect the simulated birth distributions. In Fig. 1, this effect is most apparent in the thermal eccentricity distribution with no gap. The lack of a gap allows small values Top row: Cumulative eccentricity distribution. Second row: Eccentricity probability density. Third row: Cumulative period distribution. Bottom row: Period probability density. We consider the initial distributions after applying stability criteria (see the text). Initial conditions with a gap are represented by dashed lines. In the top two panels, thermal, uniform, and f(e) ∝ e -0.5 (Sana & Evans 2011;Sana et al. 2012) are shown in green, blue, and red, respectively. Note that those simulations with no gap have fewer high eccentricity systems, a result of having smaller values of a 1i while still needing to satisfy the RL criterion. In the bottom panels, we show a uniform distribution in logarithmic period (grey lines) and f(log P/days) ∝ (log P/days) -0.5 (crimson lines) (Sana et al. 2012(Sana et al. , 2013)). Since the stability criteria tend to curtail large values of a 1 , the PU and PS gap distributions appear similar. We observe a similar albeit less pronounced effect in the distributions without a gap (solid lines).

of a 1 . Given the greater abundance of small semimajor axes, the RL criterion limits the number of high eccentricity systems. Similarly, the requirement that a 1 be small compared to a 2 for a hierarchical triple, a condition reinforced by equations ( 3) and ( 4), curtails the number of inner binaries at large periods (Fig. 1). Due to this effect, the PU and PS gap period distributions strongly resemble each other, and the no-gap PU and PS distributions converge at large periods. Figure 3. RL-crossing time-scales: Here, we show the smoothed histograms of the RL-crossing time, called t RL , using a Guassian kernel density estimator for SET-PU-e1 (solid bright red curve) and SET-PS-e1-gap (dashed maroon curve). We term these systems RL binaries.

We evolve each triple system for 10 Myr. We also include conditions which, if met, result in an early termination of the integration. We consider two stopping conditions:

(i) We terminate the simulation once the system tightens and circularizes because tide-dominated systems become numerically expensive. We consider a system to be tidally tightened and circularized when a 1 < 2.1R i,j /L R,ij and e 1 < 0.001. Fig. 2 shows the typical circularization time-scales for gap and no-gap systems.

(ii) If an inner member have crosses the other's RL, or when a 1 (1e 1 )L R, ij < R i (see equation 5), we terminate the integration. We denote the time upon which the system crosses the RL as t RL . In Fig. 3, we show typical t RL . As expected, the gap systems have systematically longer RL crossing times because they are associated The left, middle, and right columns represent an RL-crossing system, a system that undergoes weak oscillations, and a tidally circularized system, respectively. The first row shows the inner binary orbital inclination (blue) and obliquities, or spin-orbit angles, ψ 1,2 (orange, green) of the stars; the second row, eccentricities of the inner and outer binaries; and the last row, the inner semimajor axis, periapsis, and the RL. Note that in the tidally tightened and circularized system, r p and a 1 converge while e 1 approaches zero. Both simulations take a 10 M primary. with longer quadrupole time-scales. The quadrupole time-scale T quad depends on P 2 2 /P 1 , where P 1 (P 2 ) correspond to the inner (outer) orbital period. Typical gap systems have log T quad > 3.5 with a median log T quad ∼ 6.7. While the no-gap distribution has a similar median (log T quad ∼ 6.8), the distribution is wider and has more systems with low quadrupole time-scales. The lower limit of the distribution shifts to a shorter time-scale, log T quad > 1.5. Unlike Naoz & Fabrycky (2014), we do not count these systems as merged products because we expect their merger times to be a few million years (Stephan et al. 2016), the same order as the stellar evolution time-scale. The outcome of this interaction remains highly uncertain without comprehensive eccentric binary interaction physics. During this process, the stars may or may not be observed as two distinct stars. Thus, to differentiate these systems, we term them RL-crossed binaries, or RL binaries for brevity. During a mass transfer process complicated by forced eccentricity oscillations from the tertiary, an observer may detect m 1 and m 2 as a binary system. We therefore define a possible observed orbit for the RL binaries with a F ∼ R RL by assuming angular momentum conservation, a plausible assumption during the final plunge of the merger. We also artificially set the eccentricities of RL binaries to 0.01, which may inflate the number of circularized systems. We caution that the true properties of these systems are highly uncertain, and the actual observed periods may vary largely. In Fig. 7, we denote this uncertainty with arrows.

Each system has three possible general outcomes, shown in Fig. 4. The final outcome depends mainly on the eccentricity excitations and the efficiency of the tides, which act to tighten and circularize the inner orbit. The 'strength' of the EKL mechanism, parametrized by (equation 4), combined with the mutual inclination i tot of the system determine the nature of the eccentricity excitations (e.g. Naoz 2016).

(i) RL crossing: When extreme eccentricity excitations occur on a shorter time-scale than the tidal forces can circularize the system, the pericenter approach can become smaller than the RL (equation 5). This type of behaviour is depicted in the left-hand side of Fig. 4. This outcome is common in the weak-tide regime examined in this paper.

(ii) Tidal tightening and circularization: During periods of higher eccentricity induced by the third star, tidal forces can act to shrink and circularize the inner orbit. The tidal forces decrease the inner semimajor axis until the inner orbit decouples from the third star. The decoupling of the outer and inner orbits results in a conservation of their individual angular momenta. As a result, the typical final semimajor axis of the inner orbit is ∼2a 1,0 (1e 1,0 ), where the subscript '0' denotes values at the start of the circularization process. The right column of Fig. 4 shows an example of this evolution. Tidal dissipation also tends to align the spin axes of the stars with the inner orbit's angular momentum vector. The spin-orbit angle ψ 1 (ψ 2 ) is defined as the angle between the inner orbit's angular momentum and the spin axis of m 1 (m 2 ). Fig. 5 shows that this angle goes to zero in most systems. However, in some cases, the system can become locked into a resonance, called a Cassini resonance:

where ψ j denotes the spin-orbit angle and j , the spin period of the two stars (j = 1, 2) in the inner orbit (e.g. Fabrycky, Johnson & Goodman 2007;Naoz & Fabrycky 2014;Stephan et al. 2016).

(iii) Weak oscillations: In many cases, the third star induces only weak eccentricity and inclination oscillations that do not result in a dramatic change of the orbital parameters. The middle column of Fig. 4 shows an example system. While the locations of these individual systems in the eccentricity distribution shift slightly, they do not introduce a noticeable net change.

The final eccentricity distribution tends to retain the same shape and curvature as the birth distribution irrespective of other conditions. We quantify their similarity and provide predictions in Section 3.4.1. This result holds for gap and no-gap simulations. To visualize this trend, we group the results by initial eccentricity distribution in Fig. 6. The initial period distribution has no discernible effect on the shape of the final eccentricity distribution.

As expected, tidal efficiency results in more circularized systems, increasing the vertical axis intercept of the cumulative distribution as shown in Fig. 6 for no-gap simulations. Less efficient tides result in an abundance of RL binaries, which can also increase the vertical axis intercept of the cumulative distribution. Recall that we artificially set eccentricity of the RL systems to 0.01, which may inflate the number of circularized systems. Irrespective of efficiency, tides have no effect on the overall shape or functional form (e.g. thermal) of the distribution.

The final long-period distribution tends to retain the same shape and curvature as the birth distribution irrespective of other conditions. We quantify the similarity and provide predictions in Section 3.4.2.

As mentioned above, we treat the RL systems as tight binaries. Fig. 7 illustrates the strong dependence of the fraction of circularized systems on tidal efficiency. The dependence of RL system fraction on tidal efficiency is less pronounced, but generally less efficient tides result in more RL-crossed systems. Fig. 7 exemplifies this trade-off between circularized systems and RL-crossed systems, where the red curves count RL binaries. An RL crossing occurs when the EKL mechanism drives extreme eccentricity excitations before tides can shrink and circularize the system. A lower tidal efficiency therefore gives systems more chances to cross the RL. With the exception of the blue curves in Fig. 7, we include RL binaries as short-period systems unless otherwise noted.

In Fig. 7, the final period distribution for gap simulations appears bimodal and exhibits a dearth of systems at intermediate periods. The peak at short-periods consists of RL and tidally tightened systems. In tidally inefficient (IET) simulations, this peak comprises RL binaries, while the vast majority of systems remain at long periods. Even with the most efficient tides (SET), 10 Myr of EKL evolution fails to fill the gap.

Since tidal efficiency is highly sensitive to stellar radius, larger radii yield more circularized systems. Additionally, during tidal capture, angular momentum conservation (e.g. Ford & Rasio 2006) yields a relation between the final binary separation and the Roche lobe, which depends linearly on the stellar radius. Therefore, as stars leave the main sequence and inflate, we expect circularized systems to have longer periods. To illustrate this behaviour, we consider three approaches to the mass-radius relation in Fig. 8. Specifically, we consider the ZAMS mass-radius relation R = 1.01M 0.57 (blue line), as well as R = 1.33M 0.55 and the TAMS mass-radius relation, R = 1.61M 0.81 (e.g. Demircan & Kahraman 1991) (green and red lines, respectively). To highlight the differences, we include a gap and make the tides unrealistically efficient (UET, unrealistic equilibrium tides). As shown in the figure, the larger radii correspond to longer periods for the tidally tightened binaries. Additionally, larger radii result in a more filled gap at intermediate periods. However, even with unrealistically efficient tides, the TAMS distribution fails to match observations.

To quantify the degree of similarity between the simulated birth and final eccentricity distributions, we fit a function of the form cdf e = κ e e 1 αe + β e , (7)

to the initial and final cumulative distributions, where α e is the index of the power law and κ e and β e are constants. The parameters have the following limits: 0 < κ e < 3, -2 < α e < 2, and 0 < β e < 0.6. We argue that the vertical axis intercept, β e , must be greater than or equal to zero because we are fitting a cumulative distribution. We fit both the initial and final distributions to show consistency between the two. As depicted in Fig. 9, α e changes by 0.3 at most, and only for the most efficient tides. We caution that while this trend holds for hierarchical triple dynamics, we do not account for other dynamical processes that may alter the eccentricity distribution. We vary the range of eccentricity values over which we fit. As more eccentric systems circularize or cross the RL, they increase the number of low-eccentricity systems. When fitting over the full range of eccentricity values, we impose the condition that e 1 > 0.01 to absorb all circularized and RL-crossed binary systems into the vertical-axis intercept. We also fit over the range 0.1 < e 1 < 0.7. A fit of the birth distribution over this range better reflects the function used to generate values because stability criteria target the extremes.

We use a function of the same form, cdf P = κ P log P α P 1 + β P , with the same limits on α P , κ P , and β P . We fit the period distribution over the range 1 < log P 1 < 3. We select the lower limit to avoid the newly formed peak of short-period systems. This lower limit must be adjusted to accommodate the width of the short-period peak, which will depend on the stellar radii (Fig. 8) and therefore the age of the binary system when it circularizes. The reason for the upper limit is twofold: the stability criteria curb the number of longperiod systems, causing our initial distributions to converge at large period, and observational campaigns may not be sensitive to longer periods (e.g. Sana et al. 2013). The birth distribution signature proves much more difficult to discern in the period distributions than in the eccentricity distributions. We find that the orbital period distribution is most consistently preserved over this range, 1 < log P 1 < 3. As depicted in Fig. 10,the power α changes by at most ∼0.3 for PU simulations and stays ≤0.5 for PS simulations.

Several groups have studied the eccentricity distributions of massive binaries. Specifically, Almeida et al. ( 2017) have examined OB-type spectroscopic systems in the Tarantula region to find that 40 per cent of systems have small eccentricities (<0.1). Furthermore, their eccentricity distribution appears uniform. However, the cumulative distribution flattens slightly for high (>0.6) eccentricities, indicating fewer systems there. In a study of 48 massive systems from the Cygnus OB2 association, Kobulnicky et al. (2014) also show a flattening in the cumulative distribution at high eccentricity and attribute it to an observational bias. While this observed flattening may indeed reflect observational biases, we do find that our stability conditions curb the number of high eccentricity systems. As a result, our initial conditions exhibit flattening at large eccentricities (see Fig. 1), which persists in the final distribution.

Similar to Almeida et al. ( 2017), Kobulnicky et al. (2014) note an abundance of low (<0.1) eccentricity systems and conclude that the distribution is uniform for e 0.6. In their review, Duchêne & Kraus (2013) also suggest that massive binaries follow a uniform eccentricity distribution, while Sana & Evans (2011) and Sana et al. (2012) give a probability density of the form f(e) ∝ e -0.5 . Moe & Di Stefano (2017) find a thermal eccentricity distribution, f(e) ∝ e, for OB-type binaries using the catalogue from Malkov et al. (2012). However, their finding pertains to wider binaries with periods between 10 and 100 d, and many of our systems have shorter periods.

We find that the power-law index is a good indicator of the initial condition. Additionally, Geller et al. (2019) determine that star cluster dynamics has little to no effect on the shape of the eccentricity distribution for binaries with modest orbital periods. These combined results imply that the final distribution resembles the birth distribution. The birth distribution of population with uniform eccentricities is therefore also uniform.

Duchêne & Kraus (2013) estimate that 30 per cent of massive binaries have periods of less than 10 d. Sana & Evans (2011) find that 50-60 per cent of systems fall at log P < 1. Similarly, the data of Almeida et al. ( 2017) indicate an abundance (∼ 40 per cent) of short-period (<1 week) systems. Kobulnicky et al. (2014) also affirm the abundance of short-period binaries.

Counting RL binaries in the no-gap simulations, systems with log P 1 < 1 represent 45 per cent, 36 per cent, and 29 per cent of the SET, MET, and IET inner binaries, respectively, in Fig. 7. The superefficient tides (SET, first row) seem too efficient. However, counting only systems with log P 1 < 3.4 to reflect observational limits, the SET and MET simulations match the estimates of Sana & Evans (2011) and Kobulnicky et al. (2014) with 51 per cent and 42 per cent short-period (<10 d) systems, respectively. The IET simulation has 34 per cent. Additionally, observations may not be sensitive to our so-called RL binaries. In that case, the IET tides become too inefficient to account for observations. Examining only tidally tightened systems, the fractions of short-period binaries fall to 46 per cent, 34 per cent, and 24 per cent of total systems with log P 1 < 3.4 for SET, MET, and IET simulations, respectively.

A power law of the form f(log P/days) ∝ (log P/days) α is often fitted to the orbital period distribution (e.g. Kobulnicky et al. 2014). The index of this power law varies in the literature. Sana et al. (2012) suggest α = -0.55, while Sana et al. (2013) find α = -0.45. Almeida et al. (2017) conclude that an index between -0.2 and -0.5 reproduces the data well depending on the range over which they fit. Kobulnicky et al. (2014) suggest that tight binary peak has some structure: the distribution is not uniform at short periods.

In Fig. 11, we show the cumulative distributions for all of our no-gap simulations over the range 0.15 < log P 1 < 3.4. Bounded by the curves cdf P ∝ (log P 1 ) 0.8 (Kobulnicky et al. 2014;Almeida et al. 2017) and cdf P ∝ (log P 1 ) 0.5 (Sana et al. 2012(Sana et al. , 2013)), the grey region depicts the range of power laws fitted to observations. While several of our simulated results resemble the power laws fitted to the data, generally, SET and PS simulations seem to overpredict the number of short-period systems, while IET simulations underpredict the number of short-period systems.

We perform a Kolmogorov-Smirnov (KS) two-sample test to compare our final period distributions to observations. We perform this test both with and without including RL binaries in our simulated sample. Comparing to the data of Sana et al. (2012), we cannot reject the null hypothesis that the observations and simulated results share a parent probability density distribution for the following simulations: MET-PS-e05, MET-PS-e1, SET-PU-e2, IET-PS-e05 both with and without RL binaries. We also perform this test with the data from Kobulnicky et al. (2014) and find that we cannot reject the null hypothesis for MET-PS-e05 and IET-PS-e05 without RL binaries; MET-PU-e2 and MET-PU-e05 with RL binaries; and SET-PU-e2 and SET-PU-e1 both with and without RL binaries. We do not assign much significance to these results. We have a number of systems, the RL binaries, with poorly understood periods in MET and IET no-gap simulations, which hinder a statistical comparison between our results and the data. Additionally, the data have far fewer systems and are subject to observational biases: at large periods, the observed distribution likely diverges from the intrinsic distribution (Sana et al. 2012). The latter explains why the KS test favours PS and SET simulations, while the fitted power laws more often coincide with MET-PU simulations (Fig. 11).

Unlike the no-gap simulations, the gap simulations seem to be in tension with observations. We plot the cumulative distributions for the gap simulations with more efficient tides in Fig. 12. The curves all exhibit a similar behaviour indicative of a bimodal distribution.

Figure 10. Traces of the initial orbital period distribution. We show the simulated initial (left) and final (right) cumulative no-gap period distributions. We overplot a fit with of the form cdf P = κ P (log P 1 ) α P + β P . We fit the distribution over the range 1 < log P 1 < 3. We find that a signature of the birth distribution is best preserved over this range.

We again perform a KS test to compare our simulated results with the observations. The bottom panel of Fig. 12 plots simulation results with the data from Sana et al. (2012) and Kobulnicky et al. (2014) that we test. We find that we can reject the null hypothesis that the simulated distributions and data are drawn from the same parent probability density distribution. Two major differences separate our simulated results from the observations: the persistence of a substantial population of long-period systems and the lack of intermediate-period binaries, as illustrated by Fig. 7. While the true separations of RL binaries -and whether they can fall at intermediate periods -is highly uncertain, a substantial long-period population will none the less persist.

A distant stellar companion can drive the long-term evolution of a massive stellar binary. Distant companions, which may be quite common, can therefore alter the observed orbital parameter distributions for massive binaries. We characterize these effects for a large variety of birth distributions and tidal efficiencies. We find the following:

(i) Spin-orbit angle distribution: Fig. 5 plots the spin-orbit angles and periods of 3000 realizations of massive triples colour coded by eccentricity. For long-period systems, no trend in spin-orbit angle and eccentricity exists. However, as systems circularize at periods of about 10 d, the orbit and stars' angular momenta align such that the spin-orbit angle goes to zero, with a few exceptions: occasionally the spin-orbit angle becomes locked in a (non-zero) resonance, called a Cassini resonance.

(ii) Eccentricity distribution: The final eccentricity distribution is an excellent indication of the birth distribution. The cumulative distribution retains the curvature of the birth distribution. A fit over Power laws from observations compared with no-gap results. We plot the cumulative distributions for our no-gap simulations with and without RL binaries. The grey shaded region encompasses the range of power laws indicated by observations. The curves cdf P ∝ (log P 1 ) 0.8 (Kobulnicky et al. 2014;Almeida et al. 2017) and cdf P ∝ (log P 1 ) 0.5 (Sana et al. 2012(Sana et al. , 2013) ) bound the grey region. To avoid visual clutter, we deviate from our preferred presentation for the period distributions, a probability density, and instead focus on cumulative distributions here. . We plot the cumulative distributions for select gap simulations with and without RL binaries. These simulations have more efficient (SET or MET) tides and therefore the most short-period systems. The gray shaded region represents observationally constrained power laws. The curves cdf P ∝ (log P 1 ) 0.8 (Kobulnicky et al. 2014;Almeida et al. 2017) and cdf P ∝ (log P 1 ) 0.5 (Sana et al. 2012(Sana et al. , 2013) ) bound this region. Lower Panel: Probability density with data. We plot the probability densities for select gap simulations with and without RL binaries. The gray and black dashed histograms represent data from Kobulnicky et al. (2014) and Sana et al. (2012), respectively. the range 0.1 < e 1 < 0.7 of a function cdf e ∝ e αe 1 quantifies this trend (see Fig. 9). Generally, the change in α e is within 0.3, with a tendency to flatten -or render more uniform -the eccentricity distribution.

(iii) Period distribution: A signature of the birth period distribution persists at 1 < log P 1 < 3. We fit the period cumulative distribution with a power law over this range. The index changes by ∼0.3 for the uniform initial condition and remains 0.5 for cdf P ∝ (log P 1 ) 0.5 (PS) initial condition simulations (Fig. 10).

(iv) Short-period binaries: Observations indicate an abundance of short-period binaries (e.g. Duchêne & Kraus 2013). In our simulations, the tidal efficiency determines both the resulting number of short-period binaries and the dominant type of shortperiod system. Less efficient tides give systems more chances to cross the RL, while efficient tides yield more circularized, tight binaries. The former results in an artificial peak in our final period distribution because the final properties of such systems remain uncertain. Due to angular momentum conservation, the final periods of short-period systems depend on the stellar radii. We treat the stellar radii as constant. However, realistically, the radii will expand as the stars age. Systems which tidally tighten or cross the RL at later times will therefore fall at longer periods (Fig. 8).

(v) Initial period gap: Sana et al. (2017) suggest that massive binaries may form with large separations and tighten over time to match the parameter distributions of older populations. The EKL mechanism in concert with tidal dissipation represents a channel for producing hardened binaries. However, the EKL mechanism fails to build up a sufficient population of short-period binaries if we begin with a lower period cut-off of 9 months (Fig. 12). We perform a brief calculation to assess whether type-II migration may fill a 9 month gap in the inner period distribution. Following Armitage (2007), the type-II migration time-scale τ can be written as τ ∼ 2/(3α)(h/a 1 ) -2 -1 , where α is related to the viscosity, h is the scale height, and h/a 1 represents the disc aspect ratio. denotes the angular velocity. A 9-month period corresponds to a roughly 2 au semimajor axis. Taking the primary and secondary to have masses 10 and 5 M , respectively, is approximately 0.7 yr -1 . We assume that h/a 1 ≈ 0.07 (e.g. Ruge et al. 2013) and α ≈ 0.01 (Armitage 2007). We find that the time-scale is about 2 × 10 4 yr. Type-II migration therefore represents a mechanism to bridge an initial period gap, as suggested by Sana et al. (2017), with observations of older populations. Acting over a short timescale to tighten binaries, type-II migration may fill the gap in the period distribution. If a binary system undergoing migration has a third companion, the gas will mostly suppress the gravitational perturbation of the tertiary star, although the system may develop an inclination between the disc and the companion (e.g. Martin et al. 2014).

(vi) Comparison with observations: no initial period gap and moderate tides. We compare our cumulative distributions for the final inner orbital period with the observationally constrained power laws in Fig. 11. Many of our results fall in the region bounded by the power laws indicated by the literature (e.g. Almeida et al. 2017). Howerver, generally, simulations with moderately efficient tides (MET) and a uniform birth period distribution match observed distributions well. Several observations indicate a uniform eccentricity distribution except at high eccentricity (e.g. Kobulnicky et al. 2014). Only simulations that begin with a uniform eccentricity distribution produce a uniform distribution as the end result.

Unlike previous studies of EKL evolution in triple stellar systems (e.g. Naoz & Fabrycky 2014;Bataille et al. 2018;Moe & Kratter 2018), we find that the final period and eccentricity distributions carry a clear signature of the initial distributions. This behaviour is a consequence of the short time-scale of evolution (∼10 Myr) and the METs of stars with radiative envelopes. We expect that as stellar evolution increases the stellar radii and causes mass-loss, the orbital configurations will significantly alter.

SN thanks Howard and Astrid Preston for their generous support. SCR acknowledges support from the Eugene Cota-Robles Fellowship. We overplot a fit with of the form cdf P = κ P (log P 1 ) α P + β P . We fit the distribution over the range 0.15 < log P 1 < 3.4 to examine overall trends in the period distribution. We overplot a fit with of the form cdf P = κ P (log P 1 ) α P + β P . We fit the distribution over the range 0.15 < log P 1 < 3.4 to examine overall trends in the period distribution.

Example two-part fit of period distribution. We fit the most promising (PU, MET) simulated results with in two parts, 0.15 < log P 1 < 1 and 1 < log P 1 < 3. We use the e1 simulation and note that the eccentricity does not affect this result.

We apply the method described in Section 3.4.2 to all of our no-gap simulations to find traces of the initial period distribution in the final result. We show the other five simulations in Fig. C1. As noted in Section 3.4.2, the index α of the power law fitted over the range 1 < log P 1 < 3 changes by at most ∼0.3 for PU simulations and stays 0.5 for PS simulations. We overplot a fit with of the form cdf P = κ P (log P 1 ) α P + β P . We fit the distribution over the range 1 < log P 1 < 3 to find a signature of the birth distribution.

This paper has been typeset from a T E X/L A T E X file prepared by the author.