text:
# Hyades dynamics from N-body simulations: Accuracy of astrometric radial velocities from Hipparcos Based on observations by the ESA Hipparcos satellite, and on the 𝑁-body code NBODY6 by Sverre Aarseth, publicly available at ftp://ftp.ast.cam.ac.uk/pub/sverre/

## 1 Introduction

The Hyades is the nearest rich open cluster and as such has played a fundamental role in astronomy as a first step on the cosmic distance ladder and as a test case for theoretical models of stellar interiors (Lebreton lebreton00 (2000)). From the first use of the converging point method by Boss (boss08 (1908)) up to the use of pre-Hipparcos trigonometric parallaxes by van Altena et al. (altena97 (1997)), an important goal in astrometry has been the determination of an accurate distance to the cluster. With the advent of the Hipparcos Catalogue (ESA esa (1997)) the Hyades lost its unique status for distance calibration, but as the depth and internal velocity field of the cluster were well resolved by Hipparcos, focus could instead be turned to its three-dimensional structure and kinematics (Perryman et al. perryman98 (1998)). A deeper understanding of the dynamics and evolution of the cluster should now be possible through detailed comparison with $`N`$-body simulations.

Thanks to the accurate Hipparcos measurements, the Hyades has recently acquired a completely new role as a practical standard in observational astrophysics: it is one of very few objects outside the solar system for which the accurate radial motion can be determined by geometric means, i.e. without using the spectroscopic Doppler effect. From a combination of Hipparcos parallaxes and proper motions, Madsen et al. (madsen02 (2002)) obtained “astrometric radial velocities” for individual Hyades stars with a then estimated standard error of about 0.6 km s-1. Currently the Hyades is the only cluster for which astrometric radial velocities are derived with individual accuracies better than 1 km s-1, but the technique may be extended to many more objects with future space astrometry missions (Dravins et al. 1999b ).

Astrometric radial velocities are important mainly because they make it possible to determine the absolute lineshifts intrinsic to the stars, through comparison with spectroscopic measurements. Such lineshifts are caused for instance by convective motions and gravitational redshift in the stellar atmospheres (Dravins et al. 1999a ). Absolute lineshifts could previously only be observed in the solar spectrum, but are now within reach for a range of spectral types through the use of astrometric radial velocities. The present paper is part of a research programme at Lund Observatory in which absolute lineshifts are determined and used as a diagnostic tool in stellar astrophysics (Dravins et al. dravins97 (1997), 1999b ; Lindegren et al. lindegren00 (2000); Madsen et al. madsen02 (2002); Gullberg & Lindegren gullberg02 (2002)).

A major uncertainty in the astrometric radial velocities originates in the internal velocity dispersion of the cluster, which limits both the accuracy of the cluster motion as a whole, and that of the individual stars. A primary goal of the present investigation is to find out whether a better understanding of the internal velocity structure of the cluster, obtained through $`N`$-body calculations, can be used to improve the accuracy of the astrometric radial velocities.

Sect. 2 briefly recalls the kinematic information, including astrometric radial velocities, that can be derived from Hipparcos data. Sect. 3 describes the model used to simulate the evolution of the cluster up to its present state, and its subsequent observation, as well as the main properties derived from the simulations. Implications for the accuracy of the astrometric radial velocities are discussed in Sect. 4, followed by a discussion of non-modelled effects in Sect. 5, and conclusions.

## 2 Cluster kinematics derived from astrometry

Since an ultimate aim of the present programme is to confront spectroscopic measurements of line shifts in stellar spectra with independent measurements of the stellar motions, it is essential that the kinematic data, including the radial velocities, are derived without using the spectroscopic Doppler effect. Dravins et al. (1999b ) describe several methods to derive the radial motion of stars by purely geometric means, i.e. using astrometric data. Of these, the moving-cluster method has been successfully applied to several open clusters and OB associations, in particular the Hyades (Lindegren et al. lindegren00 (2000); Madsen et al. madsen02 (2002)). The principle of the moving-cluster method is very simple: let $`\theta `$ be the angular size of the cluster and $`R`$ its distance. Assuming its linear size $`R\theta `$ to be constant, we have $`\dot{R}\theta +R\dot{\theta }=0`$, where the dot signifies time derivative. Since $`R`$ is known from trigonometric parallaxes, the astrometric radial velocity of the cluster follows as $`\dot{R}=R\dot{\theta }/\theta `$.

In practice, several kinematic parameters are simultaneously estimated from the astrometric data of the cluster member stars, using the method of maximum likelihood (Lindegren et al. lindegren00 (2000)). Some features of the method, relevant for the subsequent discussion, are recalled hereafter.

The estimated parameters include the common space velocity of the cluster ($`v_0`$), the individual stellar parallaxes ($`\pi _i`$ for star $`i`$), and the internal velocity dispersion ($`\sigma _v`$). The astrometric radial velocity of an individual star $`i`$ is then calculated as $`\widehat{v}_{\mathrm{r}i}=r_i^{}\widehat{v}_0`$, where $`r_i`$ is the unit vector towards the star and the caret $`\widehat{}`$ signifies estimated quantities. As part of the procedure, improved parallaxes $`\widehat{\pi }_i`$ are obtained for the individual stars. In the Hyades, these are 2–5 times more precise than the original Hipparcos parallaxes which have errors around 1–1.5 mas. The improvement results from a combination of trigonometric and kinematic parallaxes, where the latter follow from the proper-motion components *along* the cluster motion, which are inversely proportional to distance. The kinematically improved parallaxes allow a very precise mapping of the spatial structure of the cluster. The maximum likelihood estimate of $`\sigma _v`$ is unfortunately biased. Instead the proper motions *perpendicular* to the cluster motion are used to estimate the velocity dispersion according to the method described in Lindegren et al. (lindegren00 (2000)), Appendix A.4. For each star, a goodness-of-fit statistic $`g_i`$ is also obtained from the maximum-likelihood estimation (see Lindegren et al. (lindegren00 (2000)) for a thorough discussion of $`g_i`$). The statistic is primarily used to reject stars whose astrometric data do not fit the cluster model well enough; a rejection limit of $`g_{\mathrm{lim}}=15`$ was normally used, although a stricter limit (10) or no limit at all ($`\mathrm{}`$) were also tried. For the retained stars, the $`g_i`$ values (which are then $`g_{\mathrm{lim}}`$) could be regarded as a quality index, with a lower value indicating a better fit to the cluster model.

The error in the estimated astrometric radial velocity, $`\widehat{v}_{\mathrm{r}i}`$, has two parts. The first part is due to the error in the common space motion of the cluster, $`\widehat{v}_0`$. Its uncertainty depends on global properties of the cluster such as its distance, angular extent, and richness, as well as on the accuracy of the astrometric data. The second part is due to the star’s peculiar motion relative to the cluster centroid. This part depends only on the dispersion of the peculiar motions along the line of sight, which for a uniform, isotropic velocity dispersion equals $`\sigma _v`$. In most of the clusters for which the method has been applied, the main uncertainty comes from the first part, i.e. the error in the cluster’s space motion. In the Hyades, however, the uncertainty in $`\widehat{v}_0`$ is small enough (0.36 km s-1 along the line-of-sight; Madsen et al. madsen02 (2002)) that the total uncertainty in the astrometric radial velocities is dominated by the contribution from the internal velocity dispersion (0.49 km s-1 according to the estimate in the same source).

On the other hand, the assumption of a constant and isotropic velocity dispersion throughout the cluster may be rather simplistic. Theoretically, one expects at least a variation with distance $`r`$ from the centre of the cluster, and possibly also a variation with stellar mass due to the equipartition of kinetic energy. For instance, in a simple Plummer (plummer15 (1915)) potential we have

$$\sigma _v^2(r)=\frac{GM}{6\sqrt{r_c^2+r^2}}$$
(1)
(Gunn et al. gunn88 (1988); Spitzer spitzer87 (1987)), where $`M`$ is the cluster mass and $`r_c`$ the core radius ($`3`$ pc for the Hyades). According to Eq. (1), $`\sigma _v`$ should decrease by one third as one moves two core radii away from the centre, and become even smaller further out in the cluster; but this trend is obviously broken at some distance by tidal forces. Clearly, these effects must be also reflected in the accuracy of the astrometric radial velocities. Attempts to measure the radial variation of dispersion in the Hyades from astrometry were inconclusive (Madsen et al. madsen01 (2001)), but $`N`$-body simulations could help to establish to what extent such variations exist in real clusters.

## 3 Dynamical simulation of the Hyades cluster

### 3.1 Previous N-body simulations

It is not new to use the Hyades as a comparison with $`N`$-body simulations. Aarseth (aarseth77 (1977)) discussed the dynamical relevance of the central binary 80 Tau (HIP 20995) in the context of binary formation and evolution in stellar systems as described by $`N`$-body simulations. Oort (oort79 (1979)) discussed the flattening of the Hyades parallel to the galactic plane by comparing observations with the $`N`$-body simulations by Aarseth (aarseth73 (1973)). Kroupa (1995c ) simulated the evolution of star clusters and found excellent agreement between the models and the Hyades luminosity function, concluding that the initial conditions of the cluster could to a large extent be reconstructed. An initial mass of the Hyades protocluster of some $`1300M_{}`$ was suggested. von Hippel (hippel98 (1998)) used numerical simulations of clusters and data on Hyades white dwarfs, among others, to conclude that the white-dwarf mass fraction is relatively insensitive to kinematic evolution. Portegies Zwart et al. (portegies01 (2001)) discussed the evolution of star clusters which were given initial conditions to represent open clusters, including the Hyades. A good model fit to the Hyades was obtained, thus illustrating the possibility to estimate the initial conditions for an observed star cluster.

What is new in the present study is that the three-dimensional kinematics of the Hyades is investigated through a direct comparison of the Hipparcos observations with a realistic $`N`$-body model, evolved till the present age of the cluster, as well as the objective to estimate the accuracy of the astrometric radial velocities from such a comparison.

### 3.2 Basic cluster data

Perryman et al. (perryman98 (1998)) made a detailed study on the Hyades based on Hipparcos data and a compilation of spectroscopic radial velocities from the literature. They identified 197 probable member stars, which constitute the initial Hyades sample (Hy0) used for the present study. When comparing with the simulated cluster, only stars within 20 pc from the cluster centre are considered, due to the radial limitation in the $`N`$-body code (Sect. 3.3). Adopting the cluster centre of mass in equatorial coordinates, $`(+17.36,+40.87,+13.30)`$ pc from Perryman et al. (perryman98 (1998)), and using the kinematically improved parallaxes (Sect. 2), a subset of 178 stars (Hy0r) was found within a radius of 20 pc. The cluster has a general space velocity of $`(5.90,+45.65,+5.56)`$ km s-1 in equatorial coordinates (Madsen et al. madsen02 (2002)).

Perryman et al. (perryman98 (1998)) note that a redetermination of membership with the above cited centre of mass will reduce the number of member stars outside 10 pc by 10 stars while keeping the same number of stars inside 10 pc. The true number of member stars in the Hy0r sample is then probably smaller than the 178 stars.

Hipparcos is nominally complete to $`V7.3+1.1|\mathrm{sin}b|`$ for spectral types later than G5 (or $`BV>0.8`$). However, it is known that the actual limit is somewhat fuzzy, due to photometric errors and other complications. Therefore, a conservative completeness limit of $`V7`$ mag is assumed for this study. Choosing a fainter completeness limit like e.g. $`V8`$ mag will, however, not significantly affect the outcome of the simulations as it will be shown later (Table 1). The actual number of fainter Hyades members is not known. However, at least seven single white dwarfs have been found (e.g. Reid reid96 (1996)), and this number can also be used as a constraint on the model.

