EVOLUTION OF THE BINARY FRACTION IN DENSE STELLAR SYSTEMS

Observations and recent theory strongly suggest that the initial mass function (IMF) is universal among non-zero metallicity stars (e.g., Chabrier 2003). Indeed, Bate (2009b) suggested that radiative feedback may naturally regulate the star formation process so as to produce an IMF that is only weakly dependent on the properties of the progenitor molecular cloud. Naively, one would also expect that other features of the initial stellar population—like the binary fraction—should be nearly universal. Hydrodynamical star formation simulations yield companion star frequencies and binary fractions that are largely independent of the properties of the progenitor molecular cloud (although the statistics in some cases are marginal), and are quite consistent with observations (Bate 2009a; Bate et al. 2003; Bate & Bonnell 2005).

Observations of stars in low stellar density environments where dynamics is unimportant, such as the solar neighborhood, yield a binary fraction of ∼50% among solar-type stars, with an increasing trend with primary mass (e.g., Duquennoy & Mayor 1991; Fischer & Marcy 1992). Open clusters similarly show such large binary fractions (Fan et al. 1996). However, observations of dense globular cluster cores typically yield binary fractions that are significantly smaller. Hubble Space Telescope (HST) observations of the core-collapse cluster NGC 6397 yield a binary fraction of ≈5% in the core and ≈1% beyond the half-mass radius (Davis et al. 2008). For the canonical non core-collapse cluster 47 Tuc, the binary fraction is ≈13% (Albrow et al. 2001). The core binary fraction generally ranges from a few percent to tens of percent, approaching 50% in some cases for less dense clusters (Sollima et al. 2007). Where measured, the binary fraction outside the core is always smaller (see the table in Davis et al. 2008). The question naturally arises: are the currently observed relatively small core binary fractions in globular clusters consistent with initially larger binary fractions of ∼50%?

There are many strongly coupled processes that determine the evolution of the core binary fraction in a dense stellar system. Stellar evolutionary processes alone can affect the properties of a binary greatly, causing it to expand or shrink via mass transfer or winds, circularize via dissipative effects, lose mass, receive a systemic velocity kick due to a supernova, or disrupt or merge. The properties of the binary feed into the dynamical interaction rate with other stars or binaries, causing it to interact more or less frequently depending on its semimajor axis, eccentricity, mass, and systemic velocity. A strong dynamical interaction of a binary can disrupt it, exchange one of its members for an incoming star, cause its orbit to expand or shrink, modify its eccentricity, increase its systemic velocity via gravitational recoil, or cause two or more stars to physically collide. The dynamically modified binary properties feed back into binary stellar evolution, possibly initiating or halting mass transfer, or increasing tidal effects. In contrast to stellar evolutionary processes, the dynamical interaction rate depends on the cluster density and velocity dispersion, which evolve with time. Since binaries are typically more massive than single stars, mass segregation can increase their numbers in the core at the expense of single stars. The tidal effects of the host galaxy will preferentially strip single stars from the halo of the cluster.

For a globular cluster of typical mass (∼10⁵ M_☉) and size (half-mass radius r_h ∼ 3 pc), its global evolution can be divided into three phases according to the timescales of the relevant physical processes. At early times (∼few × 10 Myr), the evolution is largely driven by stellar evolutionary mass loss from the most massive stars in the cluster. At intermediate times (∼few Gyr), as mass loss from stellar evolution has slowed, the evolution is driven primarily by two-body relaxation. At late times (possibly beyond a Hubble time; Hurley 2007), when the core has reached sufficiently high density for binaries to strongly interact dynamically and release enough energy to prevent core collapse, the properties of the cluster are determined by the makeup of the binary population in this quasi-equilibrium "binary burning" phase.

The core binary fraction is clearly a quantity that is affected by nearly all physical processes operating in a cluster, and is of obvious observational interest. Comparing observed core binary fractions with simulation results (in combination with other observables) is thus a good measure of our theoretical understanding of cluster evolution. There can be dramatic differences in definition between the observed binary fraction and what theorists call the binary fraction, however.

