In Search of the Thermal Eccentricity Distribution

TL;DR

This paper challenges the assumption that star cluster dynamics quickly thermalize binary star eccentricities, showing instead that initial conditions are likely preserved and significantly impact binary evolution predictions.

Contribution

The authors develop a rapid semi-analytic model and compare it with detailed simulations to demonstrate that dynamical encounters do not thermalize binary eccentricities within cluster lifetimes.

Findings

Dynamical encounters rarely convert uniform to thermal eccentricity distributions.

Observed thermal distributions likely reflect formation conditions, not dynamical evolution.

Models assuming initial thermal distributions may overestimate binary merger rates.

Abstract

About a century ago, Jeans (1919) discovered that if binary stars reach a state approximating energy equipartition, for example through many dynamical encounters that exchange energy, their eccentricity distribution can be described by : dN/de = 2e. This is referred to as the thermal eccentricity distribution, and has been widely used for initial conditions in theoretical investigations of binary stars. However, observations suggest that the eccentricity distributions of most observed binaries, and particularly those with masses < 5 Msun, are flatter than thermal and follow more closely to a uniform distribution. Nonetheless, it is often argued that dynamical interactions in a star cluster would quickly thermalize the binaries, which could justify imposing a thermal eccentricity distribution at birth for all binaries. In this paper we investigate the validity of this assumption. We develop our own rapid semi-analytic model for binary evolution in star clusters, and also compare with detailed N-body and Monte Carlo star cluster models. We show that, for nearly all binaries, dynamical encounters fail to convert an initially uniform eccentricity distribution to thermal within a star cluster's lifetime. Thus, if a thermal eccentricity distribution is observed, it is likely imprinted upon formation rather than through subsequent long-term dynamical processing. Theoretical investigations that initialize all binaries with a thermal distribution will make incorrect predictions for the evolution of the binary population. Such models may overpredict the merger rate for binaries with modest orbital separations by a factor of about two.

Figures

Equations

$$f(e)=2ede.$$

$$a_{\text{hs}}=\frac{G{m}_1{m}_2}{2m_3{\sigma }_0^2},$$

Full text

I. Introduction

If a population of binaries undergoes enough dynamical encounters that exchange energy, one may expect the population to reach a state approximating energy equipartition, and the distribution of energies to follow a Boltzmann distribution. This “thermalized” population of binaries was first investigated by Jeans (1919), who derived the resulting distribution of binary orbital eccentricities (and periods) for such a population in statistical equilibrium. The result, reached also by Ambartsumian (1937), Heggie (1975) and Kroupa (2008), is that the eccentricities should relax to the distribution function:

[TABLE]

All values of $e^{2}$ would be equally likely. In other words, the resulting distribution has many more high-eccentricity binaries than low-eccentricity binaries (see, e.g., Figures 2 and 3). This distribution is known as the “thermal eccentricity distribution”.

Such dynamical encounters occur far more frequently in star clusters than in the field. Therefore, one might expect that star clusters would be the best place to look for the thermal eccentricity distribution. Furthermore, observations suggest that most stars with masses $\gtrsim 0.5M_{\odot}$ were born in clusters (e.g. Lada & Lada, 2003), many of which quickly dissolve to populate the Galactic field. Therefore even binaries that are currently in the field may have undergone sufficient encounters in their birth environments to convert any primordial eccentricity distribution to thermal.

This line of reasoning, and the elegant formula, has elevated the thermal eccentricity distribution to be the distribution of choice for initial conditions in the majority of published star cluster models, population synthesis studies, and analytic investigations of binaries in both the field and in star clusters.

Ever since Jeans (1919), astronomers have searched observationally for this thermal distribution in various populations of binaries. However, in nearly all cases (including the study that Jeans 1919 compared against) the binaries are not observed to have a thermal distribution. To provide some recent catalogs, we direct the reader to the review article by Duchêne & Kraus (2013) and to Moe & Di Stefano (2017), and references therein.

In summary, current observed samples of binaries with primary stars between $\sim$0.8$M_{\odot}$ and $\sim$5$M_{\odot}$ have eccentricity distributions that are flatter than a thermal distribution. For background, the seminal Duquennoy & Mayor (1991, DM91) study divided their sample of solar-type field binaries at an orbital period of 1000 days, finding a bell-shaped distribution for the short-period sample, and, after a significant incompleteness correction, an indication of a thermal distribution for the long-period sample. More recently, the DM91 study has been superseded by Raghavan et al. (2010) and also Moe & Di Stefano (2017). The Raghavan et al. (2010) study gathered a complete volume-limited sample of solar-type stars. Using the then newly available Hipparcos data, they found that the DM91 sample was contaminated by parallax errors; specifically, that 44% of the DM91 sample actually lie outside the DM91 selection criteria, 38% of the stars that meet their criteria with current observations were not included in the DM91 sample, and several stars were erroneously included in the DM91 sample due to incorrect spectral-type assignments in Hipparcos. In short, the Raghavan et al. (2010) study clarified and improved the completeness of the solar-type binary sample in the solar neighborhood. Their analysis of the eccentricity distribution shows a uniform distribution for all binaries with periods longer than the circularization period, with no significant difference when cutting at an orbital period of 1000 days. Moe & Di Stefano (2017) confirm this result in a very careful and thorough analysis, and find that for late-type binaries (even those with orbital periods $3<\log P<5$), the eccentricity distribution is discrepant with, and flatter than, a thermal eccentricity distribution. Studies of solar-type binaries in open clusters come to a similar conclusion, and don’t observe a thermal eccentricity distribution, regardless of orbital period, out to the completeness limits (e.g. Geller & Mathieu, 2012; Geller et al., 2013). The conclusion is that late-type binaries are observed to follow more closely to a uniform distribution, for binaries with periods beyond the reach of tides.

Furthermore, the review by Duchêne & Kraus (2013) concludes that for all the samples investigated (of all spectral types, and observational methods), the eccentricity distributions are all inconsistent with thermal. The Moe & Di Stefano (2017) study supports this finding for all but the visual binaries with O5 - B5 primary stars and periods between $3.6<\log P($days$)<4.6$ (10-100 years), from the Malkov et al. (2012) catalog. Moe & Di Stefano (2017) find this sample to have an eccentricity distribution that is consistent with thermal. Conversely, Duchêne & Kraus (2013) investigate spectroscopic binaries in the SB9 catalog and the catalogs of Abt (2005) and Sana et al. (2012), and find the OB stars with $2<\log P($days$)<4$ have an eccentricity distribution inconsistent with thermal (and closer to uniform), though they note that the catalogs are likely incomplete (and, likely biased toward low eccentricities; Geller & Mathieu, 2012). Moe & Di Stefano (2017) also confirm a flatter than thermal eccentricity distribution for shorter-period ($P=10-500$ days) early-type spectroscopic binaries, after correcting for selection effects.

Keeping the various observational biases and incompleteness in mind, the general conclusion is that most observed binary samples have eccentricity distributions that are flatter than thermal, and more closely consistent with a uniform distribution.