Perryman et al. (perryman98 (1998)) estimated the cluster age to be $`625\pm 50`$ Myr, and this age is what is assumed in the following. It should be mentioned that in a more recent work by Lebreton et al. (lebreton01 (2001)), based on kinematically improved parallaxes from Dravins et al. (dravins97 (1997)), only an upper limit of 650 Myr could be estimated due to the lack of a clear turn-off point (cf. top diagram in Fig. 1). In the same work they also estimated the metallicity to $`\text{[Fe/H]}=0.14\pm 0.05`$ dex. The interstellar extinction is negligible: Taylor (taylor80 (1980)) found only a very small colour excess $`E(BV)=0.003\pm 0.002`$ mag.

From various studies, a large fraction of the stars are known to be binaries. In the compilation by Perryman et al. (perryman98 (1998)), 75 of the 197 probable member stars were either identified as binaries in the Hipparcos Catalogue or previously known as spectroscopic binaries (their Table 2). Patience et al. (patience98 (1998)) found three new binaries from a speckle imaging survey of Hyades members, plus one marked as binary in the Hipparcos Input Catalogue (HIC; Turon et al. turon92 (1992)), but not found by Hipparcos. In the Tycho Double Star Catalogue (Fabricius et al. fabricius02 (2002)), an additional 21 binaries were identified. The eclipsing binary system HIP 17962 = V471 Tau (e.g. Werner & Rauch werner97 (1997), and references therein) must also be included in the list of Hyades binaries. We thus end up with 101 known binaries in the Hy0 sample, yielding a minimum multiplicity of 0.51 companions per primary. For the Hy0r sample (within 20 pc of the cluster centre) the minimum multiplicity is 0.53. To include some more binary statistics, binaries with periods $`P<10`$ days have been taken from the compilations on the open–cluster database WEBDA111available at http://obswww.unige.ch/webda/.

The above values of the multiplicity are only lower limits to the true multiplicity, because of the difficulty to detect binaries in some intervals of separation $`\rho `$ (or period $`P`$) and magnitude difference $`\mathrm{\Delta }m`$ (or mass ratio $`q`$). In restricted intervals, the searches can however be considered complete. For instance, Hipparcos probably detected practically all binaries with $`0.2<\rho <2`$ arcsec and $`\mathrm{\Delta }m<2`$; cf. Fig. 3.2.106 in Vol. 1 of (ESA esa (1997)), where 17 are found in Hy0r. Patience et al. (patience98 (1998)) observed a high fraction of Hyades stars that were also observed by Hipparcos. The 17 binaries they found with $`0.1<\rho <1.07`$ arcsec and $`q0.4`$ must therefore also be regarded as a nearly complete sample.

Hipparcos effectively observed for about 37 months ($``$3 years) spread over a period of nearly 4 years. This means that the proper motions of binaries may be significantly affected by the orbital motion of the photocentre, which must be taken into account when simulating the Hyades proper motions (Sect. 3.4). In order to reduce this effect in the observational analysis, proper motions from the Tycho-2 catalogue (Høg et al. hog00 (2000)) have also been used, where available. In the solution for the cluster kinematics, the Tycho-2 proper motions yield slightly, but systematically smaller radial velocities ($`v_r(`$HIP$`)v_r(`$Tycho-2$`)=+`$0.9 km s-1) than do the Hipparcos data for the $`g_{\mathrm{lim}}=15`$ sample (Madsen et al. madsen02 (2002)), which can be explained by the mean difference of -0.4 mas yr-1 of the proper motions in right ascension of what was considered the best sample. In declination, the mean difference of the proper motions is 0.0 mas yr-1. Although the expected deviations between the Hipparcos and Tycho-2 Catalogues are generally under 0.5 mas (Urban et al. urban00 (2000)), the result from the Hyades might reflect some subtle bias in the Tycho-2 proper-motion system. Since the Tycho-2 system of proper motions was effectively calibrated onto the Hipparcos system, greater confidence should be put on the solution based on the Hipparcos data. The Tycho-2 data should therefore only be used to study the internal velocity structure of the cluster, where a possible bias is not important.

In addition to the Hy0r sample (which thus includes all 178 probable members within a radius of 20 pc from the cluster centroid), the following samples are also discussed: Ty0r, which is the same as Hy0r but with proper motions from Tycho-2 replacing those in the Hipparcos Catalogue; Hy1r, which is the subset of 85 stars in Hy0r for which there is no known indication of multiplicity; and Ty1r, which is the same as Hy1r but with Tycho-2 proper motions.

It has been suggested that there might be systematic errors in the Hipparcos parallaxes for at least some open clusters (Pinsonneault et al. pins98 (1998)). The discussion shall not be repeated here, but it should just be stated that there is a general consensus that the mean Hyades parallax is not affected by any correlation between positions and parallaxes (Narayanan & Gould narayanan99 (1999); van Leeuwen vanleeuwen00 (2000); Lindegren et al. lindegren00 (2000)). This problem, if it exists, has been neglected in the simulations.

### 3.3 N-body model of the Hyades cluster

The dynamical evolution of a Hyades-type open star cluster was simulated using the well-known $`N`$-body code NBODY6 (Aarseth aarseth99 (1999), aarseth00 (2000)). The code incorporates algorithms to deal with stellar (including binary) encounters (Mikkola & Aarseth mik93 (1993), mik96 (1996), mik98 (1998)) and stellar evolution (Hurley et al. hurley00 (2000)). For the present study, no modifications were made to the code. Some of the non-modelled effects are discussed in Sect. 5.

External perturbations are represented by a fixed, galactic tidal field. The cluster is assumed to move in a circular orbit at the present distance of the Sun from the galactic centre. The angular velocity is $`\mathrm{\Omega }=AB`$, where $`A=14.4`$ km s-1 kpc-1 and $`B=12.0`$ km s-1 kpc-1 are Oort’s constants. This gives rise to tidal forces plus a Coriolis force (cf. Chandrasekhar chandra42 (1942), Ch. 5.5).

To set up the initial cluster configuration, stars are randomly picked from the initial mass function (IMF) described by Kroupa et al. (kroupa93 (1993)), until the required total particle number has been reached. Binaries are included as described below. Stars are initially deployed randomly in a Plummer potential (Plummer plummer15 (1915); Spitzer spitzer87 (1987)) with virial radius $`r_v=4`$ pc. During the evolution of the cluster, stars are kept in the simulation as long as they are within two tidal radii ($`21`$$`23`$ pc). The simulation is run until the cluster reaches an age of 625 Myr.

The reason for choosing one single age was to have a fixed parameter for comparing different model realisations. The age uncertainty is not important regarding the conclusions about the current dynamics since the cluster has been relaxed for quite a while.

Binaries are generated by randomly pairing stars picked from the IMF. This gives an almost uniform distribution in the logarithm of the mass ratio ($`\mathrm{log}q`$), i.e. a strong preference for small $`q`$, similar to what has been observed for G-dwarf systems (Duquennoy & Major duquennoy91 (1991)). The semimajor axis ($`a`$) is selected from a uniform distribution in $`\mathrm{log}a`$ with an upper cut-off at 3000 AU (Quist & Lindegren quist00 (2000)). The period distribution is afterwards generated by NBODY6 based on the modelling by Kroupa (1995a , 1995b ) with minimum period 1 day, and binaries merged if $`a10R_{}`$. The initial distribution of eccentricities $`e`$ is assumed to be thermal, i.e. with a probability density function $`2e`$ (Kroupa 1995b ).

The only free model parameters are thus the total particle number and the initial binary fraction (or multiplicity). Their determination is discussed in Sect. 3.5.

### 3.4 Transformation to observables

From NBODY6, the luminosity and temperature is obtained for each star. These parameters are transformed to the observational plane $`(BV,M_V)`$ using Kurucz’s colour tables (e.g., Kurucz kurucz79 (1979) and Buser & Kurucz buser92 (1992)) for $`\text{[Fe/H]}=0.10`$. Johnson’s $`V`$ is used instead of the Hipparcos magnitude $`Hp`$, because of the lack of adequate transformations for the latter. For binaries, the combined colour and magnitude are calculated and plotted in order to get results that are directly comparable with Hipparcos data. In view of the very small interstellar reddening (Sect. 3.2), $`E_{BV}=0.0`$ is assumed.

When comparing the simulated and observed HR diagrams it should be borne in mind that the theoretical models and colour transformations may produce non-negligible errors. Observed discrepancies for the Hyades amount to some 0.05 mag in $`BV`$ or 0.3 mag in $`M_V`$ in the cool end of the main sequence (Castellani et al. castellani01 (2001)). No (empiric) corrections for this effect have, however, been made.

In order to mimic the real Hyades cluster, as observed by Hipparcos, the simulated present-day cluster is “observed” from the same distance as the real Hyades and given the same centroid velocity relative the Sun. Small errors in the “observed” $`V`$ magnitudes (standard deviation 0.0015 mag) are introduced, and parallaxes and proper motions, including observational errors, are generated following the same procedure as in Lindegren et al. (lindegren00 (2000)). The simulated sample includes all stars brighter than the completeness limit $`V=7`$, plus a random selection of the fainter stars matching the real sample in the number of stars per magnitude interval. It is assumed that the Hyades stars in the Hipparcos Catalogue with $`V>7`$ mag are not subject to any selection effects, although it cannot be ruled out due to a sometimes impenetrable selection procedure of Hipparcos objects in open clusters (Mermilliod & Turon mer89 (1989)).

Binaries receive different treatments depending on the magnitude difference ($`\mathrm{\Delta }m`$), period ($`P`$), and angular separation ($`\rho `$), in order to simulate how they were treated in the Hipparcos data analysis (see Sect. 1.4.2 in the Hipparcos Catalogue). Here, $`\mathrm{\Delta }m=M_{V2}M_{V1}`$, where subscripts “1” and “2” refer to the primary and secondary components. For certain combinations of these parameters, Hipparcos effectively observed the motion of the photocentre of the system. In the remaining cases the centre of mass were observed. The former systems include those with $`P`$ 0.1 to 20 years and $`\rho `$ 10 mas, or $`P>`$ 10 years and $`\rho `$ 100 mas; the short-period binaries ($`P<`$ few months), which may deviate from a single-star solution (the “stochastic” (X) solution), although they may have been detected as binaries by Hipparcos; and systems with Hipparcos magnitude difference $`\mathrm{\Delta }Hp>4`$ mag. For these systems, the component velocities $`v_1`$ and $`v_2`$ are combined into a single velocity of the photocentre,

$$v_{\mathrm{ph}}=(1\xi )v_1+\xi v_2,$$
(2)
where

$$\xi =\frac{I_2}{I_1+I_2}=\frac{1}{1+10^{0.4\mathrm{\Delta }m}}$$
(3)
is the fractional intensity of the secondary. $`I`$ is the intensity for each component given by $`I10^{0.4M_V}`$. For these systems, a single proper motion is derived from $`v_{\mathrm{ph}}`$.

The resulting simulated astrometric data are subject to exactly the same maximum-likelihood estimation procedure as was used for the real cluster (Lindegren et al. lindegren00 (2000)). In particular, astrometric radial velocities and kinematically improved parallaxes are derived for the individual stars or binaries. The error in the estimated parallaxes is in the range 0.2–1.0 mas (an error of 0.5 mas corresponds to approximately 1 pc in the cluster centre). The improved parallaxes are used to compute distances from the cluster centre, which allow to count the number of stars within a certain radius. Furthermore, for any subsample of the stars, the velocity dispersion can be estimated from the proper-motion residuals (Sect. 2).

### 3.5 Model fitting

In order to tune the model parameters, it is necessary to make several simulations for the same parameter values but using different initialisations of the random number generator. The average of the different random realisations is then compared with the observational data, and the input parameters adjusted accordingly. The quantities to be compared are the radial distribution of the stars, their total number above a given magnitude limit, and binary statistics. Also the number of giants (defined as $`M_V<1`$ and $`BV>0.5`$) and the number of single white dwarfs are used to constrain the model.