When measured with the common offset main-sequence (MS) method, MS–MS binaries are detected by their appearance as distinctly brighter MS objects. The observed binary fraction is defined as the ratio of the number of these "binary sequence" objects to the total number of objects in the MS and the binary sequence, corrected for the assumed number of binaries with mass ratio so small they would blend in with the MS.

The theorists' definition of the binary fraction is typically the ratio of the number of binaries to the total number of "objects" (single stars or binaries). Furthermore, computational theorists tend to consider only "hard" binaries. That is, binaries with binding energy greater than the typical particle energy, which typically become more tightly bound (harden) as a result of encounters (Heggie & Hut 2003). Soft binaries—binaries with binding energy less than the typical particle energy in a cluster, which typically become less tightly bound (soften) or dissociate completely—are less frequently considered. We consider in detail the difference between the observational and theoretical definitions of the binary fraction below, as well as the importance of soft binaries.

Recently, two very different simulation methods have been used to study the evolution of the binary fraction. Ivanova et al. (2005) have developed a simplified Monte Carlo (MC) method in which a dense, massive cluster is modeled as a constant-density core plus halo (to simulate the long-lived binary burning phase that clusters may reach late in their evolution). Binaries and stars are evolved via the population synthesis code StarTrack (Belczynski et al. 2008), and the strong dynamical interactions of binaries are integrated numerically with Fewbody (Fregeau et al. 2004). Objects move between the core and the halo due to mass segregation and systemic velocity changes resulting from dynamical encounters. In this approach, the core mass increases slowly with time, with very few stars leaving the core after mass segregating into it.

Ivanova et al. (2005) found, generally, that the core binary fraction decreases significantly with time. Even for a modest core density of 10³ pc⁻³, they found that an initial binary fraction of 100% yields a core binary fraction of 27% at 14 Gyr. For the density of 47 Tuc, they found that a 100% initial binary fraction yields an 8% core binary fraction at 14 Gyr. It should be noted, however, that these figures include substantial numbers of soft binaries—binaries that are so wide they are quickly destroyed by dynamical encounters. If only the hard binaries in these simulations are counted, an initial binary fraction of 25% in a 10³ pc⁻³ core density cluster yields a 15% core binary fraction at 14 Gyr. For a density of 10⁵ pc⁻³ the core binary fraction evolves from an initial 25% to 7%.

Hurley et al. (2007) have used a direct N-body method, coupled with the BSE single and binary stellar evolution routines (Hurley et al. 2000, 2002), to study the evolution of the binary fraction. The great benefit of this method is that it makes no simplifying assumptions about the underlying cluster evolution. On the other hand, it is computationally expensive, currently limiting its application to open clusters or globulars with low initial binary fractions. Hurley et al. (2007) found that the core (hard) binary fraction generally increases with time. For a cluster of 5 × 10⁴ stars with a central density of ∼10^3.5 pc⁻³, the core binary fraction rises from an initial 20% to 52% at 9 Gyr. For lower initial densities, the degree of increase of the core binary fraction is similar.

On the face of it, the discrepancy between the two methods appears irreconcilable. However, the two methods operate at very different core densities and cluster masses, both of which affect the half-mass relaxation time and hence the mass segregation timescale, as well as the binary dynamical interaction rate. Figure 1 shows the evolution of the various models in core number density–binary fraction space. Note that the binary fractions plotted here include only hard binaries. Each model's evolution is represented as a simple arrow, with the tip at the final properties, and the tail at the initial properties. It is clear from this figure that the two methods represent very different regions of parameter space, and could simply be displaying different aspects of the same underlying physics. The only point of concern is the N-body model starting at ∼10^3.5 pc⁻³ and evolving toward a higher binary fraction, nestled between two MC models evolving in the opposite direction.

Zoom In Zoom Out Reset image size

