# Consequences of Dynamical Disruption and Mass Segregation for the Binary Frequencies of Star Clusters

###### Abstract

The massive (13,000–26,000 M), young (15–30 Myr) Large Magellanic Cloud star cluster NGC 1818 reveals an unexpected increasing binary frequency with radius for F-type stars (1.3–2.2 M). This is in contrast to many older star clusters that show a decreasing binary frequency with radius. We study this phenomenon with sophisticated -body modeling, exploring a range of initial conditions, from smooth virialized density distributions to highly substructured and collapsing configurations. We find that many of these models can reproduce the cluster’s observed properties, although with a modest preference for substructured initial conditions. Our models produce the observed radial trend in binary frequency through disruption of soft binaries (with semi-major axes, AU), on approximately a crossing time ( Myr), preferentially in the cluster core. Mass segregation subsequently causes the binaries to sink towards the core. After roughly one initial half-mass relaxation time ( Myr) the radial binary frequency distribution becomes bimodal, the innermost binaries having already segregated towards the core, leaving a minimum in the radial binary frequency distribution that marches outwards with time. After 4–6 , the rising distribution in the halo disappears, leaving a radial distribution that rises only towards the core. Thus, both a radial binary frequency distribution that falls towards the core (as observed for NGC 1818) and one that rises towards the core (as for older star clusters) can arise naturally from the same evolutionary sequence owing to binary disruption and mass segregation in rich star clusters.

###### Subject headings:

(stars:) binaries: general - galaxies: star clusters: individual (NGC 1818) - (galaxies:) Magellanic Clouds - stars: kinematics and dynamics - methods: numerical## 1. Introduction

Observations of star-forming regions, young and old star clusters, and the Galactic field indicate that most stars reside in binary or higher-order multiple systems (e.g. Kouwenhoven et al., 2005, 2007; Raghavan et al., 2010; Geller et al., 2010; Kraus et al., 2011; Geller & Mathieu, 2012; King et al., 2012). Since most stars (with masses M) form in clusters or groups (Lada & Lada, 2003), many of which quickly dissolve to populate the Galactic field (Adams & Myers, 2001), connecting observed binary frequencies within these different environments is critical to our understanding of star formation. In dense star clusters, the frequency of binary stars can be significantly modified by close encounters with other stars, and these processes can be studied in detail through sophisticated -body simulations.

In most star clusters, the rates of binary disruption (particularly for wide binaries) greatly exceed the rates of binary creation. As explained by Heggie (1975), binaries that have low binding energies relative to the kinetic energies of stars within the cluster, known as “soft” binaries, tend to become even less bound (“softer”) on average as a result of stellar interactions, and are eventually disrupted. Conversely binaries with high binding energies relative to the kinetic energies of stars within the cluster (“hard” binaries) become more tightly bound (“harder”). Most hard binaries are unlikely to be disrupted by stellar encounters, but dynamical hardening and binary evolution processes can also destroy very hard binaries. These dynamical processes are expected to proceed most rapidly in the denser core of the cluster where stellar encounters are most frequent.

Offsetting these disruption processes is the preferential tidal stripping of low-mass single stars from the cluster by the Galactic potential and their ejection from the core due to dynamical encounters, which, after an initial stage of rapid soft binary disruption, can lead to a roughly constant global binary frequency over many Gyr (Hurley et al., 2005; Geller et al., 2013).

Locally, at different radii within a cluster, the binary frequency also evolves owing to two-body relaxation processes and dynamical friction. On average, binaries have a higher total mass than their single-star counterparts, and therefore energy exchange between these two groups tends to cause binaries to sink towards the core of the cluster. This dynamical mass segregation results in the binary frequency rising in the cluster core and falling in the halo. Indeed, many Milky Way open and globular clusters show a rising binary frequency towards the cluster core, which is interpreted to be the result of mass segregation (e.g. Mathieu & Latham, 1986; Geller & Mathieu, 2012; Milone et al., 2012).

Binary disruption and mass segregation processes compete to determine the radial distribution of the binary frequency within a dense star cluster. Observations of star clusters of different ages allow us to study empirically the timescales for these processes and to verify the predictions from -body simulations. The rich star cluster NGC 1818, located in the Large Magellanic Cloud (LMC), is the youngest rich cluster where the radial dependence of the binary frequency has been measured (Elson et al., 1998; de Grijs et al., 2013; Li et al., 2013b), and is therefore very important for our understanding of the early evolution of binaries in star clusters.

NGC 1818 has an age of 15 - 30 Myr, a total mass of 13,000 M to 26,000 M,
a central surface mass density of M pc (de Grijs
et al., 2002; Mackey &
Gilmore, 2003), and a total binary frequency^{1}^{1}1We define the
binary frequency as , where is the number of binaries, is the number of single stars, and “…” signifies
higher-order multiples. () estimated to be between
55% and 100% for F-type stars with masses between 1.3 M and 1.6 M (Hu et al., 2010).
Importantly, the binary frequency in NGC 1818 is observed to decrease towards the core of the cluster (de Grijs
et al., 2013; Li
et al., 2013b), in
stark contrast to the more typically observed radial dependence where the binary frequency rises towards the core, as also predicted from
mass segregation processes.

de Grijs et al. (2013) and Li et al. (2013b) suggest that this radial trend of the binary frequency observed in NGC 1818 may result from the processes of binary disruption. In this paper we test this hypothesis through sophisticated -body modeling of the cluster. In Section 2, we describe the simulation method and define our grid of cluster simulations. In Section 3, we discuss the evolution of the radial dependence of the binary frequency and the contribution of dynamical binary disruption and mass segregation. Then, in Section 4, we compare the simulations directly with the observations and discuss our ability to reproduce the observed trend in binary frequency within our models. In Section 5 we discuss the implications of these results for the origins of the radial binary frequency distribution in NGC 1818 and relate this result to similar observations in other rich star clusters. Finally, in Section 6 we provide our conclusions.