Nevertheless, in dynamical star cluster models, the thermal eccentricity distribution is often imposed on the primordial population. Some authors assume that encounters will thermalize the binary population so quickly that any adjustment time could be neglected. We test this hypothesis in this paper. Other authors assume that the thermal eccentricity distribution is imposed by the binary-formation process. Some go further to apply a numerical “eigenevolution” prescription (Kroupa, 1995; Belloni et al., 2017) prior to $N$-body or population synthesis modeling, that will flatten a thermal eccentricity distribution for binaries with shorter periods, and more closely approximates the results from DM91. Cluster dynamical processes that then act on this “eigenevolved” population of binaries, may then change the eccentricity distribution again into what is ultimately observed.

Eccentricity can be excited dynamically through close strong encounters or long-range flybys. One strong encounter can dramatically change the eccentricity, while flybys contribute in a more cumulative long-term manner. For strong encounters, exchanges and dynamical captures may be an efficient method to induce thermal eccentricities. For example, Kouwenhoven et al. (2010) and Perets & Kouwenhoven (2012) show that wide dynamically captured binary stars and star-planet systems (respectively), formed within $N$-body star cluster simulations, have a thermal eccentricity distribution. Fregeau et al. (2004) show that binary-single scattering experiments that result in a binary containing a merger product, may exhibit a thermal (or similar to thermal) eccentricity distribution. The common thread in these particular references is that binaries formed through dynamical exchanges or captures tend to have thermalized eccentricities.

All else being equal, binaries in wider orbits experience more encounters, and therefore, assuming encounters lead to thermal eccentricities, wide binaries may be the easiest to thermalize. However, within a star cluster, the widest binaries are “soft”; the binding energy of a soft binary is less than the typical kinetic energy of a star moving at the velocity dispersion of the cluster. Therefore soft binaries are often disrupted during encounters. For wide binaries, the thermalization of the eccentricity distribution is a race between eccentricity excitation and disruption through encounters.

In this paper, we model these processes, searching for the emergence of a thermal eccentricity distribution. We will show that if the initial eccentricity distribution is far from thermal, this hypothesis is incorrect. Kroupa (1995) and Kroupa & Burkert (2001) also challenged this hypothesis, and showed that if all binaries are born with the same eccentricity, their set of $N$-body star cluster models did not thermalize the binaries. We build upon these results using a different method and a more empirically motivated initial eccentricity distribution.

The outline of this paper is as follows. In Section II, we describe our semi-analytic model for evolving binaries within a cluster environment. In Section III, we compare and validate this semi-analytic model with more detailed $N$-body and Monte Carlo models, and also investigate the solar-type main-sequence (MS) binaries in these detailed models for a thermal eccentricity distribution. In Section IV, we explore a grid of semi-analytic models in search of the parameter space that could produce a thermal eccentricity distribution in star clusters. Then in Section V, we shift focus toward remnant binaries; this population may be the best place to look for a thermal eccentricity distribution, due to their more rich dynamical histories. In Section VI, we investigate the impact of the initial eccentricity distribution on the binary merger rate. Finally, in Sections VII and VIII, we discuss these results and provide our conclusions.

II. Semi-analytic Model for Binary Evolution

There are three primary effects that work together to change a given binary star’s eccentricity and semi-major axis over time: strong encounters, flyby encounters, and the internal evolution of the binary. Here we will define strong encounters as those where an incoming star (or binary) has a pericenter distance that is less than or equal to the semi-major axis of the binary. Complementarily, we define flybys as those where an incoming star (or binary) has a pericenter distance beyond the binary’s semi-major axis. The internal evolution of the binary can include changes based on stellar evolution, mass transfer, tides and magnetic braking, etc. In our model, we include the effects of tides and magnetic braking, but choose to exclude other effects such as stellar evolution and mass transfer (though we remove binaries that cross the Roche radius, see Section II.3).

The encounter parameters depend on the time evolution of the cluster. Specifically, we require the time-varying number of stars, total cluster mass, mean stellar mass, half-mass radius, core radius, central density, central velocity dispersion, escape velocity and binary fraction. These can be derived from, for example, a detailed $N$-body model, or a more rapid analytic model like EMACSS (Alexander & Gieles, 2012; Gieles et al., 2014; Alexander et al., 2014).

Once the cluster parameters are known, our method (described below) can enable very rapid calculations of the time-varying distributions in binary semi-major axes and eccentricities. The main bottleneck is integration of the differential equations related to tides, magnetic braking and flybys. Nonetheless, this method is “embarrassingly parallel”; each binary can be evolved on its own and in parallel with all others, in a similar manner as done for population synthesis models. This offers a significant speed up when compared with detailed $N$-body and Monte Carlo cluster models (with remarkably similar results; see Figures 2 and 3).

Previous authors have developed other analytic models for the dynamical evolution of a binary population. Notably, Marks et al. (2011) constructed a numerical transfer function, based on the first 5 Myr of $N$-body evolution to approximate the dynamical evolution of binaries. Marks & Kroupa (2011) applied this model to dissolving clusters to reproduce the field population (focusing on DM91 observations). Also Giersz et al. (2016) showed close agreement between the Marks et al. (2011) method and MOCCA models in the first 5 Gyr. Sollima (2008) also developed an analytic prescription for dynamical processing of binaries in globular clusters, though did not include changes to eccentricity. Our model focuses specifically on the evolution of the eccentricity and semi-major axis (and period) distributions, and is designed to test the hypothesis that encounters can quickly thermalize a population of binaries. We explain our method in detail below.

II.1. Strong Encounters

We use the binary-single and binary-binary encounter timescales from Leigh & Sills (2011) to estimate the time until an incoming single or binary star will pass within the semi-major axis of the binary of interest. For the encounter parameters (e.g., density, velocity dispersion, etc.), we take the values within the cluster core. If the timescale for such an encounter is less than the cluster age, we assume that a strong encounter occurs for that binary. Traditionally, these strong encounters are evolved using direct $N$-body scattering codes, such as FEWBODY (Fregeau et al., 2004). Such codes can be very efficient, but often this is a relatively time consuming calculation. There is a long history of methods and models to speed up such calculations, going back at least $\sim$30 years to the classic Fokker-Plank paper by Gao et al. (1991). (We will not attempt to summarize this body of work here.) Recently, Valtonen & Karttunen (2006) and Leigh et al. (2016, 2018) developed an analytic approximation to the statistical distribution of the outcomes from stellar encounters. We follow the same procedure as Leigh et al. (2018), which incorporates a Monte Carlo sampling of the encounter parameter space to estimate the final semi-major axis and velocity of the binary, given randomly chosen initial encounter parameters.

Specifically, our procedure is as follows. We first select a three-dimensional velocity for the binary of interest and the incoming object (either single or binary, based on the encounter time). We draw these velocities from a lowered-Maxwellian (King, 1965), which depends on the time-evolving properties of the cluster model (see above), and the mass of the incoming object. We then project the velocities of the binary of interest and the incoming object to determine the relative velocity.