To elucidate the evolution of the binary fraction, and to address the discrepancy between the existing N-body and simplified MC models, we have performed a grid of simulations with our newly upgraded MC cluster evolution code. Note that our MC code is very different from that of Ivanova et al. (2005). While their code assumes a constant core density with time and samples binary interactions using MC techniques, our code self-consistently models the global evolution of a cluster, using MC techniques to sample the stellar distribution function when applying the effects of two-body relaxation. The naming clash is the unfortunate consequence of the popularity and applicability of MC techniques in general.

Our MC code self-consistently models the evolution of star clusters due to the effects of two-body relaxation, evaporation through a Galactic tidal boundary, dynamical scattering interactions of binaries, physical stellar collisions, and now single and binary stellar evolution. The details of the method and its implementation are described in detail elsewhere (Joshi et al. 2000, 2001; Fregeau et al. 2003; Fregeau & Rasio 2007). Here we focus on the addition of stellar evolution.

For coding simplicity and for more directed comparisons with existing N-body simulations, we have incorporated the BSE single and binary stellar evolution routines in our MC code (Hurley et al. 2000, 2002). In our code stellar evolution is performed for each object (single star or binary) during a time step in step with dynamics. Since at early times a cluster can lose a lot of mass due to supernovae, we make sure to limit the time step so that no more than a small fraction of the total cluster mass is lost in one step (typically we set this fraction to 10⁻³).

To test that our inclusion of the stellar evolution routines is accurate, we have compared with the N-body results of Hurley (2007), who evolved N = 10⁵ cluster models with binary fractions ranging from 0% to 10%. The results are shown in Figure 2, which displays the evolution of the core to half-mass radius ratio (r_c/r_h) with time. The data from Hurley (2007) were extracted from that paper using ADS's Dexter Applet (Demleitner et al. 2001). For reference, we also plot the evolution of a model with stellar evolution turned off, and a model without an external tidal field. Since the model without stellar evolution reaches core collapse in under 1 Gyr, and since the model with no tide differs only minimally from the models with all physics turned on, it is clear that stellar evolution drives the evolution of the cluster. These models are thus a good test of our treatment of stellar evolution. The agreement with N-body is quite good given the vastly different methods, although the MC models tend to expand more at early times due to supernovae. The peculiar feature that the evolution does not appear to depend strongly on initial binary fraction is reproduced in our models. At late times (≳15 Gyr), our models begin to diverge with N-body. This is likely due to the fact that the clusters have lost roughly 70% of their stars by this time, and our apocenter-based treatment of the tide will tend to underestimate tidal mass loss as the number of cluster stars decreases, when an energy-based criterion is more appropriate (e.g., Giersz et al. 2008).

Zoom In Zoom Out Reset image size

A more detailed discussion of the implementation of stellar evolution in our MC code, and a suite of tests verifying the accuracy of implementation, will be presented in the next paper in the MC series (S. Chatterjee et al. 2010, in preparation). Given the vast differences between the N-body method and our MC method, we take the agreement between our models and those of Hurley (2007) in Figure 2 as a sign that our implementation of BSE in our code is at least consistent with that in N-body.

There is one aspect of our method that deserves special mention, however. It is generally believed that if a cluster avoids a collisional runaway phase (e.g., Freitag et al. 2006a, 2006b) the stellar-mass black holes (BHs) formed early in a cluster's lifetime will quickly sink to the core and dynamically decouple from the rest of the cluster, undergoing their own evolution, much like an independent small star cluster (Sigurdsson & Hernquist 1993). The BH subsystem will quickly dissolve through its own internal dynamics, ejecting all but one or two of the BHs on a timescale of ≲1 Gyr. Aside from removing nearly all BHs from the cluster, the result is a mild energy injection into the cluster core, causing it to expand somewhat at early times (Mackey et al. 2007). A typical star cluster of N = 10⁶ objects will contain a subsystem of up to ≈10⁻³N = 10³ BHs evolving independently in the core (O'Leary et al. 2006). For the N = 10⁵ clusters considered in this work, the number is ≈100. The MC method is not designed to handle subsystems of less than a few hundred objects, since they are often far from spherically symmetric, and large angle scattering dominates. (Note that we treat small-N encounters up to N = 4 via direct integration.) We therefore truncate the mass function at 18.5 M_☉ (the largest progenitor mass not resulting in a BH) for the runs presented here. The resulting discrepancy in r_c/r_h is important only at early times and as Figure 2 shows is minimal.