## 2. Simulation Method

We use the NBODY6 code (Aarseth, 2003) to model the dynamical evolution of rich star clusters, with the goal of reproducing the observations of NGC 1818 at the cluster’s age. NBODY6 includes stellar and binary evolution (Hurley et al., 2000, 2002), models dynamical encounters with binaries in detail, including those leading to the disruption of binaries, and mass segregation arises naturally in NBODY6 from two-body relaxation processes. Much of our method is identical to that of Geller et al. (2013), except that here we choose some different initial conditions for the cluster, which we describe below.

Observations of NGC 1818 suggest that the cluster has a total mass between 13,000 M and 26,000 M (de Grijs et al., 2002; Mackey & Gilmore, 2003). We choose to begin all of our models with 36,000 stars chosen from a Kroupa (2001) IMF, with masses between 0.1 M (approximately the hydrogen-burning limit and the lowest stellar mass with detailed models guiding the stellar evolution code) and 50 M, which produces an initial cluster mass of 21,800 M. Extrapolating from the results of Hurley et al. (2005) and Geller et al. (2013) suggests that this initial mass will remain within the observed mass range by an age of 30 Myr (the maximum age estimate for the cluster), despite mass loss from stellar evolution, dynamical ejections and tidal stripping from the Galactic potential, and indeed our evolved NGC 1818 models confirm this result. We do not model the embedded phase of the cluster here. Instead we begin our simulations at after gas expulsion and with all stars on the zero-age main sequence.

Mackey & Gilmore (2003) fit the observed surface brightness distribution of NGC 1818 with an EFF profile (Elson et al., 1987), defined by

(1) |

with a central surface brightness of = 3.35 0.02 L pc, = 12.50 0.78 arcsec (3.04 0.19 pc, using the canonical LMC distance modulus of 18.5, which equates to a scale of 4.116 arcsec pc, as in Mackey & Gilmore 2003) and = 2.76 0.12. Also, for reference, one can calculate the King core radius from the EFF parameters as,

(2) |

which here equates to 2.45 pc. Mackey & Gilmore (2003) adopt a mass-to-light ratio of 0.08, which implies a central mass surface density of 180 4 M pc. We aim to reproduce these observations at the age of NGC 1818 within our simulations.

The EFF model is similar to a King (1966) or Plummer (1911) model (which are more typical initial conditions for -body star cluster simulations) except in the halo, where the EFF model maintains a slightly higher surface brightness and results in a slightly more extended cluster. For simplicity, here we begin our simulations with stars distributed according to a Plummer model, and we describe the details of these models below. At the age of NGC 1818, these simulations are consistent with the EFF model to within the radial extent of the de Grijs et al. (2013) and Li et al. (2013b) studies (see Figure 5).

Observations suggest that many clusters may form with subvirial velocities (Peretto et al., 2006; André et al., 2007; Tobin et al., 2009; Proszkow et al., 2009). Therefore we investigate different simulations with initial virial ratios of = 0.5 (equilibrium) and = 0.3 and 0.1 (collapsing). In order to determine the initial length scale, we first ran simulations with a range of initial virial radii for each value and determined the initial virial radius for a given that best reproduces the observed surface density profile at the age of NGC 1818. We find that an equilibrium model reproduces the observations with an initial virial radius of 7 pc (equivalent to an initial half-mass radius of pc and a Plummer scale radius of pc ). A model also reproduces the observations with an initial virial radius of 7 pc. A model requires an initial virial radius of 10 pc ( pc, pc). For this study, we do not simulate clusters with supervirial initial conditions, although it is conceivable that some supervirial initial conditions may reproduce the observed surface density profile of NGC 1818 at the cluster age.

Many young clusters are also observed to be substructured and well represented by fractal density distributions (e.g. Larson, 1995; Kraus & Hillenbrand, 2008; Cartwright & Whitworth, 2004; Sánchez & Alfaro, 2009). Therefore, in addition to the smooth models at different virial ratios, we also explore different “clumpy” models. We impose fractal distributions on top of the Plummer models described above (see Figure 1) using the McLuster code (Küpper et al., 2011), with slight modifications to match those we made to NBODY6 for defining the initial binaries (see Geller et al., 2013), which creates initial conditions that can be easily read directly into NBODY6. We follow Goodwin & Whitworth (2004) and investigate clusters with fractal dimensions of = 1.6, 2.0, 2.6, and = 3.0 (which corresponds to the smooth models with no clumping defined above). Again, we first ran simulations with a range of initial virial radii for each combination of and to identify the initial length scale that best reproduces the observed surface density profile at the cluster age.

To define the stellar velocities for the fractal initial conditions, McLuster first draws random velocities from a normal distribution such that the mean motion of the stars in each subgroup is zero and the mean magnitude of the velocities for stars in each subgroup is unity. In addition to this normally distributed velocity, the velocity of each “parent” is added to the velocities of each of its “children”, so each generation as an ensemble moves like its parent. Finally, the velocity of each star is multiplied by the circular velocity from the Plummer model at the star’s radius. Thus, since the velocities are normally distributed around the parents’ (and grandparents’) velocities, the initial velocities are coherent but random. There is evidence that most (and perhaps all) substructured clusters form with subvirial velocities (Girichidis et al., 2012), and furthermore in order to erase the initial substructure as rapidly as is observed, subvirial initial conditions may be required (Goodwin & Whitworth, 2004). For our purposes we build models with clumpy and subvirial as well as equilibrium initial conditions to investigate the correspondence of such models with observations of NGC 1818.