The finally adopted (protocluster) model comprises 200 single stars and 1200 binaries, i.e. an initial multiplicity of 0.86 companions per primary. The total initial mass is $`1100`$$`1200M_{}`$. This is slightly less than previous estimates of $`1200`$$`1500M_{}`$ (Reid reid93 (1993)) or $`1300M_{}`$ (Kroupa 1995c ). This smaller initial mass was found necessary in order to correctly reproduce the number of observed stars with the given IMF. The true initial mass of the Hyades is probably higher due to non-modelled mass loss (Sect. 5). According to the simulations, the total current mass of the Hyades stars is $``$ 460 $`M_{}`$ with a tidal radius of $``$ 11 pc. Observationally, Reid (reid92 (1992)) made the estimation 410–480$`M_{}`$ while Perryman et al. (perryman98 (1998)) estimated $`400M_{}`$ in their Hipparcos study of the cluster.

An example of the observational HR diagram for one of the model realisations is shown in Fig. 1, together with the corresponding observed diagram for the Hyades cluster. In addition to the standard deviation introduced in $`V`$, a standard deviation of 0.01 in $`BV`$ is also introduced in the model HR diagram to make the colour distribution appear more realistic. This standard deviation includes both observational errors and the effects of peculiar stars, stellar rotation, etc. Apart from the previously mentioned possible discrepancy in the cool end of the main sequence, and the fact that the giant stars are too red in the simulations, the general agreement is reasonable. The precise colours of the giants are, however, irrelevant in the context of this study.

Table 1 shows some statistics computed from this model, after evolution to an age of 625 Myr and transformation to the observables, together with the corresponding observed numbers. From Table 1 it appears that the distribution of stars with radial distance and apparent magnitude in the Hyades is well reproduced by the model cluster. The number of stars decreases when we go from the constraints based on the true parallaxes to the constraints based on the estimated parallaxes, and the number decreases even further when we use the observed parallaxes. This is a result of observational errors affecting the parallaxes, and mostly for the smallest sphere $`r3`$. In fact, the resemblence in the three columns is so good that it shows the modelling of the errors are in accordance with reality. The underabundance of stars in the models relative to the observations in the range $`10<r<20`$ pc can be explained by an overestimation of stars outside 10 pc by Perryman et al. (perryman98 (1998)). They argued that using another centre of mass in the Hyades would lead to fewer stars in the halo (Sect. 3.2).

It has been much more difficult to reproduce the observed binary statistics (Table 2). Bright binaries with high mass ratio or small magnitude difference are underproduced. Even if every star in the protocluster were assumed to be a binary (multiplicity 1.0), the model would still predict too few binaries of these characteristics. The observed sample also has significantly more known short-period binaries ($`P<10`$ days) than obtained in the simulations. These discrepancies indicate that the model distributions in mass ratio and/or semi-major axis would need some adjustment. Alternatively, a higher initial mass leading to more binaries with the required properties could be an explanation assuming non-modelled mass loss of preferentially low mass stars. However, the discrepancies are not dramatic and for the present study it was preferred not to change the relevant code in NBODY6.

Since the initial multiplicity must be very high to fit the observed binary statistics without being in contradiction with the observed number of Hyades member stars, the degree of degeneracy between the two free input parameters (initial particle number and initial multiplicity) is small.

The simulations could in principle be “inverted” to derive an age, by for instance stopping the modelling when the realisations appear similar to observed structural or dynamical features in the Hyades. But the non-modelled effects leading to mass loss during the dynamical evolution will be a major uncertainty (Sect. 5).

### 3.6 Observed kinematics versus simulated data

#### 3.6.1 Dispersion versus cluster radius

In a Plummer potential, the velocity dispersion decreases with cluster radius according to Eq. (1). At some radius, however, the relation is expected to break down when the stars have left the cluster potential and become subject to the Galactic field. In the following this possible structure is investigated.

The various observed samples (Hy0r, Hy1r, Ty0r, Ty1r), as well as the different realisations of the adopted cluster model, are analysed by means of the maximum-likelihood method mentioned in Sect. 2. The samples are divided according to distance ($`r`$) from the cluster centroid in order to determine if there is a radial variation of the kinematics. The ranges in $`r`$ have not been chosen at random: 3 pc is approximately the core radius while 10 pc is approximately the tidal radius. Table 3 summarises the results for the number $`N`$ of stars (or systems) and the estimated velocity dispersion $`\widehat{\sigma }_v`$.

The analysis method includes the rejection procedure designed to “clean” the cluster membership described in Sect. 2 with the goodness-of-fit statistic $`g_i`$ calculated for each star. For the model simulations, no results are given for $`g_{\mathrm{lim}}=\mathrm{}`$ because of their sensitivity to run-away stars. In the observed sample such cases were already removed by Perryman et al. (perryman98 (1998)). It should be noted that the cleaning process successively reduces the estimated internal velocity dispersion, because the latter is based on the proper-motion residuals, which are also reflected in $`g_i`$. This is most clearly seen for the Hy0r sample at all radii, and for the other samples at $`r>10`$ pc. The reason that there seems to be more stars for e.g. Ty1r at $`g_{\mathrm{lim}}=15`$ than $`g_{\mathrm{lim}}=\mathrm{}`$ for certain ranges in $`r`$ is that kinematically improved parallaxes have been used to calculate the distance from cluster centre. Since it is a different solution for each $`g_{\mathrm{lim}}`$, the kinematically improved parallaxes may change slightly.

Kinematically, one cannot in general distinguish between actual non-member stars and member stars with a deviating space motion. The most probable reason for a member star not to follow the common space motion of the cluster is that it is a binary in a non-modelled orbit. As explained in Sect. 3.2, this effect should be greater for the samples based on the Hipparcos proper motions than when using the Tycho-2 data. Comparing the results for Hy0r and Ty0r as function of $`g_{\mathrm{lim}}`$ suggests that binaries are the main cause for deviating proper motions out to $`r10`$ pc, while for the greater radii they are partly caused by actual non-members.

The last two columns in Table 3 show the estimated and true dispersions from 20 realisations of the model. It appears that $`g_{\mathrm{lim}}=10`$ yields a correct estimation of the dispersion, while $`g_{\mathrm{lim}}=15`$ leads to an over-estimation of $`\sigma _v`$. Using $`g_{\mathrm{lim}}=10`$, the cluster as a whole (inside 20 pc) yields a dispersion of $`0.23\pm 0.02`$ km s-1, with no clear dependence on $`r`$. The model cluster yields a slightly larger value (0.30 km s-1) and shows a 20% decrease from the centre outwards. It should be noted that two of the 20 models yield estimated values as small as the observations ($``$ 0.23 km s-1). The dispersions $`\overline{\sigma }_v`$ characterise the stars in the simulated Hyades sample, and not the total number of stars in the cluster. Due to the limiting magnitude, stars with masses less than 0.5–0.6 $`M_{}`$ do not contribute to the velocity dispersions in the table, just as with the observations.

Madsen et al. (madsen01 (2001)) found some rather large radial variations of the velocity dispersion in the Hyades, but could not conclude whether the structure was real or not. From the present simulations it is concluded that the observed structure is probably spurious: similar variations (of either sign) can be seen in some of the model realisations, although they are absent in the average of the realisations.

Hitherto in studies of open clusters, only in the Pleiades has an indication of a relationship between $`r`$ and (the tangential component of) $`\sigma _v`$ been found (van Leeuwen vanleeuwen83 (1983)). In the globular cluster M15, however, a velocity dispersion decreasing from the centre out to 7 arcmin and then increasing was found by Drukier et al. (drukier98 (1998)). They interpreted it as an indication of heating of the outer part of the cluster by the galactic tidal field. But how the minimum at 7 arcmin was related to the tidal radius or other quantities remained unclear. Heggie (heggie01 (2001)) argued that heating might be an incorrect interpretation since the effect can also be seen in $`N`$-body simulations of star clusters moving under influence of a steady tidal field (cf. Giersz & Heggie giersz97 (1997)). In the models here, the same trend is seen, with a minimum in the $`r\sigma _v`$ relation just inside 10 pc (the mean tidal radius of the models is between 10 and 11 pc).

#### 3.6.2 Dispersion versus stellar mass

Theoretically we should also expect a decreasing velocity dispersion with higher mass, or correspondingly lower absolute magnitude, due to equipartition of kinetic energy. This should in turn lead to dynamical mass segregation, with the massive stars more concentrated to the centre of the cluster. This effect may have been seen in IC 2391 (Sagar & Bhatt sagar89 (1989)) and Praesepe (Holland et al. holland00 (2000)). Perryman et al. (perryman98 (1998)) found a clear mass segregation in the Hyades from the number density of stars in various mass groups as a function of distance from the centre. Direct searches by Lindegren et al. (lindegren00 (2000)) and Madsen et al. (madsen01 (2001)) for a relation between the observed velocity dispersion and mass (or absolute magnitude), did however prove inconclusive. Evidence of any equipartition of kinetic energy is best sought among the stars in the core of the cluster (Inagaki & Saslaw inagaki85 (1985)). For the present study, a limiting radius of 3 pc is therefore used. This is approximately the core radius of the Hyades.

In the Hipparcos Catalogue, often only the common absolute magnitude for a binary is available, and not the absolute magnitudes for both components. Since it is the mass that is interesting, only the samples without known binaries should be used, to ensure a reasonably unique correspondence between absolute magnitude and mass. In the simulated samples, binaries with a difference in absolute magnitude between the combined absolute magnitude of the two components in the binary and the the primary component of more than 0.1 mag have been removed. This is the simplest way to simulate the hy1r sample.

The remaining stars with $`r<3`$ pc in the hy1r sample are separated in four intervals of absolute magnitude, with divisions at $`M_V=2.1`$, $`3.4`$, and $`5.4`$ mag, approximately corresponding to the masses $`1.8`$, $`1.4`$, and $`1.0M_{}`$. The estimated dispersions in these intervals are $`0.17\pm 0.13`$, $`0.20\pm 0.11`$, $`0.24\pm 0.11`$ km s-1, and no solution for the last interval. The uncertainties are too large to allow any firm conclusion, although the expected trend is there. For comparison, the simulations gave an average dispersion going from 0.28 to 0.36 km s-1 in the same intervals.

#### 3.6.3 Other determinations of the dispersion

Several studies of the velocity dispersion of the Hyades have been performed during the years. In a detailed discussion by Gunn et al. (gunn88 (1988)), who performed a spectroscopic investigation of the cluster, a mean dispersion of 0.23 km s-1 was derived from a Plummer model. Their result agreed with the velocity dispersion obtained from the most precise spectroscopic radial velocities in their Hyades sample. However, it is important to note that the result of 0.23 km s-1 is dependent on the estimated $`M`$ and $`r_c`$, where the mass is the major uncertainty. Perryman et al. (perryman98 (1998)) also used a Plummer model and got 0.21 km s-1 for the central velocity dispersion. Again this value was derived by estimating the mass and the core radius. Compared to this work the values are 50% lower, but can be explained by the uncertainty in the estimation of the masses. Makarov et al. (makarov00 (2000)) used Tycho-2 proper motions to discuss the velocity dispersion of the Hyades, and found the velocity dispersion to be 0.32 km s-1 for the stars with the most precise proper motions. If known spectroscopic binaries were removed, the velocity dispersion decreased to 0.22 km s-1. The last value agrees well with the value obtained with Tycho-2 proper motions in Table 3.

## 4 Accuracy of astrometric radial velocities

From the cluster simulations and subsequent application of the maximum-likelihood method (Sect. 2) the astrometric radial velocities are estimated for the individual stars (or systems), $`\widehat{v}_{\mathrm{r}i}`$. Of course, the true radial velocities $`\overline{v}_{\mathrm{r}i}`$ are also known directly from the simulation. Thus the estimation errors $`\mathrm{\Delta }_{ij}=\widehat{v}_{\mathrm{r}i}\overline{v}_{\mathrm{r}i}`$ are known. Here, index $`j`$ is used to distinguish the different realisations of the cluster model. With $`_k`$ denoting an average over index $`k`$, the following statistics are computed:

$$\mathrm{\Delta }_j=\mathrm{\Delta }_{ij}_i$$
(4)
is the “cluster bias” in realisation $`j`$ (i.e., the common error for all stars in the cluster);