We have performed several simulations of evolving clusters for a grid in initial binary fraction and initial cluster virial radius (or equivalently, central density). All our simulations start with N = 10⁵ objects initially (an object being either a binary or a single star), and like the simulations of Hurley et al. (2007), assume a Plummer density profile with no primordial mass segregation, a "standard" Galactic tide (cluster at 8.5 kpc from Galactic center, 10¹¹ M_☉ Galactic mass enclosed), a Kroupa et al. (1993) IMF, and only hard binaries. Our IMF extends from 0.15 to 18.5 M_☉, binary secondary masses are drawn from a distribution flat in the mass ratio, the semimajor axis a is drawn from a distribution flat in log a from a minimum of a_min = 5(R₁ + R₂), where R_i are the individual stellar radii, to a maximum corresponding to an orbital velocity of the lighter member equal to the local velocity dispersion, and the eccentricity is drawn from a thermal distribution truncated at the value corresponding to contact at a_min. Note that our large a cutoff for wide binaries is equivalent to the hard–soft boundary for equal-mass stars (Fregeau et al. 2006).

Figure 3 shows the evolution of our models in core number density–binary fraction space. It is clear that for all the but the highest initial binary fraction cases, the core binary fraction increases with time. The observational data points seem to be consistent only with cluster models that started with relatively low central densities (∼10³ pc⁻³) and small hard binary fractions (∼5%). As we discuss in the following section, the core binary fraction is typically estimated observationally by measuring the fraction of MS stars belonging to the binary MS, and convolving it with an assumed binary mass ratio distribution. It is not a priori evident that this MS binary fraction reflects the underlying true binary fraction.

Zoom In Zoom Out Reset image size

Why does the core binary fraction generally tend to increase with time? As mentioned above, there are many strongly coupled processes that affect the core binary fraction. However, the general trend can be understood approximately as an interaction between mass segregation of binaries into the core, and the destruction of binaries preferentially in the core.

Figure 4 shows the evolution of our MC cluster evolution model starting with f_b = 0.05 and n_c ≈ 10^2.5 pc⁻³. As the evolution of r_c/r_h in the top panel shows, the cluster core contracts steadily until it enters a phase of binary burning at the relatively late time of ∼12 Gyr. The bottom panel shows the evolution of the half-mass radius of single stars, r_h,s (solid line), and the half-mass radius of binaries, r_h,b (dashed line), relative to the overall cluster half-mass radius, r_h. The differential mass segregation between the single and binary populations is evident, with the single stars expanding slightly relative to the bulk of the cluster, and the binaries contracting significantly. The quantity r_h,b/r_h decreases steadily until ∼11 Gyr due to mass segregation. It then begins to increase due to destruction of binaries preferentially in the cluster core by strong dynamical interactions and perturbed stellar evolution (see, e.g., Ivanova et al. 2005, for a discussion of perturbed binary evolution).

Zoom In Zoom Out Reset image size

Figure 5 shows the evolution of the number of single stars in the core, N_c,s (solid line), and number of binaries in the core, N_c,b (dashed line) for the same model. The quantity N_c,s declines steadily due to standard gravothermal evolution, in which the cluster core becomes denser and smaller in number with time (e.g., Binney & Tremaine 2008). The quantity N_c,b, on the other hand, is roughly steady until ∼11 Gyr due to mass segregation of binaries into the cluster core.

Zoom In Zoom Out Reset image size

As suggested by Figure 5, the core mass decreases with time, as expected from standard gravothermal evolution. This is in contrast to the simplified MC method of Ivanova et al. (2005), in which the core mass steadily increases with time, due primarily to mass segregation of binaries into a core of fixed density.