All of our simulations that begin with stars initially distributed in clumpy density distributions relax to smooth density distributions within the range of age estimates for NGC 1818 (see e.g., the bottom panel of Figure 1), although some of our most highly substructured and subvirial simulations retain some minor clumpiness at an age of 15 Myr (which is erased by 30 Myr).

We list in Table 1 the initial virial ratio (), initial fractal dimension () and initial virial radius () for the models that match the observed surface density distribution of NGC 1818 at the cluster’s age, as well as the number of simulations () and the percentage of lines of sight to the simulations that match the observed radial distribution of the binary frequency (, discussed in Sections 4 and 5). In the following, we will refer to a particular model using the convention of [,].

(pc) | (%) | |||
---|---|---|---|---|

0.5 | 3.0 | 7 | 20 | 28.13 |

0.5 | 2.6 | 8 | 2 | 99.3 |

0.5 | 2.0 | 9 | 2 | 85.4 |

0.5 | 1.6 | 11 | 2 | 78.0 |

0.3 | 3.0 | 7 | 10 | 17.03 |

0.3 | 2.6 | 9 | 2 | 90.7 |

0.3 | 2.0 | 10 | 2 | 53.3 |

0.3 | 1.6 | 11 | 2 | 97.4 |

0.1 | 3.0 | 10 | 10 | 22.77 |

0.1 | 2.6 | 10 | 2 | 94.7 |

0.1 | 2.0 | 11 | 2 | 50.5 |

0.1 | 1.6 | 11 | 2 | 20.3 |

^{†}

^{†}footnotetext: is the initial virial ratio. is the initial fractal dimension (with = 3.0 corresponding to non-fractal, smooth, initial conditions). is the initial virial radius of the Plummer model. is the number of simulations with the given [,,]. is the percentage of sight lines that reproduce the observed radial distribution of the binary frequency.

NGC 1818 is located at about 3.8 degrees (3.3 kpc) from the center of the LMC. We simulate the effects of the cluster orbiting within the potential of the LMC by placing the cluster at 3.3 kpc from the center of a point mass of M on a linearized circular orbit. Although the modeling of the cluster’s orbit could be more sophisticated, precise cluster orbital parameters are currently unknown and, more importantly, after only 30 Myr (the likely maximum age of NGC 1818) the effects on the cluster from the LMC tidal field are minimal.

Observations of NGC 1818 suggest a total binary frequency of up to 100% (at least for F-type stars; Hu et al., 2010). We choose to begin all of our models with a 100% binary frequency with no radial or primary mass dependence. We follow empirically defined initial distributions for the orbital parameters of the binaries in our models, as in Geller et al. (2013), which agree with observations of solar-type binaries in young open clusters (e.g. M35; Geller et al., 2010) and the field (Raghavan et al., 2010).

Specifically, we draw the initial binary orbital periods from a log normal distribution with a mean of [days]) = 5.03 and . We allow the initial binaries to populate the entire log normal distribution, even though some of these binaries are born soft, because we are interested in investigating the dynamical disruption of soft binaries within the simulations. In practice, the initial log normal distribution allows binaries with periods up to days. Consequently, a small fraction of binaries (1%) are assigned orbital separations that are larger than the distance from the binary’s center of mass to its nearest star. Such binaries would likely not form in real clusters, and are promptly disrupted dynamically at the start of the simulation, causing a drop in the binary frequency of about 0.2% to 5% inside of one and about 0.2% to 2% outside of one (depending on the initial density distribution). This decrease is similar in magnitude to the uncertainty in the binary frequency from Poisson counting statistics (of about 1% inside or outside of ), and is small compared to the 30% difference in binary frequencies of the inner- and outer-most bins observed by de Grijs et al. (2013) and Li et al. (2013b) in NGC 1818.

The NGC 1818 F-type binaries with mass ratios have a mass-ratio distribution consistent with , where or (de Grijs et al., 2013; Li et al., 2013b). In our simulations, we define the initial mass ratios by first taking two masses, and , from the Kroupa (2001) IMF (within the mass limits defined above), combining these masses ( = + ), and then choosing a new mass ratio from a uniform distribution such that the new primary and secondary masses sum to equal , and the new secondary mass is greater than 0.1 M but less than the new primary mass. For the mass range observed by de Grijs et al. (2013) and Li et al. (2013b), this procedure produces a mass-ratio distribution (over all ) that favors low values and approximately follows an distribution (as also suggested by Kouwenhoven et al. 2005, 2007 and Reggiani & Meyer 2011 for binaries of a range of spectral types). Finally, we draw the initial eccentricities from a Gaussian distribution with a mean of with . These binaries are then evolved through the Kroupa (1995) pre-main-sequence evolution prescription, with the same modifications as in Geller et al. (2013).

We produce multiple realizations of each simulation using different initial random seed values to address the stochastic effects present in -body simulations. Moreover, for all simulations with a given [,], the initial stellar and binary parameters (e.g., positions, velocities, binary periods and mass ratios, etc.) are all drawn from the same respective distributions, but each simulation randomizes the parameters to produce a unique initial stellar population. For the smooth (fractal dimension ) models, we ran 20 simulations and 10 simulations each for and 0.1. For the substructured models with a fractal dimension , we ran two simulations each. We discuss our method for combining the results from these simulations below. Our primary goal here is not to identify the specific combination of [] with which NGC 1818 most likely formed, but rather to investigate whether we can reproduce the observations of NGC 1818 with virial or subvirial and smooth or clumpy initial conditions. As we show below, this number of simulations is sufficient to answer this question.

