The Dynamical Evolution of the Pleiades

We present the results of a numerical simulation of the history and future development of the Pleiades. This study builds on our previous one that established statistically the present-day structure of this system. Our simulation begins just after molecular cloud gas has been expelled by the embedded stars. We then follow, using an N body code, the stellar dynamical evolution of the cluster to the present and beyond. Our initial state is that which evolves, over the 125 Myr age of the cluster, to a configuration most closely matching the current one. We find that the original cluster, newly stripped of gas, already had a virial radius of 4 pc. This configuration was larger than most observed, embedded clusters. Over time, the cluster expanded further and the central surface density fell by about a factor of two. We attribute both effects to the liberation of energy from tightening binaries of short period. Indeed, the original binary fraction was close to unity. The ancient Pleiades also had significant mass segregation, which persists in the cluster today. In the future, the central density of the Pleiades will continue to fall. For the first few hundred Myr, the cluster as a whole will expand because of dynamical heating by binaries. The expansion process is aided by mass loss through stellar evolution, which weakens the system's gravitational binding. At later times, the Galactic tidal field begins to heavily deplete the cluster mass. It is believed that most open clusters are eventually destroyed by close passage of a giant molecular cloud. Barring that eventuality, the density falloff will continue for as long as 1 Gyr, by which time most of the cluster mass will have been tidally stripped away by the Galactic field.


Introduction
Despite recent advances in the field of star formation, the origin of open clusters remains a mystery. It is now generally accepted that all stars are born within groups. These groups are at first heavily embedded within molecular clouds, their members obscured optically by copious interstellar dust. By the time the stars are revealed, only about 10 percent are in open clusters (Miller & Scalo 1978;Adams & Myers 2001). The remainder are in either T-or OB associations, both destined to disperse within a few Myr. In contrast, the stars within open clusters are gravitationally bound to each other, and the group can survive intact for several Gyr (Friel 1995). How do molecular clouds spawn these relatively rare but stable configurations?
One intriguing aspect of the mystery is that open clusters are intermediate in their properties between T-and OB associations. The former are relatively sparse in projected stellar density, and contain up to about 100 members (e.g., Kenyon & Hartmann 1995;Luhman 2007). The latter, as exemplified by the nearby Orion Nebula Cluster, begin with extraordinarily high density (McCaughrean & Stauffer 1994) and contain well over a thousand members (Hillenbrand 1997), far more than the eponymous O and B stars. A published compilation of Galactic open clusters (Mermilliod 1995) shows them to have from a few hundred to roughly a thousand stars, i.e., just in the middle range. Apparently, systems born with either too low or too high a population and density are fragile, while the relative minority falling in between can survive.
There is already an extensive literature on young, bound clusters, both observational and theoretical (for a review, see Elmegreen et al. 2000). Models for their origin, dating back at least to Lada et al. (1984), have focused on the need for a high star formation efficiency in the parent cloud. A standard computational technique, using N-body simulations, is to create stars in a background potential well, remove that potential through various prescriptions, and then assess the result (e.g., Goodwin & Bastian 2006;Baumgardt & Kroupa 2007). Some researchers using this approach, implemented either analytically or numerically, have hypothesized that open clusters are the bound remnants of expanding OB associations (Adams 2000;Kroupa et al. 2001). In recent years, most theoretical ideas have been motivated by fluid dynamical simulations of turbulent, collapsing clouds (Klessen et al. 2000;Vázquez-Semademi et al. 2003;Li et al. 2003). While much insight has been gained from these collective investigations, there has generally been too little contact of the theory with actual groups. A clear advance would be made if we could establish empirically the original state of one or more observed clusters. We would then be in a position to gauge how these particular systems were produced by star formation activity in their parent clouds.
As a first step in this direction, Converse & Stahler (2008, hereafter Paper I) undertook a quantitative study of the present-day structure of a well-studied, relatively nearby open cluster, the Pleiades. We derived statistically, using a maximum likelihood analysis, such key properties as the stellar density distribution, mass function, overall binary fraction, and correlation between the component masses of these binaries. The point of that study was to provide the endpoint for any calculation of the system's previous evolution.
We now take the second step. The age of the Pleiades has been determined, from observations of lithium depletion, to be 125 Myr (Stauffer et al. 1998). Using the publicly available code Starlab (Portegies-Zwart et al. 2001, Appendix B), we have run a suite of N-body calculations over just this time period, to find that initial state which evolved to the current cluster, as gauged by our previous investigation. In doing so, we also establish the detailed history of the group over that epoch, and even into the future.
A key assumption here is that the Pleiades divested itself of cloud gas relatively soon after its birth. There is currently no direct means to assess the duration of the initial, embedded phase, either in the Pleiades or any other open cluster We may take a clue from T associations, which are still surrounded (but not completely obscured) by molecular gas. No systems are observed with ages exceeding about 5 Myr, a striking fact first noted by Herbig (1978). Presumably, older groups consisting of post-T Tauri stars have already driven away their clouds and are merged observationally into the field population. If a similar embedded period held for the Pleiades, it indeed represents a small fraction of the total age. Hence, we can establish, with some confidence, the cluster's structure just after cloud dispersal. A future study will investigate, using a combination of gaseous and stellar dynamics, how this early configuration itself arose.
In Section 2 below, we describe in more detail our approach to the problem. We define the parameters characterizing both the initial configuration of the cluster and the evolved system. We then outline our strategy for finding the optimal initial state, i.e., the one whose descendent matches most closely the current Pleiades. Our actual numerical results are presented in Section 3. Here, we give the detailed properties of the inferred initial state. We also describe how the cluster changed up to the present, and how it will develop in the future. One of our key findings is that the cluster's evolution did not proceed in the classic manner associated with dynamical relaxation (Binney & Tremaine 1987, Chapter 8); we explore the origin of this discrepancy. Finally, Section 4 discusses the implications of our findings on the earlier, embedded evolution of this, and other, open clusters.