We note also that mass segregation of a binary into the core implies, by energy conservation, a preferential expansion of lighter single stars in the vicinity of the binary. (Energy conservation is roughly applicable because the mass segregation timescale is shorter than the local relaxation timescale, by a factor of M/m, where M is the mass of the segregating object and m is the mass of a background star.) This effect is not included in the code of Ivanova et al. (2005), and is likely an important factor in the discrepancy between their results and ours.

Another important factor, as suggested by Figures 2 and 4, is that the long-lived, high-density binary burning phase assumed by Ivanova et al. (2005) may not be generic for globular clusters. Instead, the "core contraction" phase may last a Hubble time, and the cluster cores we observe now may have been much less dense in the past (Fregeau 2008). Although the central density in our models increases steadily with time, the local density at the half-mass radius decreases with time, resulting in final half-mass relaxation times that are a factor of ∼3 longer than their initial values. In the cases where our models do enter the binary burning phase before a Hubble time, we find that the core binary fraction in this phase steadily decreases with time. This behavior is consistent with the results of Ivanova et al. (2005).

While the core binary fraction in the majority of our models increases with time, the overall cluster binary fraction remains roughly constant with time. This is in good agreement with the Hurley et al. (2008) N-body models inspired by NGC 6397, and supports their use of the currently observed binary fraction near the half-mass radius as a measure of the primordial binary fraction (although the validity of comparison with NGC 6397 is not obvious, since the N-body models end with a factor of 5–10 fewer stars than NGC 6397 currently contains). For the low-density f_b = 0.05 model just described, 39% of the initial binary population remains at 14 Gyr, 43% escape the cluster due to the tidal field (compared with the 60% of single stars that escape in the same fashion), 9% are destroyed via strong dynamical interactions of binaries, and 8% are destroyed via binary stellar evolutionary processes (possibly perturbed by dynamics). In other words, in this case the overall binary fraction remains roughly constant with time due to a balance between preferential tidal stripping of single stars in the outskirts and preferential destruction of binaries in the cluster core.

When using the offset MS method, what observers measure is in fact the number of MS–MS binaries with mass ratios q ≳ 0.5 relative to the total number of objects appearing in the MS (which may include apparent single MS stars, comprised of a MS star plus dim compact object companion). This fraction is then corrected to account for the low mass ratio MS–MS binaries that blend into the single MS, by adopting an assumed mass ratio distribution. This final corrected figure is what is usually quoted as the "observed binary fraction." However, there is no a priori reason to believe this quantity reflects the underlying binary fraction among stars of all types. Hurley et al. (2008) showed that, for the low binary fraction cluster models they considered (f_b ≲ 10%), the observed binary fraction is a good measure of the true binary fraction in the outer regions of a cluster, but can be a serious overestimate in the core.

Since the general nature of the relationship between the observed and true binary fraction is not obvious, we have plotted in Figure 6 the evolution of our models in core number density–observed core binary fraction space. The observed binary fraction is calculated as N_MS-MS/(N_MS-MS + N_MS + N_MS-CO), where N_MS-MS is the number of MS–MS binaries of any mass ratio in the core, N_MS is the number of single MS stars in the core, and N_MS-CO is the number of MS–compact object binaries in the core that appear near the MS. We count a MS–compact object binary as near the MS if the total luminosity of the binary is less than 10% more than that of the MS star (this corresponds to a magnitude increase of 0.1), and if the luminosity-weighted temperature of the binary is less than 10% different from that of the MS star. Like the true binary fraction plotted in Figure 3, small initial binary fraction models evolve toward larger binary fractions. However, large initial binary fraction models evolve toward drastically smaller observed binary fractions. As a relatively extreme example, the model with an initial binary fraction of 75% and initial central density of ∼10^2.5 pc⁻³ has an observed core binary fraction of just 33% at the end of the simulation. The true core binary fraction at the end of the simulation is 91%. Of the core binaries, 23% are MS–MS binaries, 32% are compact object–compact object binaries, and 44% are MS–compact object binaries (see Table 1 for more details). As expected, the discrepancy between the observed and the true binary fraction is due to compact object–compact object binaries not being counted in the observed tally, and MS–compact objects masquerading as single stars on the MS.