## 3. Dynamical Evolution of the Binary Frequency

Before we compare the -body simulations to the observations of NGC 1818, we discuss in general the evolution of the radial distribution of the binary frequency. We choose here to focus on the 20 [0.5,3.0] (equilibrium and smooth) simulations, although the general evolutionary sequence we discuss below is common to all of our models. In Figure 2 we show the binary frequency (over all primary masses and binary mass ratios) as a function of radius from the cluster center. We show six representative times in the figure, namely the crossing time, Myr, and multiples of the initial half-mass relaxation time, Myr. (Although not plotted in this figure, the initial binary frequency has no radial dependence.)

After one crossing time, the binary frequency decreases towards the core of the cluster, due to the disruption of wide binaries. This is illustrated further in Figure 3, where we show the cumulative distributions of semi-major axes for binaries inside and outside of one initial half-mass radius () at a crossing time. Compared to the initial semi-major axis distribution, both distributions at are truncated at shorter separations by disruptive stellar encounters. The very wide binaries observed in the field are soft even in the cluster halo. Importantly, the semi-major axis distribution of binaries found inside of one is shifted to even shorter separations than that of binaries outside of one . This is also seen in Figure 4, where we show the maximum separation for binaries inside and outside of as functions of time. There is a steep drop in the maximum semi-major axis very early on due to the rapid disruption of very soft binaries from the primordial population in both the inner and outer cluster regions. However, the binaries outside of are able to retain companions at much wider separations than can survive inside of , up to about 8 .

This difference between the inner and outer binaries results from the higher velocity dispersion and higher density in the cluster core relative to the halo. The higher velocity dispersion leads binaries in the core to move more rapidly, on average, relative to other stars than do binaries in the halo. Thus encounters within the core are more energetic and can disrupt tighter binaries. Also, the higher density results in a higher encounter rate in the core than in the halo (Leigh & Sills, 2011). Furthermore, early in the cluster’s evolution, the stars do not have sufficient time to mix throughout the cluster, and instead most experience the dynamical environment near where they were born. These effects combine to produce a decreasing binary frequency towards the cluster core at , which is maintained for roughly one .

As time progresses in the simulations, cluster-wide mass segregation effects begin to control the radial dependence of the binary frequency. The arrows in the five lower panels in Figure 2 mark the theoretical “” values or “zones of avoidance” from, e.g., Mapelli et al. (2004) and Ferraro et al. (2012). In a given panel, represents the radius inside of which the local dynamical friction timescale (Binney & Tremaine, 1987) for a binary with a mass equivalent to that of the mean binary mass in the cluster is shorter than the simulated time. Qualitatively, the value predicts the radius inside of which the binaries should experience the effects of dynamical friction and therefore fall towards the center of the cluster. The result of this process is to increase the binary frequency in the core at the expense of the binary frequency towards the cluster halo.

At , the radial dependence of the binary frequency is high in the core, then drops to a minimum near (at roughly one ), and then rises again towards the halo. The rate of dynamical binary (and triple-system) creation in our models is far too low to account for this increase in the binary fraction (cf. Li et al., 2013b). Instead, this phenomenon is due to mass segregation processes. Binaries inside of have already fallen towards the core as the cluster begins to become mass segregated. However binaries outside of have not experienced enough dynamical friction to migrate fully towards the core. Also, at , the binary frequency at each radial bin has decreased to well below the values at . By this time more disruptive encounters have occurred locally, and also binaries that formed in the halo have had more time to orbit through the denser central regions of the cluster where encounters are more energetic.

Over time, the value marches out towards the halo. By the binary frequency maintains the rising distribution towards the core inside of but develops a roughly flat distribution outside of . By the binary frequency increases continuously from the halo to the core of the cluster. Also at , both the binaries inside and outside of have experienced sufficient encounters to regain similar semi-major axis distributions (see Figures 3 and 4).

As a check for the dependence of these results on the long-period cutoff of the initial binary orbital period distribution, we reanalyzed the [0.5,3.0] simulations considering all binaries with initial orbital periods days as single stars. This reanalysis results in the same binary frequency in the core of the cluster at (and 30 Myr), and a progressively lower binary frequency towards the halo, as compared to the original analysis. Binaries in the core with days are already broken up by one crossing time, where the binaries in the halo are not. The decreasing trend in binary frequency towards the core of the cluster remains. After one the radial distributions of the binary frequency from both analyses are indistinguishable, since nearly all of these wide binaries have been disrupted, even in the halo.

## 4. Comparison Between Observations and Simulations of NGC 1818

Now we compare the simulations directly with observations of NGC 1818. Here we choose to investigate individual simulations rather than combining the different simulations for a given [] pair, as in the previous section. We will focus specifically on three simulations defined by [0.5,3.0], [0.3,2.0], and [0.1,2.6], and with initial conditions that lead to a radial dependence of the binary frequency that matches particularly well the observations from de Grijs et al. (2013) and Li et al. (2013b) at the age of NGC 1818.

We choose to compare the simulations at an age of 30 Myr to the observations. At 15 Myr (the minimum age estimate for the cluster), the surface density profiles for these simulations are nearly identical to those at 30 Myr, although the subvirial simulations have a slightly higher central surface density, a residual of the early collapse and subsequent restructuring that occurs on roughly a crossing time. The binary frequencies in these three simulations are also slightly higher at 15 Myr than at 30 Myr (although still consistent with the observations), since the binaries have had less time to undergo disruptive dynamical encounters with other stars.