For the mass of an incoming single star, we take the time average of the cluster’s mean stellar mass, weighted by the number of stars in the cluster, over the time spanning from the previous encounter (or the start of the model) up until the current encounter time. We weight by the (time-evolving) number of stars in an attempt to capture the most likely incoming star; in practice this does not dramatically change the mass we derive, mostly because the mean stellar mass is relatively constant over the course of the cluster evolution (at $\sim 0.3M_{\odot}-0.6M_{\odot}$). If the timescales instead predict a binary-binary encounter, we multiply this mass by $(1+q)$, where $q$ is a mass ratio drawn randomly from a uniform distribution. Note that this does not account for mass segregation (which may cause the mean mass of incoming objects in the core to be higher than that averaged over the entire cluster). Accounting for this level of detail is beyond the scope of our model, and as we show below, even with these simple assumptions, our model agrees closely with more sophisticated simulations.

Given the parameters of the target binary, incoming star or binary, the incoming velocity, and an assumed pericenter distance equal to the current binary semi-major axis, we follow the method of Leigh et al. (2018) and Valtonen & Karttunen (2006) to draw a random final velocity for the incoming object. This final velocity depends on the binary’s energy, angular momentum, semi-major axis, and eccentricity (prior to the encounter), and the masses of all objects. Assuming energy is conserved, and given this final velocity of the incoming object, we can calculate the binary’s final energy and therefore its final semi-major axis.

For the final eccentricity, we simply assume that resonant encounters will result in a random draw from the thermal distribution, and non-resonant encounters will not change the eccentricity. This assumption for resonant encounters stems from the numerical results, discussed in Section I, that dynamically formed binaries show a thermal eccentricity distribution. During resonant encounters, there are often exchanges, and the final binary will likely have similar properties to a dynamically formed binary. (Our simplified treatment does not account specifically for exchanges.) We estimate that an encounter may be resonant if the relative velocity, between the binary of interest and the incoming object, is less than the critical velocity, as defined in Fregeau et al. (2004), Equation 1. We randomly select half such encounters to result in a thermal eccentricity; the remaining half do not alter the eccentricity. This is, of course, an over-simplification, but we will show that this procedure qualitatively reproduces the results of detailed $N$-body cluster models in Section III.

Encounters with a relative velocity larger than the critical velocity are not expected to be resonant. Indeed, most such encounters disrupt the binary (e.g. Fregeau et al., 2004). For the fraction of these encounters that don’t disrupt the binary, we choose not to modify the eccentricity. In reality, these encounters likely do change the eccentricity, but there is no simple model to describe the resulting eccentricity change. Nonetheless, most such surviving binaries become wider, and therefore more susceptible to disruption from the next encounter.

This method is appropriate for evolving a distribution of binaries, though it is not expected to exactly reproduce the outcomes of detailed $N$-body calculations on a per-binary basis.

II.2. Flyby Encounters

Between the individual strong encounters, we assume that many flybys are ongoing. For flybys, we use the results from Heggie & Rasio (1996), who calculate analytic cross sections ($\Sigma$) for a given change in eccentricity resulting from the combined effects of weak perturbations to a given binary by other objects at pericenter distances beyond the binary’s semi-major axis. These cross sections depend on the component masses, semi-major axis and eccentricity of the binary, the typical mass of objects in the cluster, and the velocity at infinity. We assume the typical mass of an object in the cluster is equal to the cluster’s (time-evolving) mean stellar mass, as described above. We assume the velocity at infinity, $v$, to be equal to the time average of the central velocity dispersion of the cluster, weighted by the (time evolving) number of stars (in a similar manner as for the mean stellar mass). We estimate the central density, $n$, in the same manner as $v$. We then use a simple encounter rate estimate, $\Gamma=n\Sigma v$, along with the duration of time between encounters (or time until the first encounter, or until the desired end time, if appropriate), to solve for the expected change in eccentricity for a given binary.

We begin with all non-zero eccentricities, and therefore begin using the cross section from Equation 19 in Heggie & Rasio (1996). For a non-zero initial eccentricity, there is an equal probability that the change in eccentricity resulting from a flyby will be positive or negative. To determine the sign of this change in eccentricity, we simply draw a random number from a uniform distribution, and assign a negative change in eccentricity if the random draw is $>0.5$. If the evolution of a given binary produces a circular orbit, we then switch to Equation 25 in Heggie & Rasio (1996), with the appropriate time duration remaining in the encounter window. (For simplicity here, we only study binaries with equal masses in the semi-analytic model, and therefore only the exponential regime in Heggie & Rasio 1996 can produce a non-zero change in eccentricity for a circular binary). For an initially circular binary, the change in eccentricity resulting from encounters can only be positive. Rasio & Heggie 1995 show that this regime can be very important for compact-object binaries.

Unlike the strong encounters, we assume that this change in eccentricity occurs smoothly over the full time duration between strong encounters. We therefore divide the expected change in eccentricity resulting from these calculations by the encounter duration to estimate an $\dot{e}$ from flybys, which we use along with the $\dot{e}$ from tides (see below), to integrate the total change in eccentricity between encounters.

II.3. Internal Evolution: Tides and Magnetic Braking

We integrate the differential equations from the Hut (1981) weak friction tide model coupled to the differential equations from Hurley et al. (2002) for magnetic braking, and the $\dot{e}$ from flybys (Section II.2) to estimate the changes in eccentricity, semi-major axis and spin between encounters. We follow a similar method to Hurley et al. (2002) to calculate the parameters for the Hut (1981) equations. However, as we are not including stellar evolution, we make some additional order-unity simplifications (which we can consider as contributing to the overall uncertainty in the strength of tides from this model on stars with convective envelopes); namely, we set the envelope mass and radius equal to the stellar mass and radius, respectively, the initial stellar spins to $2\pi/p_{0}$ (where $p_{0}$ is the initial binary orbital period), and the radius of gyration to that of a solid sphere ($r_{g}=2/5$). We also allow for a multiplicative factor that can be applied to the $\dot{a}$ and $\dot{e}$ differential equations for tides; as suggested by, e.g., Belczynski et al. (2008) and Geller et al. (2013), this factor may be as high as 50 - 100 in order to match observed solar-type binary systems. By default, and for simplicity in the majority of the paper (unless otherwise noted), we set this value to unity.

In principle, we could also include the various changes to the binary orbital parameters, stellar spins, and stellar masses, resulting from stellar evolution by making use of a rapid binary evolution code, such as BSE (Hurley et al., 2002). That is beyond the scope of this work. Furthermore, for the majority of this paper, we will focus on binaries with component masses of $\lesssim$1 $M_{\odot}$, which (for solar metallicity) will remain on the MS for $\gtrsim$11 Gyr.

II.4. Collisions, Coalescence and Disruptions

Changes to the binary orbital elements resulting from strong encounters, flybys and tides (as described above) can eliminate a binary from the population by direct collisions (due to high eccentricity), coalescence and disruptions. In our model, this can happen in a variety of ways.

Flybys can cause the eccentricity to increase enough that the two stars in the binary physically collide or orbit inside the respective Roche radius ($r_{L}$, defined here in the same units as the semi-major axis; Eggleton, 1983).

The semi-major axis, $a$, of a binary can be reduced both by strong encounters and through tides with magnetic braking. We assume the components of the binary coalesce (merge) if $a(1-e)<r_{L}$, and that the components would collide if $a(1-e)<r_{*}$, (where $r_{*}$ is the star’s radius).

