Issue 
A&A
Volume 536, December 2011



Article Number  A64  
Number of page(s)  12  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201118021  
Published online  08 December 2011 
Simulations of the Hyades
^{1} Astronomisches RechenInstitut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstrasse 1214, 69120 Heidelberg, Germany
email: aernst@ari.uniheidelberg.de
^{2} MaxPlanckInstitut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
^{3} National Astronomical Observatories of China, Chinese Academy of Sciences, Datun Lu 20A, Chaoyang District, 100012 Beijing, PR China
^{4} The Kavli Institute for Astronomy and Astrophysics at Peking University, Yi He Yuan Lu 5, Hai Dian Qu, 100871 Beijing, PR China
^{5} Main Astronomical Observatory, National Academy of Sciences of Ukraine, Akademika Zabolotnoho 27, 03680 Kyiv, Ukraine
Received: 5 September 2011
Accepted: 4 October 2011
Context. Using recent observational data we present Nbody simulations of the Hyades open cluster.
Aims. We make an attempt to determine the initial conditions of the Hyades cluster at the time of its formation to reproduce the presentday cumulative mass profile, stellar mass and luminosity function (LF).
Methods. We performed direct Nbody simulations of the Hyades in an analytic Milky Way potential. They account for stellar evolution and include primordial binaries in a few models. Furthermore, we applied a Kroupa initial mass function and used extensive ensembleaveraging.
Results. We find that evolved singlestar King initial models with King parameters W_{0} = 6−9 and initial particle numbers N_{0} = 3000 provide good fits to the presentday observational cumulative mass profile within the Jacobi radius. The bestfit King model has an initial mass of 1721 M_{⊙} and an average mass loss rate of −2.2 M_{⊙}/Myr. The Kband LFs of models and observations show reasonable agreement. Mass segregation is detected in both observations and models. If 33% primordial binaries are included, the initial particle number is reduced by 5% as compared to the model without primordial binaries.
Conclusions. The presentday properties of the Hyades can be reproduced well by a standard King or Plummer initial model when choosing appropriate initial conditions. The degeneracy of goodfitting models can be quite high due to the large dimension of the parameter space. More simulations with different Rochelobe filling factors and primordial binary fractions are required to explore this degeneracy in more detail.
Key words: methods: data analysis / stars: luminosity function, mass function / open clusters and associations: individual: Hyades / methods: numerical
© ESO, 2011
1. Introduction
The Hyades (“rain stars”) are one of the bestknown open star clusters on the northern sky. They are located in the constellation Taurus, close (on the celestial sphere) to the Pleiades or “Seven Sisters”, the most prominent open star cluster on the northern night sky. The Hyades and the Pleiades form the “golden gate of the ecliptic”, i.e. the ecliptic separates the two clusters on the celestial sphere. The Hyades cluster is markedly older than the Pleiades cluster. It has an age of 625 ± 50 Myr derived from the helium abundance in combination with isochrone modeling, which includes convective overshooting (Perryman et al. 1998). The shape of the Hyades HertzsprungRussell (HR) diagram with a narrow and welldefined main sequence has been recently studied by von Leeuwen (2009). Previous Nbody simulations of the Hyades have been performed by Portegies Zwart et al. (2001), Madsen (2003, see also references therein) and Chumak et al. (2005). These studies are in many respects similar to the present study (although special emphasis is placed on different aspects of the modeling of the Hyades). The study by Portegies Zwart et al. (2001) differs from the present study by the fact that they did not include a kick for the white dwarfs but argued that they could be hidden in binary systems with considerably brighter companion stars. Madsen (2003) concentrated on investigating the accuracy of astrometric radial velocities and compared the HR diagrams and internal velocity dipersions of observations and models. On the other hand, Chumak et al. (2005) studied the Hyades orbit and provided first models of the tidal tails of the Hyades.
The starting point for the present work is an observational data file (Röser et al. 2011) from which masses, positions, and velocities for 724 probable member stars of the Hyades can be derived. The PPMXL catalog (Röser et al. 2010) contains positions and proper motions of the Hyades members down to 0.116 M_{⊙}^{1}. The membership has been determined with the convergentpoint method. Given the right ascension, declination, and the respective proper motions, the convergentpoint method predicts a heliocentric distance (secular parallax) and radial velocity for each candidate member. Thus the full 6D phase space information is given. Röser et al. determined the individual masses from the masstoluminosity relation by Pinsonneault et al. (2004), from Dartmouth isochrones (Dotter et al. 1998), and from BCAH isochrones (Baraffe et al. 1998). Binarity is known only for a minor portion of the sample (Röser et al. 2011).
Table 1 gives a few parameters of the Hyades cluster as determined from the observational data file. We determined the tidal (or Jacobi) radius (King 1962) iteratively from the observational data according to the procedure described in Sect. 4 of Ernst et al. (2010) under the assumption that the solar radius is at R_{0} = 8 kpc (values for R_{0} = 8.5 kpc are given in brackets) and that β = κ/Ω = 1.37, where κ and Ω are the epicyclic and circular frequencies of a nearcircular orbit at the solar radius. The incompleteness in the stellar mass function sets in somewhere between 0.25 and 0.17 M_{⊙} (Röser et al. 2011).
Hyades parameters as determined from the observational data file.
This paper is organized as follows: Sect. 2 describes the numerical methods we use. Section 3 describes our integration of the Hyades orbit. In Sect. 4 we discuss our parameter space in detail. Section 5 summarizes the results. Section 6 contains the discussion and conclusions.
2. Numerical method
List of galaxy component parameters.
We solve the Nbody problem of the Hyades in an analytic background potential of the Milky Way. For the background Milky Way potential, we use an axisymmetric threecomponent model, where the bulge, disk, and halo are described by PlummerKuzmin models (Miyamoto & Nagai 1975) with the potential (1)The parameters a,b, and M of the Milky Way model are given in Table 2 for the three components. In the limit a → 0 this model reduces to a Plummer model (Plummer 1911). On the other hand, in the limit b → 0 the model reduces to a Kuzmin model (Kuzmin 1956).
Figure 1 shows the rotation curve of the threecomponent model of the Milky Way. As in our previous works, the parameters of the threecomponent model are chosen such that the rotation curve matches that of the Milky Way (Dauphole & Colin 1995). At the solar radius R_{0} = 8.0 kpc, which we have assumed in this study, the value of the circular velocity is V_{0} = 234.2 km s^{1}. The values of Oort’s constants A and B are consistent with the observed values (A,B) = (14.5 ± 0.8, −13.0 ± 1.1) km s^{1} kpc^{1} derived by Piskunov et al. (2006).
Figure 2 shows an “artistic” RGBA image of the (logarithmic) density structure of the threecomponent PlummerKuzmin model. The three components are overblended with an alpha channel with 90% transparency. The halo is plotted in blue, the disk in green, and the bulge in red.
Fig. 1 Rotation curve (at z = 0) of the threecomponent PlummerKuzmin model of the Milky Way. 
Fig. 2 “Artistic” RGBA image of the (logarithmic) density structure of the threecomponent PlummerKuzmin model. The three components are overblended with an alpha channel with 90% transparency. The halo is plotted in blue, the disk in green, and the bulge in red. 
For the solution of the Nbody problem in the analytic background potential of the Milky Way, the new direct Nbody code nbody6tidgpu (see Appendix A for a description) is used. Since this program is based on the code nbody6 by Aarseth (see Aarseth 1999, 2003), it is possible to follow stellar evolution within the HertzsprungRussell diagram. The implementation of stellar evolution for individual stars of the Nbody system is based on analytic formulae from which radius, luminosity, and stellar type can be derived as functions of the initial mass and age (Aarseth 2003; for the stellar evolution recipes see Hurley et al. 2000, and 2001). In particular, this approach includes the formation of stellar remnants like white dwarfs (WDs), neutron stars (NSs) and black holes (BHs). In addition to modeling standard evolutionary mass loss by stellar winds (which is implemented in many Nbody codes available today), Aarseth’s code and its variants include routines for treating large instantaneous mass loss due to special events like the occurrence of supernovae and the formation of planetary nebulae. (It is not very easy to correct energy, all forces, and first derivatives within the Nbody system for large instantaneous mass loss of one particle in the case of, say, a supernova.) Thus it is possible to model a star cluster realistically with nbody6 (and its variants) even in the first 55 Myr of evolution in which the supernova events occur (e.g., Hansen & Kawaler 1994, their Eq. (1.88)).
The code nbody6 and its variants optionally apply velocity kicks to stellar remnants when a supernova occurs or a planetary nebula forms. The random kick velocities for supernova kicks are drawn from a Maxwellian distribution with a 1D velocity dispersion of ≈ 190 km s^{1} corresponding to a mean of ≈ 300 km s^{1} (Hansen & Phinney 1997). For all WDs, the 1D dispersion of the Maxwellian is taken to be 5 km s^{1} which corresponds to a mean of ≈ 8 km s^{1} (e.g., Fellhauer et al. 2003, for a lower limit on that dispersion). We stress that the final number of WDs in our models depends critically on this dispersion. The escape velocity from the center of the cluster is ≈ 6 − 7 km s^{1} for an open cluster with a mass of 2000 M_{⊙} and a halfmass radius of 3.5 pc, assuming that it has a Plummer profile (Plummer 1911).
Finally, we note that in nbody6tidgpu, all stars are kept in the simulation forever, i.e. there is no optional removal of escapers as in nbody6 and its other variants. The reason is that we would like to follow the stellar orbits in the tidal tails and investigate their properties.
3. Orbit
Fig. 3 3D orbit of the Hyades in the threecomponent PlummerKuzmin model of the Milky Way. The integration time is 625 Myr. 
We integrated the center of mass orbit of the Hyades backwards in time in the analytic Milky Way potential (i.e. the threecomponent PlummerKuzmin model described in Sect. 2). We neglected encounters with giant molecular clouds and the effect of spiral arm passages and disk shocking. The integration time (625 Myr) was equal to the most probable age of the Hyades (Perryman et al. 1998). The full 3D orbit is shown in Fig. 3. The orbit integration with a simple integrator yielded the initial position of the Hyades at the time of its formation. From the initial position of the Hyades we ran full Nbody models up to the present time.
The presentday position of the Hyades was taken to be The first terms in Eqs. (2) − (4) correspond to the center of mass of the observed Hyades spatial distribution with respect to the position of the Sun. The second terms in Eqs. (2) and (4) represent the solar position in the Milky Way. We used (R_{0},z_{0}) = (8.00,0.02) kpc (consistent with Piskunov et al. 2006).
The presentday velocity was given by The first terms in Eqs. (5) − (7) represent the velocity of the center of mass of the observed Hyades distribution with respect to the Sun. The second terms in Eqs. (5) − (7) correspond to the peculiar motion (U,V,W)_{⊙} of the Sun with respect to the local standard of rest (LSR) according to Schönrich et al. (2010). The third term in Eq. (6) is the orbital speed of the LSR at R_{0} = 8 kpc in the threecomponent PlummerKuzmin model of Sect. 2.
With these parameters, the Hyades orbit around the Galactic center is oriented in a clockwise direction as seen from the Galactic north pole. The orbit has a considerable vertical oscillation in the zdirection between ≈ ± 74 pc. In the radial direction the orbit is mildly eccentric and oscillates between 7087 pc and 8638 kpc. Roughly three orbits are completed in 625 Myr.
4. Parameter space
We were looking for the bestfitting King models (King 1966) and Plummer models (Plummer 1911) to the observations of Röser et al. after 625 Myr of evolution with respect to the final cumulative mass profile M(r) within 0 < r < 9 pc (Jacobi radius). We compare the cumulative mass profiles of models and the observational data and not the density profiles because the cumulative mass profiles are of a better quality with respect to noise. Second, we note that the observational data of Röser et al. contains only 724 stars, but the observational presentday mass function (PDMF) drops below the completeness limit, which lies between 0.25 and 0.17 M_{⊙}. If we assume that the PDMF is shallow at the lowmass end according to Eq. (9) below, ≈ 100 stars would be missing in the lowmass bins above 0.116 M_{⊙} (see Röser et al. 2011, their Fig. 10). We took that into account for the comparison of the particle numbers and total masses of observation and model.
All our models are initially Rochelobe filling, i.e. the initial 99% Lagrange radius r_{99%} was set to be equal to the tidal (or Jacobi) radius r_{J}. The 99% Lagrange radius is defined as the radius that contains 99% of the cluster mass. The Jacobi radius is given by (King 1962) (8)where G, M, Ω, and κ are the gravitational constant, the cluster mass, the circular, and the epicyclic frequencies of a nearcircular orbit, respectively.
We defined the initial mass function (IMF) as the number of stars per unit logarithmic mass interval. It is given by (9)where c_{1} − c_{4} are four constants (determined by the normalization and three conditions for the transitions between the mass regimes). The exponents in our Eq. (9) are the exponents α_{0} − α_{3} given in Kroupa (2001), his Eq. (2), adjusted for logarithmic mass bins. Thus our IMF is a Kroupa (2001) IMF. We forced the lowest and the highest mass to take on the values m_{l} = 0.08 M_{⊙} and m_{h} = 100.00 M_{⊙}, respectively. Thus the first power law in (9) for the brown dwarf regime is obsolete. In nbody6tidgpu, the individual masses were randomly drawn from the above distribution (9) using a routine by Kroupa with a correction by Weidner.
For all models we switched on stellar evolution in nbody6tidgpu. For the metallicity we used the observed value for the Hyades Z = 0.024 (Perryman et al. 1998).
More than 400 individual models have been computed. The run time for one model on a dualXeon 3.2 GHz with a GeForce 8800 GTS 512 graphics card was about 45 min.
In a first step, a grid of 30 King models (King 1966) in the parameter space (W_{0}, N), where W_{0} is the King parameter and N the initial particle number, was integrated for an overview of the evolution of different models. We used for the King parameters W_{0} = 3,4,5,6,7,8,9,10,11,12 and particle numbers N = 3000,4000,5000. Thus the grid cell size in the parameter space is given by (ΔW_{0},ΔN) = (1,1000). We did not include primordial binaries in these models.
Using this first grid, we were able to notice trends in the choices of N and W_{0}. The variation in W_{0} turned out to be weak. Also, the grid was too coarse to reach the desired accuracy in N.
In the second step, we therefore continued with an iterative search in the parameter space (W_{0},N) for the bestfitting models to the observations without primordial binaries. For this purpose we refined our first grid in the parameter space with respect to N to a grid cell length of ΔN = 125 particles. On the other hand, it turned out to be useful to coarsen the grid with respect to W_{0} to a grid cell height of ΔW_{0} = 3, i.e. we restricted the investigation to the models with W_{0} = 3,6,9,12. Thus the grid cell size of the second grid in the parameter space was given by (ΔW_{0},ΔN) = (3,125).
Our criterion for the best fit is that the inner cumulative mass profiles within r = 9 pc (Jacobi radius) of observations (Röser et al. 2011, their Fig. 6) and models (ensemble mean, see below) show the best agreement.
We did not run all models in this second grid. Rather, we started from the bestfitting model in the first grid for a given W_{0} and varied the particle number N in steps of ΔN = 125 particles in order to find the bestfitting model in the second grid.
At the same time, we performed n = 15 runs with different random number seeds for each initial model in the second grid. From these 15 runs we obtained physical quantities by ensemble averaging. For example, we compared the observations with the mean cumulative mass profile of 15 models to find the best fit acording to the criterion mentioned above. The random number generator is invoked mainly for the following purposes:

drawing of initial masses, positions, and velocities from the givenprobability distributions (Kroupa 2001; King 1966);

assignment of kick velocities for WDs, NSs, and BHs;

assignment of perihelion, node, and inclination of primordial binaries.
We notice that we did not simply ensembleaverage the initial model from 15 initial models, which differed by their random number seeds. In particular, we did not perform only one single run from such an ensembleaveraged initial model but really performed 15 individual runs, because

we are looking for standard deviations in the physical quantities;

the modeling of the population of stellar remnants should be as realistic as possible.
In order to verify that we found indeed the bestfitting ensembles, we looked at the cumulative mass profile of at least one of the neighboring ensembles with respect to N.
The ensembles in the second parameter grid are named, for example, en3000W6 (N_{0} = 3000,W_{0} = 6), i.e. the letters “en” for “ensemble” are followed by the initial particle number N_{0} and the King parameter W_{0}.
We also looked for the bestfitting ensemble of singlestar Plummer models (Plummer 1911) with respect to the inner cumulative mass profile (within r = 9 pc) also using steps of ΔN = 125. The choice of r_{99%}/r_{J} fixes the Plummer radius, while the King models still have the free parameter W_{0}. The Plummer model runs are discussed in Sect. 5.2.4.
We ran the bestfitting ensemble of singlestar Plummer models with primordial binaries in order to see their effect on the evolution. For models with primordial binaries we defined the primordial binary fraction f_{b} and the fraction of stars in binaries f_{sb} as (10)where N_{s} and N_{b} are the initial numbers of binaries and single stars, respectively (e.g., Kroupa 1995; Jahreiß & Wielen 2000; Heggie et al. 2006). The periodgenerating function for the primordial binaries was given by (Kroupa 1995, Eq. (11b)) (11)with a uniform random variate X ∈ [0,1] . The parameters are η = 2.5, δ = 45.0, and P_{min} = 5 days. For the binary eccentricities a thermal distribution was adopted with the eccentricitygenerating function e^{2} = X where e is the eccentricity and X ∈ [0,1] a uniform random variate. We did not assume any correlation between the masses of the binary components (e.g., Eggleton et al. 1989).
5. Results
For the statistical analysis, we applied the quantities (12)where n = 15 is the number of models in an ensemble, x_{i} are the individual measurements, is their mean, σ the standard deviation (i.e. the mean error of a measurement in a single model), and σ_{m} the standard error (i.e., the mean error of the mean).
5.1. Cumulative mass profiles
The contamination with field stars in the observational data file is negligible for r < 9 pc, 7.5% at 9 pc < r < 18 pc, and 30% at 18 pc < r < 30 pc (with respect to the number of stars; Röser et al. 2011), where r is measured from the center of the cluster.
Three comparisons of the cumulative mass profiles are shown in Fig. 4. The top panel shows the comparison for the ensemble of models en3000W6 (N_{0} = 3000,W_{0} = 6). Shown is the mean cumulative mass profile of the ensemble of 15 models, i.e., we averaged over the 15 single cumulative mass profiles of the individual models. Compared with the observations, the ensemble en3000W6 is the bestfitting of all our models according to the criterion stated in Sect. 4. Assuming that the model is perfect, the field star contamination is M_{model} − M_{obs}/M_{obs} = 2% (3 pc), 2% (9 pc), 12% (18 pc), 21% (30 pc), where M_{mod} and M_{obs} are theoretical and observed cumulative masses.
Fig. 4 Comparison of the presentday cumulative mass profiles of observations and models. The thick solid (red) lines in the middle of the models represent the mean, which is ensembleaveraged over 15 runs. The standard error is shown as a filled area. The dashed (blue) line is calculated from the observational data. Top: the ensemble en3000W6. Middle: the ensemble en3000W9. Bottom: the ensemble en3875W12. 
The secondbest fitting ensemble of King models (en3000W9; N_{0} = 3000,W_{0} = 9) is shown in the middle panel of Fig. 4. The ensemble mean of the cumulative mass profile also fits the inner part of the observational profile fairly well. The field star contamination here is M_{model} − M_{obs}/M_{obs} = 12% (3 pc), 5% (9 pc), 18% (18 pc), 26% (30 pc).
The bottom panel shows the comparison of the cumulative mass profiles for the observations and the ensemble en3875W12 (N_{0} = 3875,W_{0} = 12). This model has been chosen because the agreement with the observations is very good in the core (r < 3 pc) and for the final total mass contained within a radius of r = 30 pc from the cluster center. In the range 3 pc < r < 18 pc, the fit is fairly bad, i.e. inconsistent within 1σ_{m}. A better fit of the core region cannot be obtained for W_{0} = 12. For W_{0} = 3 the models dissolve so fast that a good fit with respect to the inner cumulative mass profile or the final total mass has not been obtained at all. In Appendix B we show the best fits for the presentday total mass within r = 30 pc from the center, which, however, do not agree with the inner cumulative mass profile.
5.2. Past evolution and presentday state
Parameters of the best fitting ensemble of models en3000W6.
Fig. 5 Particle number within r = 30 pc from the cluster center as a function of time for all 15 runs of the ensemble en3000W6. One can see the nonMarkovian process of star cluster dissolution. The mean final particle number and the standard deviation for a single run is shown as a (red) dot with vertical errorbars at t = 625 Myr. The future evolution of the mean particle number within r = 30 pc is shown in the inset. There the times when the cluster has reached 20, 10 and 5% of its initial particle number are shown as (red) dots. 
The bestfitting ensemble of models (en3000W6) for the inner cumulative mass profile is an ensemble average over 15 runs with different random number seeds. The initial and final parameters of the ensemble en3000W6 are given in Table 3. We note that (1) the final parameters are obtained after 625 Myr of evolution, i.e. at the present time; (2) the final “observed” particle number is the number of stars with m > 0.116 M_{⊙} and that (3) the highest mass is obtained by excluding BHs.
In the following discussion, we look at stars within distances of r = 30 pc of the cluster center, within r = 9 pc (Jacobi radius) and at stars within r = 3 pc (core). Figure 5 shows the evolution of the particle number within r = 30 pc from the cluster center. It can be seen how the 15 cluster models of the ensemble en3000W6 dissolve as time proceeds. The mean final particle number and the standard deviation for a single run is shown with errorbars at t = 625 Myr. The standard deviation is rather high due to the nonMarkovian evolution of the models, which produces rms scatter between the individual models of an ensemble. The average particle loss rate until t = 625 Myr is dN/dt = −3.6/Myr. The average mass loss rate until t = 625 Myr is dM/dt = −2.2 M_{⊙}/Myr. The inset of Fig. 5 is explained in Sect. 5.3.
Figure 6 shows the presentday state of the cluster with its tidal tails for one run of the ensemble en3000W6. The tidal tails have reached a length of 800 pc after 625 Myr of evolution. However, the kickedout NSs and BHs (and a few WDs) hurried ahead along the cluster orbit during the evolution and reside at t = 625 Myr spread out along the orbit well beyond the tips of the tidal tails at 800 pc.
Table 4 shows the values of the 3D velocity dispersion σ_{3D} for the ensemble en3000W6 within 30, 9, and 3 pc of the cluster center. The three values are consistent with the fact that the velocity dispersion in the core is higher than in the halo within the Jacobi radius and rises again outside of the Jacobi radius. The values can be compared with those given in Röser et al. (2011).
Fig. 6 The presentday state of one run of the ensemble of models en3000W6 after 625 Myr of evolution. The tidal tails have reached a length of 800 pc. 
3D velocity dispersions at t = 625 Myr for the ensemble en3000W6 for stars 30, 9 and 3 pc of the cluster center.
5.2.1. Stellar mass functions
Figure 7 shows the average realization of the Kroupa (2001) IMF of the ensemble en3000W6 with lowest mass m_{l} = 0.08 M_{⊙} and initial particle number N_{0} = 3000. Shown is the mean value of log (N) from the ensemble averaging over 15 runs. The errors are 2σ_{m} standard errors. The analytical slopes and the transition at 0.5 M_{⊙} are also shown as lines. Note that the depletion in the lowest mass bin is due to the binning.
Figure 8 shows a comparison of the PDMFs between observations and ensemble en3000W6. The top panel shows the stars within a radius of r = 30 pc from the center; the middle panel shows the stars within a radius of r = 9 pc from the center, i.e. within the Jacobi radius; the bottom panel shows the stars within a radius of r = 3 pc from the center (core). Shown is the mean value of log (N) from the ensemble averaging over 15 runs. The errors are 1σ_{m} standard errors. These two particle numbers N of model and observations are indeed realized in the top panel. This cannot be seen easily because of an “optical illusion”: the deviations between ensemble and observation carry more weight at high values of log (N). Furthermore, the completeness limit of the observations is around 0.25 M_{⊙} (Röser et al. 2011).
Fig. 7 The average realization of the Kroupa (2001) IMF of the ensemble en3000W6 with lowest mass m_{l} = 0.08 M_{⊙}, highest mass m_{h} = 100.00 M_{⊙}, and initial particle number N_{0} = 3000. Shown is the mean value of log (N) from the ensemble averaging over 15 runs. The 2σ_{m} standard errors are shown. The analytical slopes and the transition at 0.5 M_{⊙} are also shown as lines. 
Fig. 8 Comparison of the PDMFs between observations and ensemble of models en3000W6. Top panel: stars within a radius of 30 pc of the center; middle panel: stars within a radius of 9 pc of the center (Jacobi radius); bottom panel: stars within a radius of 3 pc of the center (core). Shown is the mean value of log (N) from the ensemble averaging over 15 runs. The errors are 1σ_{m} standard errors. The blue bars are derived from the observational data with N = 724 stars and lowest mass m_{l} = 0.116 M_{⊙}, while the red bars are derived from the ensemble of models with ⟨ N_{f} ⟩ = 736 stars and lowest mass m_{l} = 0.08 M_{⊙}. The onset of incompleteness of the observations is shown as a vertical line in the top panel. 
The following details can be seen in Fig. 8: in all panels the upper main sequence (m > 1 M_{⊙}) is overabundant in the observations. The reason may be simply that the IMF of our models is wrong. Also, unresolved binaries could be the reason (see Sect. 5.2.4). Binaries are not resolved below 3′′, i.e. 133 astronomical units at a distance of 46 pc (Röser, priv. comm.). In the upper two panels, the PDMF of the ensemble is nearly flat at the lowmass end while the observed PDMF decreases below the completeness limit of the observations. A dip in the PDMFs between m = 0.5 − 0.6 M_{⊙} can be seen in the upper two panels. This dip in the model coincides with the location of the Wielen dip (Kroupa et al. 1990) which is related to the shape of the massluminosity relation (Kroupa 2002). The stellar evolution in our models is different for m < 0.7 M_{⊙} and m > 0.7 M_{⊙} (Hurley et al. 2000).
Fig. 9 Presentday cumulative mass function for the ensemble en3000W6 for r < 9 pc from the cluster center. Shown is the sum of all masses which are higher than a given mass. The errors are 1σ_{m} standard errors. The completeness limit of the observations is shown as a vertical solid line. 
Figure 9 shows the cumulative stellar mass function for the ensemble en3000W6 for stars within a radius of r = 9 pc of the cluster center, i.e. within the Jacobi radius. Shown is the sum of all masses that are higher than a given mass. The histogram includes main sequence stars, giants, and WDs. Stellar mass BHs are excluded and NSs do not occur. Up to the completeness limit of the observations, there is a discrepancy of 10%, which leaves room for further optimization in future models.
5.2.2. Stellar evolution
Number counts of stellar types for the bestfitting ensemble of models (en3000W6).
Table 5 shows the number counts of stellar types for the bestfitting ensemble of models (en3000W6) for all 3000 stars (initial and final), for stars within a radius of 30 pc (final), 9 pc (final inside Jacobi radius), and 3 pc (final inside core). The stellar number counts are averaged over the 15 runs of the ensemble. The stellar types are the same as in the classification in the beginning of Sect. 4 of Hurley et al. (2000). The standard errors from the ensembleaveraging over 15 runs are also given. We find 1 ± 0 stellar mass BHs within r = 3 pc. Almost all stellar mass BHs and all NSs have been kicked out because of supernova kicks and are spread out along the cluster orbit. We find 27 ± 2, 22 ± 2, and 9 ± 1 carbonoxygen WDs within r = 30, 9, and 3 pc, respectively. These numbers can be compared to the observations (Schilbach & Röser 2011).
Figure 10 shows a presentday synthetic HertzsprungRussell diagram for one run of the evolved ensemble en3000W6. Shown are stars within a radius of r = 30 pc from the cluster center. The synthetic diagram in Madsen (2003) is more suited to direct comparison with observations since he transforms effective temperatures to B − V colors using color tables and bolometric magnitudes to V magnitudes. In addition, he introduces a small scatter in both V and B − V. We cannot directly compare to the observations since the data in B − V are missing for the data set used in the present study.
Figure 11 show the comparison of the presentday luminosity functions (PDLFs) for the ensemble en3000W6 and the observations by Röser et al. (2011). In the top, middle and bottom panels stars within a radius of 30, 9 (Jacobi radius) and 3 pc (core) from the cluster center have been used in the statistics. Shown is the mean value of log (N) as a function of the 2MASS K_{s} band magnitude (observations) and K band magnitude (model) from the ensemble averaging over 15 runs. The errors are 1σ_{m} standard errors.
Fig. 10 Presentday synthetic HertzsprungRussell diagram for one run of the evolved model en3000W6. Shown are stars within a radius of r = 30 pc from the cluster center. The final particle number within r = 30 pc for this run is N_{f} = 836. 
Fig. 11 PDLFs in the K band for the ensemble en3000W6 and the observations (Röser et al. 2011). For explanations see the text. Shown is the mean/measured value of log (N) as a function of the absolute K_{s} (observations)/K (model) band magnitude. The errors are 1σ_{m} standard errors from the ensemble averaging. 
To obtain the absolute K band magnitudes from the theoretical luminosities given by nbody6, we applied the compilation of infrared bolometric corrections given in Appendix A of Muench (2002) with the calibration m_{bol, ⊙ } = 4.75. The applied conversion formula is (13)where M_{K}, L/L_{⊙}, M_{bol, ⊙ }, and BC_{K} are the absolute K band magnitudes, the luminosity in units of the solar luminosity, the bolometric magnitude of the Sun, and the K band bolometric correction, respectively.
In Fig. 11 there are giants in the bin centered at M_{K} = −1.75 of the observational LF. For these luminous objects it is hard to obtain the correct observational M_{K} owing to overexposure (Roeser & Schilbach, priv. comm.). On the other hand, in the theoretical LF the WDs with log _{10}(T_{eff}) = 4.0−4.4 populate the faint end (M_{K} > 10). Between 1 < M_{K} < 5 it can be seen that the observational M_{K}’s are overabundant as compared to the model M_{K}’s.
5.2.3. Mass segregation
A first unmistakable sign of mass segregation is the increased mean mass within the core in Table 3 as compared to the mean mass within the Jacobi radius. Figure 12 shows the average cumulative mass for observation and ensemble en3000W6. The average cumulative mass is calculated only for stars with mass m > 0.116 M_{⊙}. Also shown is the line for the ensemble enb2125PL with 33% primordial binaries (see Sect. 5.2.4 below). The figure clearly shows the segregation of masses within the Jacobi radius.
Fig. 12 Average cumulative mass for observations and ensembles en3000W6 and enb2125PL. The average cumulative mass is calculated only for stars with mass m ≥ 0.116 M_{⊙}. 
Figure 13 shows the result of a detailed analysis of mass segregation based on the minimum spanning tree (MST) method developed by Olczak et al. (2011). They define a measure of mass segregation Γ_{MST}, where e_{i} are the lengths of the n MST edges. The superscript “ref” refers to a sample of n + 1 random stars from the entire population while the superscript “mass” refers to an equalsized sample of the most massive stars. Note that γ_{MST} has the dimension of length while Γ_{MST} is a dimensionless measure. Γ_{MST}·(ΔΓ_{MST})^{ − k} > 1 implies mass segregation of the n + 1 most massive stars with significance kσ.
Fig. 13 Comparison of presentday mass segregation between models and observations. Top figure: analysis of the evolved (for 625 Myr) ensemble en3000W6. Bottom figure: analysis of the observational data of Röser et al. (2011). Each single plot contains on the lefthand side the “cumulative” Γ_{MST} (points) for the 5 (red), 10 (green), 20 (blue), 50 (magenta), 100 (cyan), 200 (brown), 500 (orange), and all most massive stars (black). The righthand side shows the “differential” Γ_{MST} (points) for the 5 (red), 6 − 10 (green), 11 − 20 (blue), 21 − 50 (magenta), 51 − 100 (cyan), 101 − 200 (brown), 201 − 500 (orange), and 501all most massive stars (black). The error bars mark the 1, 2, and 3σ uncertainties. A value of one marks the unsegregated state. 
The top panel of Fig. 13 shows the analysis of the evolved (for 625 Myr) ensemble en3000W6, while the bottom panel shows an analysis of the observational data of Röser et al. (2011). Each single plot contains on the lefthand side the “cumulative” Γ_{MST} (points) for the 5, 10, 20, 50, 100, 200, 500, and all most massive stars (black). On the righthand side the “differential” Γ_{MST} (points) is shown for the 5, 6 − 10, 11 − 20, 21 − 50, 51 − 100, 101 − 200, 201 − 500, and 501all most massive stars. The error bars mark the 1σ, 2σ, and 3σ uncertainties. A value of one marks the unsegregated state. The higher the values, the more segregated the given mass group is (compared to a set of random groups). The results of this analysis are the following.