At 30 Myr, the total masses for the [0.5,3.0], [0.3,2.0] and [0.1,2.6] simulations are 19,518.6 M, 19,064.4 M and 17,664.8 M, respectively. For reference, the mean masses for all simulations with these respective [] pairings are 19,550 140 M, 19,160 100 M and 18,046 380 M (showing uncertainties of one standard deviation for the simulations and half of the range in masses for the two simulations with ). As the initial virial ratio decreases, the simulated clusters lose more mass early on, due to both violent relaxation processes that can eject more stars from the core of the cluster and the more loosely bound halo stars that are more vulnerable to tidal stripping (as the =0.1 and =0.3 simulations began with an initial virial radius of 10 pc, compared to 7 pc for the =0.5 simulation). These simulation masses agree well with the range of estimates from de Grijs et al. (2002) and Mackey & Gilmore (2003).

At 30 Myr, the total binary frequencies in the [0.5,3.0], [0.3,2.0], and [0.1,2.6] simulations are 73%, 78%, and 80%, respectively. For reference, the mean binary frequencies for all simulations with these respective [] pairings are 73% 1%, 77.9% 0.1%, and 79% 1% (with uncertainties defined as above for the masses). Again, these values agree well with the observational estimates for NGC 1818.

Next we investigate each simulation in projection, as an observer of NGC 1818 would only have access to two dimensions of spatial data. Strictly, the orbit of the simulated cluster within the LMC potential defines a preferred line-of-sight projection for the cluster relative to an observer on Earth. However, given the uncertainties in the true cluster orbital parameters and our simplistic modeling of the LMC potential, we choose not to define a particular projection to the cluster. Furthermore, asymmetries within the simulations beginning with cool and clumpy conditions are not expected to exactly reproduce the initial cluster conditions, and are only meant as statistical approximations. Therefore we select 1000 different lines of sight randomly distributed over the sphere and project each simulation along these sight lines. In the following we combine these results to show the range in possible observed projections for a cluster like NGC 1818.

In Figure 5, we show the 30 Myr mass surface density profile for these three simulations compared with the EFF model fit to the observations by Mackey & Gilmore (2003). Although these simulations began with Plummer profiles, all simulations reproduce the observed surface density profile to within the range from our 1000 sight lines out to the 20 pc limit of the de Grijs et al. (2013) and Li et al. (2013b) observations. These results are not particularly sensitive to the chosen projection, nor are they sensitive to the specific simulation chosen from these sets. Importantly, reproducing the observed surface density profile, total mass, and binary frequency indicates that we have correctly modeled the dynamical environment for binaries in the real cluster within our simulations.

In Figure 6 we investigate the cumulative radial distribution of the binary frequency at 30 Myr in the three simulations and compare these results to the observations in the bottom panel. Note that here we plot the results in cumulative form rather than binned as in Figure 2, to match the analysis of the observations. Moreover, here each point plots the binary frequency inside the given radius. For comparison, we plot in the top panel mean results from all simulations in the respective [] pairs in three-dimensional shells without limiting the simulations by primary mass or mass ratio. The [0.5,3.0] and [0.1,2.6] simulations show a relatively large range in binary frequencies at each bin between the individual simulations with the given [,] combination. We will return to this below.

The remaining three panels show only the individual simulations we focus on here (and which are also shown in Figure 5). In the second panel from the top, we show cumulative binary frequency distributions from these three simulations in projection for binaries of all primary masses and mass ratios. A visual comparison of the uncertainties plotted in this panel (which show the range within which lie 95% of the projections) reveals that the [0.5,3.0] simulation is most symmetric at 30 Myr, followed by the [0.3,2.0] simulation, and finally the [0.1,2.6] simulation.

In the bottom two panels, we only include stars with primary masses between 1.3 M and 2.2 M,
matching the mass range from the Li
et al. (2013b)
observations^{2}^{2}2Because the observations were selected to be within a range in magnitude, some near equal-mass binaries with
primary masses M are included in the observational sample that are not included in the analysis of the simulations,
and the opposite is true at the bright end of the observed magnitude range for primary masses M.
We assume that the contributions to the observed binary frequency from binaries gained and lost at the faint and bright ends of the magnitude range,
respectively, approximately cancel each other out..
Here, in general, the uncertainties
become larger towards the cluster core
relative to the distributions in projection using all primary masses and mass ratios.
The increased uncertainties are primarily due to Poisson noise owing to the decreased number of objects when selecting only the observed regime
in primary mass (and also mass ratio in the bottom panel).

Finally in the bottom panel we further limit the simulations to only include binaries with mass ratios (and count binaries with as single stars in our calculations of the binary frequencies). In this panel we compare the simulations directly to the observations from Li et al. (2013b). All simulations closely reproduce the observations.

To test the relative difference between the observations and a given simulation, we performed tests comparing all 1000 lines of sight to the three simulations shown in Figure 6 against the observations (bottom panel). Specifically, for each line of sight, we calculate

(3) |

where and are the simulated and observed cumulative binary frequency inside the given radius, respectively, and is the error on the observed cumulative binary frequency at the same radius. We have not constrained the simulation results, and therefore we have eight degrees of freedom (one for each bin in Figure 6). We find that, respectively, 97.4%, 99.9%, and 100% of the [0.5,3.0], [0.3,2.0], and [0.1,2.6] projections are indistinguishable from the observations (at the level). Thus the specific line of sight for these simulations does not change the correspondence with the observations.