$$ϵ_{\mathrm{int}}=(\mathrm{\Delta }_{ij}\mathrm{\Delta }_j)^2_{ij}^{1/2}$$
(5)
is the “internal standard error” of the astrometric radial velocities (i.e., the dispersion of the individual values around the cluster bias); and

$$ϵ_{\mathrm{tot}}=\mathrm{\Delta }_{ij}^2_{ij}^{1/2}$$
(6)
is the “total standard error” of the astrometric radial velocities (i.e., including the cluster bias). Clearly $`ϵ_{\mathrm{int}}`$ is the relevant statistic for the precision of *relative* astrometric radial velocities within a given cluster, while $`ϵ_{\mathrm{tot}}`$ is relevant for the accuracy of *absolute* astrometric radial velocities. Both $`ϵ_{\mathrm{int}}`$ and $`ϵ_{\mathrm{tot}}`$ can be computed for various subsets depending on observable quantities such as the goodness-of-fit measure $`g_i`$, radial distance $`r`$, and mass or absolute magnitude. An interesting question is whether it is possible to observationally define subsets with reduced $`ϵ_{\mathrm{int}}`$ or $`ϵ_{\mathrm{tot}}`$.

The results presented below are based on solutions using the rejection limit $`g_{\mathrm{lim}}=15`$, although the results for $`g_{\mathrm{lim}}=10`$ are very similar. Any conclusions from these simulations are also applicable to the astrometric radial velocities published in Madsen et al. (madsen02 (2002)).

### 4.1 Standard errors versus goodness-of-fit

In Fig. 2 (top) the internal and total standard errors of the astrometric radial velocities are shown versus the goodness-of-fit $`g_i`$. The absence of any significant trend shows that $`g_i`$ is not a useful criterion for selecting “good” astrometric radial velocities. Even stars with $`g_i>10`$ are not worse than the rest in terms of radial-velocity precision. This somewhat counter-intuitive result can be understood if the line-of-sight component of the peculiar velocities is statistically independent of the tangential component. This is obviously the case for truly random motions, but one might expect that large proper-motion errors caused by photocentric motion in binaries should be correlated with large errors in the radial component.

### 4.2 Standard errors versus radius

The bottom part of Fig. 2 shows the internal and total standard errors of the astrometric radial velocities versus the distance $`r`$ from the cluster centre. In this case the standard errors clearly decrease from the centre out to 7–8 pc radius, after which they seem to increase again.

The initial decrease (for $`r<8`$ pc) is roughly in agreement with the Plummer model in Eq. (1) for $`M460M_{}`$ and $`r_c2.7`$ pc.

### 4.3 Standard errors versus mass and absolute magnitude

In Fig. 3, the internal standard errors of the astrometric radial velocities are plotted versus the true masses of the stars or systems (top) and versus the absolute magnitudes (bottom). The sample is divided at 3 pc (see Sect. 3.6.2). Inside 3 pc there is a clear difference in the velocity dispersion between the highest masses and 1 $`M_{}`$, although not as much as for a full equipartition of kinetic energy ($`\sigma _vM^{1/2}`$). The effect is much smaller outside of 3 pc. The velocity dispersion also seems to decline again for stars with masses less than 1 $`M_{}`$.

The effect can still be seen when the dispersion is plotted versus absolute magnitude instead of mass (Fig. 3, bottom), although the trend is less clear because of the many binary systems, for which there is no unique correspondence between system mass and total luminosity.

Together with the results of the previous section we can conclude that the practical minimum for the internal error of the astrometric radial velocities in the Hyades is around 0.20 km s-1, which is achieved for stars at an intermediate distance ($`2`$–3 core radii) from the cluster centre. At that distance there is little equipartition of kinetic energy, so it does not matter much if more or less massive stars are selected.

## 5 Non-modelled effects

The validity of the conclusions above depends critically on the realism of the $`N`$-body simulations. A number of non-modelled effects, and their possible impact on the results, are briefly considered below.

*Time-dependent tidal field:* When star clusters move through the galactic disk, they are subject to tidal shocks, and shock heating from the bulge. These effects are important to consider here since they increase the random motion of the stars. For globular clusters it has been found that tidal shocks accelerate significantly both core collapse and evaporation (Gnedin et al. gnedin99 (1999)).

In the case of open clusters, Bergond et al. (bergond01 (2001)) estimated that those with high-$`z`$ oscillations lose some 10–20% of the mass integrated over the lifetime of the cluster, mainly in low-mass stars, through disk-shocking. The Hyades have a low vertical velocity ($`W=6`$ km s-1 relative to the LSR), and therefore only oscillates with an amplitude of about 50 pc in $`z`$. Since this is small compared with the thickness of the disk, the disk-crossings should not cause much additional heating. The radial oscillations in the galactic plane, having an amplitude of 2 kpc, may be more important. The present $`N`$-body model assumes that the cluster moves in a circular galactic orbit. Thus it cannot be excluded that it underestimates the mass loss by perhaps some 5–10% of the initial mass. Preferentially, the lowest-mass stars leave the cluster, forming tidal tails (Combes et al. combes99 (1999)). Although this would slightly affect the estimation of the velocity dispersion, it would have only a very small effect on the number of observed stars of spectral type earlier than M0.

*Molecular clouds:* Terlevich (terlevich87 (1987)) studied open cluster $`N`$-body models with initially 1000 particles and moving in a circular orbit at 10 kpc from the galactic centre (i.e., assumptions comparable with this work). She concluded that the timescale for encounters with giant molecular clouds is of the same order of magnitude as the present age of the Hyades. Since such an encounter would probably be catastrophic, it can be assumed that the Hyades have not been exposed to such a meeting. More abundant are encounters with smaller interstellar clouds. They will not shorten the lifetime of open clusters significantly but may contribute to the tidal heating of the outer regions in a given cluster. Wielen (wielen75 (1975)) stated that gravitational shocks due to interstellar clouds will produce a significant flattening (up to 1:2) of the halo of the cluster perpendicular to the galactic plane. For the Hyades the flattening is 1:1.5 (Perryman et al. perryman98 (1998)). Since the galactic tidal field is also contributing to the flattening, it is doubtful if the Hyades have had any but minor interactions with interstellar clouds.

Perryman et al. (perryman98 (1998)) examined the possibility that the Hyades recently experienced an encounter with a massive object causing a tidal shear in the outer regions of the cluster, but excluded it based on the impulsive approximation (Spitzer spitzer58 (1958); Binney & Tremaine binney87 (1987)). Lindegren et al. (lindegren00 (2000)) included more velocity components in their model to test for non-isotropic dilation, and concluded that if such an effect existed it had to be higher than 0.01 km s-1 pc-1 to be detected with Hipparcos data. Effects from a tidal heating are thus not detectable in the Hyades with current astrometric precision.

*Brown dwarfs:* Despite extensive searches, no single-star brown dwarf (BD) candidate has been found in the Hyades (Reid & Hawley reid99 (1999); Gizis et al. gizis99 (1999); Dobbie et al. dobbie02 (2002)). Reid & Hawley (reid99 (1999)) found that the lowest-mass Hyades candidate star (LH 0418+13) has a mass of 0.083 $`M_{}`$, placing it very close to the hydrogen-burning limit. The only promising candidate brown dwarf in the Hyades is the unresolved companion in the short-period system RHy403 (Reid & Mahoney reid00 (2000)). Of course, the faintness of these substellar objects make them hard to observe, but still, the conclusion seems to be that the number today is quite small.

Adams et al. (adams02 (2002)) performed extensive simulations with a modified version of NBODY6 to model the brown dwarf population in open clusters, and concluded that the effects of different brown-dwarf populations were minimal, leaving the dynamics of the cluster largely unchanged.

The IMF in the version of NBODY6 used here cannot produce brown dwarfs, so this must be considered when defining the initial binary fraction. The IMF for brown dwarfs, or substellar masses, is very uncertain. Kroupa (kroupa01 (2001)) argues that a power-law value of $`\alpha =0.3\pm 0.7`$ is the most reasonable. Since stellar masses with $`M<0.08M_{}`$ are not produced in the code, one must represent the star–BD binary systems either as single stars or by overproducing binaries with secondary components slightly above the BD limit. Thus an initial binary fraction of 86% was assumed, which corresponds approximately to unity if brown dwarfs had been included. Based on the investigations of Adams et al. (adams02 (2002)), and considering that Hipparcos did not observe stars less massive than M0 stars in the Hyades, the above approximation should be sufficient for the present purpose.

*Cluster rotation:* Gunn et al. (gunn88 (1988)) did a comprehensive study of the rotation of the Hyades, but had to conclude that it was at most of the same size as their statistical error. Nonetheless they stated that their results suggested a cluster rotation, but not higher than 0.015 km s-1 pc-1.

Perryman et al. (perryman98 (1998)) did a thorough study of the velocity residuals and concluded that they were consistent with a non-rotating system and the given observational errors. Lindegren et al. (lindegren00 (2000)) tested the Hyades for rotation by assuming solid-body rotation parameters, but found that it was too small to be detected, setting an upper limit of 0.01-0.02 km s-1 pc-1. If this upper limit should equal the true rotation of the Hyades, then the effect is non-negligible at 10 pc compared to the internal error. But there seems to be nothing in the present study suggesting such a rotation.

But probably the solid-body assumption is too simple. In the globular cluster $`\omega `$ Centauri, Merritt et al. (merritt97 (1997)) found that only at small radii could the rotation be approximated by a solid-body. Beyond that the rotation falls off. Einsel & Spurzem (einsel99 (1999)) did theoretical investigations on the influence of rotation on the dynamical evolution of collisional stellar systems, that could explain the findings by Merritt et al. (merritt97 (1997)). In fact, it seems that only inside the half-mass radius is it reasonable to talk about a solid-body rotation (cf. Kim et al. kim02 (2002)).

Although it is unlikely that the cloud in which the Hyades formed had zero angular momentum, there currently exists no certain measure of the rotation. In the model, it is instead assumed that the effects are sufficiently small and can be ignored.

*Expansion:* During the evolution of a cluster parts of it expand and parts of it contract. Under the assumption that the relative expansion rate equals the inverse age of the cluster, Dravins et al. (1999b ) estimated that an isotropic expansion of the Hyades would lead to a bias in the astrometric radial velocity of 0.07 km s-1 of the centroid velocity. This is completely negligible and any expansion effects have been ignored.

To summarise, it appears that none of these non-modelled effects would affect the results very significantly. While the modelling of tidal fields and brown dwarfs could be improved, the possible effect of cloud encounters remains an uncertainty which cannot easily be included in the modelling of a specific cluster such as the Hyades. Although NBODY6 allows encounters with interstellar clouds, the option has not been used in the present study.

## 6 Conclusions

A dynamical model of the Hyades cluster, based on $`N`$-body simulations using the NBODY6 code, has been fitted to the astrometric information available in the Hipparcos and Tycho-2 catalogues in order to study the accuracy of astrometric radial velocities. The number of stars as function of magnitude, their three-dimensional distribution, and the distribution of proper motions have been adequately reproduced by the model, as well as basic binary statistics. No spectroscopic radial velocities have been used in the present study (except for the initial membership determination by Perryman et al. perryman98 (1998)) meaning that the results should be directly comparable with the astrometrically determined radial velocities of Hyades stars by Lindegren et al. (lindegren00 (2000)) and Madsen et al. (madsen02 (2002)).

From the simulations it is concluded that the velocity dispersion of the Hyades decreases from $`\sigma _v0.35`$ km s-1 at the centre of the cluster to nearly $`0.2`$ km s-1 at 7–8 pc from the centre. Outside the tidal radius of 10–11 pc, the dispersion slightly increases again. Compared with previous studies of the velocity dispersion in the centre of the Hyades, the results here indicate a somewhat higher value.

