Star Clusters with Primordial Binaries:I. Dynamical Evolution of Isolated Models
Abstract
In order to interpret the results of complex realistic star cluster simulations, which rely on many simplifying approximations and assumptions, it is essential to study the behavior of even more idealized models, which can highlight the essential physical effects and are amenable to more exact methods. With this aim, we present the results of body calculations of the evolution of equalmass models, starting with primordial binary fractions of 0  100 %, with values of ranging from 256 to 16384. This allows us to extrapolate the main features of the evolution to systems comparable in particle number with globular clusters.
In this range, we find that the steadystate ‘deuterium main sequence’ is characterized by a ratio of the core radius to halfmass radius that follows qualitatively the analytical estimate by Vesperini & Chernoff (1994), although the dependence is steeper than expected. Interestingly, for an initial binary fraction greater than , the binary heating in the core during the post collapse phase almost saturates (becoming nearly independent of ), and so little variation in the structural properties is observed. Thus, although we observe a significantly lower binary abundance in the core with respect to the FokkerPlanck simulations by Gao et al. (1991), this is of little dynamical consequence.
At variance with the study of Gao et al. (1991), we see no sign of gravothermal oscillations before halfmass relaxation times. At later times, however, oscillations become prominent. We demonstrate the gravothermal nature of these oscillations.
keywords:
globular clusters – methods: Nbody simulations – stellar dynamics.1 Introduction
The evidence for primordial binary stars in old globular star clusters continues to accumulate steadily. Even though no comprehensive review has appeared since Hut et al. (1992), recent additions include Pulone et al. (2003), Bellazzini et al. (2002), Albrow et al. (2001).
The body modeling of such clusters has also developed rapidly, with studies devoted to exotic species associated with binaries, blue stragglers, intermediate mass black holes, hierarchical systems, and so on (see, for example, Hurley et al., 2001; Portegies Zwart & McMillan, 2002; Aarseth, 2001; Portegies Zwart et al., 2004; Hurley et al., 2005).
There have been relatively few body studies devoted to the effects of binaries on the structural evolution of model clusters. McMillan et al. (1990), McMillan et al. (1991), McMillan & Hut (1994), Heggie & Aarseth (1992), Aarseth & Heggie (1993), Kroupa (1995) and de La Fuente Marcos (1996) described a wide range of results with values of which are generally rather modest by modern standards (though not inappropriate for application to many open clusters). More recent studies include that of Wilkinson et al. (2003), whose focus was on the evolution of the core radius with age. Generally speaking, however, in recent years the effect of primordial binaries on the structural evolution of rich systems has been studied using approximate methods based on the FokkerPlanck equation or related approaches. The study by Gao et al. (1991) [hereafter Gao et al.], using the FokkerPlanck model with finite differences, established an interesting set of parameters which has been taken up by Giersz & Spurzem (2000, 2003), using their hybrid gas/Monte Carlo code, and by Fregeau et al. (2003, 2005), whose code was pure Monte Carlo (with computed binary interactions in the later paper).
For all the usual reasons, it is advisable to check the results of approximate methods by means of body simulations, and it is the main purpose of the present paper to do so, based on the standard model defined by Gao et al. We present the results of an extensive set of direct body simulations, performed with up to particles and different ratios for the number of initial binaries to single stars, ranging from to . Following an outline of the physical picture in Sec. 2, the specification of our simulations is given in Sec. 3.
The subsequent two sections contain the main presentation of our results, with emphasis on (i) a comparison with the results of Gao et al. (Sec. 4), and (ii) an analysis in terms of the dependence on and on the primordial binary fraction (Sec. 5). In fact these topics cannot be cleanly separated. For example, the main systematic differences observed with respect to Gao et al. are in the fraction of binaries in the core, and possibly in the size of the core, in the postcollapse phase. After extrapolating our results to the number of particles used by Gao et al., we clearly observe a lower fraction of binaries in the core (by around a factor ). Extrapolation of the core radius is more problematic, however (see Sec. 5.1.2). In this section we also compare our results with those of other approximate numerical methods, and with the analytical estimate of the core radius due to Vesperini & Chernoff (1994). The last section of the paper provides an extended summary.
2 The Physical Picture
For an isolated cluster, the most important characteristics describing the mass distribution are the core radius and the halfmass radius . Fig. 1 shows how both evolve for the case of a Plummer model distribution of single stars, all of which have equal mass, with no primordial binaries present. The halfmass radius remains roughly constant during the core collapse phase, while the core shrinks at an increasing rate until the first core bounce takes place. Around this time, one or more binaries are formed dynamically, through simultaneous threebody encounters. The energy production in those binaries halts core collapse, and also initiates an expansion of the halfmass radius. In the postcollapse phase, the core expands and contracts dramatically, while new binaries are formed and older binaries are first hardened and then ejected.
Fig. 2 shows the striking difference that primordial binaries make to this simple picture. The presence of primordial binaries causes the core to shrink roughly twice as fast as in the singlestar case. This effect is largely due to mass segregation: the binaries, between twice as massive as single stars, tend to precipitate toward the center. In contrast, the core contraction is much less deep in the case of primordial binaries: the core contraction levels off after the core shrinks by less than a factor ten, whereas in the single star case core collapse continues till the core is one hundredth of its original size. Another difference is that in the primordial binary case, the halfmass radius starts increasing very early on, a factor four earlier than in the single star case. Since there is very little mass loss in these early stages of cluster evolution, an increase in halfmass radius indicates an injection of energy into the bulk motion of the stars and the center of masses of the binaries. This points to an efficient energy mechanism at work, which then also explains the early halting of the core contraction.
This is indeed the physical explanation. Binaries can give off energy in encounters between binaries and single stars^{1}^{1}1Though such encounters can result in the ejection of single stars and binaries, the cluster still gains energy, by the process known as “indirect heating”. (and also in binarybinary encounters), the probability of which is given by the product of the densities of binaries and single stars, whereas the formation rate of new binaries is proportional to . The relevant scales for these densities are set by the ninety degree turnaround distance , which is of order in standard units, far smaller than typical interparticle distances, even in the core. If we define the interaction densities as the number of particles per unit of volume , and denote them by a hat, we have and . In a model with primordial binaries, and are comparable, at least in order of magnitude, and it follows that . In other words, at a given density the rate of energy generation by primordial binaries far exceeds the rate of energy generation by binaries that first have to be formed at that density (Hut, 1989).
Another way of describing the difference between these two forms of evolution is that Fig. 1 shows a collapse to the ‘hydrogen main sequence’, whereas Fig. 2 shows a collapse to the ‘deuterium main sequence’, to borrow metaphors from stellar evolution. Let us take a gas cloud of a solar mass, consisting of only hydrogen and helium. When this protostar contracts sufficiently, nuclear reactions in the core will at some point balance the radiation losses at the photosphere. At this point, by definition the star has landed on the (hydrogen) main sequence. If the star would be made of pure helium, it would shrink further, until it would reach the helium main sequence, where the central temperature and pressure are high enough to generate the same energy at the rate at which it is lost at the surface. In practice, however, protostars never contain pure hydrogen, but typical have a very small admixture of deuterium. This deuterium can burn to form helium at lower densities and temperatures than hydrogen can, for a similar reason that primordial binaries can ‘burn’ gravitationally at lower densities than single stars. As a result, a contracting protostar first gets hung up at the deuterium main sequence, where it has a radius significantly larger than it will have soon afterwards, when its central deuterium reservoir is burned up, and the star settles on the (hydrogen) main sequence.
Here we witness something similar: our contracting cluster core settles onto an almost stable extended core configuration, while burning the primordial binaries. Gradually, when this primordial fuel is being exhausted, the core shrinks since the burning process becomes less efficient when the binary fraction becomes less. This slow shrinking of the core continues until (almost) all primordial binaries are burned up, after which the new binaries have to be formed by the process of simultaneous threebody encounters between single stars. Note that the value in Fig. 2 toward the end is dropping toward the values seen in Fig. 1 from the moment of the first core collapse. This behavior is more evident in Fig.5.
3 Simulations: Setup and Analysis
3.1 Initial conditions and definitions
The models shown in Figs. 1 and 2, like all those studied in this paper, are isolated, with stars of equal mass. The initial distribution is a Plummer model. In the standard model of Gao et al., initially 10% of the objects are binaries, as in Fig. 2, with energies in the range –, where is the central temperature of the system. For a Plummer model , where is the mean kinetic energy per particle of the system (the binaries being replaced by their barycenter). Gao et al. used FokkerPlanck models with parameters corresponding to particles to study the evolution. They also carried out runs with the initial fraction of binaries equal to and .
In this paper we report on the result of body calculations, which we performed using Aarseth’s NBODY6 code (Aarseth, 2003), which has been slightly modified to provide additional diagnostics on the spatial distributions of single and binary stars (see, e.g., Fig. 10). Starting from a minimum of 256 particles, we have made runs with no primordial binaries, 2% primordial binaries, the three Gao et al. choices of 5%, 10%, and 20% primordial binaries, as well as runs with 50% and 100% primordial binaries (Table 1). As in the case of Gao et al., we have defined here the primordial binary fraction as the fraction of objects that are binaries:
(1) 
where and are the initial number of binaries and single stars, respectively. This implies that the fraction of the total mass in binary stars, in the case of equal masses, is larger in the following way:
For example, for a run with 10% primordial binaries, whereas .
Note that denotes the number of original objects, i.e. . When we discuss a run with and 50% primordial binaries, we are dealing with stars.
All our results are presented using standard units (Heggie & Mathieu, 1986) in which
where is the gravitational constant, the total mass, and the total energy of the system of bound objects. In other words, does not include the internal binding energy of the binaries, only the kinetic energy of their centerofmass motion and the potential energy contribution where each binary is considered to be a point mass. The corresponding unit of time is . For the relaxation time, we use the following expression
A quantity of great interest is the core radius, which we here define as in NBODY6, i.e. essentially
where is the distance of the th star from the density centre, and the density around each particle is computed from the distance to the fifth nearest neighbour (Casertano & Hut, 1985). In practice outlying stars are omitted from the sum.
The runs were performed using idle time on many workstations through the Condor system^{2}^{2}2http://www.cs.wisc.edu/condor/: this allowed us to obtain in a few months what would have taken a few years on a single dedicated workstation.
:  0%  2%  5%  10%  20%  50%  100% 