In Table 1 we show all [] pairs that reproduce the observed surface density profile of NGC 1818 at 30 Myr, and also the combined percentages of sight lines to these simulations that result in a trend in binary frequency with radius that is indistinguishable from the observations (). Specifically, we perform tests comparing the Li et al. (2013b) observations to all 1000 sight lines towards a given simulation as described above. Then is defined as the percentage of sight lines to all simulations for a given [] pairing which result in a -value of less than .

## 5. Discussion

All of the simulations in Table 1 reproduce the observed surface density profile of NGC 1818 at 30 Myr, and certainly there are more combinations of initial [,,] that will also reproduce these observations. About 60% of the sight lines to these simulations (the mean of the 12 values in Table 1) also reproduce the observed radial trend in binary frequency. Furthermore, if we take the average value at each radius of the cumulative binary frequency at all 1000 projections for all simulations of a given initial [,,], we find that two thirds of these sets of initial conditions evolve to be indistinguishable from the observations at 30 Myr. Moreover, the available observations do not constrain the -body simulations to one particular set of initial conditions. Instead, it appears that a wide range of initial virial ratios and fractal dimensions can reproduce the observations at the age of NGC 1818.

Although our primary goal here is not to define the precise primordial conditions with which NGC 1818 was likely born, the values we calculate here for our simulations do indicate a preference towards clumpy initial conditions. The mean values for the four respective = 0.1, 0.3, and 0.5 models at all values are 47% 17%, 65% 19%, and 73% 16%. Given the relatively large range in values for the simulations at each value of , we do not detect a preferred virial ratio for reproducing the observations. However, the mean values for the three respective models with = 3.0, 2.6, 2.0, and 1.6 at all values are 23% 3%, 95% 3%, 63% 11%, and 65% 23%. A clumpy primordial population reproduces the observations more closely than a smooth primordial population.

As noted above, in some cases we see significant variations between simulations of the same [] pairing (for example, see the uncertainties in top panel of Figure 6). Furthermore, although two simulations of the same [] pairing may both reproduce the observed surface density profile, only one may reproduce the observed radial dependence of the binary frequency, and in many cases this difference cannot easily be attributed to projection effects (at least given our 1000 sight lines). Note that we include stars of all masses in the surface density profile, but only stars within the primary-mass and mass-ratio range of Li et al. (2013b) in the radial dependence of the binary frequency, which comprise 3% of the total stellar population of the cluster at 30 Myr. This may indicate that small asymmetries in the primordial population and Poisson noise are important for determining the radial dependence of the binary frequency at a few , particularly when examining a magnitude- or mass-limited sample. The larger number of simulations beginning with smooth () initial density distributions alleviates much of this uncertainty, but additional simulations may be required to determine very robust values for [] pairings.

In summary, we confirm here the hypothesis of de Grijs et al. (2013) and Li et al. (2013b) that the observed decreasing radial trend in binary frequency towards the core of NGC 1818 can be reproduced by the early dynamical disruption of binaries by stellar encounters within the cluster. This trend can be produced within clusters born both in equilibrium or collapsing as well as with smooth or clumpy primordial density distributions, although our analysis indicates a modest preference for clumpy primordial conditions for NGC 1818.

For a cluster that is born with a binary population that has similar distributions of orbital parameters to those of the solar-type binaries in the Galactic field (with no radial dependence), the disruption of soft binaries naturally produces a decreasing trend in binary frequency towards the cluster core in roughly a crossing time. As described in Section 3 and illustrated in Figures 2 to 4, the higher density and velocity dispersion in the core result in more frequent and more energetic encounters that ionize wide binaries, while binaries of similar binding energies can survive in the less dense halo (see also Sollima, 2008). The formation of this type of early distribution relies on primordial soft binaries, and likely will not arise if, for example, the primordial binaries are cut off at the initial hard-soft boundary of the core.

At the age of NGC 1818, we find that our simulations that start with more clumpy and subvirial initial conditions (lower values of both and ) have the largest differences between the binary frequencies inside and outside of . For example, the percent difference between the mean binary frequencies inside and outside of for the [0.1,1.6] simulations is 17% at 30 Myr, while for the [0.5,3.0] simulations we find only a 7% difference at the same age. This can also be seen in the top two panels of Figure 6, where the distribution for the [0.5,3.0] model is flatter than for those of the initially subvirial and clumpy simulations. Clusters with more highly substructured and subvirial initial conditions undergo a more rapid relaxation process (e.g. McMillan et al., 2007) and, in turn, disrupt (wide) binaries more efficiently.

After roughly six initial half-mass relaxation times (about 2040 Myr) both the core and halo binary populations have experienced sufficient encounters to regain similar semi-major axis distributions. Over this time, mass segregation processes become most important in determining the radial dependence of the binary frequency in the cluster. Figure 2 illustrates this turn-over from early dynamical disruption to mass segregation.

Notably, the radial dependence of the binary frequency develops a minimum value after about 1 (340 Myr), which marches outwards towards the halo until about 4 (1360 Myr), when the outer cluster no longer maintains a second peak in binary frequency. In our simulations this bimodal radial distribution in the binary frequency manifests because of both early binary disruption and subsequent mass segregation processes. Furthermore, because the radial distribution of the binaries changes predictably with time (moving from decreasing towards the core, to a bimodal distribution, to rising towards the core), observations of the radial dependence of the binary frequency in a cluster can help to constrain the cluster’s dynamical age (e.g., the number of that the cluster has lived through).