Zoom In Zoom Out Reset image size

An initial population of binaries that contains a substantial soft component can be a significant cluster energy sink, since the soft binaries are destroyed in dynamical scattering interactions. The result is that the core of a cluster born with many soft binaries will quickly contract as those binaries are ionized. Could soft binaries increase the concentration of a cluster so much that it would become core collapsed?

The total energy in soft binaries, for a distribution flat in the log of the semimajor axis, is simply

$$\begin{equation} E_{\rm b,s} = \frac{N_b (E_{\rm b,hs} - E_{\rm b,amax})}{\ln (a_{\rm max}/a_{\rm min})} \approx \frac{N_b E_{\rm b,hs}}{\ln (a_{\rm max}/a_{\rm min})} \,, \end{equation} \tag{ 1 }$$

where a_min and a_max are the limits on the semimajor axis distribution, E_b,hs is the energy of a binary at the hard–soft boundary, E_b,amax is the energy of the least-bound binary, N_b is the total number of binaries, and we have assumed that a_max ≫ a_hs. Note that this can also be written in the physically more insightful form:

$$\begin{equation} E_{\rm b,s} \approx \frac{N_{b,s} E_{\rm b,hs}}{\ln (a_{\rm max}/a_{\rm hs})} \,, \end{equation} \tag{ 2 }$$

where N_b,s is the number of soft binaries, and a_hs is the semimajor axis at the hard–soft boundary. Assuming for simplicity a cluster of equal-mass objects (binaries in this case) of mass m_ave with mean one-dimensional velocity dispersion σ, this becomes

$$\begin{equation} E_{\rm b,s} \approx \frac{\frac{3}{2} N_b m_{\rm ave} \sigma ^2}{\ln (a_{\rm max}/a_{\rm min})} \,. \end{equation} \tag{ 3 }$$

From the virial theorem, the total mechanical energy of a cluster is simply $E_{\rm clus} = -\frac{3}{2}Nm_{\rm ave}\sigma ^2$, where N is the number of cluster objects, so

$$\begin{equation} \frac{E_{\rm b,s}}{|E_{\rm clus}|} \approx \frac{N_b}{N \ln (a_{\rm max}/a_{\rm min})} \,. \end{equation} \tag{ 4 }$$

For an admittedly optimistic binary fraction of 1 (N_b = N), and realistic binary semimajor axis limits of a_min = 5 × 10⁻² AU (corresponding to a contact binary during the pre-MS phase) and a_max = 10³ AU (corresponding to a 10⁷ day orbital period), the energy in soft binaries is ≈10% of the total cluster mechanical energy.