The internal velocity dispersion contributes to the random errors of the astrometric radial velocities with the same magnitude. This is significantly less than the $`\sigma _v=0.49`$ km s-1 estimated in Madsen et al. (madsen02 (2002)) directly from the Hipparcos observations. This discrepancy can be understood with reference to Table 3 as an *overestimation* from the observed data when the less strict rejection limit $`g_{\mathrm{lim}}=15`$ was used. Thus the previous estimate of the internal standard error (due to the dispersion) can now be almost halved.

In fact, stars with an expected velocity dispersion as low as 0.20 km s-1 can be selected for studies that compare astrometric and spectroscopic radial velocities in order to disclose astrophysical phenomena causing spectroscopic line shifts. However, it should be remembered that the total standard error, including the uncertainty of the motion of the cluster centroid, is still of order 0.55–0.65 km s-1 (Fig. 2, bottom), in agreement with the previous estimate.

Attempts to see a radial dependence of the velocity dispersion with Hipparcos and Tycho-2 astrometry have been inconclusive. The observed relation is essentially flat for the most optimal sample. Given the uncertainty of the estimated velocity dispersions when the stars are divided into radial shells, this result is not surprising. Similar examples can be found in the simulations. Only when the mean relation is computed from several realisations of the cluster model do the variations become clear. In particular, it appears that the structure of dispersion/radius relation reported by Madsen et al. (madsen01 (2001)) does not reflect typical dynamical properties of the cluster, but could result by chance or from some (unknown) mechanism related to the photocentric motions of undetected binaries.

The fit has yielded an estimate of the initial cluster mass of 1100–1200 $`M_{}`$ and of the initial multiplicity, which appears to be very high (possibly near 100%, if brown-dwarf companions are included). The current cluster mass is estimated to be $``$ 460 $`M_{}`$ with a tidal radius of $``$ 11 pc and a mean velocity dispersion within $`r<3`$ pc of 0.32 km s-1.

Some of the differences between observations and simulations could be due to some of the non-modelled features discussed in Sect. 5, which would lead to a higher initial particle number in the model and which might also solve some of the discrepancies noted in the binary statistics. The development of numerical tools such as NBODY6 to include e.g. a time-dependent tidal field would allow an improved realism of the Hyades model, and to study the effect on the accuracy of astrometric radial velocities from assumed negligible contributions to the velocity field with respect to the Hipparcos precisions.

The method used to estimate astrometric radial velocities discussed in Sect. 2 cannot eliminate of the error contribution from the internal dynamics of the cluster, no matter how precise the astrometry might be. The velocity dispersion therefore sets a fundamental limit on the accuracy of astrometric radial velocities, and as a consequence the results from the simulations presented here also apply to planned astrometric space missions such as GAIA (Perryman et al. perryman01 (2001)), even though it has been Hipparcos observations of the Hyades that have been simulated.

The Hipparcos and Tycho-2 catalogues contain the best available astrometry to study the internal velocity structure of the nearest open cluster, the Hyades. To study it in greater detail, even better astrometry is needed. The GAIA mission, in combination with improved $`N`$-body simulations, will make it possible to observe directly the internal velocity field of the Hyades, and give us insight in the kinematics of the Hyades in particular and open clusters in general.

###### Acknowledgements.

I thank Sverre Aarseth for making NBODY6 freely available, Tim Adams for helping me with the code, and Lennart Lindegren, Melvyn B. Davies, and Dainis Dravins for useful comments and valuable suggestions.


format: text

text:
# A search for X-ray variability in the open cluster NGC 2516

## 1 Introduction

NGC 2516 is a rich, well studied, galactic open cluster. It has a low metallicity (eg Cameron 1985, Jeffries, James & Thurston 1998), is nearby ($``$387pc; Jeffries, Thurston & Pye 1997), is young ($`10^8`$ years; eg Feinstein, Marraco, Mirabel et al 1973) and has a relatively low line of sight extinction. Because of its low metallicity, it is an excellent target to test how coronal activity and stellar structure depend on metal content.

NGC 2516 has been the subject of several studies in the X-ray band. It was observed using ROSAT by Dachs & Hummel (1996), Jeffries, Thurston & Pye (1997) and Micela et al (2000). More recently it was the subject of a study using XMM-Newton (Sciortino et al 2001) and also Chandra (Harnden et al 2001). These studies have concentrated on determining the cluster X-ray luminosity function: Sciortino et al (2001) found that dwarf G and K stars in NGC 2516 were less luminous than those in the Pleiades (which has a similar age but solar metallicity ratio). In contrast, the luminosity function of dwarf M stars in NGC 2516 and the Pleiades were indistinguishable.

In this paper, we show the results of a search for X-ray variable objects in NGC 2516 using XMM-Newton data. With the relatively large number of X-ray objects it is an excellent target to determine the frequency of stellar flaring activity in a young, metal poor cluster. Further, NGC 2516 has been observed on 8 occasions with XMM-Newton since it has been used as a calibration target.

## 2 Observations

The satellite XMM-Newton was launched in Dec 1999 by the European Space Agency. It has the largest effective area of any imaging X-ray satellite. The EPIC instruments have imaging detectors covering the energy range 0.15–10keV with moderate spectral resolution. The observation log is shown in Table 1. The first two observations were centered at 7h 58m 20s, –60 52$`^{^{}}`$ 13$`^{^{\prime \prime }}`$ (2000), and the rest at 7h 58m 22s, -60 45$`^{^{}}`$ 36$`^{^{\prime \prime }}`$ (2000). (We do not discuss the data taken using the Reflection Grating Spectrometer since the sources were comparatively faint and source confusion would be a problem. The Optical Monitor was in blocked mode).

The data were processed using the XMM-Newton Science Analysis Software (SAS) v5.3.3. An image was generated from each dataset in the energy range 0.2-5.0keV. We removed data at the edge of the detector furthest from on-axis since these rows had increased noise. We also selected those events which were tagged as FLAG=0: this excludes events close to, for instance, hot pixels and chip edges. Time intervals which showed high levels of particle background were excluded. The resulting images were then source searched using standard SAS tasks. Events were then extracted around a specified radius from the source positions. The events for each source were then subject to a Discrete Fourier Transform in the time interval of 1min to 0.7$`\times `$Exposure length. All power spectra which showed a peak above a certain criteria above the mean power were flagged and inspected by eye. We also had a second flagging procedure which indicated which objects had a count rate which varied over a certain value.

After the completion of this process, 4 sources were found to show significant variability in their X-ray light curves. Their X-ray positions were compared to the list of X-ray sources in NGC 2516 and their corresponding optical counterparts (Jeffries, Thurston & Pye 1997). We found counterparts for all 4 sources: the offset between the XMM-Newton and optical counterpart was typically better than 6$`^{^{\prime \prime }}`$ (which is the FWHM of sources in the EPIC pn). We show their X-ray positions along with their counterpart of the X-ray catalogue of Jeffries et al (1997) in Table 2. For each source we examined the other datasets to determine if the source was also detected at other epochs. In each case the source was detected but at low count rates (typically 0.01-0.02 ct/s). However, they showed no significant variation in these other datasets.

We show the (background subtracted) X-ray light curves for each source in Figure 1. Three of the sources show evidence for a rapid flare, although in one of those (JTP80) the observation ends shortly after the start of the flare. The other source (JTP122) shows a more gradual rise in brightness. The most dramatic source is JTP138 which shows a gradual rise in flux lasting $``$12 min followed by a sharp increase in flux taking place in less than 2 min. After reaching maximum, it reaches half its maximum intensity after $``$20 mins. The duration of the flare lasts around 90 min. We show in Figure 2 the light curves in the 1–3keV and 3–10keV energy bands for JTP 138 in time bins of 500 sec together with its softness ratio. Although the hardness ratio curve shows some evidence for a softening over the duration of the flare, it has large errors. However, there is some evidence for a small increase in flux in the 1–3keV light curve just before the sharp rise that is not seen in the 3–10keV light curve, although again this may just be due to the low signal to noise.

## 3 Source Spectra

We extracted spectra from each of the 4 variable sources: only the time interval of the flare duration and single and double events were used. We also extracted background spectra from a nearby source free area off the detector. We used the ‘canned’ response files appropriate to the location of the source on the detector and generated an auxiliary file using the SAS task arfgen. We fitted the spectra using the X-ray fitting package XSPEC (Dorman & Arnaud 2001).

All spectra could be fitted using an absorbed one or two temperature plasma model (we used the MEKAL model in XSPEC). We show the spectrum of the brightest flare (JTP138) in Figure 3: the fit is good with no evidence for significant residuals at particular energies. This source also showed evidence for an Fe line at 7.0$`\pm 0.2`$keV (the other sources did not). We show the spectral fits in Table 3.

All spectra show a moderate amount of absorption, typically $`5\times 10^{20}`$ cm-2, and temperatures of a few keV. In determining their luminosity we assume that all 4 sources are members of NGC 2516 and its distance is 387 pc (Jeffries, Thurston & Pye 1997). In the case of JTP138, where we observe the total duration of the flare, the energy released was $`8\times 10^{34}`$ ergs. In the case of JTP80, where we observed only the start of the flare, the energy released was greater than $`1\times 10^{34}`$ ergs. There was a similar energy release observed in JTP114 where most of the flare was observed.

## 4 The nature of the sources

We have identified each of the 4 X-ray variable sources with an optical counterpart (Table 2). Jeffries et al (1997) identify JTP122 as having an AOV spectral type: the other 3 sources have not been classified spectroscopically. Further, JTP138 and JTP122 were both classed as cluster members by King (1978) from a proper motion survey. To determine the likelihood of the other two sources being cluster members we plotted the $`V,BV`$, data of Sung et al (2002) and traced the position of the main sequence. We then over-plotted the optical colours of our X-ray variable sources (Figure 4). All 4 sources lie close to the main sequence of NGC 2516: we conclude that each source is likely to be a cluster member.

Stellar coronal activity is common in stars which have a convective layer and a relatively rapid rotation period. Enhanced levels of activity are often seen in close binary systems, for example the RS CVn systems which typically consist of a giant or sub-giant primary with a main-sequence or sub-giant companion. In many RS CVn systems the rotation of one of the components is synchronised with the orbital motion, leading to enhanced rotation and high levels of activity, seen as large spot coverage and flares in X-ray, optical and radio bands. Additionally X-ray activity has also been seen on early-type stars, such as Castor A and B. Here it is thought the flaring component may originate from low-mass, late-type companions (Güdel at al. 2001). X-ray flares have also been reported in pre-main sequence stars (eg Schmidt 1994) and single solar type stars (eg Landini et al 1986).

The energy released in the X-ray flares we observed in NGC 2516 were $`\times 10^{3435}`$ ergs and lasted around 1–2 hrs. This contrasts to the energy released on a flare seen on Algol (7$`\times 10^{36}`$ ergs) which lasted over a day (Ottmann & Schmitt 1996). In dwarf M stars flares have been seen in the $`U`$ band (eg Byrne, Doyle & Butler 1984) although they are relatively rare and last several 10’s of minutes. In X-rays, dwarf M stars show flares of energies typically $`10^{3034}`$ ergs with a duration less than 1 hour (eg Cheng & Pallavicini 1991). We consider that it is unlikely that the flares originate from pre-main sequence stars since NGC 2516 is not a star forming region. Solar type stars have also been found to show flares in various energy bands with energies $`10^{3338}`$ ergs, although they are very rare (eg Schaefer et al 2000).

In contrast, flares have been seen on RS CVn systems which share similar characteristics to the flares we have observed. van den Oord, Mewe & Brinkman (1988) estimated the energy released in a flare on the RS CVn system $`\sigma ^2`$ CrB was 2.4$`\times 10^{34}`$ ergs and lasted $``$ 90 mins. Further, they found clear evidence for a softening of the light curve after the sharp rise in flux. This is comparable to the RS CVn system II Peg which released $`>5\times 10^{34}`$ ergs in a flare lasting more than several hours (Doyle et al 1991, Doyle, van der Oord & Kellett 1992). Their X-ray spectra are well modelled using a thermal plasma of a few keV. Further, there are several instances where Fe K$`\alpha `$ emission has been seen in stars with active stellar coronae (eg Tsuru et al 1989, Güdel et al 2001).