256  43  2  2  46  39  16  2 
512  49  2  3  46  45  24  1 
1024  51  3  4  45  44  46  1 
2048  47  3  2  36  31  21  2 
4096  9  3  2  13  8  8  2 
8192  1  1  1  3  1  
16384  1  1 
Note: This table gives the number of realisations for each value of and .
3.2 Fluctuations
One significant difference between the body models we are using and the models of Gao et al is the presence, in body models, of fluctuations. These have both a cosmetic and a scientific importance, and we consider these issues together here.
In presenting the results from body simulations, there is always the question of how much smoothing to apply. In Fig. 3, we present the results for one run, for the change in time of the core radius, for a variety of different amounts of smoothing. In some figures in this paper we apply smoothing in order to highlight the particular points made in each figure. Fig. 3 can be used as a gauge by which to judge the amount of ‘noise’ present in various figures. Generally we have used a smoothing interval of .
These fluctuations also have a scientific interest. In the context of threebody binaries (i.e. those formed in threebody encounters, without a population of primordial binaries) their role was considered in simplified models by Inagaki & Hut (1988) and Takahashi & Inagaki (1991). The role of fluctuations in an body model, though in the absence of binaries, was considered by Bettwieser & Sugimoto (1985). They showed that values of the velocity dispersion, sufficiently well defined to determine its spatial gradient, could be determined in a 1000body system, if averaged over a shell containing at least 90 particles and over at least one crossing time. With our present models we have an opportunity to consider similar issues for the generation of energy by primordial binaries.
We restricted attention to our largest simulation () and considered the change in internal energy of all the bound binaries. Within one body time unit this is usually zero, but may be negative (if a binary hardens, for example) or positive (if a binary escapes). Perhaps surprisingly, it is necessary to smooth over a long time interval (at least ) in order for the result to have a consistent sign (Fig.4).
We have also considered how these fluctuations affect the rate of interactions, in the following way. If the proportion of binaries in the core is very high, the rate of binarybinary interactions is , where is the central density. Since it follows that . Therefore a simplified model, such as a FokkerPlanck model, in which binaries interact at the correct average rate, will yield a core with parameters such that body model, however, will lead to values of the mean core radius and central velocity dispersion such that . Fluctuations in an
(2) 
and so one need not expect these mean values to agree with those of a FokkerPlanck model, even if the rate of energy generation by binaries in the FokkerPlanck model is correct. However, we have checked that the difference in eq.(2), with smoothing as in Fig.4, is very small. In a core supported by primordial binaries, the core contains many particles, and so fluctuations are unimportant by comparison, say, with a core supported only by threebody binaries.
4 Comparison with Gao et al.
The main purpose of this section is to present results of body simulations which may be compared rather directly with those in Gao et al. Therefore the sequence of presentation closely follows the figures in that paper. The main difference between our runs and those of Gao et al. is in the particle number. While Gao et al. used parameters appropriate to , even our largest run is considerably smaller. Our discussion of the dependence of our results is mainly given later, in Sec. 5. The runs discussed by Gao et al. have also been taken as a test case by Giersz & Spurzem (2000) and by Fregeau et al. (2003, 2005), and we also make occasional brief comparisons with their results.
4.1 Core and HalfMass Radius
The structural evolution is best captured by the halfmass and core radii. Fig. 2 shows the results for a run with and 10% primordial binaries, and can be compared with Fig. 1 in Gao et al., except that they scale radii by the initial scale radius of the Plummer model (i.e. about 0.59 in our body units). The expansion of is similar (about 0.7 dex in both calculations), and the main differences are in the evolution of the core radius. These are