1.
The ensemble average of the numerical models resembles theobservational data very closely. Except for a larger deviation ofthe fourth bin the agreement is excellent.

2.
Both observations and simulations show a significant degree of mass segregation, mostly above the 3σ level for the cumulative analysis.

3.
There is a clear signature of “inverse mass segregation”: Γ_{MST} < 1 for the last bin in the differential plot (i.e. the ~ 200 least massive stars). This demonstrates the removal of the lowest mass members from the inner cluster parts due to dynamical mass segregation.

4.
The signature of mass segregation agrees fairly well with a moderately (S = 0.3: cf. Šubr et al. 2008) masssegregated model star cluster of 1000 particles as shown in the middle panel of Fig. 4 in Olczak et al. (2011).
5.2.4. Primordial binaries
Fig. 14 Comparison of the presentday cumulative mass profiles of observations and models with primordial binaries. The thick solid (red) lines are derived from the ensemble en2250PL without binaries, the dashed (green) line is derived from the ensemble enb2250PL with 33% primordial binaries, and the dotted (magenta) line is derived from the bestfitting ensemble with 33% primordial binaries enb2125PL. The dashed (blue) line represents the observations. The standard errors are shown as filled areas. 
Fig. 15 PDLFs in the K band for the ensemble enb2125PL where all binaries are either resolved or unresolved. Top panel: stars within a radius of 30 pc of the center; middle panel: stars within a radius of 9 pc of the center (Jacobi radius); bottom panel: stars within a radius of 3 pc of the center (core). Shown is the mean value of log (N) as a function of the absolute K band magnitude. 
Fig. 16 Comparison of the PDMFs between ensembles enb2250PL, en2250PL, and enb2125PL. Top panel: stars within a radius of 30 pc of the center (Jacobi radius); middle panel: stars within a radius of 9 pc of the center (Jacobi radius); bottom panel: stars within a radius of 3 pc of the center (core). Shown is the mean value of log (N) from the ensemble averaging over 15 runs. The errors are 1σ_{m} standard errors. 
The recent work of Sollima et al. (2009) showed that open clusters may contain between f_{b} ≈ 35 − 70% binaries. We have run the bestfitting ensemble en2250PL of Plummer models again with a primordial binary fraction f_{b} = 33% (f_{sb} = 50% from Eq. (10)) in order to investigate the influence of primordial binaries. The ensemble has been called enb2250PL. The third letter “b” denotes primordial binaries and the following number is the number of stars (i.e. not the number of systems). Since this ensemble does not fit well to the inner observational cumulative mass profile, we reduced the initial particle number further and finally found a bestfitting model enb2125PL (also with 33% primordial binaries).
Figure 14 shows the comparison of the presentday cumulative mass profiles of observations and models with 33% primordial binaries. The standard errors are shown as filled areas. It can be seen that the profile of the ensemble enb2250PL is steeper than that of the ensemble en2250PL at small radii. The ratio of halfmass to Jacobi radius for the model enb2125PL is given by r_{h}/r_{J} = 0.18 ± 0.01. The field star contamination for the model enb2125PL is M_{model} − M_{obs}/M_{obs} = 12% (3 pc), 0% (9 pc), 13% (18 pc), and 20% (30 pc).
Figure 15 shows the PDLFs in the K band for the ensemble enb2125PL where all binaries are either resolved or unresolved. For the unresolved binaries we transformed for each star the bolometric luminosities given by nbody6 to absolute K band magnitudes using the bolometric corrections of Muench (2002) and then back to K band luminosities using the absolute Kband magnitude of the Sun M_{K, ⊙ } = 3.33. Then we added the K band luminosities of both binary components and transformed back to absolute K band magnitudes, which are shown in Fig. 15. With the present parameters of the distribution of binaries, unresolved binaries produce a difference at the faint end of the LF. However, our statistics of the binary parameters may be not fully realistic.
Figure 16 shows a comparison of the PDMFs between ensembles enb2250PL, en2250PL, and enb2125PL. It can be seen that the ensembles en2250PL and enb2125PL are depleted of highmass stars as compared to the ensemble enb2250PL for r < 30 pc and r < 9 pc. Furthermore, the comparison of ensembles en2250PL and enb2125PL shows that the ensemble enb2125PL is depleted of lowmass stars compared to the ensemble en2250PL, while highmass stars are overabundant in the ensemble enb2125PL as compared with the ensemble en2250PL.
The initial and final parameters of the ensemble enb2125PL are given in Table 6. We note that (1) the final parameters are obtained after 625 Myr of evolution, i.e. at the present time; (2) the final “observed” particle number is the number of stars with m > 0.116 M_{⊙} and that (3) the highest mass is obtained by excluding BHs.
Parameters of the ensemble enb2125PL with 33% primordial binaries.
Number counts of stellar types that are the component of a binary for the ensemble enb2125PL.
Table 7 shows the number counts of stellar types that are the component of a binary for the ensemble enb2125PL for all 2125 stars (final), for stars within a radius of 30 pc (final), 9 pc (final inside Jacobi radius), and 3 pc (final inside core). The stellar number counts are averaged over the 15 runs of the ensemble. The stellar types are the same as in the classification in the beginning of Sect. 4 of Hurley et al. (2000). There are no BHs or NSs in binaries. Also, there are no WDs in binaries within 30 pc radius of the cluster center. They are all kicked out by binary kicks. Since WDs in binaries are observed in the Hyades (Schilbach & Röser 2011) they may not be subject to strong kicks.
5.3. Future evolution
The future evolution of the ensemble en3000W6 is shown in the inset of Fig. 5. There the times when the cluster has reached 20, 10 and 5% of its initial particle number are shown. The times are t_{20} = 695 Myr, t_{10} = 875 Myr, and t_{5} = 1029 Myr. The last time may be taken as an approximate dissolution time.
6. Discussion
The presentday properties of the Hyades can be reproduced well by a standard King or Plummer initial model when choosing appropriate initial conditions. We have found a model (en3000W6), which imitates the cumulative mass profile within the Jacobi radius of the observed Hyades spatial distribution very well. Also, the derived PDMF and the K band LF show a good agreement between observations and model.
From our models we can make the following statements about the Hyades:

The tidal tails of the Hyades have a length of ≈ 800 pc today if they are not destroyed by passing giant molecular clouds, spiral arm passages, disk shocking or other effects.

Rochelobe filling W_{0} = 6−9 singlestar King initial models (King 1966) with the initial particle number N_{0} = 3000 provide a good fit to the presentday inner cumulative mass profile.

The bestfit, singlestar, Rochelobe filling King model has an initial mass of M_{0} = 1721 M_{⊙}, an initial Jacobi radius of r_{J} = 16.2 pc, and an average mass loss rate of dM/dt = 2.2 M_{⊙}/Myr.

A reasonably good fit is obtained with a Plummer model with 33% primordial binaries and a ratio of halfmass to Jacobi radius of r_{h}/r_{J} = 0.18. Here the initial particle number is N_{0} = 2125 and mass are is M_{0} = 1230 M_{⊙}. The average mass loss rate is dM/dt = 1.4 M_{⊙}/Myr. The observed average cumulative mass is reproduced very well.

Mass segregation is clearly detected in both the observations and models of the Hyades cluster.

The Hyades contain only 1 ± 0 stellar mass black holes. The probability is very high that nearly all NSs and BHs are kicked out by supernova kicks.

The number of WDs critically depends on the kick velocity that is adopted for WD kicks. This kick velocity has not yet been well constrained. Assuming a mean kick velocity of ≈ 8 km s^{1} (comparable to the escape velocity from the cluster center), we obtain 27 ± 2, 22 ± 2, and 9 ± 1 carbonoxygen WDs within a radius r < 30, r < 9 (Jacobi radius), and r < 3 pc (core) from the cluster center, respectively. For the comparison to the observations we refer to the article by Schilbach & Röser (2011). We do not find WDs in binaries within r < 30 pc in our model. They are all kicked out in the form of binary kicks. However, some WDs in binaries are observed in the Hyades (Schilbach & Röser 2011). This suggests that WDs in binary systems are not subject to strong kicks.