The semi-major axis of a binary can be increased by strong encounters. This is particularly relevant for soft binaries, as described in Section I. Soft binaries have low binding energies relative to the typical kinetic energies of the cluster stars. Therefore the typical outcome for an encounter involving a soft binary is for the incoming star to donate energy to the binary, widening the binary’s orbit. If the incoming star donates enough energy, it will disrupt the binary. Within our strong encounter approximation, we assume that a binary is disrupted if the final binary energy is positive. (In principle, tides can also increase the semi-major axis with appropriate spins, but in practice this is not relevant in our model.)

II.5. Outcomes for Individual Binaries

As a demonstration of the full semi-analytic model described above, we run 100 draws for binaries with three different initial semi-major axes, of 1, 5 and 10 AU respectively, all with an initial eccentricity of 0.5 and with component masses of $m_{1}=m_{2}=1M_{\odot}$, and show the results in Figure 1. We evolve these binaries within a cluster with initially 105 stars. For the strong encounter rates, we assume a binary fraction of 50%.

For reference, we can estimate the hard-soft boundary, $a_{\text{hs}}$, of a cluster using the virial theorem,

[TABLE]

where $m_{3}$ is the mass of the incoming object (either a single star or the combined mass of a binary; a typical mass for an evolved star cluster is $\sim 0.5M_{\odot}$), and $\sigma_{0}$ is the initial (three-dimensional) central velocity dispersion. Our use of the initial velocity dispersion is supported by detailed $N$-body models that show that the hard-soft boundary is imposed early on (Geller et al., 2013, and see also Section III in this paper). A cluster of this size has a hard-soft boundary at $\sim 10$AU, and is predicted by EMACSS to dissolve at slightly beyond 8 Gyr.

The strong encounter times here are identical for all binaries of a given semi-major axis, but the parameters of each encounter are chosen randomly (from the appropriate distributions, as described in Section II.1). By running many realizations of these binaries, we can study the distribution of final semi-major axes and eccentricities, shown in the horizontal histograms on the right side of the figure.

For tighter binaries (e.g., the 1 AU binaries here), encounters are too infrequent to significantly change the binary orbital parameters. As the semi-major axis approaches the hard-soft boundary, encounters will often change the binary orbital parameters. For the wider binaries, Figure 1 shows that the eccentricity tends to smoothly increase or decrease due to flybys and tides, with the occasional jump in semi-major axis, and possibly also eccentricity (depending on the relative velocity of the stars), at the times of strong encounters.

The 5 AU binaries in this cluster are hard binaries, and, as expected, encounters tend to decrease the semi-major axes. Conversely, the 10 AU binaries here are nearly soft, and, as also expected, the encounters tend to increase the semi-major axes. For both the 5 and 10 AU binaries, encounters spread out the eccentricity distribution. Many of the wider binaries are also removed from the population through collisions, coalescence and disruptions, as described in Section II.4.

Though we focus this paper on the distributions of binary orbital parameters, we remark here on an interesting and relatively common outcome of the interplay between flybys and tides. Flybys can gradually increase a binary’s eccentricity. The “spikes” seen in Figure 1, when a binary changes from an increasingly high eccentricity to a decreasing eccentricity, are times when tides begin to dominate over flybys. Wide binaries in particular may be driven quickly to very high eccentricity, and therefore very small pericenter distance, by flybys. As the pericenter distance decreases, tides become more important, and when combined with magnetic braking, can drive the binary to coalescence. This is clearly very relevant for the initially 10 AU binaries shows in Figure 1. In general, wide binaries may be particularly susceptible to high-eccentricity driven mergers (or collisions) due to eccentricity pumping from flybys in star clusters, which may lead to the production of exotic stars like blue stragglers and sub-subgiants (Leonard, 1989; Leigh & Sills, 2011; Giersz et al., 2013; Kaib & Raymond, 2014; Leiner et al., 2017; Geller et al., 2017a, b).

III. Validation by Comparisons to $N$-body and Monte Carlo Star Cluster Models

In the following, we test our semi-analytic model by comparing to more detailed star cluster simulations. In Section III.1, we compare against an $N$-body open cluster model from the literature created using the nbody6 code. In Section III.2, we compare against a Monte Carlo globular cluster model created using the MOCCA code.

For both comparisons, we use the time-evolving cluster parameters directly from the $N$-body and Monte Carlo models, respectively, as input to our semi-analytic models, to, e.g., define the encounter timescales. As we show below, our semi-analytic model does a remarkable job of reproducing the dynamical evolution of binaries in more detailed models. These comparisons show that we can indeed use this simpler (and faster) approach to make confident statements about the dynamical generation of the thermal eccentricity distribution in star clusters.

III.1. NGC 188 $N$-body Model

Geller et al. (2013) produced a detailed $N$-body model of the old (7 Gyr) open star cluster NGC 188, using nbody6 (Aarseth, 2003). The initial conditions for the model were based on empirical data for the cluster, and particularly for the binary population. The $N$-body model was initialized with 39,000 stars, a half-mass radius of 4.6 pc, and a galactocentric radius of 8.5 kpc. (We refer the reader to Geller et al. 2013 for details on other initial conditions and parameters of their model.) The NGC 188 $N$-body model matches the observed solar-type binaries in the real cluster quite well, for binaries with periods $\lesssim 3000$ days (the completeness limit of the Geller & Mathieu 2012 survey).

We use the cluster structure, velocity and mass parameters from the NGC 188 $N$-body model as input to our semi-analytic model, and show the results in Figure 2. For the semi-analytic model, we create $10^{5}$ binaries, with $m_{1}=m_{2}=1M_{\odot}$, and orbital periods and eccentricities drawn from the same distributions functions as used for the $N$-body model. We also increase the tidal circularization rate (specifically, the differential equations for $\dot{a}$ and $\dot{e}$) by a factor of 50 (Belczynski et al., 2008; Geller et al., 2013), to approximate the similar increase imposed in the $N$-body model. For the $N$-body model, we combine all 20 unique simulations in the NGC 188 model, and select only binaries where the combined mass of the two components is $>=1.5M_{\odot}$ and both are MS stars. (The cluster turnoff at the age of NGC 188, and at the cluster’s roughly solar metallicity, is $\sim$1.1$M_{\odot}$.) This mass limit is meant to include binaries with similar masses, and therefore similar dynamical histories, to the equal-mass $m_{1}=m_{2}=1M_{\odot}$ binaries in our semi-analytic model, while also allowing for a large enough sample size. We exclude binaries from the $N$-body model that were dynamically formed (through exchanges, three-body formation, etc.), by only including binaries that were paired at the start of the simulation.

On the left-hand side of Figure 2, we plot the distributions in eccentricity and period for the $N$-body model, in red, and our semi-analytic model in black and white. The final distributions are shown in solid lines (while the initial conditions are shown in dashed lines). The agreement between the $N$-body and semi-analytic models is encouraging. In particular we note the truncation of the binary period distribution at the hard-soft boundary. Note that neither the $N$-body nor the semi-analytic model, converts the initially Gaussian eccentricity distribution to thermal. Indeed the overall eccentricity distribution hardly changes throughout the entire cluster lifetime.