in the body model the initial decrease of is smaller, and remains considerably larger than in the FokkerPlanck model; and

in the body models (even for ), there is no sudden onset of deep core collapses, as observed in the FokkerPlanck model at about .
We now discuss these two issues in turn.
The fact that the contraction of the core is less deep than in the model of Gao et al. could partly be an effect of the dependence which we observe in our body models (Fig. 18). Our conclusion, however, depends on how the extrapolation to (the value used by Gao et al.) is carried out. If we adopted the theory of Vesperini & Chernoff 1994 (see Equation 3) we would expect , which still exceeds by about 60% the value of about 0.025 found by Gao et al. The data from our simulations, however, suggest a steeper decrease of the core size with , and so the number given by Gao et al. could well be consistent with our direct simulations.
The difference in values of could be caused by other factors, including the relative abundance of single and binary stars (Fig. 11), which is also dependent (Fig. 21). A further factor to be borne in mind is that the core radius should depend on the efficiency at which binaries produce energy in three and fourbody collisions, and Gao et al. were well aware of the uncertainties in the analytical cross sections which they adopted. Finally, the definition of core radius in the two types of simulation may be different, though the initial values are quite similar.
Comparison with other results in the literature is beset by similar uncertainties. We report simply that Fregeau et al. (2003) found an initial contraction of the core fairly comparable to our result (despite the different values of ), but then Fregeau et al. (2005) obtained a result comparable with that of Gao et al. (1991) when binary interactions were computed directly. In the hybrid model of Giersz & Spurzem (2000), on the other hand, the initial decrease of the core radius was only about 0.2 dex (Table 2).
Now we consider the onset of gravothermal oscillations in our simulations. While gravothermal oscillations would not be expected to occur in models as small as that shown in Fig. 2 (Goodman, 1987), for deep gravothermal oscillations are observed in singlecomponent body systems (cf. Makino 1996). Indeed we observe gravothermal oscillations in our run with . However their onset occurs later than in the simulations of Gao et al. (see Figs. 5, 6), though again we do not know how the time of onset should vary with ; this difference may or may not be significant.
Further evidence, though it is contradictory, is provided by the results of other authors (Table 2). By comparison of their two papers, Giersz & Spurzem (2003) concluded that the time of onset of gravothermal oscillations is greatly affected by the treatment of binary interactions. The treatment of Fregeau et al. (2003) more closely resembles that used by Giersz & Spurzem in their earlier paper, but they found (Fregeau et al. (2005)) that the onset of gravothermal oscillations is delayed when binary interactions are computed directly.
Source  