The including of 33% primordial binaries reduces the initial number of Hyades stars N_{0} by ≈ 5% as compared to the model without primordial binaries.

Under the assumption that the singlestar ensemble en3000W6 is a good approximation of reality, we find that the Hyades will be dissolved except for 5% of its initial number of stars in ≈ 400 Myr from now.
The degeneracy of goodfitting models can be quite high due to the large dimension of the parameter space, which is spanned by the initial particle number N_{0}, the King parameter W_{0}, and two extra dimensions: the primordial binary fraction f_{b} and the Rochelobe filling factor r_{99%}/r_{J}. In this work we did not explore the latter two degrees of freedom in detail. More simulations with different Rochelobe filling factors and primordial binary fractions are required to explore this degeneracy in more detail.
Acknowledgments
This work was supported by the German Research Foundation (DFG) under the Collaborative Research Center SFB 881 (“The Milky Way System”) at the University of Heidelberg. Simulations were performed on the GRACE supercomputer (grants I/80 041043 and I/84 678680 of the Volkswagen Foundation and 823.219439/30 and /36 of the Ministry of Science, Research and the Arts of BadenWürttemberg). C.O. appreciates funding by the German Research Foundation (DFG), grant OL 350/11. P.B. acknowledges the special support by the NAS Ukraine under the Main Astronomical Observatory GRAPE/GPU/GRID computing cluster project. P.B.’s studies are also partially supported by the program Cosmomicrophysics of NAS Ukraine. The authors would like to thank Siegfried Röser and Elena Schilbach for providing their observational data and the reading of the manuscript, Keigo Nitadori and Sverre Aarseth for making nbody6gpu available and Jarrod Hurley for developing of the stellar evolution routines in nbody6. Furthermore, useful conversations with Siegfried Röser, Elena Schilbach and Sverre Aarseth are also acknowledged. A.E. would like to thank his father, Reinhard Ernst, for a discussion.
References
 Aarseth, S. J. 1999, PASP, 111, 1333 [NASA ADS] [CrossRef] [Google Scholar]
 Aarseth, S. J. 2003, Gravitational Nbody simulations – Tools and Algorithms (Cambridge: Cambridge Univ. Press) [Google Scholar]
 Baraffe, I., Chabrier, G., Allard, F., & Hauschild, P. H. 1998, A&A, 337, 403 [NASA ADS] [Google Scholar]
 Carraro, G., Ng, Y. K., & Portinari, L. 1998, MNRAS, 296, 1045 [NASA ADS] [CrossRef] [Google Scholar]
 Chumak, Ya. O., Rastorguev, A. S., & Aarseth, S. J. 2005, Astron. Lett., 31, 308 [NASA ADS] [CrossRef] [Google Scholar]
 Copenhagen Univ. Obs., Inst. of Astronomy, Cambridge, UK & Real Instituto Y Observatorio de LaArmada, F. E. S. 2006, VizieR Online Data Catalog, 1304, [Google Scholar]
 Dehnen, W. 1993, MNRAS, 265, 250 [NASA ADS] [CrossRef] [Google Scholar]
 Dotter, A., Chaboyer, B., Jevremović, D., et al. 2008, ApJS, 178, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Eggleton, P. P., Fitchett, M. J., & Tout, C. A. 1989, ApJ, 347, 998 [NASA ADS] [CrossRef] [Google Scholar]
 Ernst, A. 2009, Ph.D. Thesis, University of Heidelberg, http://www.ub.uniheidelberg.de/archiv/9375 [Google Scholar]
 Ernst, A., Just, A., & Spurzem, R. 2009, MNRAS, 399, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Ernst, A., Just, A., Berczik, P., & Petrov, M. I. 2010, A&A, 524, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fellhauer, M., Lin, D. N. C., Bolte, M., Aarseth, S. J., & Williams, K. A. 2003, ApJ, 595, 53 [Google Scholar]
 Hansen, C. J., & Kawaler, S. D. 1994, Stellar interiors: Physical Principles, Structure and Evolution (New York: Springer) [Google Scholar]
 Heggie, D. C., Trenti, M., & Hut, P. 2006, MNRAS, 368, 677 [Google Scholar]
 Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Tout, C. A., Aarseth, S. J., & Pols, O. R. 2001, MNRAS, 323, 630 [NASA ADS] [CrossRef] [Google Scholar]
 Janes, K. A., Friel, E., Montgomery, K., Phelps, R., & Marschall, L. 1992, Mem. Soc. Astron. Ital., 63, 283 [NASA ADS] [Google Scholar]
 Just, A., & Peñarrubia, J. 2005, A&A, 431, 861 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 King, I. R. 1962, AJ, 67, 471 [NASA ADS] [CrossRef] [Google Scholar]
 King, I. R. 1966, AJ, 71, 64 [Google Scholar]
 Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P. 2002, Science, 295, 82 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kroupa, P., Tout, C. A., & Gilmore, G. 1990, MNRAS, 244, 76 [NASA ADS] [Google Scholar]
 Kustaanheimo, P., & Stiefel, E. L. 1965, J. Reine Angewandte Mathematik, 218, 204 [CrossRef] [Google Scholar]
 Kuzmin, G. G. 1956, Astron. Zh., 33, 27 [Google Scholar]
 Madsen, S. 2003, A&A, 401, 565 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 McLachlan, R. 1995, SIAM J. Sci. Comp., 16, 151 [CrossRef] [Google Scholar]
 Mikkola, S., & Aarseth, S. J. 2002, Cel. Mech. Dyn. Astron., 84, 343 [Google Scholar]
 Mikkola, S., & Tanikawa, K. 1999, Cel. Mech. Dyn. Astron., 74, 287 [NASA ADS] [CrossRef] [Google Scholar]
 Muench, A. A. 2002, Ph.D. Thesis, University of Florida [Google Scholar]
 Nitadori, K., & Aarseth, S. J. 2011, MNRAS, submitted [Google Scholar]
 Olczak, C., Spurzem, R., & Henning, T. 2011, A&A, 532, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perryman, M. A. C., Brown, A. G. A., Lebreton, Y., et al. 1998, A&A, 331, 81 [NASA ADS] [Google Scholar]
 Portegies Zwart, S. F., McMillan, S. L. W., Hut, P., & Makino, J. 2001, MNRAS, 321, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Preto, M., & Tremaine, S. 1999, AJ, 118, 2532 [NASA ADS] [CrossRef] [Google Scholar]
 Pinsonneault, M. H., Terndrup, D. M., Hanson, R. B., & Stauffer, J. R. 2004, ApJ, 600, 946 [NASA ADS] [CrossRef] [Google Scholar]
 Piskunov, A. E., Kharchenko, N. V., Röser, S., Schilbach, E., & Scholz, R.D. 2006, A&A, 445, 545 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Plummer, H. C. 1911, MNRAS, 71, 460 [CrossRef] [Google Scholar]
 Röser, S., Demleitner, M., & Schilbach, E. 2010, AJ, 139, 2440 [NASA ADS] [CrossRef] [Google Scholar]
 Röser, S., Schilbach, E., Piskunov, A. E., Kharchenko, N. V., & Scholz, R.D. 2011, A&A, 531, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schilbach, E., & Röser, S. 2011, A&A, in press, DOI: 10.1051/00046361/201117688 [Google Scholar]
 Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829 [NASA ADS] [CrossRef] [Google Scholar]
 Sollima, A., Bellazzini, M., Smart, R. L., et al. 2009, MNRAS, 401, 577 [Google Scholar]
 Šubr, L., Kroupa, P., & Baumgardt, H. 2008, MNRAS, 385, 1673 [NASA ADS] [CrossRef] [Google Scholar]
 Sundman, K. F. 1912, Acta Math., 36, 105 [CrossRef] [Google Scholar]
 Tremaine, S., Richstone, D. O., Byun, Y.I., et al. 1994, AJ, 107, 634 [NASA ADS] [CrossRef] [Google Scholar]
 van Leeuwen, F. 2009, A&A, 497, 209 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yoshida, H. 1990, Phys. Lett. A, 150, 262 [Google Scholar]
 Zacharias, N., Finch, C., Girard, T., et al. 2010, AJ, 139, 2184 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: nbody6tidgpu