Density and Velocity Distribution
As in Paper I, we model the Pleiades as a perfectly spherical system, although the cluster is observed to be slightly elongated (Raboud & Mermilliod 1998). This elongation seems to have been created by the tidal gravity of the Galaxy (Wielen 1974), which would have exerted influence throughout the cluster's dynamical history. We assume that this modest tidal stretching had negligible effect on the internal evolution, and that relatively few stars were lost by tidal stripping over the Pleiades age. Thus we can safely ignore the associated Galactic potential. Similarly, we ignore mass loss through stellar evolution, which is negligible for our adopted mass function, over the 125 Myr age of the Pleiades. In Sections 3.2 and 3.3, we present simulations that include both effects Returning to our standard runs, we further assume that the cluster was in virial equilibrium following expulsion of the gas. Any significant departure from equilibrium would be erased on a dynamical time scale, about 10 Myr for our input parameters and therefore much shorter than the evolutionary span of interest. Two popular choices for spherical equilibria are King (1966) models and polytropes. As will be explained in Section 2.3, King models do not include low enough density contrasts for a full exploration of initial states. We therefore used polytropes, which are more versatile in this regard. In polytropes, the stellar distribution function f , i.e., the number of stars per volume in configuration and velocity space, is given by Here, A is a normalization constant, while E, the relative energy per unit mass, is In this last equation, v denotes the stellar speed. Thus, E is the negative of the physical energy, and also has an offset in its zero point, as conventionally defined. This offset is embedded in the relative potential Ψ(r), which is related to the usual gravitational potential Φ(r) by The tidal radius r t marks the outer boundary reached by cluster stars. By construction, the relative potential Ψ is positive inside the cluster and falls to zero at r = r t .
The number of stars per unit volume is found by integrating the distribution function over velocity space. The manipulations here are standard (Binney & Tremaine 1987, Chapter 4), so we give only the essential results. We let m denote the stellar mass, assumed provisionally to be identical for all cluster members. Then ρ, the mass density of stars, is Here, v max (r) ≡ 2 Ψ(r) is the maximum speed for a star at radius r. For such a star, E = 0. The total physical energy per unit mass, Φ(r) + v 2 /2, is Φ(r t ), so the star can just reach r t . Using equation (2), equation (4) becomes We define a new variable θ ≡ arcsin (v/ √ 2Ψ), so that ρ(r) = 2 5/2 π m A Ψ n n π/2 0 cos 2n−2 θ dθ (6) = (2 π) 3/2 m A Γ(n − 1/2) Γ(n + 1) Ψ n (r) .
To calculate the relative potential Ψ(r), we use Poisson's equation. For our spherical system, this is 1 where ρ 0 and Ψ 0 are the central values of ρ(r) and Ψ(r), respectively. We define a dimensionless potential as ψ ≡ Ψ/Ψ 0 , and a dimensionless radius as ξ ≡ r/r 0 , where the scale radius r 0 is Since ρ = ρ o ψ n , the new potential obeys the Lane-Emden equation: with boundary conditions ψ(0) = 1 and ψ ′ (0) = 0. The nondimensional tidal radius ξ t ≡ r t /r 0 is the point where ψ falls to zero. 1 For any chosen polytropic index n, equation (10) can readily be solved numerically. Our task is to translate this nondimensional solution into a physical model of the initial cluster. Given n, the basic quantities characterizing the cluster are: N tot , the total number of stars; r v , the virial radius; and m, the mean stellar mass. 2 The virial radius is defined by Here, M ≡ N tot m is the total cluster mass, while W is the gravitational potential energy: In Appendix A, we show how to obtain the dimensional quantities r 0 , ρ 0 , and Ψ 0 from our three input parameters and the solution ψ(ξ).
Given the scale factors r 0 and ρ 0 , we know the dimensional mass density ρ(r). We then populate space with stars according to a normalized distribution p 1 (r) such that p 1 (r) dr is the probability a star is located between r and r + dr. This probability is simply The actual position vector of each star is then distributed isotropically within each radial shell.
Finally, we require an analogous distribution for stellar speeds. At a given radius, the speed must be consistent with the prescribed energy distribution. Let p 2 (v|r) dv be the probability that the speed lies between v and v + dv given that its radius is r. Clearly, Replacing d 3 r by 4 π r 2 dr and d 3 v by 4 π v 2 dv, we have, after using equation (13), We take f = f (E) from equation (1) and use the definition of E from equation (2), finding Here, the relative potential is calculated at each r from Ψ = Ψ 0 ψ(ξ), where we recall that ξ is the nondimensional radius. Given the stellar speed, i.e., the magnitude of the velocity vector, the direction of that vector is again distributed isotropically in space.