On the right-hand side of Figure 2, we show cumulative eccentricity distributions in bins of final orbital period (in the log). Here we see that the eccentricity distributions trend toward thermal for the wider binaries. (And also that the short-period binaries become more circular, due to tides). However, even the widest binaries do not achieve a thermal distribution.

III.2. MOCCA Globular Cluster Model

Next we compare our semi-analytic model with a MOCCA Monte Carlo globular cluster model (Hypki & Giersz, 2013; Giersz et al., 2013, 2015; Askar et al., 2017). This MOCCA model was initialized with $1.2\times 10^{6}$ total objects111One object is either a single star or a binary., ($9.13\times 10^{5}M_{\odot}$), a half-mass radius of 1.2 pc, and a 30% binary fraction. Masses of single stars and the primary stars of binaries are drawn from a Kroupa (2001) initial mass function in the range of $0.08-150M_{\odot}$. The binaries were initialized with a log-uniform initial semi-major axis distribution (up to 200 AU and such that the initial pericenter distance is more than 4 times the radius of a 0.08 $M_{\odot}$ star) and a uniform initial eccentricity distribution (see also Figure 3). The secondary masses in binaries are chosen from a uniform mass-ratio distribution, such that the IMF is preserved (Oh & Kroupa, 2016). Internal binary evolution uses an upgraded version of BSE (Hurley et al., 2002; Belloni et al., 2018; Giacobbo et al., 2018). The cluster was evolved in a standard tidal field for the solar neighborhood for 12 Gyr. The global properties of this model at 12 Gyr (total mass of $\sim 4\times 10^{5}M_{\odot}$ and half-light radius of $\sim 3$ pc) are broadly representative of the global properties of relatively massive globular clusters in the Milky Way (Harris, 2010).

We use the cluster structure, velocity and mass parameters from the MOCCA model as input into our semi-analytic model and show the results in Figure 3, in the same format as for the NGC 188 $N$-body model in Figure 2. Again, for simplicity in the semi-analytic model, and for more easy comparison to Figure 4, we include only equal mass binaries with both components at 1$M_{\odot}$. (We do not increase the tidal circularization rate, as we did in the comparison to the NGC 188 $N$-body model.) From the MOCCA model, we include binaries with two MS stars that have a combined mass between 1$M_{\odot}$ and 2$M_{\odot}$, and again we only include binaries that were paired in the initial population (to exclude dynamically formed binaries), for Figure 3. The turnoff in the MOCCA model at 12 Gyr (and at $Z=0.001$) is $\sim 0.85M_{\odot}$; therefore the most massive MS binaries in the MOCCA model at 12 Gyr are still slightly less massive than those in our semi-analytic model. This will introduce some inconsistencies between the encounter times in the MOCCA and semi-analytic models. Nonetheless, as is clear from Figure 3, the agreement between these models is very close.

As in the comparison to the NGC 188 $N$-body model, our semi-analytic approach reproduces the shape of the final period distribution remarkably well (and note that the $N$-body and MOCCA models began with different initial period distributions). The disruption of binaries, and the definition of the hard-soft boundary is captured faithfully in the semi-analytic model. The overall eccentricity distributions are also very close.

There is a slight excess of circular binaries in the semi-analytic model; this is primarily due to tides, which are apparently more efficient in our model than in the MOCCA model (perhaps due, at least in part, to the simplifying assumptions we’ve made for the stellar structure in the semi-analytic model). Here, we are more concerned with the long-period binaries, beyond the reach of tides.

The overall result from this comparison (and the comparison to the NGC 188 $N$-body model) is that dynamical evolution in the star cluster will not convert a uniform (or Gaussian, as in the NGC 188 model) eccentricity distribution into a thermal distribution within the cluster lifetime.

IV. Parameter Space Search for Thermalization

In the previous sections, we describe our semi-analytic model and validate the results by detailed comparisons against $N$-body and Monte Carlo cluster models. Here, in order to efficiently cover a larger parameter space in cluster initial conditions, we switch to the EMACSS code (Alexander & Gieles, 2012; Gieles et al., 2014; Alexander et al., 2014). EMACSS can rapidly evolve a cluster, given an initial number of stars, half-mass radius and galactocentric distance, and provides an estimate for the time evolution of the number of stars, total stellar mass, and half-mass radius (among other parameters).

We estimate additional cluster parameters based on a Plummer model (Plummer, 1911). We also use the output from EMACSS with the results from Webb & Leigh (2015), to estimate an appropriate dynamically modified mass function at each EMACSS output time, which then provides an estimate of the mean stellar mass in the cluster. Finally we estimate the (assumed constant) binary fraction starting by first taking a 50% total binary fraction, as appropriate for field solar-type binaries (Raghavan et al., 2010), then truncating the period distribution at the hard-soft boundary (Equation 2), which thereby defines a new binary fraction dependent on the central cluster velocity dispersion (in a similar manner as in Geller & Leigh, 2015). These values allow us to calculate the encounter timescales needed for our semi-analytic model. (Note, we verified this EMACSS-based approach by comparing to the MOCCA and $N$-body models from Section III, using the same initial conditions. The resulting models are nearly identical to those shown in Figures 2 and 3.)

We construct a grid of semi-analytic models covering a range in initial number of stars, half-mass radius (and therefore density), and binary component masses. We show results from this grid in Figure 4. In the top panel we show clusters of a range in initial number of stars, all with an initial half-mass radius ($r_{\textrm{hm,0}}$) of 2pc, and at a galactocentric radius ($r_{\textrm{g}}$) of 8.5 kpc. In the middle panel we show a cluster with initially 105 stars, $r_{\textrm{g}}=8.5$kpc, and with a range in densities (and velocity dispersions). In the top two panels, we only consider binaries with component masses each equal to 1$M_{\odot}$. In the bottom panel, we again show a cluster with initially 105 stars, $r_{\textrm{g}}=8.5$kpc, and $r_{\textrm{hm,0}}=2$pc , but we consider different component masses for the binary of interest. All cluster models are run using EMACSS and are evolved until either the cluster dissolves or for a Hubble time. (A cluster born with $\gtrsim 10^{5}$ stars, $r_{\textrm{hm,0}}=2$pc, and $r_{\textrm{g}}=8.5$kpc, will not dissolve in a Hubble time, which is responsible for the upturn in the solid and dashed lines in the top panel toward larger $N_{\mathrm{stars,0}}$.) We initialize the binaries in each cluster with a log-normal period distribution, and a uniform eccentricity distribution (as in Raghavan et al. 2010).

The results of this parameter space search are shown in Figure 4. On the left, the colored regions show the respective exponents, $\eta$, resulting from power-law fits ($p_{e}\propto e^{\eta}$) to the eccentricity distribution in bins of the semi-major axis. For reference, a thermal eccentricity distribution has $\eta=1$, and a uniform distribution has $\eta=0$. As is clear from the figure, when we include both strong encounters and flybys (left-most panel), no model in our grid converts a uniform eccentricity distribution to thermal. For some cluster and binary parameters the power-law exponent to the eccentricity distribution becomes larger than for a uniform distribution, but never reaches $\eta=1$. Interestingly, if we exclude flybys, there are a few regions of parameter space that do reach thermal (due to our assumption of imposing thermal eccentricities after sufficient strong encounters, see Section II.1). We return to this in Section VII below.

