© ESO, 2011

1. Introduction

The Hyades (“rain stars”) are one of the best-known 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 Hertzsprung-Russell (HR) diagram with a narrow and well-defined main sequence has been recently studied by von Leeuwen (2009). Previous N-body 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_⊙¹. The membership has been determined with the convergent-point method. Given the right ascension, declination, and the respective proper motions, the convergent-point 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 mass-to-luminosity 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₀ = 8 kpc (values for R₀ = 8.5 kpc are given in brackets) and that β = κ/Ω = 1.37, where κ and Ω are the epicyclic and circular frequencies of a near-circular 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).

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

We solve the N-body problem of the Hyades in an analytic background potential of the Milky Way. For the background Milky Way potential, we use an axisymmetric three-component model, where the bulge, disk, and halo are described by Plummer-Kuzmin 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 three-component model of the Milky Way. As in our previous works, the parameters of the three-component model are chosen such that the rotation curve matches that of the Milky Way (Dauphole & Colin 1995). At the solar radius R₀ = 8.0 kpc, which we have assumed in this study, the value of the circular velocity is V₀ = 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 three-component Plummer-Kuzmin 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 N-body problem in the analytic background potential of the Milky Way, the new direct N-body 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 Hertzsprung-Russell diagram. The implementation of stellar evolution for individual stars of the N-body 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 N-body 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 N-body 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 half-mass 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

We integrated the center of mass orbit of the Hyades backwards in time in the analytic Milky Way potential (i.e. the three-component Plummer-Kuzmin 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 N-body models up to the present time.

The present-day 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₀,z₀) = (8.00,0.02) kpc (consistent with Piskunov et al. 2006).

The present-day 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₀ = 8 kpc in the three-component Plummer-Kuzmin 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 z-direction 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 best-fitting 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 present-day 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 low-mass end according to Eq. (9) below,  ≈ 100 stars would be missing in the low-mass 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 Roche-lobe 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 near-circular 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₁ − c₄ 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 α₀ − α₃ 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 dual-Xeon 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₀, N), where W₀ 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₀ = 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₀,Δ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₀. The variation in W₀ 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₀,N) for the best-fitting 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₀ to a grid cell height of ΔW₀ = 3, i.e. we restricted the investigation to the models with W₀ = 3,6,9,12. Thus the grid cell size of the second grid in the parameter space was given by (ΔW₀,Δ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 best-fitting model in the first grid for a given W₀ and varied the particle number N in steps of ΔN = 125 particles in order to find the best-fitting 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 ensemble-average 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 ensemble-averaged 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 best-fitting 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₀ = 3000,W₀ = 6), i.e. the letters “en” for “ensemble” are followed by the initial particle number N₀ and the King parameter W₀.

We also looked for the best-fitting ensemble of single-star 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₀. The Plummer model runs are discussed in Sect. 5.2.4.

We ran the best-fitting ensemble of single-star 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 period-generating 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 eccentricity-generating function e² = 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₀ = 3000,W₀ = 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 best-fitting 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.

The second-best fitting ensemble of King models (en3000W9; N₀ = 3000,W₀ = 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₀ = 3875,W₀ = 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₀ = 12. For W₀ = 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 present-day 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 present-day state

The best-fitting 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 non-Markovian 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 present-day 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 kicked-out 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).

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₀ = 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).

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 low-mass 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 mass-luminosity 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).

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

Table 5 shows the number counts of stellar types for the best-fitting 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 ensemble-averaging 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 carbon-oxygen 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 present-day synthetic Hertzsprung-Russell 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 present-day 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.

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 ₁₀(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.

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 equal-sized 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σ.

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 left-hand side the “cumulative” Γ_MST (points) for the 5, 10, 20, 50, 100, 200, 500, and all most massive stars (black). On the right-hand side the “differential” Γ_MST (points) is shown for the 5, 6 − 10, 11 − 20, 21 − 50, 51 − 100, 101 − 200, 201 − 500, and 501-all 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) mass-segregated 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

The recent work of Sollima et al. (2009) showed that open clusters may contain between f_b ≈ 35 − 70% binaries. We have run the best-fitting 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 best-fitting model enb2125PL (also with 33% primordial binaries).

Figure 14 shows the comparison of the present-day 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 half-mass 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 K-band 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 high-mass 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 low-mass stars compared to the ensemble en2250PL, while high-mass 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.

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₂₀ = 695 Myr, t₁₀ = 875 Myr, and t₅ = 1029 Myr. The last time may be taken as an approximate dissolution time.

6. Discussion

The present-day 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.

• Roche-lobe filling W₀ = 6−9 single-star King initial models (King 1966) with the initial particle number N₀ = 3000 provide a good fit to the present-day inner cumulative mass profile.

• The best-fit, single-star, Roche-lobe filling King model has an initial mass of M₀ = 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 half-mass to Jacobi radius of r_h/r_J = 0.18. Here the initial particle number is N₀ = 2125 and mass are is M₀ = 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 carbon-oxygen 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₀ by  ≈ 5% as compared to the model without primordial binaries.

• Under the assumption that the single-star 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 good-fitting models can be quite high due to the large dimension of the parameter space, which is spanned by the initial particle number N₀, the King parameter W₀, and two extra dimensions: the primordial binary fraction f_b and the Roche-lobe 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 Roche-lobe 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 041-043 and I/84 678-680 of the Volkswagen Foundation and 823.219-439/30 and /36 of the Ministry of Science, Research and the Arts of Baden-Württemberg). C.O. appreciates funding by the German Research Foundation (DFG), grant OL 350/1-1. 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

Appendix A: nbody6tidgpu

The program nbody6tidgpu is a new direct N-body code for integrating the N-body 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 three-component Plummer-Kuzmin model of theMilky Way described in Sect. 2(Miyamoto & Nagai 1975);

• power-law models with and without a central black hole;

• Dehnen models (Dehnen 1993; Tremaine et al. 1994);

• Plummer model (Plummer 1911).

In addition to the equations of motion of the N-body problem, the program solves the equations of motion of the star cluster orbit around the Galactic center with an eighth-order composition scheme (Yoshida 1990; for the coefficients see McLachlan 1995). The tidal force of the background potential acts on all particles in the N-body 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 leapfrog-like 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 three-component Plummer-Kuzmin model described in Sect. 2 is used as a background potential. The acceleration for the Plummer-Kuzmin model is given by the expression (A.1)For the tidal correction to the equations of motion of the N-body 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

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₀ = 3125,W₀ = 6) while the bottom panel shows the ensemble of models en3375W9 (N₀ = 3375,W₀ = 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

All Figures