Stellar Masses: Single and Binary
Thus far, we have described a cluster that is composed of members with identical mass. In actual practice, we assign masses to the stars according to a realistic distribution. The parameters of this mass function for the initial cluster are among those we vary to obtain an optimal match between the evolved system and the present-day Pleiades. In the course of evolution, some stars will be given enough energy, through three-body interactions, to escape the cluster. The most massive ones die out over 125 Myr. It is therefore not obvious that the initial mass function is identical to that found today.
We suppose that the distribution of stellar masses in the young Pleiades was similar in form to the initial mass function for the field population. In recent years, large-scale surveys of lowluminosity objects, combined with spectroscopy, have established an accurate initial mass function down to the brown dwarf limit (e.g., Covey et al. 2008). The consensus is that the original power law of Salpeter (1955) for masses above solar is joined at the lower end by a lognormal function. This basic form appears to hold in diverse environments, including young clusters (Chabrier 2005).
Let φ(m) dm be the probability that a star's mass (in solar units) is between m and m + dm. We posit that this probability is where B, C, α, and the joining mass µ are all constants. We set m min = 0.08 and m max = 10 (although we also tested a higher mass limit; see §3.2 below). The variable y is given by Here, m 0 is the centroid of the lognormal function, and σ m its width. For input parameters α, m 0 , and σ m , the constants B, C, and µ are determined by the normalization condition and by requiring that φ(m) and its first derivative be continuous at m = µ. Analytic expressions may be found for the three constants, which we do not display here.
Most stars are not single objects, but have binary companions. Indeed, the Pleiades today is especially rich in binaries (see §3.4 of Paper I). Such pairing must have been present in the initial cluster. We therefore view the positions and velocities established in the previous subsection as pertaining to N tot stellar systems, rather than individual stars. Similarly, the symbol m used, e.g., in equation (13), actually denotes the average system mass, after accounting for binaries. We specify the global binary fraction as a parameter b, which gives the probability that a system actually consists of two stars. Conversely, a fraction 1 − b of the systems are indeed single stars. Their mass is distributed according to the probability φ(m).
Our analysis of the present Pleiades in Paper I showed that the masses of the component stars within binaries are correlated. Such a correlation must also have been present at early times. Accordingly, we include the effect in our initial state. Within the fraction b of systems that are binaries, we first independently assign masses to each component, using the probability distribution φ(m). After identifying the primary mass m p and the secondary mass m s , we then alter the latter Here, γ is an input parameter that measures the degree of mass correlation within binaries (see also eq. (42) of Paper I). If γ = 0, the component masses are uncorrelated, while γ = 1 corresponds to perfect matching.
We give our binaries randomly inclined orbital planes, and a period and eccentricity distribution characteristic of present, solar-type binaries both in the field (Duquennoy & Mayor 1991) and in the Pleiades itself (Bouvier et al. 1997). If p 3 (P) dP is the fraction of systems with periods between P and P + dP, then p 3 (P) is lognormal: where z ≡ log P − log P 0 √ 2 σ P .
We set the centroid period to log P 0 = 4.8 and the width to σ P = 2.3, where the period is measured in days. The eccentricity distribution has a thermal distribution: as motivated both by observations (Duquennoy & Mayor 1991) and theory (Heggie 1975).
The initial cluster described thus far is homogenous, in the sense that any volume containing an appreciable number of systems has the same average system mass. However, there have long been claims of observed mass segregation in young clusters, i.e., an increase of average stellar mass toward the center (Sagar et al. 1988;Jones & Stauffer 1991;Moitinko et al. 1997;Stolte et al. 2006). The present-day Pleiades also exhibits this phenomenon to a striking degree (see §4.2 of Paper I). We want to see if this property developed on its own or was inherited from an earlier epoch. We accordingly include a quantitative prescription for mass segregation in our initial state.
One system has a higher probability of being near the cluster center than another, in a timeaveraged sense, if its relative energy E is greater. Mass segregation therefore manifests itself as a correlation between the system mass m and E. This fact was noted by Baumgardt et al. (2008), who used it to implement a specific procedure for mass segregation. Here we have adopted a variant of their method that allows us to include the effect to a variable degree. We first assign E-and m-values to all member systems according to equation (1) and our prescription for binary masses. We then place the systems in two lists -the first ordered by increasing m, and the second by increasing E. When we first construct these lists, the ranking of the system in the first is unrelated to its ranking in the second. This is the case of zero mass segregation. There would be perfect mass segregation if the two rankings were identical.
Let us quantify the intermediate case. For a star of given mass, we find its index in the massordered list. To assign an energy to that star, we choose the second (energy-ordered) index from a Gaussian distribution centered on the mass index. The width of this distribution, denoted σ E , can be infinite (no mass segregation) or zero (perfect segregation). More generally, we define a parameter β, the degree of mass segregation, which varies between 0 and 1. After some trial and error, we adopted the following prescription relating the width σ E to β: The logarithmic dependence on β ensures that σ E has the desired behavior in the extreme limits. The proportionality with N tot ensures that our algorithm gives the same degree of biasing in clusters of any population.
In summary, β becomes another input parameter that we vary within the initial configuration. As we will see, having a non-zero β is critical to obtaining a proper match between the evolved cluster and the Pleiades today. Mass segregation was therefore present at a relatively early epoch.
Since relatively massive stars preferentially reside near the center, our imposition of mass segregation alters the shape of the gravitational potential from that of a single-mass polytrope.
Relative to the total gravitational energy, the total kinetic energy has a value slightly below that for virial equilibrium. We rescaled all stellar velocities by a uniform factor to restore exact equilibrium. In practice, this factor was typically about 1.05.
For convenient reference, Table 1 lists the full set of our input parameters for the starting state. Anticipating the results detailed in §3 below, the last column gives the numerical value of each parameter in the optimal configuration. We also list the associated uncertainties. The meaning of these uncertainties and how they were assessed, will also be discussed presently.