The question, of course, is if this amount of energy is sufficient to make a cluster concentrated enough to appear to be core collapsed. For our working definition of core collapse we assume that a cluster core can be resolved with HST if its radius is at least 1 arcsec in size. At a typical cluster distance of 10 kpc, this corresponds to ∼0.05 pc. Starting with a King model of a given mass, binary fraction, central concentration W₀, and half-mass radius r_h, we calculate the total mechanical energy of the cluster within the half-mass radius, E_h (Binney & Tremaine 2008). We then calculate the energy of the soft binaries, E_b,s. This energy will be absorbed from the cluster when those binaries are destroyed in dynamical interactions in and around the cluster core. Keeping r_h fixed (since the timescale for destruction of soft binaries is shorter than the half-mass relaxation time), we then calculate a new King model with half-mass energy E'_h = E_h − E_b,s (note that E_b,s > 0 by construction). For the new King model we calculate the new central velocity dispersion and hence the new hard–soft boundary (which has moved to a smaller binary semimajor axis), calculate the energy available in the newly soft binaries, and iterate until the solution converges. For a 5 × 10⁵ M_☉ cluster with half-mass radius r_h = 5 pc and initial core radius r_c = 1.9 pc (W₀ = 6, concentration c = log₁₀(r_t/r_c) = 1.25), an initial binary fraction of 100% with the semimajor axis distributed flat in log a from 5 × 10⁻² AU to 10³ AU is sufficient to drive the cluster to a W₀ = 10, c = 2.3 King model with core radius r_c = 0.16 pc. (A W₀ = 10 King model has maximal binding energy within r_h for fixed r_h and mass.) This is quite close to core collapsed, and may even be classified as such if viewed with a ground-based telescope. In fact, Wiyanto et al. (1985) showed that clusters enter the self-similar stage of evolution (the "onset" of core collapse) when W₀ > 7.4, so such a model would reach core collapse quickly. We have repeated this calculation with a binary distribution that is lognormal in orbital period, as in Duquennoy & Mayor (1991) or Fischer & Marcy (1992), with 〈log₁₀P_d〉 = 4.8 and $\sigma _{\log _{10} P_d}=2.3$, where P_d is the period in days, and with the same limits on the semimajor axis as above. The results are unchanged with this binary distribution, largely because its peak lies at wider orbits than the hard–soft boundary for globular clusters.

To test this scenario numerically, we have run models with a binary distribution extending well beyond the hard–soft boundary, to P = 10⁷ days. Our "high density" model (cluster mass 9 × 10⁴ M_☉, standard wide mass spectrum, r_h = 1.0 pc initially) with f_b = 0.9 (including soft binaries), evolves from f_b,c = 0.9 and r_c = 0.6 pc to f_b,c = 0.4 and r_c = 0.05 pc in just 3 Myr (see Figure 7). This is in striking agreement with the energy argument above, which predicts rapid evolution to r_c = 0.1 pc for this model. Note that the energy argument assumes all soft binaries will be destroyed on a short timescale. To achieve this in practice requires efficient mass segregation of binaries into the core, which has been aided in this case by a wide mass spectrum, at the expense of inaccuracy in calculating E_b,s. After the rapid initial contraction of the core, the cluster quickly (after a few Myr) enters into a long-lived binary burning phase.

Zoom In Zoom Out Reset image size

From the preceding discussion, it is evident that the dynamical importance of soft binaries should not be ignored. If a cluster is born with significant numbers of soft binaries, its evolution may be vastly different from a similar cluster containing only hard binaries. First, the rate of binary destruction is greatly enhanced in clusters containing soft binaries, yielding a starkly decreasing binary fraction with time. Second, the binary burning phase is reached quickly (within a few Myr) due to soft binary destruction. When only hard binaries are present, the binary burning phase may not be reached within a Hubble time, as shown, for example, by Figure 2. The implications for our understanding of the current dynamical states of Galactic globular clusters are profound, as certain properties of clusters can be explained by the majority of clusters currently being in the initial "core contraction" phase, and not yet in binary burning (Fregeau 2008). We have provided here just a cursory analysis of the effects of soft binaries. A more detailed study should certainly be undertaken in the future.

Independent of the details, it seems clear that the hard binary fraction in the core of a dense stellar system will generally increase with time (with the exception of an initial hard binary fraction ≳90%). Yet there is no compelling evidence that clusters should be born with binary fractions smaller than the typical field value of ∼50%, and observations yield core binary fractions of just ∼10% in Galactic globular clusters. If the observations are to be taken at face value, how then can they be consistent with large initial binary fractions? One possibility, as pointed out by Davis et al. (2008), is that the binary fraction is a strong function of primary mass (Lada 2006), with the single star fraction increasing to ∼75% for M dwarfs and lighter stars. A Kroupa et al. (1993) IMF with a 25% binary fraction from 0.1 to 0.5 M_☉ and a 50% binary fraction from 0.5 to 100 M_☉ yields an overall binary fraction of just 32%.