We suggest that JTP 138, JTP 80 and JTP 114 are good candidates for RS CVn systems. The fact that JTP 122 has been classified as AO V spectral type may indicate that it has a dwarf M star as a companion. Followup spectroscopy will be able to test these assertions.

## 5 The frequency of flaring activity

NGC 2516 has been observed using XMM-Newton for 97.4 ksec (cf Table 1). We detected 4 significantly variable X-ray sources over this time implying the detection of one X-ray variable source every $``$25ksec. The first two observations were offset from the other observations by $`7^{^{}}`$ (cf Table 1). In this first field we detected 113 X-ray sources and in the second 151 (using the longest exposures in these fields). We detected 1 variable source in field 1 and 3 in field 2. This gives 2.7$`\times 10^4`$ variable X-ray sources per ksec per X-ray source in field 1 and 3.0$`\times 10^4`$ in field 2. Although simplistic, if every source in the field was a potential X-ray variable source, then these rates would imply that each source would show a flare every $``$40 days.

## 6 Comparison with solar flares

Flaring activity has been seen on the Sun at many wavelengths and timescales. We have searched for solar flares that have a similar temporal profile to ‘JTP 138’ using the Geostationary Operational Environmental Satellites (GOES). The satellites are in geostationary orbit, their purpose being to monitor the Sun’s X-ray emission continuously.

We have identified several solar flares which have similar timescales to that seen in JTP 138. Figure 5 (curve a) shows the light curve of a ‘small’ solar flare which has an intensity several orders of magnitude lower than the largest flares in the current solar cycle. Like JTP 138 it shows a small increase in flux (lasting $``$10 mins) before a rapid rise which takes around 10 mins to reach maximum flux. We have estimated its temperature using ratios of the low energy channel on GOES (1-8Å, 1.5–12.4keV) with the high energy channel (0.5-4Å, 3–25keV). The temperature reached a maximum of 0.6 keV, with a luminosity of 7$`\times `$1025 ergs s-1. The second example solar flare is also shown in Figure 5 (curve b). This flare is brighter, and reaches a temperature of 0.8 keV, with a luminosity of 2$`\times `$1026 ergs s-1. In contrast to the small solar flare and JTP 138, there is no evidence for a pre-flare increase in flux. Both solar flares last around 1–1.5 hrs, similar to JTP 138. In Figure 6 we show the energy resolved light curve of the small solar flare. Similar to the flare seen in JTP 138, the solar flare reaches a maximum in hard X-rays earlier than in soft X-rays.

Although the shape and duration of these solar flares are similar to that of JTP 138, they have lower temperatures and much lower luminosities. The luminosity may simply be a reflection of the volume of the emitting plasma. In the case of RS CVn systems, the flare size can extend to the size of the binary system (or indeed, originate in the intra-binary region as seen in the RS CVn system CF Tuc, Gunn et al 1997). For reasonable system parameters and an orbital period of 2 days, the separation is $`5\times 10^{11}`$ cm. Moreover, the Suns rotation rate has been slowing down since it joined the main sequence. It therefore has lower magnetic activity compared to stars which have recently joined the main sequence (as in the case of stars in NGC 2516).

Numerous explanations for the triggering mechanism of solar flares have been developed over the years. The most popular current theory for converting magnetic energy into kinetic and thermal energies is that of magnetic reconnection (eg Priest & Forbes 1999). The site of the reconnection takes place high in the solar corona above magnetic loops, or between at least two loop systems that interact. At the reconnection site, the electrons are accelerated and move along the loops with a high velocity ($`v0.5c`$). Some of these electrons will bombard the chromosphere which is denser and cause this plasma to be heated. This plasma will then evaporate into the loops causing the soft X-ray emission such as is shown in Figure 5 (cf Harra 2000 for more details). Although the details are still uncertain, observations support this view.

Feldman, Laming & Doschek (1995) found a correlation between solar and stellar flare temperatures and emission measures. Shibata & Yokoyama (1999) proposed a model based on a magnetic reconnection model with heat conduction and chromospheric evaporation to account for this correlation. Based on our spectral fits to JTP 138, we find an emission measure of 7$`\times 10^{53}`$ cm-3 for JTP 138. For a temperature of 6$`\times 10^7`$ K for the hot component, Figure 1 of Shibata & Yokoyama (1999) suggests a magnetic field strength of $``$100 G for this source. Measuring stellar magnetic field strengths is difficult. However, observations of various RS CVn systems show evidence for magnetic field strengths of the order $``$100–1000 G (eg Donati 1999), indicating that a field strength of 100 G for JTP 138 is reasonable.

## 7 Conclusions

We have performed a X-ray time variability survey of the open cluster NGC 2516. We have found 4 systems which have shown significant variability. Three of these systems show light curves which are characteristic of stellar flaring activity. Based on the flare duration and their spectral properties we suggest that they are RS CVn systems. The fourth system has been identified with a AO V type star and propose that it may have an unseen dwarf M star which shows flaring activity. An optical spectroscopic study of these objects is required. Finally, this study shows that open clusters are excellent objects to search for stellar flaring activity and may be relatively common.

## 8 Acknowledgments

This paper is Based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). We thank Mat Page for the use of some of his software. We acknowledge the use of the WEBDA database of galactic clusters.


text:
# Condition for unambiguous state discrimination using local operations and classical communication

## Acknowledgements

The author is grateful to the University of Hertfordshire for the award of a Postdoctoral Research Fellowship.


text:
# A method to include the spin-fluctuation in the ab-initio electronic-structure calculation.

## Abstract

We present a new method of ab-initio electronic-structure calculation including the spin-fluctuation (SF) self-consistently. We start from the Luttinger-Ward functional given as the sum of the LDA functional plus the temperature-dependent part of the SF energy functional. The size of interactions used in it are determined in a similar manner to the self-consistent renormalization theory by Moriya and Kawabata. Obtained paramagnetic susceptibilities on Pd, Ni, Fe, and fcc-Co above $`T_c`$ show rather good agreements with experiments.

In magnetic metals like Fe, LDA gives quite good description for ground state properties at zero-temperature $`T=0`$; LDA reproduces the so-called Slater-Pauling curve obtained by experiments very well , also does the spin waves . But neither LDA nor its crude extension to finite temperature contains the contributions due to the spin-fluctuation (SF). SF for magnetic materials should have strong temperature dependence, and enforce them paramagnetic above their critical temperatures $`T_\mathrm{c}`$. In Refs.  and , methods including SF had developed with adiabatic approximation in an ab-initio scheme; they first calculated the adiabatic surface of all the configuration of static SF, then calculated the partition function assuming spin variables as classical. It works rather well, though these approximations could be problematic in paramagnetic phases at high $`T`$ where the Stonar-like excitations and spin-wave like excitation could mix well. In addition, the ground state, where they calculated the adiabatic surface, can be changed at higher $`T`$. Especially in Pd, which is ”almost-ferromagnetic”, these methods have not been applied to, presumably because the change of ground state would cause problems. As the opposite limit where the Heisenberg-like model works well, we know ”metallic weak-ferro-magnet”, which are well described by the self-consistent renormalization (SCR) theory  given by Moriya and Kawabata. They evaluate SF in a kind of random-phase-approximation (RPA), where we should determine the size of interaction in a self-consistent (SC) manner. Here we present a new ab-initio SC method including SF, where the size of interaction is determined self-consistently as in the SCR theory. As results compared with experiments, we show paramagnetic susceptibilities for Pd, Ni, Fe, and fcc-Co above $`T_c`$. Our calculations are based on a finite-temperature version of AkaiKKR , in the atomic sphere approximation (ASA).

In the treatment like LDA+$`U`$ and their extensions like LDA+’dynamical mean field theory (DMFT)’, a key problem is in the determination of the effective interaction $`U`$. Lichtenstein, Katsnelson, and Kotliar presented such a scheme to describe magnetic susceptibilities of Ni and Fe at Hight $`T`$. Their method can take into account higher-level of correlation effects than our treatment, though their $`U`$ is given externally . Very recently, Biermann, Aryasetiawan, and Georges suggested $`GW`$+DMFT to determine them self-consistently . Our method can be identified as a simplified version of such a method.

We start from the Luttinger-Ward (LW) functional $`\mathrm{\Omega }[g]`$ at finite temperature. It can be written as

$`\mathrm{\Omega }[g]=\mathrm{Tr}\left[\mathrm{log}(g)+{\displaystyle \frac{g}{\tau }}\right]+\mathrm{Tr}\left[{\displaystyle \frac{^2}{2m}}g\right]+E_{\mathrm{ext}}[n]`$ (1)
$`+E_{\mathrm{Cou}}[n]+E_{\mathrm{xc}}[g]`$ (2)
Here $`g_\sigma (1,2)`$ is the finite-temperature Green function, and $`1𝐫_1\tau _1`$. $`n_\sigma (𝐫)`$ denotes the electron density and is treated as the functional of $`g_\sigma (1,2)`$ in Eq.(2). We may identify the first term in right-had side (RHS) in Eq.(2) as ”the population-entropy” term because it is reduced to be $`T`$ times the entropy of the system in the non-interacting case. The second term is the kinetic energy term, and $`E_{\mathrm{ext}}[n]`$ is the external-potential term containing the chemical potential. Terms $`E_{\mathrm{Cou}}`$ and $`E_{\mathrm{xc}}`$ are related to the Coulomb interaction between electrons. Sum of them are referred to as $`\mathrm{\Phi }[g]`$ in literatures . An advantage starting from the LW functional is in the conservation laws; they apparently hold if our approximated $`\mathrm{\Omega }[g]`$ keeps the corresponding symmetries; such approximations are called as the Baym-Kadanoff approximation . Our approximation for $`E_{\mathrm{xc}}[g]`$ is

$$E_{\mathrm{xc}}[g]E_{\mathrm{xc}}^{\mathrm{LDA}}[n]+\overline{\mathrm{\Omega }^{\mathrm{SF}}}[g],$$
(3)
where the SF functional $`\overline{\mathrm{\Omega }^{\mathrm{SF}}}`$ contains just the temperature-dependent part of the SF energy written as $`\overline{\mathrm{\Omega }^{\mathrm{SF}}}=\mathrm{\Omega }^{\mathrm{SF}}\mathrm{\Omega }^{\mathrm{SF0}}`$, where $`\mathrm{\Omega }^{\mathrm{SF0}}`$ is $`\mathrm{\Omega }^{\mathrm{SF}}`$ at $`T=0`$ so as to remove the temperature-independent part. We use $`E_{\mathrm{xc}}^{\mathrm{LDA}}[n]`$ given by . This $`\overline{\mathrm{\Omega }^{\mathrm{SF}}}`$ should have strong temperature-dependence. The stationary-point condition $`\delta \mathrm{\Omega }[g]/\delta g=0`$ gives $`g_\sigma (1,2)`$; it is written as

$`V_{\mathrm{LDA}}^\sigma (𝐫)=V_{\mathrm{ext}}^\sigma (𝐫)+V_{\mathrm{Cou}}^\sigma (𝐫)+V_{\mathrm{LDA}\mathrm{\_}\mathrm{xc}}^\sigma (𝐫),`$ (4)
$`V_{\mathrm{eff}}^\sigma (𝐫_1𝐫_2,\omega _p)=V_{\mathrm{LDA}}^\sigma (𝐫_1)\delta (𝐫_1𝐫_2)+{\displaystyle \frac{\delta \overline{\mathrm{\Omega }^{\mathrm{SF}}}}{\delta g^\sigma (𝐫_1𝐫_2,\omega _p)}},`$ (5)
$`\left[i\omega _p\left({\displaystyle \frac{_{𝐫_1}^2}{2m}}+V_{\mathrm{eff}}\right)\right]g^\sigma (𝐫_1𝐫_2,\omega _p)=\delta (𝐫_1𝐫_2),`$ (6)
where $`\sigma =\pm 1`$ for $``$ and $``$. $`\delta \overline{\mathrm{\Omega }^{\mathrm{SF}}}/\delta g`$ is the contribution to the self-energy from SF. $`\omega _p`$ denotes the fermion Matsubara frequencies.