Interestingly, the similarly aged LMC star cluster NGC 1805 shows evidence for a minimum in the radial distribution of the binary frequency near a radius of about 15–20 arcsec (about 3.6–4.9 pc; Li et al., 2013b). This may indicate that NGC 1805 is more dynamically evolved than NGC 1818 despite having the same chronological age. Indeed, the total masses of NGC 1805 and NGC 1818 of [M]) = 3.45 and 4.01 and the numbers of stars of 6000 and 9000 (Li et al., 2013b), and the core radii of 1.33 pc and 2.45 pc (Mackey & Gilmore, 2003), imply present-day half-mass relaxation times for NGC 1805 and NGC 1818 of about 115 Myr and 215 Myr, respectively. Detailed models of NGC 1805 will be important to understand if such a bimodal distribution can develop from an initial binary population born with no radial dependence in frequency (or distributions of orbital parameters and masses) after only about one fourth of a half-mass relaxation time.

Bimodal radial distributions have also been observed previously for the frequency of blue stragglers in many globular clusters. Mapelli et al. (2004) and Ferraro et al. (2012) argue that such distributions are the result of dynamical friction and mass segregation processes acting on the binary progenitors of the blue stragglers, much like we see here for the binaries in our NGC 1818 simulations. Beccari et al. (2013) find that the blue stragglers and possibly also the binaries in the globular cluster NGC 5466 show bimodal radial distributions with minimum values located at approximately the same radial distance from the cluster center, a significant step towards confirming this hypothesis. As noted by Ferraro et al. (2012), the blue stragglers can also serve as an indicator of a cluster’s dynamical age, assuming that these same mass segregation processes govern the evolution of the blue straggler radial distribution in a cluster.

Alternatively these observed bimodal blue straggler distributions may be the result of different formation processes, owing to the different densities in the different cluster regions (Ferraro et al., 1997; Li et al., 2013a). Dynamical ejections of blue stragglers from the cluster core, perhaps during their dynamical formation, may also contribute to forming such bimodal blue straggler radial distributions (Sigurdsson et al., 1994).

Finally, we note that Elson et al. (1998) found the opposite radial trend for the binary frequency of NGC 1818 in their observations of binaries with primary masses between 2 M and 5.5 M and . Specifically they found a binary frequency of 35% 5% inside the core (using a core radius of 2 pc) which decreases to 20% 5% outside of about 3 core radii. de Grijs et al. (2013) attribute this discrepant result to blending and the near-vertical morphology of the stellar main sequence within the magnitude range observed by Elson et al. (1998), which may have compromised their analysis of the observations. When we examine the observed mass range of Elson et al. (1998) in our [0.5,3.0] simulations, all have significantly lower binary frequencies (for ) in the core than the Elson et al. (1998) value of 35% 5% (and note that we began our simulations with %), which adds further support to the conclusions of de Grijs et al. (2013) regarding the analysis of the observations.

Elson et al. (1998) also ran one -body simulation to model the cluster binary population. They were unable to simulate the full cluster or the high binary frequency due to computational limitations. Instead they simulated a cluster with roughly half the expected initial mass of NGC 1818, an initial binary frequency of about 17%, and soft binaries were not included (4000 AU was the upper limit to the binary separations). They drew stellar positions and velocities from a smooth Plummer model in virial equilibrium. In their simulation, they found marginal evidence for an increase in the binary frequency within the core as compared to the halo, but they note that given the relatively small sample sizes the error bars on the binary frequency at different radial bins were sufficient to allow for a flat distribution as well. These early models showed that mass segregation could be effective in the mass range of the Elson et al. (1998) sample, potentially explaining a rise in the binary frequency towards the core. However, given the indirect nature of the models and the exclusion of the soft end of the binary distribution, the ability to reproduce the Elson et al. (1998) observations with -body models should be revisited, and particularly to understand if both the Elson et al. (1998) and Li et al. (2013b) results can be accommodated within one cohesive direct -body simulation of NGC 1818 (which we will address in a future paper).

## 6. Conclusions

NGC 1818 is a very valuable star cluster for our understanding of the early evolution of the binaries in a dense stellar environment. The cluster exhibits a rare radial dependence of the binary frequency that decreases towards the core of the cluster, whereas many other (generally older) star clusters show the opposite trend (e.g. Mathieu & Latham, 1986; Geller & Mathieu, 2012; Milone et al., 2012). Here we show, through a grid of sophisticated -body simulations, that the observed surface density profile of NGC 1818 (Mackey & Gilmore, 2003) and the radial dependence of the binary frequency (de Grijs et al., 2013; Li et al., 2013b) can be reproduced simultaneously at the cluster’s age by -body simulations with both initially smooth or substructured and equilibrium or collapsing stellar populations, with a modest preference for substructured initial conditions. The radial distribution of the binary frequency in a rich star cluster can transition smoothly over time from a uniform primordial radial distribution, to one that decreases towards the core at early ages, to one that rises towards the core at later ages. Thus both rising and falling radial distributions in binary frequency can arise naturally from the evolution of a binary population within the same rich star cluster as a consequence of both dynamical disruption and mass segregation of the binaries.

## References