Another possibility is that most binaries born in clusters are soft relative the cluster velocity dispersion, in which case they will be destroyed very quickly by dynamics. If the binary period distribution is uniform in log P from 0.1 to 10⁷ days as in Ivanova et al. (2005), the 32% overall binary fraction just suggested corresponds to a hard binary fraction of merely ∼10% for a cluster with central density 10⁶ pc⁻³. As demonstrated above, the early, rapid destruction of soft binaries may lead to a binary burning phase within a short time (≲5 Myr, depending on initial conditions).

Aside from the initial binary properties, could it be that observations are under-counting the binary fraction significantly? When we measure the binary fraction using an offset MC method similar to what observers use, we find that clusters with large initial binary fractions (f_b ≳ 0.5) evolve toward smaller observed core binary fractions (f_b ≲ 0.5). The discrepancy between the observed and true core binary fractions is caused by compact object–compact object binaries not being counted in the sample, and MS–compact object binaries masquerading as single stars.

Could a binary be sufficiently wide to be resolved as two single stars and hence missed as a binary? For the wide-field camera on HST, one requires two turnoff mass stars in a binary to be separated by roughly 4 pixels for the binary to be resolved. For a cluster at a distance of 10 kpc, this corresponds to a binary separation of ∼4 × 10³ AU. For a cluster with a velocity dispersion of 10 km s⁻¹, this corresponds to a binary hardness of Gm/av²_σ ≈ 2 × 10⁻³, which is too soft to survive dynamically for even a short time.

We have briefly described our inclusion of the BSE single and binary stellar evolution routines in our MC globular cluster evolution code (Hurley et al. 2000, 2002). We have compared with the results of direct N-body simulations and found good agreement, suggesting that our implementation of BSE in our code is consistent with that in N-body.

We have used our newly upgrade MC code to study the evolution of the core hard binary fraction in star clusters, and in particular attempt to settle the apparent disagreement between direct N-body and simplified MC techniques on its evolution. We find that the core binary fraction generally increases with time, even for low initial core density models (n_c ≈ 10^2.5 pc⁻³), with only very small initial binary fraction models (f_b ≲ 0.05) producing the core binary fractions of ∼10% observed today. The increase in the core binary fraction with time can be understood as an imbalance between mass segregation of binaries into the core (and single stars out of the core) and the destruction of binaries in the core directly via strong dynamical encounters, and indirectly via dynamical perturbation of binary stellar evolution processes. The overall cluster binary fraction remains roughly steady with time, due to the additional effect of preferential tidal stripping of single stars from the cluster outskirts.

This evolution, however, refers to the true binary fraction. When measuring the core binary fraction using an offset MC method analogous to what observers use, we find that the observed core binary fraction can seriously underestimate the true core binary fraction. This results from compact object–compact object binaries not being counted and MS–compact object binaries masquerading as single stars in the observed tally. In the course of creating more detailed models of M 67, 47 Tuc, M 4, and NGC 6397 to be compared with observations, we are now developing a data reduction pipeline that includes simulations of spectra for every star. Among our near future plans is the creation of a cluster sky map in different bands, to which we can apply the MS binary detection method to determine more accurately how many binaries are missed by the method.

Most of our discussion concerned hard binaries. However, we also considered the effects of a substantial population of soft binaries. We found that the energy sink represented by soft binaries (for a typical binary period distribution) is sufficient to cause the core of a typical globular cluster to contract significantly. The result is not only a rapid, efficient destruction of a significant number of binaries at early times, but also a much earlier onset of the binary burning phase, resulting in enhanced binary destruction in the core with time.

The authors thank C. Heinke and J. Hurley for data and helpful discussions, and the anonymous referee for constructive comments that improved the draft and for pointing out the insightful form in Equation (2). J.M.F. acknowledges support from Chandra/Einstein Postdoctoral Fellowship Award PF7-80047. F.A.R. acknowledges support from NASA Grant NNG06GI62G at Northwestern University. This research was completed at KITP while the authors participated in the spring 2009 program on "Formation and Evolution of Globular Clusters," and was supported in part by the NSF under Grant PHY05-51164.