It is difficult to include the full non-locality and $`\omega `$-dependence of $`\delta \overline{\mathrm{\Omega }^{\mathrm{SF}}}/\delta g`$. So we restrict the variational space as follows and try to find the optimum solution within it. $`V_{\mathrm{LDA}}`$ gives exact solutions of $`\mathrm{\Omega }[g]`$ if we omit $`\overline{\mathrm{\Omega }^{\mathrm{SF}}}`$. Therefore we take a variational space of $`g_\sigma (1,2)`$ so that it is generated by $`V_{\mathrm{LDA}}`$ plus the onsite magnetic fields $`B^R`$ just acting on only $`d`$ electrons in each atomic-sphere (AS) $`R`$. Thus $`g_\sigma (1,2)`$ within the variational space is generated by solving a set of SC Eqs. for given $`B^R`$. The set consists of Eq.(6) and

$`V_{\mathrm{eff}}^\sigma (𝐫_1𝐫_2)=V_{\mathrm{LDA}}^\sigma (𝐫_1)\delta (𝐫_1𝐫_2)+\sigma \mu _\mathrm{B}{\displaystyle \underset{R}{}}B_R\widehat{P}_R.`$ (7)
Here $`\widehat{P}_R=\delta (r_1r_2)_mY_{2m}(\theta _1,\varphi _1)Y_{2m}^{}(\theta _2,\varphi _2)`$ is the projection operator to the onsite $`d`$ channel. Note that the set $`\{B_R\}`$ for all $`R`$ is only the variational parameters, which determine $`g_\sigma (1,2)`$. Therefore, our problem is reduced to searching stationary point of $`\mathrm{\Omega }[g]`$ as a functional $`\{B_R\}`$. (As in our static approximation $`E_{\mathrm{xc}}[g]E_{\mathrm{xc}}^{\mathrm{LDA}}[n]`$, the stationary point means minimum.) This means that we have to solve $`0=\delta \mathrm{\Omega }/\delta g\times \delta g/\delta B_R=\left(V_{\mathrm{eff}}+V_{\mathrm{LDA}}+\delta \overline{\mathrm{\Omega }^{\mathrm{SF}}}/\delta g\right)\times \delta g/\delta B_R`$. This determines $`B_R`$ as the solution of an inversion equation

$`{\displaystyle \frac{\delta \overline{\mathrm{\Omega }^{\mathrm{SF}}}}{\delta g}}{\displaystyle \frac{\delta g}{\delta B_R}}={\displaystyle \underset{R^{}}{}}{\displaystyle \frac{M_R^{}}{B_R}}B_R^{}^{\mathrm{SF}},`$ (8)
where $`M_R=\mu _\mathrm{B}_Rd^3r(n_{}^dn_{}^d)`$. $`B_R=B_R^{\mathrm{SF}}+B_R^{\mathrm{ext}}`$ if we add external magnetic field $`B_R^{\mathrm{ext}}`$ as a probe. $`M_R`$ is the spin magnetic moment of $`d`$ electrons in $`R`$. This kind of inversion equation like Eq.(8) is also used in the optimized-effective-potential method . Equations (6),(7), and (8) give a set of SC equations. With the restriction of the variational space, all the quantities in RHS in Eq.(2) are treated as the functional of $`B_R`$ in Eq.(7), or corresponding magnetic field $`M_R`$.

As for $`\mathrm{\Omega }^{\mathrm{SF}}[g]`$, we give it in a simple model in the onsite approximation; it consists of contributions from each sites as $`\mathrm{\Omega }^{\mathrm{SF}}[g]=_R\mathrm{\Omega }_R^{\mathrm{SF}}[g_R]`$. Here $`g_R`$ denotes the $`d`$-channel onsite part of $`g_\sigma (1,2)`$ as

$`g_R^\sigma (𝐫,𝐫^{},\omega _p)={\displaystyle \underset{m}{}}g_{mR}^\sigma (\omega _p)\mathrm{\Phi }(r,i\omega _p)`$ (9)
$`\times Y_{2m}(\theta ,\varphi )Y_{2m}^{}(\theta ^{},\varphi ^{})\mathrm{\Phi }(r^{},i\omega _p).`$ (10)
Here we now take just the diagonal part with respect to the magnetic quantum number $`m`$. $`\mathrm{\Phi }(r,i\omega _p)`$ denotes solutions of the radial Schrödinger equation. $`\mathrm{\Omega }_R^{\mathrm{SF}}[g_R]`$ is given as

$`\mathrm{\Omega }_R^{\mathrm{SF}}[g_R]={\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\nu _q}{}}Q_0^{\mathrm{RPA}}(I_R^\mathrm{L}\chi _R^{0\mathrm{L}}(\nu _q))`$ (11)
$`+{\displaystyle \frac{2}{\beta }}{\displaystyle \underset{\nu _q}{}}Q_0^{\mathrm{RPA}}(I_R^\mathrm{T}\chi _R^{0\mathrm{T}}(\nu _q)),`$ (13)
$`Q_0^{\mathrm{RPA}}(s)\mathrm{log}(1s)+s`$
$`\nu _q`$ denotes the bosonic Matsubara frequencies. The bare longitudinal and transversal onsite spin-polarization functions, $`\chi _R^{0\mathrm{L}}(\nu _q)=\widehat{S}_{zR}\widehat{S}_{zR}_0`$ and $`\chi _R^{0\mathrm{T}}(\nu _q)=\frac{1}{2}\widehat{S}_{xR}\widehat{S}_{xR}+\widehat{S}_{yR}\widehat{S}_{yR}_0`$, are given as

$`\chi _R^{0\mathrm{L}}(\nu _q)={\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{m}{}}{\displaystyle \underset{\omega _p}{}}\{g_{mR}^{}(\omega _p)g_{mR}^{}(\nu _q+\omega _p)`$ (14)
$`+g_{mR}^{}(\omega _p)g_{mR}^{}(\nu _q+\omega _p)\}W_{mR}(i\omega _p,i\nu _q),`$ (15)
$`\chi _R^{0\mathrm{T}}(\nu _q)={\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{m}{}}{\displaystyle \underset{\omega _p}{}}\{g_{mR}^{}(\omega _p)g_{mR}^{}(\nu _n+\omega _p)`$ (16)
$`+g_{mR}^{}(\omega _p)g_{mR}^{}(\nu _n+\omega _p)\}W_{mR}(i\omega _p,i\nu _q),`$ (17)
$`W_{mR}(i\omega _p,i\nu _q)=\left\{{\displaystyle _0^{R_{\mathrm{max}}}}r^2𝑑r\mathrm{\Phi }_{mR}(r,i\omega _p)\mathrm{\Phi }_{mR}(r,i\nu _q+i\omega _p)\right\}^2,`$ (18)
where $`W_{mR}`$ is the square of the radial part of integrals. $`R_{\mathrm{max}}`$ denotes the radius of AS $`R`$. $`\mathrm{\Omega }_R^{\mathrm{SF}}[g_R]`$ is the SF energy given by the onsite-only RPA, though we treat the effective interactions $`I_R^\mathrm{L}`$ and $`I_R^\mathrm{T}`$ as the functionals of $`B_R`$. They are determined in a SC manner explained later. Together with $`\mathrm{\Omega }_R^{\mathrm{SF0}}`$ calculated in the same procedure with taking the $`T0`$ limit in Eq.(13), we can obtain $`\overline{\mathrm{\Omega }^{\mathrm{SF}}}[g]=_R\overline{\mathrm{\Omega }_R^{\mathrm{SF}}}`$, where $`\overline{\mathrm{\Omega }_R^{\mathrm{SF}}}\mathrm{\Omega }_R^{\mathrm{SF}}\mathrm{\Omega }_R^{\mathrm{SF0}}`$. Corresponding to our onsite approximation on $`\overline{\mathrm{\Omega }^{\mathrm{SF}}}[g]`$, we only take diagonal terms on $`R`$ in the inversion Eq.(8). That is, we just take $`\delta g_R/\delta B_R`$, and $`M_R/B_R`$, neglecting the off site contributions. Then Eq.(8) is simplified as

$`B_R^{\mathrm{SF}}=\left({\displaystyle \frac{\delta \overline{\mathrm{\Omega }_R^{\mathrm{SF}}}}{\delta g_R}}{\displaystyle \frac{\delta g_R}{\delta B_R}}\right)\left({\displaystyle \frac{M_R}{B_R}}\right)^1.`$ (19)
In practice, $`\chi _R^0(\nu _q)`$, which denotes each term contained in $`\chi _R^{0\mathrm{L}}`$ and $`\chi _R^{0\mathrm{T}}`$ of Eqs.(15,17), is calculated through these formulas symbolically written as

$`\chi _R^0(\nu _q)={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\mathrm{\Gamma }^0(\nu )d\nu }{i\nu _q\nu }},`$ (20)
$`\mathrm{\Gamma }^0(\nu )={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega \left\{f(\omega )f(\nu +\omega )\right\}X(\nu +\omega )Y(\omega )W_{mR}(\omega ,\nu ).`$ (21)
Here $`f(\omega )`$ in this convolution Eq.(21) denotes the Fermi distribution function; $`X(\omega ),Y(\omega )`$ denotes the partial density of states (DOS) $`D_{mR}^\sigma (\omega )=\mathrm{Im}g_{mR}^\sigma (\omega )/\pi `$. The factor $`W_{mR}(\omega ,\nu )`$ introduces a natural cutoff for the convolution of DOS in Eq.(21).

It is necessary to include contributions due to $`\widehat{S}_{zR}^m\widehat{S}_{zR}^m_0`$ which carries SF with other magnetic angular momentum $`m0`$ in addition to the contribution of $`m=0`$ ($`\widehat{S}_{zR}^{m=0}`$ is $`\widehat{S}_{zR}`$ given above). If we assume spherical symmetry, we can evaluate it as $`\widehat{S}_{zR}^m\widehat{S}_{zR}^m_0\frac{5|m|}{5}\widehat{S}_{zR}\widehat{S}_{zR}_0`$ by counting the number of allowed combinations of $`m_1`$ and $`m_2`$ for $`m=m_1+m_2`$; Then we can include the contributions by using

$`Q^{\mathrm{RPA}}(s){\displaystyle \underset{m=4,\mathrm{}4}{}}Q_0^{\mathrm{RPA}}\left({\displaystyle \frac{5|m|}{5}}s\right)`$ (22)
instead of $`Q_0^{\mathrm{RPA}}(s)`$ in Eq.(13). Our system is not spherically symmetric but we utilize this just as a convenient method in order to evaluate $`\mathrm{\Omega }_R^{\mathrm{SF}}[g_R]`$ only through quantities $`\chi _R^{0\mathrm{L}}`$ and $`\chi _R^{0\mathrm{T}}`$ in addition to $`I_R^\mathrm{L}`$ and $`I_R^\mathrm{T}`$.