The program nbody6tidgpu is a new direct Nbody code for integrating the Nbody problem in an analytic background potential using Graphics Processing Units (GPUs). The program is essentially based on nbody6 by Aarseth (Aarseth 1999, 2003). nbody6tidgpu is not parallelized with MPI routines but instead runs on multicore processors (via OpenMP) with one or two GPU(s). For this purpose, special routines in CUDA have been added by Nitadori in collaboration with Aarseth. The variant of nbody6 with GPU support is termed nbody6gpu (Nitadori & Aarseth 2011).
Implementing different galactic tidal fields into nbody6gpu is the contribution of some of the present authors and the special feature of nbody6tidgpu. This implementation is based on the implementation in the (parallel) code nbody6gc which is described in detail in Ernst et al. (2009) and Ernst (2009). Several analytic models for the background potential are part of the code, for example,

the threecomponent PlummerKuzmin model of theMilky Way described in Sect. 2(Miyamoto & Nagai 1975);

powerlaw models with and without a central black hole;

Plummer model (Plummer 1911).
In addition to the equations of motion of the Nbody problem, the program solves the equations of motion of the star cluster orbit around the Galactic center with an eighthorder composition scheme (Yoshida 1990; for the coefficients see McLachlan 1995). The tidal force of the background potential acts on all particles in the Nbody system and is added as a perturbation to the KS regularization (Kustaanheimo & Stiefel 1965) of nbody6.
The program is able to handle the integration of star cluster orbits with almost all eccentricities. For nearly radial orbits a time transformation can be switched on to guarantee an exact orbit integration during the pericenter passages. The leapfroglike symplectic integration schemes allow for adaptive time steps if one applies a Sundman transformation to the time (Sundman 1912; Mikkola & Tanikawa 1999; Preto & Tremaine 1999; Mikkola & Aarseth 2002).
It is also possible to switch on dynamical friction with different χ functions and different realizations of the Coulomb logarithm, i.e. fixed or variable according to Just & Peñarrubia (2005). Since the symplectic composition schemes are by construction suited to Hamiltonian systems, the dissipative dynamical friction force requires special attention. It is implemented in nbody6tidgpu with an iterative implicit midpoint method (see, e.g., Mikkola & Aarseth 2002).
For the current study, the threecomponent PlummerKuzmin model described in Sect. 2 is used as a background potential. The acceleration for the PlummerKuzmin model is given by the expression (A.1)For the tidal correction to the equations of motion of the Nbody problem, the time derivative of acceleration (jerk) for the analytic background potential is needed to guarantee energy conservation and exactness of the orbit integration on a high level. The jerk is given by lengthy expressions that depend also on the velocity components V_{x},V_{y}, and V_{z} and that we do not reproduce here. However, they are also implemented in nbody6tidgpu.
Appendix B: Bad fits
Fig. B.1 Comparison of the cumulative mass profiles of observations and ensemble of models. The thick solid (red) lines in the middle for the model represent the mean that is ensembleaveraged over 15 runs. The standard error of the mean for the ensemble is shown as a filled area. The dashed (blue) line is calculated from the observational data. Top: ensemble of models en3125W6. Bottom: ensemble of models en3375W9. 
Figure B.1 shows the two best fits for the total mass of Hyades members within 30 pc from the cluster center for the models within the second parameter grid of Sect. 4, i.e. the grid with ΔN = 125. The top panel shows the ensemble of models en3125W6 (N_{0} = 3125,W_{0} = 6) while the bottom panel shows the ensemble of models en3375W9 (N_{0} = 3375,W_{0} = 9). It can be seen that the agreement between observations and ensemble of models is not even consistent within 2σ_{m} for the largest part of the radial range.
All Tables
3D velocity dispersions at t = 625 Myr for the ensemble en3000W6 for stars 30, 9 and 3 pc of the cluster center.
Number counts of stellar types for the bestfitting ensemble of models (en3000W6).
Number counts of stellar types that are the component of a binary for the ensemble enb2125PL.
All Figures
Fig. 1 Rotation curve (at z = 0) of the threecomponent PlummerKuzmin model of the Milky Way. 