This paper  
Giersz & Spurzem (2000)  
Giersz & Spurzem (2003)  
Gao et al. (1991)  
Fregeau et al. (2003)  
Fregeau et al. (2005) 
Notes: is the change in the logarithm of from until the end of the initial contraction onto the “deuterium main sequence”; is the time (in units of ) of onset of gravothermal oscillations. Most of these values have been estimated quite approximately from the figures in the quoted papers.
4.2 Total Mass in Singles and Binaries
The destruction of the binary population is shown in Fig. 7, which is comparable with Fig. 2a in Gao et al. The destruction of binaries takes place at a quite similar rate in the two models, but escape of single stars is much larger in the body simulations. Indeed in the body models this somewhat diminishes the early increase in the numbers of single stars which is so prominent in the FokkerPlanck model, and which is caused by the destruction of binaries; nevertheless this feature is still visible in our body results. In the hybrid model of Giersz & Spurzem (2000) the number of single stars decreases monotonically, though the residual mass in binaries at is very similar to the value in our Fig. 7.
The cumulative mass lost by escape (separately in binary and single stars) is shown in Fig. 8, which can be compared with Fig. 27 in Giersz & Spurzem (2000). The mass of single escapers at time is well matched, but we find that the mass of binary escapers is about 50% larger than in the hybrid model.
In our model with , the decrease of the total mass by time is rather similar to that found by Fregeau et al (2003, Fig.4; 2005, Fig.2), but is about 50% larger in our model than theirs by time . We find that the total potential and kinetic energy of the bound members varies with in a manner very similar quantitatively to the result depicted in their Fig. 2 (Fregeau et al. (2003)).
4.3 Energy Distribution of Binaries
Fig. 9 shows the evolution of the distribution of internal binding energies of the binaries, and can be compared with Fig. 2b in Gao et al. Both diagrams use a logarithmic energy scale (abscissa), but with different definitions; for Gao et al. it is based on the discrete energy levels they adopted, whereas in our case it is determined by the binning in the output of NBODY6. The range of energies shown in the figure of Gao et al. (which is slightly smaller than the range of energies in the initial conditions) corresponds in our diagram to the range , where is defined in the caption.
Harder binaries were not plotted by Gao et al., as they were treated as being dynamically inert. By the end of their run this group of binaries accounted for about one third of the remaining binaries. In our body simulations, very few such extremely hard binaries are found. Simple theoretical ideas show that a binary is liable to be ejected (in a binarysingle interaction) if (Goodman, 1984), where is the central potential. For a Plummer model this gives , which corresponds to a value of about in the units of the horizontal axis of Fig. 9, in good agreement with our results. The central potential deepens by a factor of only about 2 (Fig. 15).
We consider that the escape of very hard binaries in the body simulations accounts for the main differences between our result (Fig. 9) and that of Gao et al. While we find that the distribution of (at fixed time) declines at the highest binding energies, in the results of Gao et al. the distribution of binding energy is monotonically increasing. Both kinds of simulation agree in the depletion of lesshard pairs, though in Gao et al. the depletion extends to somewhat harder pairs. To quantify these effects briefly, we remark that, for , the lower 10 percentile increases by a factor of about 3.5, and the mean increases by a factor of about 2.3, which is smaller because of the preferential destruction of softer binaries.
4.4 Spatial Distribution of Binaries
The evolution of the spatial distribution of binaries is indicated in two ways in Figs. 10 and 11, which can be compared with Fig. 3 in Gao et al.
Our results for the evolution of the halfmass radius of the single and binary stars (Fig. 10) are very similar to those of Gao et al. and Giersz & Spurzem (2000), and rather at variance (at least in later stages of the evolution) with those of Fregeau et al. (2003); their Fig. 4 indicates the halfmass radius of the binaries exceeding its initial value by about , and exceeding the current halfmass radius of the single stars after about 120. In Fregeau et al. (2005) the evolution of the halfmass radius is, if anything, still faster.
Within the core the agreement with Gao et al. is poorer. Our values for the density ratio (in binaries and singles: Fig. 11) are considerably smaller than the values found by Gao et al.; those authors find towards the end of the steady binary burning phase. This quantity is dependent (see Sec. 5.2.2), but we do not consider this sufficient to explain the disagreement. (Extrapolation from our data suggests for the Gao et al. number of particles).
From the dynamical point of view the early core contraction and increase in can be seen as manifestations of mass segregation. Spitzer (1969) showed that this could lead to equipartition between two species (here binaries and single stars) if the initial mass in the heavier species was sufficiently small. His criterion corresponds, in the present situation, to . Therefore all of the evolution we observe takes place in the regime corresponding to the mass segregation instability.
It is worth considering whether the topics discussed in this subsection and the last are interconnected, i.e. whether the distribution of binding energies is different at different radii. Fig. 24 of Giersz & Spurzem (2000) implies that the removal of soft pairs mainly affects binaries in the core. They also find a gap in the numbers of binaries at intermediate radii, which was predicted by Hut et al. (1992). These trends seem to be confirmed in our simulations (see Fig. 12 bottom right panel), but the number of surviving binaries is too low to draw firm conclusions.
4.5 Lagrangian Radii
We have already discussed the evolution of the halfmass radius. The evolution of other Lagrangian radii is displayed in Fig. 13, which can be compared with Fig. 5 in Gao et al., except that the choice of mass fractions differs at some points from theirs. The differences in the Lagrangian radii are at the level of quantitative detail. For example, we find that the 10% Lagrangian radius grows by about 0.1 dex up to , whereas Gao et al. find an increase by about 0.3 dex.
4.6 Structural Parameters
Gao et al. point out that much of the evolution of their model can be described approximately as homological. Part of the evidence for this is their finding that some dimensionless measures of their model vary little over the evolution, even though individual quantities, such as the halfmass radius, may change by large factors. We find, as do Gao et al., that the structural parameter is almost constant (Fig. 14; note the magnified vertical scale).
The value of is often adopted for this scaled energy in simple models of bound stellar systems (Spitzer 1987, equation (110)).
Another important quantity is the central potential, a scaled form of which is presented in Fig. 15. Though our scaling differs from that adopted by Gao et al., the fractional change between the initial value and the value in the postcollapse phase is comparable.
5 Dependence on N and binary fraction
5.1 Core evolution
5.1.1 The time of core collapse
Determination of the time of core collapse requires first a careful operational definition of this term, which is much less clear in these models than it is in the absence of primordial binaries. If is large enough, when there is a deep and sharp core collapse (Fig. 1), and the determination of the time of core collapse is relatively unambiguous. It does, however, vary from one simulation to another (Spurzem & Aarseth, 1996). In the presence of primordial binaries the problem is illustrated in Fig. 16. The initial collapse of the core is almost linear, and we may determine the point at which this ends () by finding the intersection of two linear fits to the data, for and respectively. This has been done iteratively. Alternatively, one may attempt to determine the end of core collapse by finding the minimum of . This requires smoothing of the data, and in general the result will depend on the smoothing interval.
Both methods present advantages and disadvantages. The linear fitting method is very robust for the determination of the first line (i.e. for ), which plays the major role for setting the time of core collapse. In fact this first line is much steeper than the second (that fits the core radius in the phase of post collapse). This second fit is more sensitive to the noise and to the extent of the fitting window (and thus in particular to the end time of the simulation), so that the uncertainty associated with the determination of the core radius at the time of core collapse, defined as the intersection of the two linear fits, is larger than the uncertainty in the time. It may be also argued that, with this definition, the core radius determination is somewhat artificial, since it does not correspond to a realized (or averaged) value in the simulation. In fact, at the time of core collapse the core radius in the simulation is bigger than the intersection of the two fits (see Fig. 16).
The second method, i.e. the determination of the time of core collapse as the time at which the minimum value of is attained, may be significantly affected by the intrinsic random fluctuations in the measure, due to its sensitivity to the local details of the curve. After some experimentations, we found that a good choice for minimizing the errors is given by adopting a smoothing window with a width of for and slightly bigger () for . We note, however, that the local nature of this estimate can play also a welcome role, since to determine the time of core collapse, and the associated core radius, the simulations may be stopped just shortly after core collapse. In contrast, the double linear method requires data for several tens of relaxation times during the postcollapse phase.
To evaluate the relative performance of the two estimators we have applied both methods to a large number of simulations (a major subsample of the simulations with primordial binaries; see Table 1). As can be guessed from Fig. 16, we found a systematic bias of approximately between the two methods in the measure of the time of core collapse, with the linear fit giving smaller times. The numerical experiments also confirm that the variance of the measure is marginally smaller ( less) for the linear fit method; the results are only very weakly Ndependent. The values for the core radius at core collapse are, on the other side, similar for the two methods.
In what follows we present results from the second method (minimum of ) for determining and , and the possible limitations of its significance need to be borne in mind. The main reason for choosing this procedure is the fact that this has allowed us to consider as valid simulations all those that reached at least after core collapse. In fact, while formally all the runs should have lasted for at least , some of them were stopped prematurely due to numerical difficulties. Given the large number of experiments (see Table 1), it was impractical for us to manually restart all the runs that required assistance. For each entry in Table 1 we have, however, at least one completed simulation, and on average the completeness of the sample is above .
Our results are shown in Fig. 17. From this figure it may be concluded that, at a given value of , the time of core collapse is essentially independent of the initial binary fraction , if %, and smaller than the collapse time of a model without primordial binaries by a factor of roughly . This compares well with the expectation based on a two components model with single stars only (and mass ratio 1:2) studied by Inagaki & Wiyanto (1984), that would predict a relative ratio of the order . While this lack of dependence on may seem surprising, it should be observed that, at late times in core collapse, the binary fraction in the core has increased considerably (Fig. 21), to 50% by number (or more) even if % initially. Therefore the core is practically saturated with binaries, approximately independent of the initial value of . Under these circumstances it is not surprising that the behavior of the core varies little with . Further evidence of this “saturation” effect is visible below in Figs.18 and 19.
Fregeau et al. (2003) give a plot showing that the “corecollapse time” exhibits a roughly linear increase with . This is quite at odds with our result, but the issue is a semantic one. Fregeau et al. (2003) do not call the initial contraction “core collapse”, preferring to reserve this term for the deep collapses akin to gravothermal oscillations; these occur much later, when the binaries are almost fully depleted. Thus the time they discuss might correspond to the deeper collapses in Fig. 5 after about , as discussed in Sec. 4.1, whereas our figure refers to the end of the initial phase of contraction at about .
5.1.2 Core radius in the phase of steady binary burning
The information from our simulations on this topic is summarised in Fig. 18. We notice again a saturation effect for initial binary fractions , and a significant variation with .
To explore these dependences further, we now compare the ratio of the core to half mass radius in the postcollapse phase with the theoretical estimate given by Vesperini & Chernoff (1994). By balancing the energy production in the core with the rate of expansion at the half mass radius, they get the following estimate for this ratio (their eq. 9):
(3) 
where the various quantities are defined as follows. The quantity is a parameter depending on and , where is the ratio of the expansion timescale to , and and are the velocity dispersion in the core and at the halfmass radius, respectively; we have assumed typical values for these parameters: and (Heggie & Aarseth, 1992; Goodman & Hut, 1989). The quantity is the binary fraction in the core defined as
Finally, and are coefficients for the efficiency of binarysingle and binarybinary burning and depend on the distribution of binding energy of the binaries. These coefficients have been computed by a numerical implementation of Equation 9 in Vesperini & Chernoff (1994), which has been checked against the results given by those authors. For simplicity, to compute the interaction crosssections we have assumed a flat distribution in between and kT, where the upper and lower limits have been estimated from the measured binding energy distribution (see Fig. 9).
For a fixed initial binary fraction this model is compared with the results of our simulation in Fig. LABEL:fig:comparison ( is assumed to be constant). Clearly the dependence is different. This result appears indeed rather puzzling. In fact it is unlikely that differences in the current fraction of binaries in the core, , are responsible, because of the saturation effect. The and coefficients are also independent, since we have checked that the distribution of binary binding energy, that is the only relevant input, evolves remarkably similarly in simulations with from to . In the worst case, our simplified approach (i.e. assuming a uniform distribution in for computing the analytical model) thus introduces only a constant multiplicative bias in eq. 3. In addition, this bias should be limited, since Figs. 5 and 6 in Vesperini & Chernoff (1994) suggest that the core radius is not particularly sensitive to the upper and lower limits of this distribution, at least within the limits of variation we observe (Fig. 9). Other factors which may influence the efficiency of binary heating are the total mass and central potential, but varies with by at most about up to and does not show any systematic trend with (with an average value, after the core collapse, of in NBODY6 units for our runs with ).
An interesting issue, raised by the referee, is the time scale on which binaries release energy. The time scale for close interactions between a binary of semimajor axis and single stars of space density may be estimated from the cross sections of Hut (1984) as approximately , where is the relative velocity of a binary and third body when far apart, and is the stellar mass. Comparing with the local relaxation time (Spitzer (1987)), we find that , where is the root mean square threedimensional velocity of single stars and is the Coulomb logarithm. Since this depends on there is an dependent regime of binaries which are effectively inert on the timescale of a given simulation, and this potentially introduces an dependence in the efficiency of binary interactions. On the other hand the time scale for interactions involving our hardest primordial binaries is approximately . Since the central relaxation time is much smaller than the halfmass relaxation time, we conclude that all primordial binaries have enough time to interact within the time scale we have considered (up to ).
We have also measured the value for in the simulations, confirming, immediately after core collapse, the values reported by Heggie & Aarseth (1992), i.e. , and observing basically no dependence. In this regard we note that Gao et al. measured . From our set of simulations this value appears to be unexpectedly low, even when account is taken of the fact that their model is isotropic (cf. Fig.2 in Takahashi 1996). Though we are unable to explain the discrepancy between the values of it does offer an alternative explanation of the fact that our values of appear to disagree with those of Gao et al (cf. our Sec.4.1). According to the theory of Vesperini & Chernoff a change in the value of leads to a change in the value of consistent quantitatively with what is found.
Finally, the differences in the value of at core collapse could be due to variations of the central velocity dispersion , which could have a significant impact, as . From our simulations we do not observe a significant dependence for this quantity, which appears to be close to the values quoted by Goodman & Hut (1989), . To be sure a limited dependence is present, but this goes in the opposite direction with respect to the one required to obtain a better match in Fig. 18, so that if this trend is taken into account the comparison with the analytical estimate for becomes worse: for primordial binaries we measure with , for and for . In contrast Gao et al. report a lower value, , so that this low central velocity dispersion could account for the smaller core in their Monte Carlo run, although it is not clear what is the origin of such a low velocity dispersion. Further studies would be required to clarify this point.
To summarize, from the data obtained in our simulations each term that contributes to the Vesperini & Chernoff (1994) model cannot explain why exhibits a steeper variation in than expected. Since the number of particles for the lower tail of the curve in Fig. 18 is rather modest, one possible explanation could be due to the importance of granularity effects so that a mean field approximation is not fully justified (see also the concluding discussion in Sec. 6).
At fixed a comparison with eq.(3) is shown in Fig. 19. Here is kept as a free parameter. In this case, there is satisfactory quantitative agreement within the scatter of different simulations at each value of , except possibly for the highest values. (Note that we here use the current value of the binary fraction, , averaged over the same time interval, and not the initial value. The two groups on the right, i.e. at and correspond to initial values of and , respectively.) For (not shown) we obtain comparable agreement for all , except for an indication that the theory predicts slightly too small a core radius for the very highest and the very lowest values of . It must be borne in mind, however, that the core radius in a multicomponent system is not a well defined concept. Our data for do not exhibit the slightly puzzling clumping at in Fig. 19, or the systematic offset at . Our dataset for is considerably larger, and we conclude that these features in Fig. 19 are results of the more modest sample size.
5.2 The binary and single populations
5.2.1 Total mass in singles and binaries
When the initial binary fraction is 10%, the binaries are mostly destroyed in the first 100. By the end of the simulation they account for only about 5% of the mass (Fig. 7). The situation is quite different for a simulation which contains only binaries initially (Fig. 20). By the same point in this simulation the binaries account for about 80% of the remaining mass.
5.2.2 Spatial Distribution of Binaries
The concentration of binaries in the core, after the end of core collapse, is presented as a function of and in Fig. 21. Some comparison is possible with the results of Gao et al. at fixed (their Fig. 3b), who found that the ratio increased by dex from to . An increase by this factor is consistent with what we find for all that we studied, but we repeat that the actual values we find are much smaller than those displayed by Gao et al. (Sec. 4.4).
While Fig. 21 confirms that the concentration of binaries in the core increases with the original binary fraction, , their overall spatial distribution differs less from that of the single stars, as we consider models with larger . This result is illustrated in Fig. 22, which should be compared with Fig. 10 (see caption to Fig. 22). In this respect our result differs from that of Fregeau et al. (2003): as we have already remarked (Sec. 4.4), when % they find that the halfmass radius of binaries exceeds that of singles for . This does not happen in their run with 20% binaries initially (their Fig. 5), where at all times up to the end of their run, at .
6 Conclusion
This paper has had several objectives. Our first has been to provide evidence from body simulations of the evolution of highly idealised stellar systems containing an initial population of binaries. The FokkerPlanck models of Gao et al. have become an important touchstone for subsequent modelling, but until now it has not been clear whether any disagreement with their results is due to the simplifying assumptions they made by adopting an approximate method. What we have tried to do with our simulations is to establish what actually happens, so that subequent research can rest on fewer assumptions and approximations. In this respect our work can be considered as a stepping stone to improve the physical understanding of the results from more complex and realistic numerical simulations, such as those recently performed by Portegies Zwart & McMillan (2004).
In line with this objective we have presented detailed results following closely the exposition of Gao et al., which covers the evolution of the structure of the system (core and halfmass radius, and structural parameters such as the scaled central potential), and the evolution of the binaries: their total numbers, their abundance in the core, their energies, and so on.
Our second objective has been to begin a discussion of the differences between the various models which have appeared in the literature so far, with particular reference to those of Gao et al., Giersz & Spurzem (2000) and Fregeau et al. (2003, 2005). Here our discussion has covered several of the above issues, and also the gravothermal oscillations which have been observed in almost all simulations, but with somewhat contradictory results for some characteristics, such as the time of onset of the oscillations. We have also compared our results with some of the small amount of quantitative theoretical predictions in the literature, mainly those of Vesperini & Chernoff (1994). With respect to this analytical estimate for the core to half mass radius in the postcollapse phase we find good quantitative agreement (at any fixed value of ) for the dependence on the initial binary ratio. On the other hand the observed dependence of this quantity is much steeper than one would expect on the basis of the theoretical arguments given by Vesperini & Chernoff (1994). This could be simply a low effect, but possibly some of the assumptions underlying the model could be inaccurate.
This issue is relevant to the comparison of the core radius between our models and FokkerPlanck or Monte Carlo models. With the present data we cannot give a clear answer to this point, since either extrapolation or direct simulations with would be required. A different way of gaining some insight into this problem could be to run some comparison simulations with Monte Carlo methods with in the range considered in this paper. If our results would be reproduced, they could be used to indirectly validate Monte Carlo simulations with high enough to answer this open question.
Our third objective has been to attempt to add some new insight and information on the nature of the collisional evolution of stellar systems with an initial population of binaries. In particular we have considered the time of core collapse (the definition of which we have attempted to clarify). For the period of steady binary burning which follows this initial collapse, we have determined the manner in which the core properties (core radius and binary fraction) depend on the initial binary fraction and particle number.
It is not clear which of these conclusions are applicable, even qualitatively, to real star clusters, with an initial mass function, a largely unknown initial distribution of binary parameters, a tidal field, and stellar and binary evolution. For example, Ivanova et al. (2005) have argued (on the basis of a simplified treatment of the core) that the binary fraction in the core of a dense cluster may be as little as 5%, even if initially. Now Fregeau et al. (2003) find (with a Monte Carlo code) that a system in a tidal field with initially disrupts after at most 45. If we were to guess from our Fig. 20, we would conclude that the binary fraction at this time would actually exceed 70%. Of course this argument ignores the essential physical issues of stellar evolution, collisions, etc., which underly the conclusion of Ivanova et al. (2005). On the other hand, all these models neglect factors which may be essential. Our next paper on this topic will provide detailed information on one of these which we have ignored so far – the influence of the tide.
7 Acknowledgments
We are indebted to Sverre Aarseth for the provision of NBODY6 and we are grateful to Enrico Vesperini for helpful discussions and suggestions. The input from the referee, David Chernoff, was both helpful and interesting. DCH and MT thank the Institute for Advanced Study for its hospitality while part of this work was being carried out. PH and MT similarly thank the School of Mathematics at Edinburgh University for its hospitality during later stages of this project. MT thanks Prof. Mineshige for his kind hospitality at the Yukawa Institute at Kyoto University, through the GrantinAid 14079205 of the Ministry of Education, Culture, Sports, Science and Technology, Japan. PH thanks Prof. Ninomiya for his kind hospitality at the Yukawa Institute at Kyoto University, through the GrantinAid for Scientific Research on Priority Areas, number 763, ”Dynamics of Strings and Fields”, from the Ministry of Education, Culture, Sports, Science and Technology, Japan. The numerical simulations have been performed on the Condor Cluster of the Institute for Advanced Study.
References
 Aarseth (2001) Aarseth S., 2001, in ASP Conf. Ser. 228: Dynamics of Star Clusters and the Milky Way The Formation of Hierarchical Systems. p 111
 Aarseth (2003) Aarseth S., 2003, Gravitational Nbody Simulations. Cambridge University Press
 Aarseth & Heggie (1993) Aarseth S. J., Heggie D. C., 1993, in ASP Conf. Ser. 48: The Globular ClusterGalaxy Connection A 6000Body Simulation with Primordial Binaries. p 701
 Albrow et al. (2001) Albrow M. D., Gilliland R. L., Brown T. M., Edmonds P. D., Guhathakurta P., Sarajedini A., 2001, Ap.J., 559, 1060
 Bellazzini et al. (2002) Bellazzini M., Fusi Pecci F., Messineo M., Monaco L., Rood R. T., 2002, A.J., 123, 1509
 Bettwieser & Sugimoto (1985) Bettwieser E., Sugimoto D., 1985, MNRAS, 212, 189
 Casertano & Hut (1985) Casertano, S. and Hut, P. 1985, Ap.J., 298, 80
 de La Fuente Marcos (1996) de La Fuente Marcos R., 1996, A&A, 314, 453
 Fregeau et al. (2003) Fregeau J. M., Gürkan M. A., Joshi K. J., Rasio F. A., 2003, Ap.J., 593, 772
 Fregeau et al. (2005) Fregeau J.M., Gurkan M.A.,. Rasio F.A., 2005, preprint (http://arxiv.org/abs/astroph/0512032)
 Gao et al. (1991) Gao B., Goodman J., Cohn H., Murphy B., 1991, Ap.J., 370, 567
 Giersz & Spurzem (2000) Giersz M., Spurzem R., 2000, MNRAS, 317, 581
 Giersz & Spurzem (2003) Giersz M., Spurzem R., 2003, MNRAS, 343, 781
 Goodman (1984) Goodman, J. 1984, Ap.J., 280, 298
 Goodman (1987) Goodman, J. 1987, Ap.J., 313, 576
 Goodman & Hut (1989) Goodman, J. & Hut, P. 1989, Nature, 319, 40
 Heggie & Aarseth (1992) Heggie D. C., Aarseth S. J., 1992, MNRAS, 257, 513
 Heggie & Mathieu (1986) Heggie D. C., Mathieu R. D., 1986, in LNP 267: The Use of Supercomputers in Stellar Dynamics Standardised Units and Time Scales. p 233
 Heggie & Hut (2003) Heggie D., Hut P., 2003, The Gravitational MillionBody Problem: A Multidisciplinary Approach to Star Cluster Dynamics. Cambridge University Press
 Hurley et al. (2001) Hurley J. R., Tout C. A., Aarseth S. J., Pols O. R., 2001, MNRAS, 323, 630
 Hurley et al. (2005) Hurley, J.R., Pols, O.R., Aarseth, S.J. and Tout, C.A. 2005, preprint (http://arxiv.org/abs/astroph/0507239)
 Hut (1984) Hut P., 1984, ApJS, 55, 301
 Hut (1989) Hut, P., 1989, in Dynamics of Dense Stellar Systems. ed. David Merritt (Cambridge University Press), pp. 229236
 Hut et al. (1992) Hut P., McMillan S., Goodman J., Mateo M., Phinney E. S., Pryor C., Richer H. B., Verbunt F., Weinberg M., 1992, PASP, 104, 981
 Hut et al. (1992) Hut, P., McMillan, S., Romani, R., 1992, Ap.J., 389, 527
 Inagaki & Hut (1988) Inagaki S., Hut P., 1988, fbp..coll, 319
 Inagaki & Wiyanto (1984) Inagaki, S., Wiyanto, P., 1984, PASJ, 36, 391
 Ivanova et al. (2005) Ivanova, N., Belczynski, K., Fregeau, J. M., & Rasio, F. A. 2005, MNRAS, 358, 572
 Kroupa (1995) Kroupa P., 1995, MNRAS, 277, 1522
 McMillan & Hut (1994) McMillan S., Hut P., 1994, Ap.J., 427, 793
 McMillan et al. (1990) McMillan S., Hut P., Makino J., 1990, Ap.J., 362, 522
 McMillan et al. (1991) McMillan S., Hut P., Makino J., 1991, Ap.J., 372, 111
 Makino (1996) Makino, J. 1996, Ap.J., 471, 796
 Portegies Zwart et al. (2004) Portegies Zwart S., Hut, P., McMillan, S. L. W., Makino, J., 2004, MNRAS, 351, 473
 Portegies Zwart & McMillan (2002) Portegies Zwart S. F., McMillan S. L. W., 2002, Ap.J., 576, 899
 Portegies Zwart & McMillan (2004) Portegies Zwart S. F., McMillan S. L. W., 2004, ASP Conference Series, astroph/0411188
 Pulone et al. (2003) Pulone L., De Marchi G., Covino S., Paresce F., 2003, A&A, 399, 121
 Spitzer (1969) Spitzer, L. J. 1969, Ap.J., 158, L139
 Spitzer (1987) Spitzer, L. 1987, Dynamical Evolution of Globular Clusters. Princeton, NJ, Princeton University Press, 1987
 Spurzem & Aarseth (1996) Spurzem, R., & Aarseth, S. J. 1996, MNRAS, 282, 19
 Takahashi (1996) Takahashi K., 1996, PASJ, 48, 691
 Takahashi & Inagaki (1991) Takahashi K., Inagaki S., 1991, PASJ
 Vesperini & Chernoff (1994) Vesperini, E., & Chernoff, D. F. 1994, Ap.J., 431, 231
 Wilkinson et al. (2003) Wilkinson M. I., Hurley J. R., Mackey A. D., Gilmore G. F., Tout C. A., 2003, MNRAS, 343, 1025