As an additional check, we can also use our semi-analytic approach to estimate timescales to convert a uniform eccentricity distribution into a thermal eccentricity distribution. The mean eccentricity in a thermal distribution is 2/3, while for a uniform distribution the mean eccentricity is 1/2. Therefore, the mean change in eccentricity required to convert the uniform to thermal distribution is 1/6. We can estimate the timescale for producing this change in eccentricity from flybys alone using the Heggie & Rasio (1996) cross section, averaged over a uniform distribution in initial eccentricity, along with our $\Gamma=n\Sigma v$ approximation. Since positive and negative changes to eccentricity are equally likely, here we simply double the time required to increase the eccentricity by 1/6. We then solve for the semi-major axis beyond which encounters are frequent enough to increase the eccentricity by 1/6 within a given timescale. In Figure 4, we set this timescale to the dissolution time of the cluster (as predicted by EMACCS, and defined at the time when there is less than 100 stars in the cluster), or a Hubble time if that is shorter, and show the result of this calculation with the dashed lines.

With the solid lines, in Figure 4, we show the semi-major axis beyond which we expect to have at least one strong encounter (where the incoming star has a pericenter passage within the binary semi-major axis). For our purposes here, we will simply assume that half of the strong encounters (regardless of whether they are expected to be resonant) result in a thermal eccentricity. Finally, the dotted lines in Figure 4 show an estimate of the semi-major axis at the hard-soft boundary (Equation 2).

For all other cluster parameters used in these timescale estimates we take the time average, weighted by the (time varying) number of stars from our EMACSS plus Plummer (1911) approximation, described above. We weight by the number of stars in an attempt to best approximate the cluster parameters when most encounters will take place and therefore when the majority of the changes to the eccentricity distribution will occur.

These timescale estimates suggest that a cluster could indeed thermalize a population of binaries (albeit a small population), near to the hard-soft boundary. Models including only strong encounters (middle panel in Figure 4) indeed trend toward thermal in this region of parameter space. However, the more detailed semi-analytic calculations, including both strong encounters and flybys refute this. For these models, near the hard-soft boundary we see a turnover in the power-law exponent extending to $\eta<0$. This is because higher eccentricity binaries are preferentially depleted, especially at wider semi-major axes. Take for example flybys involving a binary with an initial eccentricity of 0.5. If the eccentricity is increased by $\gtrsim 0.5$, the binary components collide. However if the eccentricity is decreased by the same amount (as is also possible through flybys), the binary will circularize, and then enter the (less efficient) flyby regime starting from zero eccentricity. This preferentially removes wide high eccentricity binaries from the population. Also recall that this region near the hard-soft boundary is expected to be severely depleted of binaries by strong encounters (e.g., see Figures 2 and 3).

On the right side of Figure 4, the color scale denotes the mean number of strong encounters ($N_{\mathrm{enc}}$) within each bin in semi-major axis, for models that include both strong encounters and flybys. (The numbers are nearly identical for models that exclude flybys). For much of the parameter space, $N_{\mathrm{enc}}<1$. For binaries with semi-major axes near the dashed lines in the figure, $N_{\mathrm{enc}}\sim 1$ (as expected). For wider binaries, and particularly those with massive components, the number of encounters can reach of order 10. If flybys are ignored, these binaries achieve a thermal eccentricity distribution. However, the contribution from flybys, which also increases toward wider and higher-mass systems, is sufficient to remove many of the high-eccentricity binaries from the distribution (through collisions and mergers).

V. Comparing to Remnant Binaries in a CMC Globular Cluster Model

Above, we show that it is difficult, if not impossible, to convert a uniform eccentricity distribution to thermal for MS binaries in the lifetime of a cluster through dynamics. Here we also compare to remnant binaries using a Cluster Monte Carlo (CMC) globular cluster model from the Northwestern group (Joshi et al., 2000; Joshi et al., 2001; Fregeau et al., 2003; Fregeau & Rasio, 2007; Chatterjee et al., 2010; Pattabiraman et al., 2013; Chatterjee et al., 2013; Rodriguez et al., 2018). This cluster was initialized with $8\times 10^{5}$ total objects ($4.85\times 10^{5}M_{\odot}$), a half-mass radius of 0.81 pc, and a binary fraction of 5%. Masses of single stars and the primary stars of binaries are drawn from a Kroupa (2001) initial mass function in the range of $0.08-150M_{\odot}$. For a given binary, the secondary mass is assigned by drawing from a uniform distribution in mass ratio. Binary orbital periods are drawn from a log-uniform distribution with cut-offs at five times the distance to physical contact of the components and the hard-soft boundary. Binary eccentricities are drawn from a uniform distribution (between zero and one). We evolve the cluster well past a Hubble time to see if the binaries will eventually thermalize. In Figure 5, we show both the solar-type MS and black-hole (BH) binary eccentricity distributions at different snapshot times from the model.

In the top panel, we show only MS binaries that were paired primordially, excluding binaries formed dynamically (as we also do for Figures 2 and 3). In the second panel from the top we include all MS - MS binaries (within our mass and period range, as defined in the Figure 5 caption). In the third panel from the top, we include only the MS - MS binaries (within our mass and period range) that have at least one component which went through a dynamical encounter that was strong enough to be sent through the direct $N$-body integrator, FEWBODY (see references above for the criteria to initiate a direct $N$-body calculation of an encounter in CMC).

In all the top three panels, we see that the eccentricity distribution trends further toward thermal as time progresses. Comparing the top two panels shows that the dynamically formed binaries are generally formed at higher eccentricities than those in the primordial population (at any given time). Likewise, those that went through a FEWBODY encounter in CMC have eccentricity distributions shifted to higher values than for the primordial population. Indeed, as we discussed in Section I, dynamically formed binaries are expected to form with thermalized eccentricities. However, even after 36 Gyr of dynamical evolution, the solar-type MS binaries do not reach a thermal eccentricity distribution.

In the bottom panel of Figure 5, we show binaries containing at least one BH. For compact object formation, CMC uses a modified prescription from that implemented in SSE (Hurley et al., 2000) and BSE (Hurley et al., 2002) by using the results of Fryer & Kalogera (2001) and Belczynski et al. (2002). Natal kicks for core-collapse neutron stars are drawn from a Maxwellian with dispersion width $\sigma=265\rm{km\,s}^{-1}$ (Hobbs et al., 2005). BHs are assumed to form with significant fallback and BH natal kicks are calculated by sampling from the same kick distribution as used for neutron stars, but reduced in magnitude according to the fractional mass of fallback material (see Morscher et al., 2015, for more details).

For the sample shown in Figure 5, we attempt to exclude any BH binary whose current eccentricity was significantly altered by mass-transfer or tides. Specifically, we exclude any binaries that have eccentricities of uniquely zero (as this results from BSE for binaries undergoing mass transfer), BH - MS binaries with periods $\leq 15$ days (roughly the tidal circularization period for old stellar populations; Meibom & Mathieu, 2005), and all BH - giant binaries. Most of the BH binaries form early, and are ejected quickly from the cluster. The one early time step shown in Figure 5 is the only time step that contains sufficient BH binaries for which to construct an eccentricity distribution. A Kolmogorov-Smirnov test shows that we cannot distinguish the eccentricity distribution for these BH binaries from a thermal eccentricity distribution.