A key point in our method is in the SC determination of $`I_R^\mathrm{L}`$ and $`I_R^\mathrm{T}`$. Here we only treat the paramagnetic ground state at finite-temperature. At first, let us consider the case without $`\overline{\mathrm{\Omega }^{\mathrm{SF}}}[g]`$; this means that $`\mathrm{\Omega }[g]`$ in Eq.(2) is reduced to be the crude finite-temperature LDA. Then $`D_{mR}^\sigma (\omega )`$ in each site $`R`$ is given as a functional of $`B_R^{\mathrm{ext}}`$ as the solution of an impurity problem; adding the perturbation $`B_R^{\mathrm{ext}}`$ only at a site $`R`$. This problem is easily treated in AkaiKKR. Then $`\chi _R^{0\mathrm{L}}(0)`$ is calculated through Eqs.(20,21) from $`D_{mR}^\sigma `$ as a function of $`B_R^{\mathrm{ext}}`$. On the other hand, we can also calculate the longitudinal local static spin susceptibility at the site $`R`$ from the numerical 2nd derivative of $`\mathrm{\Omega }[g]`$ with respect to $`B_R^{\mathrm{ext}}`$ as $`\chi _R^\mathrm{L}(0)=\widehat{S}_{zR}\widehat{S}_{zR}={\displaystyle \frac{1}{4\mu _\mathrm{B}^2}}{\displaystyle \frac{M_R}{B_R^{\mathrm{ext}}}}`$. Note that these $`\chi _R^{0\mathrm{L}}(0)`$ and $`\chi _R^\mathrm{L}(0)`$ are calculated as functions of $`B_R^{\mathrm{ext}}`$. At the same time, we also calculate $`\chi _R^\mathrm{T}(0)={\displaystyle \frac{1}{4\mu _\mathrm{B}^2}}{\displaystyle \frac{M_R}{B_R^{\mathrm{ext}}}}`$ and $`\chi _R^{0\mathrm{T}}(0)`$; we evaluate $`\chi _R^\mathrm{T}(0)`$ as if our system is spherical symmetric for simplicity. Then we can calculate $`I_R^{0\mathrm{L}}`$ so that $`\left\{\chi _R^\mathrm{L}(0)\right\}^1=\left\{\chi _R^{0\mathrm{L}}(0)\right\}^1I_R^{0\mathrm{L}}`$. With this $`I_R^{0\mathrm{L}}`$, we can calculate the longitudinal part of $`\mathrm{\Omega }_R^{\mathrm{SF}}`$. The same method is used also for $`I_R^{0\mathrm{T}}`$ so as to be $`\left\{\chi _R^\mathrm{T}(0)\right\}^1=\left\{\chi _R^{0\mathrm{T}}(0)\right\}^1I_R^{0\mathrm{T}}`$. As for $`\mathrm{\Omega }_R^{\mathrm{SF0}}`$, we can calculate it in the same manner above from the same $`D_{mR}^\sigma (\omega )`$, $`\chi _R^\mathrm{L}(0)`$, and $`\chi _R^\mathrm{T}(0)`$, but taking the $`T0`$ limit in Eq.(13) through the Matsubara frequencies. As a result, we can obtain $`\overline{\mathrm{\Omega }_R^{\mathrm{SF}}}(M_R)`$ as a function of $`M_R`$. This $`\overline{\mathrm{\Omega }_R^{\mathrm{SF}}}(M_R)`$ should give an effect to make the system paramagnetic above $`T_\mathrm{c}`$. We should have to include the effect from the beginning, that is, we rather have to start from $`\mathrm{\Omega }[g]`$ including $`\overline{\mathrm{\Omega }^{\mathrm{SF}}}[g]=_R\overline{\mathrm{\Omega }_R^{\mathrm{SF}}}(M_R)`$. This means that we can obtain a new function $`\overline{\mathrm{\Omega }_R^{\mathrm{SF}}}(M_R)`$ if we start from $`\mathrm{\Omega }[g]`$ with a trial function $`\overline{\mathrm{\Omega }_R^{\mathrm{SF}}}(M_R)`$. We impose a self-consistency condition, namely the agreement between these functions. In practice, we assume a simple form $`\mathrm{\Omega }_R^{\mathrm{SF}}=\mathrm{\Omega }_0^{\mathrm{SF}}+\frac{1}{2}\alpha ^{\mathrm{SF}}(M_R/\mu _\mathrm{B})^2`$. Then our self-consistency is just for the parameter $`\alpha ^{\mathrm{SF}}`$, as the trial $`\alpha _{\mathrm{in}}^{\mathrm{SF}}`$ and the calculated $`\alpha _{\mathrm{out}}^{\mathrm{SF}}`$ are the same. This assumption is necessary so as to make the SC equation determining $`\mathrm{\Omega }_R^{\mathrm{SF}}(M_R)`$ be closed. We can take the self-consistency used here as a simplified version of the general self-consistency scheme in the construction of the LW functional; the 2nd derivative of $`\mathrm{\Omega }[g]`$ with respect to $`g_\sigma (1,2)`$ can give the fluctuation (or 2-body propagators); on the other hand, the fluctuation can be used to construct $`\mathrm{\Omega }_{\mathrm{xc}}[g]`$ through the coupling-constant-integral formula . In the STLS theory  for homogeneous electron gas and the SCR theory , other simplified versions of the general scheme were used.

In Fig.1, we give an example of plots of $`\alpha _{\mathrm{in}}^{\mathrm{SF}}`$ vs. $`\alpha _{\mathrm{out}}^{\mathrm{SF}}`$, determining $`\alpha ^{\mathrm{SF}}`$. In our formalism, the interaction $`I_R^{0\mathrm{L}}`$ has a part that is linearly-dependent on $`\alpha _{\mathrm{in}}^{\mathrm{SF}}`$. The most divergent part in Eq.(13) occurs at $`\nu _0=0`$ in $`\mathrm{log}(1I_R^{0\mathrm{L}}\chi _R^{0\mathrm{L}}(0))=\mathrm{log}\left(\left\{\chi _R^\mathrm{L}(0)\right\}^1\chi _R^{0\mathrm{L}}(0)\right)`$; $`\mathrm{log}\left(\chi _R^\mathrm{L}(0)\right)=\mathrm{log}\left(\overline{\chi _R^\mathrm{L}(0)}+4\alpha _{\mathrm{in}}^{\mathrm{SF}}\right)`$, where $`\overline{\chi _R^\mathrm{L}(0)}`$ denotes $`\chi _R^\mathrm{L}(0)`$ at $`\alpha _{\mathrm{in}}^{\mathrm{SF}}=0`$. It diverges at $`\alpha _{\mathrm{in}}^{\mathrm{SF}}\frac{1}{4}\left\{\overline{\chi _R^\mathrm{L}(0)}\right\}^1`$. This occurs at $`\alpha _{\mathrm{in}}^{\mathrm{SF}}=0.0139`$ Ry denoted with the arrow in Fig.1. In the case of Fe shown in Fig.1 and fcc-Co, $`\left\{\overline{\chi _R^\mathrm{L}(0)}\right\}^1`$ is positive even at $`T=0`$ because LDA gives a local static magnetic moment . Therefore our treatment gives $`\alpha ^{\mathrm{SF}}\frac{1}{4}\left\{\overline{\chi _R^\mathrm{L}(0)}\right\}^1>0`$ at $`T0`$, thus $`\chi _R^\mathrm{L}(0)\mathrm{}`$. Thus $`\mathrm{\Omega }^{\mathrm{SF}}`$ keeps $`\chi _R^\mathrm{L}(0)>0`$ and forbids the existence of the static local magnetic moment at $`T=0`$. This is physically reasonable as a ground state character, though we still have an unphysical divergence at $`T0`$. In these cases, the results does not coincide with the LDA limit. On the other hand, as for Ni and Pd, which show negative $`\left\{\overline{\chi _R^\mathrm{L}(0)}\right\}^1`$, it coincides with the LDA limit with $`\alpha ^{\mathrm{SF}}0`$ at $`T0`$.

The partial DOS $`D_{mR}^\sigma (\omega )`$ are calculated with $``$ 8000 $`k`$ point in the irreducible Brillouin zone (BZ), where each poles of eigenvalues are effectively smeared with Lorentzian functions $`1/(\omega ^2+\mathrm{\Delta }_\mathrm{L}^2)`$. We needed to perform calculations with $`\mathrm{\Delta }_\mathrm{L}=0.009,0.012,0.016,0.022`$Ry cases and to take an extrapolation $`\mathrm{\Delta }_\mathrm{L}0`$. After we obtain the self-consistently determined $`\alpha ^{\mathrm{SF}}`$, we can calculate the bulk susceptibility $`\chi `$ for the whole systems. In Fig.2, we summarise our results. Used lattice constants $`a=`$ 5.47, 6.70, and 6.80 a.u. for bcc Fe, fcc Ni, and fcc-Co (experimental values at $`T_c`$); $`a=`$7.35 a.u. for fcc Pd (experimental value at $`T=0`$). They are shown with solid lines. In addition, we show the cases for 3% smaller lattice constants by the broken lines; the differences from solid lines indicate errors due to LDA because the ordinary LDA calculations give rise to errors of this level in the lattice constants. The bottom panels show $`\alpha ^{\mathrm{SF}}`$. The top panels are for the inverse of the bulk magnetic susceptibilities $`\chi ^1`$ calculated from these $`\alpha ^{\mathrm{SF}}`$, where we also show experimental results taken from Ref.. The middle panels are for $`\left\{\chi _R^{0\mathrm{L}}(0)\right\}^1`$,$`\left\{\chi _R^\mathrm{L}(0)\right\}^1`$ and $`I=\left\{\chi _R^{0\mathrm{L}}(0)\right\}^1\left\{\chi _R^\mathrm{L}(0)\right\}^1`$.

As for Pd, our LDA calculations without $`\overline{\mathrm{\Omega }^{\mathrm{SF}}}`$ in ASA at $`T=0`$ give ferromagnetism for $`a`$=7.35 a.u. and paramagnetism for $`a=7.13`$ a.u. In Ref., the LDA calculation with $`a=7.35`$ gives paramagnetic phase. This discrepancy might arise from the difference in the formalism and details of implementation . However, the difference seems not important for the overall discussions here. The temperature dependence of $`\chi ^1`$, especially slope at a high temperature region, is rather well described in spite of no tunable parameters. It fails to describe $`\chi ^1`$ around $`T=0`$ where we expect only SF around $`𝐤0`$ dominates the behaviour. It is reasonable because our formalism could work well only in the case that SF in each site behaves independently; in other words, all the modes of SF are excited rather homogeneously in BZ. As $`T`$ becomes higher, $`\alpha ^{\mathrm{SF}}`$ increases smoothly, which suppresses SF at $`𝐤=0`$ ($`\chi ^1`$ gets larger). On the other hand, the local SF $`\chi _R^{0\mathrm{L}}(0)`$ does not change much. The enhancement factor on local SF $`\chi _R^\mathrm{L}(0)/\chi _R^{0\mathrm{L}}(0)`$ due to the interaction $`I_R^\mathrm{L}`$, is not so large $`1.5`$ at $`T=0`$ and gradually decreases for larger $`T`$. In Ni, our method can also describe experimental $`T_c`$ and the slope of $`\chi ^1`$ reasonably. $`\alpha ^{\mathrm{SF}}`$ decreases above $`T1000`$K. $`\chi ^1`$ is not so linear as experiments and contains negative-quadratic dependence as $`\chi ^1T^2`$ reflecting the behaviour of $`\alpha ^{\mathrm{SF}}`$. In both cases of Pd and Ni, it seems that we have to add positive-quadratic dependence $`\chi ^1T^2`$, so as to have better agreement with experiments. We might expect such an effect if we include $`𝐤`$-dependent SF, because then we will have larger contributions even for lower temperatures.

As for Fe and fcc-Co, calculated $`\chi ^1`$ is rather linear. Slopes of $`\chi ^1`$ are too large compared with experiments, though calculated $`T_\mathrm{c}`$’s are still reasonable. The factors $`\chi _R^\mathrm{L}(0)/\chi _R^{0\mathrm{L}}(0)`$ are rather large $`5`$ in these cases. The relation $`\left\{\chi _R^\mathrm{L}(0)\right\}^10`$ at $`T0`$ should be hold in agreement with the previous discussion. This condition keeps the temperature dependence of $`\left\{\chi _R^\mathrm{L}(0)\right\}^1`$ too large, corresponding to too large a slope of $`\chi ^1`$. As is same the as cases in Ni and Pd, we expect better agreements by taking account of $`𝐤`$-dependent SF which will remove the above condition from our formalism.

We have shown results of paramagnetic susceptibilities for some metals. They give reasonable agreements with experiments, indicating the possibilities of treatments along the present line including SF self-consistently. In particular, the results seems to indicate necessity of including the $`𝐤`$-dependence of SF for the further improvement.

We thank Dr.Aryasetiawan for valuable discussions and comments on the manuscript.