In the text 
Fig. 2 “Artistic” RGBA image of the (logarithmic) density structure of the threecomponent PlummerKuzmin model. The three components are overblended with an alpha channel with 90% transparency. The halo is plotted in blue, the disk in green, and the bulge in red. 

In the text 
Fig. 3 3D orbit of the Hyades in the threecomponent PlummerKuzmin model of the Milky Way. The integration time is 625 Myr. 

In the text 
Fig. 4 Comparison of the presentday cumulative mass profiles of observations and models. The thick solid (red) lines in the middle of the models represent the mean, which is ensembleaveraged over 15 runs. The standard error is shown as a filled area. The dashed (blue) line is calculated from the observational data. Top: the ensemble en3000W6. Middle: the ensemble en3000W9. Bottom: the ensemble en3875W12. 

In the text 
Fig. 5 Particle number within r = 30 pc from the cluster center as a function of time for all 15 runs of the ensemble en3000W6. One can see the nonMarkovian process of star cluster dissolution. The mean final particle number and the standard deviation for a single run is shown as a (red) dot with vertical errorbars at t = 625 Myr. The future evolution of the mean particle number within r = 30 pc is shown in the inset. There the times when the cluster has reached 20, 10 and 5% of its initial particle number are shown as (red) dots. 

In the text 
Fig. 6 The presentday state of one run of the ensemble of models en3000W6 after 625 Myr of evolution. The tidal tails have reached a length of 800 pc. 