Characterizing the Evolved Cluster
After evolving a particular initial state for 125 Myr, we compare the outcome with the actual Pleiades. In making this comparison, it is important to "observe" the simulated cluster under the same conditions as the real one. Thus, we project the three-dimensional distribution of stars onto a two-dimensional plane, assumed to lie at the mean Pleiades distance of 133 pc (Soderblom et al. 2005). The angular separation ∆θ between each pair of stars is then determined. If ∆θ res denotes the telescope resolution, then any pair with ∆θ < ∆θ res is taken to be an unresolved point source. For the near-infrared catalog of the Pleiades analyzed in Paper I (Stauffer et al. 2007), an appropriate value of ∆θ res is 10 arcsec. Note that our unresolved sources include a small fraction (less than 0.5 percent) of triples and high-order systems, as well as a few unrelated pairs observed to be close in projection. We denote as N s the total number of point sources out to a radius from the cluster center of 12.3 pc, corresponding in angle to 5.3 • . This was the radius enclosing the catalog of sources used in Paper 1. For each simulated evolution, we compare the final N s -value with the observed Pleiades figure.
The vast majority of stars observed in the Pleiades today are on the main sequence (see Fig. 1 of Paper 1). The number of post-main-sequence objects, while relatively small, is sensitive to the shape of the stellar mass function. Hence, it is important that we reproduce, as closely as possible, the number inferred for the present-day cluster. At an age of 125 Myr, the main-sequence turnoff is about 4 M ⊙ . If N 4 denotes the number of stars (singles or primary stars in binaries) whose mass exceeds 4 M ⊙ , then this quantity, as calculated, may also be compared directly with that in the Pleiades. Similarly, we compare M tot , the total mass of all stars in the evolved cluster, with the Pleiades mass obtained through the statistical analysis of Paper 1.
One striking result of Paper 1 was the prevalence of binaries. Specifically, we found that the near-infrared fluxes of the catalogued point sources demanded that the fraction b unres = 0.68 were unresolved binaries. (Any resolved binaries were listed in the catalog as separate sources.) For our assumed resolution limit ∆θ res , we could also assess b unres computationally for each evolutionary run. Again, this is the fraction of point sources representing two or more unresolved stars. Note that b unres is less than the initially imposed binary fraction b, both because some pairs are wide enough to be resolved, and because others are torn apart in the course of evolution.
The distribution of stellar masses is, of course, another property that should be compared with the actual Pleiades. As just described, we set the form of the distribution within the initial configuration as a lognormal function with a power-law tail. The apportionment of masses within the evolved state could in principle differ, due to the escape of some stars from the cluster and the death of others with sufficiently high mass. Our procedure is first to find, within the output state, the normalized distribution of single stars. Included in this distribution are both isolated stars and the components of resolved binaries. We then peer within unresolved binaries and find the analogous distributions of primary mass m p and secondary mass m s . Finally, we record the distribution of the binary mass ratio, q ≡ m s /m p .
In Paper I, we statistically determined the stellar masses from the photometric data by assuming that the single-star distribution was a pure lognormal, with centroid m 0 and width σ m . For consistency, we characterize the evolved cluster in our simulations in a similar fashion. Now for a given single-star function and binary correlation parameter γ, the primary, secondary, and q-distributions are all uniquely determined. Appendix B outlines the mathematical derivation. The task is to vary γ, as well as m 0 and σ m , for the presumed lognormal single-star distribution until this function, as well as the primary, secondary, and q-distributions, best fit those we find directly in the numerical output. We then compare γ, m 0 , and σ m to these same quantities derived in a similar way for the observed Pleiades.
We next consider the projected density profile. We divide the cluster into radial bins that match those used in the analysis of the Pleiades. The resulting surface density of stellar systems is then fit to the empirical prescription of King (1962): Here, R is the projected radius, k is a constant with the dimensions of a surface density, and R c and R t are the core and tidal radii, respectively. We determine the values of k, R c and R t which best match the data, i.e., the same parameters determined for the real Pleiades by an analogous fitting procedure. However, only R c is used in the final optimization routine (see below). We also determine the King concentration parameter c K ≡ log (R t /R c ), both for each evolved simulation and in the real cluster. From equation (25), the central surface density Σ 0 is This is also compared to the Pleiades value.
Finally, we measure the degree of mass segregation. For the evolved cluster, we compute the cumulative fraction of systems contained within a projected radius R, both by number (f N (R)) and mass (f M (R)). As in Paper I ( §4.2), the Gini coefficient is computed as and then compared to that found in the Pleiades. Table 2 gives the full list of quantities evaluated for each evolved cluster. We also display the values found when the optimal initial state is used, as well as the corresponding figures in the actual Pleiades. Notice that m 0 and σ m for the mass function do not match those in the initial state, as given in Table 1. This discrepancy arises partly from real changes of the stellar masses, but even more from our adoption of a simple lognormal when fitting the evolved cluster. The tabulated errors for the calculated quantities were obtained by running 25 simulations, all with identical input parameters. Since we populated the cluster stochastically, according to probability distributions (e.g., φ(m) in eq. (17)), initial states differed from one another in detail. The errors represent the standard deviations for each quantity in the evolved cluster, due solely to differing realizations of the initial state. The tabulated errors for the observed Pleiades are from the calculation of Paper I. In addition, N s and N 4 are assumed to be Poisson-distributed, so that the errors are √ N s and √ N 4 , respectively.