About two thirds of the BH binaries in Figure 5 contain one BH and one non-BH (i.e., MS star or other remnant). The remainder are BH - BH binaries. A primary method for forming BH - BH binaries in the CMC models (and also in the MOCCA models) is through “three-body binary formation” (Morscher et al., 2015). In these models an analytic approximation is made to account for encounters involving three initially single stars that experience a close encounter resulting in one bound binary. In CMC models this mechanism is only applied to remnants. The resulting binary is assumed to have an eccentricity drawn from a thermal distribution. About half of the BH binaries (regardless of the stellar type of the partner) in the timestep shown in the bottom panel of Figure 5 formed through the three-body binary mechanism and retain the same partner. About 10% retained the primordial partner. The remaining $\sim$40% formed through dynamical exchanges.

To investigate further for the thermal eccentricity distribution, we take all the eccentricities resulting from FEWBODY encounters in CMC, over the entire duration of the simulation ($\sim 40$ Gyr), shown in Figure 6. In this figure, we divide the sample into MS - MS binaries, and remnant - remnant binaries. (Remnants include white dwarfs, neutron stars, and BHs.) If an individual binary undergoes multiple encounters, each resulting eccentricity is included in the figure. Here again we see that the eccentricities of MS binaries, even those that have just undergone an encounter strong enough to require FEWBODY, do not obtain a thermal distribution.

However, the remnant binaries that go through encounters do come out with a thermal eccentricity distribution. $>$99% of the remnant binary encounter products shown here were formed dynamically. About one quarter were initially bound through the three-body binary mechanism (which assumes a thermal eccentricity); those born through the three-body binary mechanism underwent, on average, $\sim$6 encounters after formation and throughout the remainder of the simulation, which may be sufficient to erase the initial thermal assumption. The remaining roughly three quarters were formed by dynamical exchanges. We return to this in Section VII.

We also investigated a similar CMC model, but with the modification that BH kicks were drawn from the same Maxwellian as for neutron stars, with a dispersion of $\sigma=265\rm{km\,s}^{-1}$. In this model, nearly all BHs are ejected quickly, and the cluster undergoes core collapse much earlier in time, which creates a smaller core radius and presumably more frequent encounters for the MS stars. However, even in this model, the MS binaries’ eccentricity distribution (for stars in our selected mass range) does not reach thermal, even after 30 Gyr. Most of the MS stars in this mass range spend the majority of their lives outside of the cluster core, where the densities are similar between models with and without many BHs.

VI. Merger Rates

The assumption of thermalization dramatically overpopulates the high-eccentricity end of the binary distribution relative to a uniform distribution, and the consequences of choosing the wrong eccentricity distribution are not negligible. Take, for example the merger rate of solar-type binaries due to slow angular momentum loss from tides coupled to magnetic braking. We can use our model, and exclude encounters, to estimate the fraction of binaries with $m_{1}=m_{2}=1M_{\odot}$ and a range in initial orbital separations, that would merge in a Hubble time. In the top panel of Figure 7 we show the ratio of the number of such mergers when we draw eccentricities from the thermal distribution ($N_{t}$) over the number when eccentricities are drawn from a uniform distribution ($N_{u}$). At very small initial semi-major axes, nearly all binaries are expected to merge in a Hubble time, but toward wider binaries, only the most initially eccentric binaries, with the smaller pericenter distances, can merge. At these modest initial orbital separations, a population of binaries with eccentricities drawn from a thermal distribution has twice the merger rate to at population with eccentricities drawn from a uniform distribution.

Next we perform a similar investigation for BH - BH binaries merging due to gravitational wave radiation. Though we have found that BH - BH binaries formed dynamically in clusters are likely to have a thermal eccentricity distribution (due to their dynamical formation pathway), it is not immediately clear what is the expected eccentricity distribution for BH - BH binaries that form directly from primordial binaries in the field, or in star clusters. Perhaps it is closer to a uniform distribution. In the bottom panel of Figure 7 we use the Peters (1964) formulae to compare the merger rate for BH - BH binaries with initial eccentricities drawn from either a thermal or uniform distribution, over a range in initial semi-major axis. We show equal-mass BH - BH binaries, with masses of 10, 20 and 30 $M_{\odot}$. Again, for the tightest binaries, all are expected to merge in a Hubble time. However toward wider binaries, the merger rate for BH - BH binaries with eccentricities drawn from a thermal distribution is twice the rate of BH -BH binaries with eccentricities drawn from a uniform distribution.

The conclusion from this investigation is that choosing an initially thermal eccentricity distribution instead of a uniform distribution boosts the merger rate by a factor of about two, at modest orbital separations, for both MS - MS and BH - BH binaries.

VII. Discussion

Although a thermal eccentricity distribution may be the expected equilibrium outcome for a population of binaries that has undergone many energy-exchanging interactions, in real star clusters it is nearly impossible to reach this state. Naively, and based on encounter-time arguments alone, one might expect that the most likely binaries that could become thermalized are those with the most massive binary component masses, in wide orbits, and within the oldest, densest and most massive clusters. However, these assumptions break down. While it is true that wider binaries experience more frequent encounters, and therefore their eccentricity distribution is modified most rapidly, the widest binaries are soft, and soft binaries are easily disrupted. Encounters are most frequent in the most massive and dense clusters. However, the more massive and more dense the cluster, the tighter the hard-soft boundary becomes. Encounters occur more frequently for the most massive binaries. However, the oldest clusters (where sufficient time has past for many encounters to occur) don’t have massive MS binaries (due to stellar evolution), and therefore this regime is not easily observable.

Figure 4 shows this tension graphically. Under optimistic assumptions, simple timescales estimates (e.g., the lines in Figure 4) suggest that only binaries with periods ranging from the hard-soft boundary down to about one decade shorter, in typical clusters could achieve a thermal distribution within a cluster’s lifetime (if born with a uniform eccentricity distribution). The result is the same when considering strong encounters and/or flybys. However, as is clear from Figures 2 and 3, this period range just below the hard-soft boundary is expected to be significantly depleted of binaries, due to encounters that lead to disruptions, collisions, or mergers, and the remaining binaries are not thermalized. Indeed, our grid of semi-analytic models shows that this region near the hard-soft boundary is difficult to thermalize, due to the preferential removal of wide high-eccentricity binaries from the distribution by flybys. Furthermore, even within a model that evolves for many Hubble times, an initially uniform eccentricity distribution is not converted to thermal (Figure 5).

Thus the conclusion is that, for MS binaries, it is very difficult, and perhaps impossible, to convert a uniform eccentricity distribution to a thermal distribution through encounters. This finding is consistent with the observations showing solar-type binaries have a uniform (not thermal) eccentricity distribution (Duchêne & Kraus, 2013; Geller & Mathieu, 2012; Geller et al., 2013; Moe & Di Stefano, 2017). Indeed, for the majority of the solar-type binaries in our models, the birth eccentricity distribution is maintained throughout the cluster lifetime.