In the text 
Fig. 7 The average realization of the Kroupa (2001) IMF of the ensemble en3000W6 with lowest mass m_{l} = 0.08 M_{⊙}, highest mass m_{h} = 100.00 M_{⊙}, and initial particle number N_{0} = 3000. Shown is the mean value of log (N) from the ensemble averaging over 15 runs. The 2σ_{m} standard errors are shown. The analytical slopes and the transition at 0.5 M_{⊙} are also shown as lines. 

In the text 
Fig. 8 Comparison of the PDMFs between observations and ensemble of models en3000W6. Top panel: stars within a radius of 30 pc of the center; middle panel: stars within a radius of 9 pc of the center (Jacobi radius); bottom panel: stars within a radius of 3 pc of the center (core). Shown is the mean value of log (N) from the ensemble averaging over 15 runs. The errors are 1σ_{m} standard errors. The blue bars are derived from the observational data with N = 724 stars and lowest mass m_{l} = 0.116 M_{⊙}, while the red bars are derived from the ensemble of models with ⟨ N_{f} ⟩ = 736 stars and lowest mass m_{l} = 0.08 M_{⊙}. The onset of incompleteness of the observations is shown as a vertical line in the top panel. 

In the text 
Fig. 9 Presentday cumulative mass function for the ensemble en3000W6 for r < 9 pc from the cluster center. Shown is the sum of all masses which are higher than a given mass. The errors are 1σ_{m} standard errors. The completeness limit of the observations is shown as a vertical solid line. 

In the text 
Fig. 10 Presentday synthetic HertzsprungRussell diagram for one run of the evolved model en3000W6. Shown are stars within a radius of r = 30 pc from the cluster center. The final particle number within r = 30 pc for this run is N_{f} = 836. 

In the text 
Fig. 11 PDLFs in the K band for the ensemble en3000W6 and the observations (Röser et al. 2011). For explanations see the text. Shown is the mean/measured value of log (N) as a function of the absolute K_{s} (observations)/K (model) band magnitude. The errors are 1σ_{m} standard errors from the ensemble averaging. 

In the text 
Fig. 12 Average cumulative mass for observations and ensembles en3000W6 and enb2125PL. The average cumulative mass is calculated only for stars with mass m ≥ 0.116 M_{⊙}. 

In the text 
Fig. 13 Comparison of presentday mass segregation between models and observations. Top figure: analysis of the evolved (for 625 Myr) ensemble en3000W6. Bottom figure: analysis of the observational data of Röser et al. (2011). Each single plot contains on the lefthand side the “cumulative” Γ_{MST} (points) for the 5 (red), 10 (green), 20 (blue), 50 (magenta), 100 (cyan), 200 (brown), 500 (orange), and all most massive stars (black). The righthand side shows the “differential” Γ_{MST} (points) for the 5 (red), 6 − 10 (green), 11 − 20 (blue), 21 − 50 (magenta), 51 − 100 (cyan), 101 − 200 (brown), 201 − 500 (orange), and 501all most massive stars (black). The error bars mark the 1, 2, and 3σ uncertainties. A value of one marks the unsegregated state. 

In the text 
Fig. 14 Comparison of the presentday cumulative mass profiles of observations and models with primordial binaries. The thick solid (red) lines are derived from the ensemble en2250PL without binaries, the dashed (green) line is derived from the ensemble enb2250PL with 33% primordial binaries, and the dotted (magenta) line is derived from the bestfitting ensemble with 33% primordial binaries enb2125PL. The dashed (blue) line represents the observations. The standard errors are shown as filled areas. 

In the text 
Fig. 15 PDLFs in the K band for the ensemble enb2125PL where all binaries are either resolved or unresolved. Top panel: stars within a radius of 30 pc of the center; middle panel: stars within a radius of 9 pc of the center (Jacobi radius); bottom panel: stars within a radius of 3 pc of the center (core). Shown is the mean value of log (N) as a function of the absolute K band magnitude. 

In the text 
Fig. 16 Comparison of the PDMFs between ensembles enb2250PL, en2250PL, and enb2125PL. Top panel: stars within a radius of 30 pc of the center (Jacobi radius); middle panel: stars within a radius of 9 pc of the center (Jacobi radius); bottom panel: stars within a radius of 3 pc of the center (core). Shown is the mean value of log (N) from the ensemble averaging over 15 runs. The errors are 1σ_{m} standard errors. 

In the text 
Fig. B.1 Comparison of the cumulative mass profiles of observations and ensemble of models. The thick solid (red) lines in the middle for the model represent the mean that is ensembleaveraged over 15 runs. The standard error of the mean for the ensemble is shown as a filled area. The dashed (blue) line is calculated from the observational data. Top: ensemble of models en3125W6. Bottom: ensemble of models en3375W9. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.