Optimization Procedure
As a first guess, we set all the input parameters equal to values appropriate for the Pleiades today. In Paper I, we characterized the present-day cluster as a King model, so we initially adopted this prescription. The analytic models of King (1966) have the distribution function Here Ψ 0 , the central value of the relative potential, is again set by our basic input quantities. The dimensionless parameter W 0 , like c K , characterizes the degree of central concentration. 3 Running a King model for 125 Myr, we found that it invariably became more centrally concentrated than the actual Pleiades. We therefore tried successively lower W 0 -values for the initial state. Now equation (28) shows that f (E) is proportional to E in the limit of small W 0 . Comparison with equation (1) reveals that such a model is equivalent to a polytrope of n = 5/2. Our search for low-concentration initial states therefore led us naturally to the polytropic models described in §2.1. Within the regime of polytropes with small n, we varied other input parameters, such as r v , as necessary.
Once our computed cluster began to resemble the Pleiades, we changed to a more systematic gradient method for refining the initial state. Let x be the vector whose elements are the 9 input parameters listed in Table 1. Similarly, let y represent the 11 evolved cluster properties of Table 2. This latter vector is, of course, a function of x. To move y toward the values characterizing today's Pleiades, we need to evaluate, in some sense, the gradient of this function.
A practical complication is one to which we alluded earlier. Even among evolutionary runs assuming an identical input vector x, the resulting y differs because of the stochastic sampling of the various assumed distribution functions. In computing the gradient, we need to take a step size h large enough that the resulting change in y exceeds that due to this realization variance. We found that the prescription h = 0.5 x sufficed for this purpose (see §5.7 of Press et al. 2002, for a more rigorous justification).
For each x, we first do 9 runs and average the result to obtain y(x). We then decrease, in turn, each element x j to x j − h j , and find the average output of two runs at each decreased x j -value. Similarly, we find the average result of two runs at each x j + h j . We thus establish the 11 × 9 matrix of derivatives D, whose elements are The change in outputs for any subsequent input change ∆x may then be approximated by Here, the vector ∆y is taken to be the difference between the current y-vector and that for the Pleiades. We may evaluate the 9 elements ∆x j by solving the 11 linear equations summarized in (29). Since the system is overdetermined, we did a least-squares fit to find that set of ∆x j which best satisfied the equations.
As we took a step in x, we evaluated how close the resulting y was to y p , the aggregate properties of the observed Pleiades. We did a χ 2 -test, where Here, each < y i > is the average y i value, established by doing 9 runs with identical input values. The standard deviation σ i includes errors in both the inferred Pleiades properties and those generated by different statistical realizations of the input state: The first righthand term is the Pleiades variance whose square root is the error given in the last column of Table 2. The quantity σ 2 <y i > is the error in the mean y i . This error in the mean is related to σ y,i , the variance in each individual y i , by For the first few x-steps, χ 2 declined, but then stalled. Beyond this point, the gradient method itself was clearly failing, as it indicated initial states which evolved to configurations less resembling the Pleiades. The difficulty was that the numerical derivatives of equation (29) were too crude to refine the initial state further. Refinements are possible in principle, but prohibitive computationally. After pushing the method to its limit, we were forced to stop the search before χ 2 reached a true minimum. We took the last state in the sequence where χ 2 declined to be the best-fit initial configuration.
Our final task was to assess the errors in all input parameters for this state. These should reflect uncertainties in properties of the actual Pleiades, as well as the variation in output parameters among different runs using identical inputs. This latter effect is quantified by the covariance matrix Y, whose elements are The averaging here refers to different realizations using identical input parameters. Standard error propagation (Cowan 1998, Section 1.6) dictates that the known Y is related to X , the desired covariance matrix of input parameters, through the derivative matrix and its transpose: We need to invert this equation to obtain X . As noted, the input errors should also reflect the observational uncertainties in the Pleiades itself. We do not know the correlation of these observational uncertainties. Thus, we use on the lefthand side of equation (35) a matrix Y ′ , formed by adding σ 2 p,i to each diagonal element Y ii .
Since D is not a square matrix, a standard inverse cannot be defined. However, the product D T D is square, and so has an inverse, provided it is not singular. As discussed in Graybill (1983), this fact allows us to define the pseudo-inverse of D: The term "pseudo-inverse" is appropriate since where I is the identity matrix. Taking the transpose of this last equation, we also find By employing equations (37) and (38), inversion of the modified equation (35) is straightforward: The errors in the initial cluster parameters of Table 1 are then the standard deviations obtained from the diagonal elements of X . Table 1 lists the optimal values for the parameters characterizing the initial state. The polytropic index n is about 3, corresponding to a volumetric, center-to-average, number density contrast of 54. 4 This particular polytrope closely resembles a King (1966) model with W 0 ≈ 1.4. Note the relatively large uncertainty in the optimal n, reflecting the fact that a range in initial density contrasts relaxes to a similar state after 125 Myr. There is much less uncertainty in the virial radius r v , which is surprisingly large compared to observed embedded clusters (see §4 below). Smaller assumed r v -values, however, evolved to systems with too high a density contrast. Figure 1 shows, as the dashed curve, the initial surface density as a function of projected radius. Also plotted (solid curve) is the evolved surface density, along with observed data from the Pleiades. Notice how the surface density decreases with time, a result of the inflation experienced by the entire cluster. This behavior contrasts with expectations from the standard acount of dynamical relaxation (e.g. Binney & Tremaine 1987, Chapter 8). The swelling of the central region that we find is consistent, however, with previous simulations of binary-rich clusters with relatively low populations (Portegies-Zwart et al. 2001). We explore further the underlying physical mechanism in §3.3 below.

Global Properties of the Cluster
Note from Table 1 that N tot , the initial number of stellar systems, is determined to within about 5 percent uncertainty. The main constraint here is the need to match N s , the final, observed number of point sources. Note also that N tot < N s throughout the evolution. Almost all the stellar systems are binaries. Some of these are wide enough that they could be resolved observationally. Thus, the total number of point-like (i.e., unresolved) sources is always higher than N tot , the number of stellar systems (resolved or unresolved). By the same token, the unresolved binary fraction, b unres = 0.68, is significantly less than the full initial binary fraction, b = 0.95. Indeed, we were forced to choose a b-value close to unity in order to make b unres close to the observationally inferred figure (see Table 2). Figure 2 quantifies the degree of mass segregation in the evolved cluster, Following the technique introduced in Paper I, we plot f M , the fractional cumulative mass at any projected radius, against f N , the fractional cumulative number. The fact that this curve rises above the dashed diagonal (f M = f N ) indicates the existence of mass segregation. The empirical f M − f N relation for the Pleiades, shown by the points with error bars, is well matched by the simulation. We were able to obtain this match only by adopting a non-zero value of β, the initial degree of mass segregation defined in equation (24). 5 Table 1, m 0 , σ m , α, and γ, concern the mass function. The number of stars escaping the cluster during its evolution is relatively small (on average, 280 of the 2400 stars present initially). Because of this small loss, and because few members evolve off the main sequence, the initial and final mass functions are essentially identical, and all four parameters are highly constrained by the observations. Note, in particular, that the exponent α of the power-law tail directly influences N 4 , the observed number of massive stars. The binary correlation parameter γ is independent of the single-star mass function, but influences the primary, secondary, and qdistributions, as described in §2.2. Any substantial variation in γ would alter the corresponding parameter obtained statistically for the observed cluster (see Fig. 5 of Paper I).  Table 2. The data points, along with error bars, represent the inferred single-star mass function for the Pleiades, obtained through the maximum likelihood analysis of Paper I. The agreement with the simulation is naturally poorest at the highest masses, since we modeled the output as a pure lognormal, in order to be consistent with the procedure adopted in Paper I.

Four quantities in
Finally, Figure 4 shows the initial distribution of the binary mass ratio q. We see how nearly equal-mass systems are strongly favored for our best-fit γ of 0.73. In the simulations, this distribution evolves almost unchanged, and closely matches the one inferred for the Pleiades today. The figure also displays the q-distribution for the hypothetical case of γ = 0. Such random pairing of stellar masses does not result in a flat curve, as one might expect. Instead, it reflects the character of the single-star mass function, which here is lognormal. As seen in Figure 4, the q-distribution for γ = 0 peaks at q = 0.34 and still vanishes as q approaches 0.