In contrast, the eccentricity distribution for binaries containing BHs in the CMC model we investigated (Section V), is consistent with thermal. A similar result was also found by Tanikawa (2013) in their $N$-body models. The high-mass regime investigated in Figure 4, is somewhat analogous to BHs in globular clusters. However, these models are evolved for the full cluster lifetime, while many high-mass objects are not retained in the cluster nearly this long, either due to rapid stellar evolution or dynamical ejections from encounters with other massive objects (neither of which are included in our semi-analytic models). For instance, in the CMC model (Figure 5), only the first snapshot contains sufficient BHs to construct an eccentricity distribution. This regime in parameter space may therefore be difficult to observe. Furthermore, given sufficient time, flybys can remove the high eccentricity binaries from the distribution; which may produce gravitational wave sources from high-eccentricity mergers or collisions of BHs, but will erase the thermal eccentricity distribution (e.g., compare the bottom-left and bottom-middle panels in Figure 4.

The post-encounter eccentricity distribution for binaries containing white dwarfs, neutron stars, and/or BHs in the CMC model is thermal (while a similar analysis of the post-encounter MS - MS binaries is inconsistent with thermal). Nearly all of these remnant binaries formed dynamically (through exchanges, three-body encounters, etc.), in a similar manner to the MS binaries in the models from Fregeau et al. (2004) and Kouwenhoven et al. (2010), which also showed thermal eccentricity distributions. For dynamically formed binaries, a thermal eccentricity distribution appears to be an appropriate choice (though this has not been verified observationally). However, we note that the CMC models do not explicitly account for eccentricity changes due to long-range flybys in the manner that we use in our semi-analytic model. These flybys may be important for high-mass objects that reside for Gyrs in dense regions of a cluster, and is worthy of further investigation.

It is also not immediately clear what is the expected eccentricity distribution for BH - BH binaries that form directly from primordial binaries in the field, or in star clusters. Furthermore, BH - BH binaries that form through the three-body mechanism in Monte Carlo star cluster models are assumed to be born thermalized; this assumption requires further verification.

One consequence for choosing a thermal eccentricity distribution, when a uniform eccentricity distribution is more appropriate, is to artificially boost the merger rate (for both MS - MS and BH - BH binaries) by roughly a factor of two (Section VI). This is most relevant for binaries at modest orbital separations. Very tight binaries are expected to merge regardless of the eccentricity. However, wider binaries may only merge if they are highly eccentric, and therefore have a small pericenter distance.

Triples can also affect the eccentricities of binaries caught, even briefly, in a secular resonance (Fabrycky & Tremaine, 2007; Naoz & Fabrycky, 2014), and in principle similar ”Kozai-Lidov type” secular effects could be produced with a binary coupled to a central IMBH, or perhaps simply to the rest of the cluster. Though not included here, the full evolution of the binary eccentricity distribution should account for triples and related secular effects.

Finally, our model highlights an additional path for stellar mergers, through wide binaries driven to high eccentricity by flybys that either results in a physical collision, or a coalescence through tides and magnetic braking (Figure 1). This mechanism is similar to that discussed by Kaib & Raymond (2014) for field binaries, and may contribute to the population of exotic stars like blue stragglers and sub-subgiants in star clusters (Leonard, 1989; Leigh & Sills, 2011; Giersz et al., 2013; Leiner et al., 2017; Geller et al., 2017a, b). We save a more in depth investigation into the rate of such mergers for a separate paper.

VIII. Conclusions

In this paper we develop and utilize a semi-analytic model for the dynamical evolution of binaries in star clusters. The model can be considered a “population synthesis” approach for cluster dynamics, in the sense that each binary is evolved in parallel, and results can be obtained significantly more rapidly than for more detailed $N$-body or Monte Carlo star cluster models. We use this semi-analytic model to investigate changes to the semi-major axis and eccentricity distributions for binaries evolved within the dynamical environment of a star cluster, to test if a thermal eccentricity distribution will emerge within a cluster due to dynamical encounters.

The thermal eccentricity distribution is predicted by theory for a population of binaries that has achieved energy equipartition, with energies following a Boltzmann distribution. It has an elegant form, is easy to implement in numerical codes, and is therefore extremely popular in theoretical investigations of populations of binaries.

However, nearly all observed binary populations show eccentricity distributions that are flatter than thermal and more closely consistent with a uniform distribution (Duchêne & Kraus, 2013; Raghavan et al., 2010; Moe & Di Stefano, 2017). Still, most empirical determinations of binary eccentricity are limited to the shortest-period binaries. Our analysis of solar-type binaries with both our semi-analytic model and more sophisticated $N$-body and Monte Carlo star cluster simulations, shows that it is difficult, if not impossible, for cluster dynamics to convert a uniform eccentricity distribution to a thermal distribution within a star cluster, even for the widest binaries. For most solar-type MS cluster binaries, the eccentricity distribution is maintained, with only minimal perturbations, throughout the cluster lifetime.

Moe & Di Stefano (2017) argue that for wide visual O5 - B5 binaries, the eccentricity distribution is consistent with thermal. These wide binaries may be analogous to the captured binaries from the Kouwenhoven et al. (2010) $N$-body simulations. If not, they likely formed with a thermal distribution at birth (simply because the timescales for encounters to thermalize the eccentricities is long compared to the ages of high-mass stars).

Indeed our models suggest that if a solar-type MS binary population is observed to have a thermal eccentricity distribution, this was likely imprinted upon birth, or perhaps shortly thereafter during some highly dynamic epoch in the cluster formation process (which is not included in our models).

The populations of binaries in numerical models that achieve a thermal distribution (without imposing the thermal distribution at birth) all formed dynamically. We see this in the BH and remnant binaries from the CMC model studied here (Figures 5 and 6), and also in the literature (Fregeau et al., 2004; Kouwenhoven et al., 2010; Perets & Kouwenhoven, 2012). Indeed, binaries that form through dynamical exchanges may be the best case where one can justifiably assume an initially thermal eccentricity distribution. Those that form through the “three-body binary formation” mechanism in Monte Carlo model are typically assumed to have thermal eccentricities (e.g. Morscher et al., 2015). Detailed scattering simulations of this mechanism are desirable to verify this assumption.

Theoretical investigations that choose to initialize all binaries with a thermal distribution will make incorrect predictions for the evolution of the binary population. For example the stellar merger rate may be overpredicted by a factor of about two (see Section VI), if an initially thermal eccentricity distribution is chosen when a uniform distribution is more appropriate.

In closing, we suggest careful consideration when choosing an initial distribution of eccentricities for binaries in star cluster simulations and population synthesis models. Though the thermal eccentricity distribution is a popular choice, and often the default, in most cases it may not be appropriate, and if so can lead to incorrect results.

We thank Douglas Heggie for his input and suggestions in constructing our semi-analytic model. This research was supported in part through the computational resources and staff contributions provided for the Quest high performance computing facility at Northwestern University which is jointly supported by the Office of the Provost, the Office for Research, and Northwestern University Information Technology. K.K. acknowledges support by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-1324585. K.K. and F.A.R. acknowledge support from NASA ATP Grant NNX14AP92G and NSF Grant AST-1716762. M.G. was partially supported by NCN, Poland, through the grant UMO-2016/23/B/ST9/0273.