- Aarseth (2003) Aarseth, S. J. 2003, Gravitational N-Body Simulations (Cambridge, UK: Cambridge University Press)
- Adams & Myers (2001) Adams, F. C., & Myers, P. C. 2001, ApJ, 553, 744
- André et al. (2007) André, P., Belloche, A., Motte, F., & Peretto, N. 2007, A&A, 472, 519
- Beccari et al. (2013) Beccari, G., Dalessandro, E., Lanzoni, B., Ferraro, F. R., Sollima, A., Bellazzini, M., & Miocchi, P. 2013, ArXiv e-prints
- Binney & Tremaine (1987) Binney, J., & Tremaine, S. 1987, Galactic dynamics (Princeton, NJ: Princeton University Press)
- Cartwright & Whitworth (2004) Cartwright, A., & Whitworth, A. P. 2004, MNRAS, 348, 589
- de Grijs et al. (2002) de Grijs, R., Gilmore, G. F., Johnson, R. A., & Mackey, A. D. 2002, MNRAS, 331, 245
- de Grijs et al. (2013) de Grijs, R., Li, C., Zheng, Y., Deng, L., Hu, Y., Kouwenhoven, M. B. N., & Wicker, J. E. 2013, ApJ, 765, 4
- Elson et al. (1987) Elson, R. A. W., Fall, S. M., & Freeman, K. C. 1987, ApJ, 323, 54
- Elson et al. (1998) Elson, R. A. W., Sigurdsson, S., Davies, M., Hurley, J., & Gilmore, G. 1998, MNRAS, 300, 857
- Ferraro et al. (2012) Ferraro, F. R., et al. 2012, Nature, 492, 393
- Ferraro et al. (1997) Ferraro, F. R., et al. 1997, A&A, 324, 915
- Geller et al. (2013) Geller, A. M., Hurley, J. R., & Mathieu, R. D. 2013, AJ, 145, 8
- Geller & Mathieu (2012) Geller, A. M., & Mathieu, R. D. 2012, AJ, 144, 54
- Geller et al. (2010) Geller, A. M., Mathieu, R. D., Braden, E. K., Meibom, S., Platais, I., & Dolan, C. J. 2010, AJ, 139, 1383
- Girichidis et al. (2012) Girichidis, P., Federrath, C., Allison, R., Banerjee, R., & Klessen, R. S. 2012, MNRAS, 420, 3264
- Goodwin & Whitworth (2004) Goodwin, S. P., & Whitworth, A. P. 2004, A&A, 413, 929
- Heggie (1975) Heggie, D. C. 1975, MNRAS, 173, 729
- Hu et al. (2010) Hu, Y., Deng, L., de Grijs, R., Liu, Q., & Goodwin, S. P. 2010, ApJ, 724, 649
- Hurley et al. (2005) Hurley, J. R., Pols, O. R., Aarseth, S. J., & Tout, C. A. 2005, MNRAS, 363, 293
- Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543
- Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897
- King (1966) King, I. R. 1966, AJ, 71, 64
- King et al. (2012) King, R. R., Parker, R. J., Patience, J., & Goodwin, S. P. 2012, MNRAS, 421, 2025
- Kouwenhoven et al. (2007) Kouwenhoven, M. B. N., Brown, A. G. A., Portegies Zwart, S. F., & Kaper, L. 2007, A&A, 474, 77
- Kouwenhoven et al. (2005) Kouwenhoven, M. B. N., Brown, A. G. A., Zinnecker, H., Kaper, L., & Portegies Zwart, S. F. 2005, A&A, 430, 137
- Kraus & Hillenbrand (2008) Kraus, A. L., & Hillenbrand, L. A. 2008, ApJ, 686, L111
- Kraus et al. (2011) Kraus, A. L., Ireland, M. J., Martinache, F., & Hillenbrand, L. A. 2011, ApJ, 731, 8
- Kroupa (1995) Kroupa, P. 1995, MNRAS, 277, 1507
- Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
- Küpper et al. (2011) Küpper, A. H. W., Maschberger, T., Kroupa, P., & Baumgardt, H. 2011, MNRAS, 417, 2300
- Lada & Lada (2003) Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57
- Larson (1995) Larson, R. B. 1995, MNRAS, 272, 213
- Leigh & Sills (2011) Leigh, N., & Sills, A. 2011, MNRAS, 410, 2370
- Li et al. (2013b) Li, C., de Grijs, R., & Deng, L. 2013b, MNRAS in press, arXiv: 1309.0929
- Li et al. (2013a) Li, C., de Grijs, R., Deng, L., & Liu, X. 2013a, ApJ, 770, L7
- Mackey & Gilmore (2003) Mackey, A. D., & Gilmore, G. F. 2003, MNRAS, 338, 85
- Mapelli et al. (2004) Mapelli, M., Sigurdsson, S., Colpi, M., Ferraro, F. R., Possenti, A., Rood, R. T., Sills, A., & Beccari, G. 2004, ApJ, 605, L29
- Mathieu & Latham (1986) Mathieu, R. D., & Latham, D. W. 1986, AJ, 92, 1364
- McMillan et al. (2007) McMillan, S. L. W., Vesperini, E., & Portegies Zwart, S. F. 2007, ApJ, 655, L45
- Milone et al. (2012) Milone, A. P., et al. 2012, A&A, 540, A16
- Peretto et al. (2006) Peretto, N., André, P., & Belloche, A. 2006, A&A, 445, 979
- Plummer (1911) Plummer, H. C. 1911, MNRAS, 71, 460
- Proszkow et al. (2009) Proszkow, E.-M., Adams, F. C., Hartmann, L. W., & Tobin, J. J. 2009, ApJ, 697, 1020
- Raghavan et al. (2010) Raghavan, D., et al. 2010, ApJS, 190, 1
- Reggiani & Meyer (2011) Reggiani, M. M., & Meyer, M. R. 2011, ApJ, 738, 60
- Sánchez & Alfaro (2009) Sánchez, N., & Alfaro, E. J. 2009, ApJ, 696, 2086
- Sigurdsson et al. (1994) Sigurdsson, S., Davies, M. B., & Bolte, M. 1994, ApJ, 431, L115
- Sollima (2008) Sollima, A. 2008, MNRAS, 388, 307
- Tobin et al. (2009) Tobin, J. J., Hartmann, L., Furesz, G., Mateo, M., & Megeath, S. T. 2009, ApJ, 697, 1103