Past Evolution
We can now describe, based on our suite of simulations, the evolution of the cluster from its initial state to the present epoch. The main trend is an overall expansion of the system. This tendency is clear in Figure 5, which shows the variation in time of the virial radius, r v . After an initial drop, lasting about two crossing times (t cross = 10 Myr) the radius steadily swells, increasing by about 40 percent to the present. From the definition of r v in equation (11), we infer that the gravitational potential energy W is decreasing in absolute magnitude, i.e., the cluster is gradually becoming less bound. Note that we do not obtain r v by calculating W directly, but through fitting the cluster at each time to a King model, and then finding the appropriate r v for the best-fit model parameters. Figure 5 shows that R c , the projected core radius, displays similar behavior to r v . After the transient phase which again lasts about two crossing times, R c also swells, albeit more slowly. Analogous early adjustments are evident in other global quantities (see . This transient results from our implementation of mass segregation, which alters slightly the gravitational potential (recall Section 2.1.2). Although the initial cluster is in virial equilibrium, the stellar distribution function is no longer a steady-state solution to the collsionless Boltzmann equation. Within the first two crossing times, the distribution readjusts to become such a solution. The core radius bounces, before settling to a value that subsequently evolves more gradually.
Expansion of a cluster's outer halo is one manifestation of dynamical relaxation. However, application of equation (4-9) of Binney & Tremaine (1987), with ln Λ = ln (0.4 N ), reveals that the relaxation time is 250 Myr, or about twice the Pleiades age. In addition, the inner cores of relaxing systems shrink, giving energy to the halo. The secular expansion of R c further indicates that we are not witnessing the usual effects of dynamical relaxation. Figure 6 provides yet another illustration of this point. Here, we see that the King concentration parameter c K remains virtually constant, again following an initial adjustment. Recall that c K ≡ log (R t /R c ), where R t is the projected tidal radius. Thus, R t and R c swell at about the same pace.
The projected surface number density, Σ(R), currently peaks strongly at R = 0 (Fig. 1). The actually central value, however, previously declined from an even higher level. Figure 7 shows this gradual decline, which is consistent with the previously noted rise in R c . Thus, R c increases from 1.6 to 2.2 pc over the period from t = 30 Myr to t = 125 Myr. Over the same interval, Σ 0 falls by a factor of 0.50, which is close to (1.6/2.2) 2 . The number of systems in the core therefore remains virtually constant as the core itself expands. The volumetric number density similarly falls in the central region.
In Paper I, we documented a strong degree of mass segregation in the current Pleiades, quantifying this property through the Gini coefficient. Another result of the current study is that G did not attain its current value through purely stellar dynamical evolution. As seen in the top curve of Figure 8, G(t) rose only slightly at first, and then remained nearly constant, even declining somewhat in the recent past. Initial states in which the parameter β was too low never attained the requisite degree of mass segregation. As an illustration, Figure 8 shows also the result from a single simulation using β = 0 initially. The Gini coefficient does grow, but not by enough to match observations. We note, parenthetically, that G(t) exhibits oscillatory behavior over the a period that roughly matches the crossing time. These oscillations (unlike the initial readjustment) were washed out in the averaging procedure that produced the top curve in the figure. Finally, we remark that G(t) appears to saturate in time. We will return to this interesting phenomenon shortly.
All the simulations we have described thus far ignored any effects of stellar evolution. We could afford this simplification because of the relatively small number of cluster members that would have evolved significantly over 125 Myr. However, the code Starlab does have the capability of tracking stellar evolution, including mass loss, through fitting formulae. As a check, we retained our usual maximum mass of 10 M ⊙ and ran 25 simulations using the best-fit initial cluster parameters, but with stellar evolution included. The mass loss from relatively massive cluster members did not have a significant dynamical effect, and the endstate of the cluster was essentially identical. With reference to Table 2, the only parameter that changed appreciably was N 4 , which fell.
In more detail, the few stars above 7 M ⊙ usually evolved to white dwarfs of approximately solar mass. On average, about 10 white dwarfs formed, of which 7 were the secondaries within binaries. Even the few that were primaries were faint relative to their main-sequence companions, and thus would be difficult to detect. Our findings are thus consistent with the observation of Fellhauer et al. (2003) that white dwarfs are generally rarer in open clusters than might be expected statistically from the initial mass function.
Finally, we relaxed the upper mass limit in the single-star mass function and allowed the maximum mass to be arbitrarily large, according to the power law in equation (17). Choosing stars stochastically from this distribution yielded a few members with masses as high as 40 M ⊙ . If we again allowed for stellar evolution and used our standard initial cluster parameters, the evolution did take a different turn. The very massive stars represented a significant fraction of the total cluster mass, and their death had a quantitative impact. As before, the cluster went through an initial adjustment, partially from the heavier stellar mass loss. The system then smoothly expanded, but at a faster pace. At 125 Myr, the projected core radius R c was 3.1 pc, or 1.5 times larger than that of the present-day Pleiades. Similarly, the central surface density Σ 0 was a factor 0.58 lower. Had we begun with very massive stars in an initial state a factor of 1.5 smaller than our standard one, a closer match would have resulted.
These results were instructive, if somewhat academic. In reality, stars more massive than about 10 M ⊙ would have inflated HII regions so quickly as to ionize and disperse the parent molecular cloud forming the Pleiades. In order to retain even a remnant, gravitationally bound cluster, the initial membership must have been very large, about 10,000 stars in the simulations of Kroupa et al. (2001). We stress that even this figure is a lower bound, as Kroupa et al. (2001) assumed a star formation efficiency in the parent cloud of 33 percent by mass. Such an efficiency is plausible within individual dense cores (Alves, Lombardi, & Lada 2007), but significantly higher than observational and theoretical estimates in cluster-forming clouds (e.g. Duerr, Imhoff, & Lada 1982;Huff & Stahler 2007).
Suppose we nevertheless adopt this scenario as a limiting case, and assume provisionally that the Pleiades progenitor contained at least 10,000 individual stars. Such groups are rare. Equation (39) of McKee & Williams (1997) gives the birthrate of OB associations based on their population of supernova progenitors (m > 8). If we use our adopted initial mass function to estimate this population, then the birthrate of relevant OB associations is 0.09 Myr −1 kpc −2 . This is a factor between 5 and 8 smaller than the total formation rate of open clusters (Adams & Myers 2001;Miller & Scalo 1978). It is unlikely, therefore, that formation through dispersing OB associations dominates, and we continue to use an upper mass limit of about 10 M ⊙ for the Pleiades.

Future Evolution
The same calculations that reconstruct the past history of the Pleiades may also be used to predict its development far into the future. It is still believed, following the original proposal by Spitzer (1958), that most open clusters are eventually destroyed by the tidal gravitational field of passing interstellar clouds, now identified as giant molecular complexes. In this project, we do not attempt to model encounters with such external bodies. However, Starlab can follow the effects of the Galactic tidal field, both through imposition of the appropriate external potential and by adding a Coriolis force to individual systems. We switched on the Galactic field, in addition to stellar evolution, and followed the cluster from its initial state for a total of 1 Gyr. While most open clusters do not survive this long, some do last up to several Gyr (Friel 1995). Our simulation thus models at least a portion of the Pleiades' future evolution. 6 Up to the present cluster age of 125 Myr, adding the Galactic tidal field and stellar mass loss made very little difference in the evolution. Beyond this point, the cluster will continue the overall expansion that characterized it in the past. As seen in Figure 9, the central density Σ • keeps declining. The falloff is roughly exponential, with an e-folding time of 400 Myr. This figure, along with Figures 10 and 11, show average results from the 4 runs we conducted. Even after averaging, the calculated Σ • displays increasing scatter for t > 700 Myr. By this time, the total population has also fallen to the point that numerical determination of Σ • (through a fitted King model) becomes problematic.
The decline in the cluster population, which was modest until the present, accelerates as stars are tidally stripped by the Galactic gravitational field. It is the lighter stars that populate the cluster's outer halo and that preferentially escape. Consequently, the average mass of the remaining cluster members rises. Figure 10 shows both these trends. Displayed here is N s , the number of systems contained within the initial Jacobi radius of 14.3 pc. 7 By 1 Gyr, the total membership has fallen to a few dozen systems. Meanwhile, the average system mass m rises, almost doubling by the end. Careful inspection of Figure 10 shows that m initially fell slightly to its present-day value. This falloff reflects the loss, through stellar evolution, of the most massive members, an effect which is eventually overwhelmed by the escape of the lightest systems. Notice again the jitter in the m curve at later times.
Despite this qualitative change in the cluster's internal constitution, the degree of mass segregation remains essentially constant until very late times. Figure 11 displays the Gini coefficient. In detail, G(t) exhibits oscillations qualitatively similar to what we saw in lower portion of Figure 8. Nevertheless, its average magnitude does not appreciable change until t ∼ 700 Myr, when it be-gins a steep descent. The very large scatter during this late epoch again reflects the diminishing population.
What accounts for these trends? During most of the evolution, some process is able to enforce mass segregation, despite the continual depletion of the lightest members. This process evidently loses its efficacy at late times, when the total population falls too low. Earlier, we showed that dynamical relaxation did not establish the present-day level of mass segregation. Although we are now spanning a period well in excess of the initial relaxation time (250 Myr), we still do not see the classic behavior -monotonic shrinking of the core that feeds halo expansion. Up to t ∼ 200 Myr, the projected core radius R c continues the increase noted earlier. Between 200 and 600 Myr, when 75 percent of cluster members escape, R c does decline slightly, from 2.2 to 1.2 pc. Thereafter, the core swells once more.
We believe that the system's overall expansion is due principally to the release of energy in three-body encounters, specifically, close passages of binaries and single stars. Over the Gyr time span, mass loss through stellar evolution and tidal stripping weakens the cluster's gravitational binding, rendering it increasingly responsive to such internal heating. We ascribe the maintenance of mass segregation, i.e., the inward drift of more massive stars, to dynamical friction with the background population. We shall revisit these key processes, binary heating and dynamical friction, momentarily. Figure 12 shows graphically how the cluster will appear far in the future. Here we show positions of the member systems projected into the Galactic plane, both at the present time and at t = 700 Myr. One sees at present a slight elongation along the direction toward the Galactic Center. This tidal stretching is well documented observationally (Raboud & Mermilliod 1998). Stars that leave tend to do so along that direction. But any appreciable excursion leads, because of the Galaxy's differential rotation, to a change in angular speed. As a result, two tidal streams develop that are orthogonal to the Galactocentric radius. These streams are present in both panels of Figure 12, but are especially noticeable in the diminished cluster shown at the right.
Since we suspected that binaries were important in the gross dynamics of the cluster, we recalculated the entire 1 Gyr evolution after effectively removing all primordial binaries from the system. We began with the same, best-fit initial state as before, which had the usual distribution of single-star masses and initial binary fraction b = 0.95. However, we replaced every binary by a single star, located at the system's center of mass and comprising the total of the component masses. In addition, we turned off both stellar evolution and the Galactic tidal field, in order to explore the evolution under the simplest conditions possible. Note that our procedure for fusing binaries into single stars preserved the total cluster mass, number of stellar systems, and mass distribution of those systems. In other words, all two-body interactions between cluster members were the same as before. The important difference is that we eliminated the source or sink of energy associated with the internal motion of binaries.
The results were both surprising and illuminating. The cluster still undergoes overall expan-sion. Figure 13 shows that the central density again falls steadily. The nominal e-folding time is again 400 Myr, but the actual decline is not well fit by an exponential. The root cause of the cluster expansion is that new binaries continually form and interact with other stars. This process occurs principally near the relatively dense cluster center, where the most massive stars reside, along with other, more representative members. The component masses in the new binaries are high, typically 8 m . Prompt formation of binaries is a well-documented occurrence in systems initially containing only single stars (Aarseth 1971), and the formation rate is greatly enhanced at higher stellar mass (Heggie 1975). In our simulations, only 3 or 4 of these systems exist at any time. Nevertheless, they are significant dynamically, because of the cluster's relatively low gravitational binding.
Any such massive binary with a separation less than 8 × 10 4 AU = 0.4 pc is hard, i.e., has a gravitational potential energy exceeding the initial mean kinetic energy of all cluster members. Thus, even the relatively wide binaries formed in these simulations, with initial separations of order 10 3 AU, are capable of heating the cluster dynamically. As has long been appreciated (Heggie 1975;Hut 1983), the encounter of a hard binary with a third star usually results in a harder (tightened) binary. Both this pair and the isolated star have more translational kinetic energy than before. The extra energy, which comes at the expense of the binary's tightening, is quickly transferred to other cluster members.
The same dynamical heating operates, of course, in all stellar groups containing binaries. However, very populous systems, such as globular clusters, have such high gravitational binding that almost all newly formed binaries are soft. In this case, energy exchange via three-body encounters has a minor effect, and the classical picture of dynamical relaxation via two-body encounters applies. In relatively sparse systems like open clusters, both primordial and dynamically formed binaries inject so much energy that they impulsively change the velocity distribution function and qualitatively influence the course of evolution. This stochastic resetting of the velocities, which was emphasized in the classic study of Terlevich (1987), is a conspicuous feature of the Pleiades evolution, both past and future.

Discussion
While undertaken primarily to reconstruct the history of the Pleiades, our study has shed light on a well-documented, but still poorly understood, feature of open clusters generally -mass segregation. We demonstrated that the currrent, rather high degree of segregation in the Pleiades could not have been the result of dynamical relaxation from a pristine state with homogeneous mass distribution. First, the cluster has only been evolving for about half its initial relaxation time. Second, a hypothetical cluster starting with no mass segregation cannot reach the present level. Quantitatively, the Gini coefficient rises, but not enough (recall Figure 9). Two conclusions may be drawn. The ancient Pleiades must have already had substantial mass segregation before it drove off the gas. Some other process, unrelated to dynamical relaxation, must drive this effect, and continues to do so long into the future (Figure 11). The most obvious candidate is dynamical friction. A relatively massive star moving through a lower-mass population, experiences a drag force, causing it sink toward the cluster center. The associated time scale for braking, t DF , can be substantially smaller than the dynamical relaxation time t relax (Spitzer 1969). According to Portegies-Zwart & McMillan (2002), the quantitative relation is for a heavy star of mass m in a background of average mass m . (2002) and other researchers have invoked dynamical friction to explain mass segregation, focusing on very populous clusters in which massive star infall leads to the runaway growth of a central black hole (see also Gürken, Freitag, & Rasio 2004). Our work reveals a further, curious aspect of the phenomenon. Figures 9 and 11 suggest, and further calculations confirm, that G(t) saturates, regardless of its initial value. Why does the degree of mass segregation level off? A possible explanation is that, as the most massive stars sink to the center, the population there becomes increasingly homogeneous. Since m /m rises, so does t DF . In other words, mass segregation through dynamical friction may be a self-limiting process.

Portegies-Zwart & McMillan
Returning to the prehistory of the Pleiades, another significant finding is the relatively large size of the initial state. The virial radius r v began at 4 pc, while the projected half-mass radius of the initial cluster was about 2 pc. For comparison, the observed half-light radii of embedded clusters, as seen in the near infrared, range from about 0.5 to 1.0 pc, with some outliers on either side (Lada & Lada 2003). Thus, the initial, gas-free Pleiades had a radius 2 to 4 times larger than typical embedded systems. It may plausibly be argued that the Pleiades is an especially populous open cluster, and therefore began as a larger configuration, far outside the typical range. With this caveat in mind, our result suggests that the system expanded during its earliest, embedded phase. This swelling, which was accompanied, or even preceded, by mass segregation, could have been due to the loss of ambient gas during the formation process. Interestingly, observations of extra-Galactic clusters appear to show a similar, early expansion phase (see Bastian et al. 2008, and references therein).
We have stressed the importance of binary heating to explain the global evolution of the Pleiades, both past and future. This is a three-body effect, not considered in classical studies of dynamical relaxation. As we indicated, binary heating is more effective in less populous systems, including open clusters. In the near future, we hope to explore further this general issue of stellar dynamics, i.e., the demarcation between systems that do and do not undergo classical, dynamical relaxation. This study will necessarily delve further into the role of binaries. We also intend to repeat our Pleiades analysis with another, relatively nearby system of comparable age, to ensure that the Pleiades results are representative for the entire class of open clusters.
Steve McMillan, one of the authors of Starlab, provided crucial assistance throughout this project. Not only did he instruct us in the workings of the code, but he even debugged portions of it at our request. Simon Portegies-Zwart, another Starlab author, also gave valued advice. Finally, we thank James Graham and Chris McKee for their continued interest and provocative questions. S. S. was partially supported by NSF grant AST-0908573.
where we have used equations (A7) and (9) to make the last transformation. Solving this equation for W ′ and substituting into equation (A3), we find This last relation gives us more information about the virial radius. We now see that the nondimensional version, ξ v ≡ r v /r 0 , obeys Thus, ξ v can be obtained at once from the nondimensional solution. Since r v itself is an input, the dimensional scale radius, r 0 , can be obtained from Similarly, the central density ρ 0 is while Ψ 0 is found from Ψ 0 = 4 π G ρ 0 r 2 0 (A15) = 4 π ξ v µ G N tot m r v .
(A16) Fig. 1.-Surface number density as a function of projected radius. The dashed curve represents our initial configuration, an n = 3 polytrope. The solid curve is a King (1962) model fit to our simulation results for the evolved cluster. Table 2 lists the parameters for this optimal model. The numerical results displayed are an average of 25 simulation runs. Also shown are Pleiades data with error bars, taken from Paper I.  -Initial distribution of the mass ratio within binaries, q ≡ m s /m p . The solid curve was obtained using a lognormal fit to the calculated single-star mass function, including the proper binary mass correlation parameter γ. The dashed curve is the hypothetical distribution obtained with the same single-star mass function, but with no mass correlation (γ = 0). The upper curve shows the three-dimensional virial radius r v , and is an average over simulation runs. The lower curve shows the projected core radius R c , and is also an average. The data point in the lower right is the observed Pleiades value for R c , along with error bars.  To the right of this curve is the observed Pleiades value. The lower curve shows the result from a single simulation run in which the mass segregation parameter β was artificially set to zero.