format: text

text:
# CERN-TH/2000-196 IFIC-00/41 FTUV-00-0712 A note on the practical feasibility of domain-wall fermions 11footnote 1Talk given by P. Hernández at the workshop on Current theoretical problems in lattice field theory, 2 – 8 April 2000, Ringberg, Germany.

## Abstract

Domain-wall fermions preserve chiral symmetry up to terms that decrease exponentially when the lattice size in the fifth dimension is taken to infinity. The associated rates of convergence are given by the low-lying eigenvalues of a simple local operator in four dimensions. These can be computed using the Ritz functional technique and it turns out that the convergence tends to be extremely slow in the range of lattice spacings relevant to large-volume numerical simulations of lattice QCD. Two methods to improve on this situation are discussed.

Introduction

The idea to realize 4-dimensional (4D) chiral fermions on the lattice by coupling 5D fermions to a 4D domain wall  has attracted a lot of attention in the lattice community (for a review see ). Although originally designed to construct chiral gauge theories, the idea can also be applied to lattice QCD in order to preserve the global chiral symmetry at zero quark mass . In this case, a 5D Wilson–Dirac operator is chosen with $`N`$ slices in the extra dimension and appropriate Dirichlet boundary conditions in the fifth dimension. In the limit $`N\mathrm{}`$, chiral zero modes exist as surface modes on the 4D boundaries, even at finite lattice spacing.

The bulk fermionic degrees of freedom are massive and can be shown to decouple in the continuum limit: the action of the 5D system is equivalent to the one corresponding to a 4D Dirac operator describing the boundary chiral modes ; similarly, the propagator of the boundary fields can be obtained from the inverse of the same 4D operator . The chirality of the surface modes in the limit $`N\mathrm{}`$ then follows  from the fact that, in this limit, this 4D Dirac operator satisfies the Ginsparg–Wilson (GW) relation , which implies an exact lattice chiral symmetry . Thus the 5D domain-wall construction in the limit $`N\mathrm{}`$ is completely equivalent to a 4D lattice formulation of Ginsparg–Wilson fermions, and satisfies all the properties that follow from the exact chiral symmetry . Moreover, if the continuum limit is taken in the extra dimension, this 4D formulation coincides with that using Neuberger’s fermion operator .

The introduction of an extra dimension makes domain-wall fermions more demanding numerically than standard Wilson fermions (the equivalent 4D formulation is similarly more demanding, owing to the non-ultralocality of the action). Nevertheless the advantage of preserving an exact chiral symmetry might compensate for the higher cost in some cases. An analysis in the free theory showed that the convergence to the exact operator as a function of $`N`$ is rapid . This gave rise to the hope that also in the interacting case domain-wall fermions could be used without too much computational overhead. However, in realistic simulations, there are indications that the convergence rate deteriorates rapidly at large values of the gauge coupling, and much larger values of $`N`$ are indeed needed , leading to a substantial computational cost.

In a recent paper , the problems found in practical simulations were traced back to the appearance of very small eigenvalues of a certain 4D operator, which control the rate of convergence in $`N`$. We have performed an independent study and confirm the analysis in . In addition, we discuss a new method to improve the domain-wall fermion operator, which differs from the one proposed in  and proves to work better numerically.

Five dimensional theory

In this section we establish our notation and collect some useful formulae, the derivation of which can be found in . The 5D domain wall operator we consider here is defined as

$$𝒟=\frac{1}{2}\left\{\gamma _5(_s^{}+_s)a_s_s^{}_s\right\}+M,$$
(1)
where $`s`$ denotes a lattice site in the fifth direction ($`a_ssa_sN`$), $`a_s`$ is the corresponding lattice spacing, and $`_s^{}`$ and $`_s`$ are the usual forward and backward derivatives. The operator $`M`$ in eq. (1) is obtained from the standard 4D Wilson operator by

$$M=D_\mathrm{w}m_0$$
(2)
with

$$D_\mathrm{w}=\frac{1}{2}\left\{\gamma _\mu (_\mu ^{}+_\mu )a_\mu ^{}_\mu \right\}.$$
(3)
Here $`_\mu ^{}`$ and $`_\mu `$ are the gauge covariant forward and backward derivatives and $`a`$ is the lattice spacing in the four physical dimensions $`\mu =1,\mathrm{},4`$. The mass parameter $`m_0`$ obeys

$$m_0>0,a_sm_0<2,am_0<2.$$
(4)
Note that the lattice spacings $`a_s`$ and $`a`$ can be different. The boundary conditions are fixed through

$`P_+\psi (0,x)=P_{}\psi (a_sN+a_s,x)=0,`$ (5)
where $`P_\pm \frac{1}{2}(1\pm \gamma _5)`$.

Appropriate interpolating fields for the quarks constructed out of the left and right boundary modes are

$`q(x)`$ $`=`$ $`P_{}\psi (a_s,x)+P_+\psi (Na_s,x),`$ (6)
$`\overline{q}(x)`$ $`=`$ $`\overline{\psi }(a_s,x)P_++\overline{\psi }(Na_s,x)P_{}.`$ (7)
A mass term can then be introduced by adding to eq. (1) the term

$$\frac{1}{2}m\left\{\overline{\psi }(a_s,x)P_+\psi (Na_s,x)+\overline{\psi }(Na_s,x)P_{}\psi (a_s,x)\right\}.$$
(8)
The two-point function of the quark fields is related to an effective 4D operator $`D_N`$ by

$$q(x)\overline{q}(x)=\frac{2aD_N}{aD_{m,N}},$$
(9)
with

$$D_{m,N}=(1\frac{1}{2}am)D_N+m.$$
(10)
In terms of the operators $`K_\pm `$,

$$K_\pm \frac{1}{2}\pm \frac{1}{2}\gamma _5\frac{a_sM}{2+a_sM},$$
(11)
$`D_N`$ is given by

$$aD_N=1+\gamma _5\frac{K_+^NK_{}^N}{K_+^N+K_{}^N}.$$
(12)
From eq. (12) it is straightforward to show that

$$aD\underset{N\mathrm{}}{lim}aD_N=1+\gamma _5ϵ(K_+K_{}),$$
(13)
which can also be written as

$`aD`$ $`=`$ $`1A(A^{}A)^{1/2}`$ (14)
$`A`$ $`=`$ $`a_sM(2+a_sM)^1.`$ (15)
It follows easily from this expression that $`D`$ satisfies the Ginsparg–Wilson relation, the only difference to Neuberger’s operator being the different definition of $`A`$. Indeed, Neuberger’s operator is readily obtained from eqs. (14) and (15) by taking the limit $`a_s0`$.

Similarly, in the limit $`N\mathrm{}`$, the fermion determinant of the 5D formulation can be written in terms of the determinant of $`D_{m,N}`$, up to local subtractions. The 5D formulation is thus completely equivalent to a 4D lattice formulation of Ginsparg–Wilson fermions satisfying an exact chiral symmetry.

A final necessary condition for this formulation to be an acceptable regularization of QCD is that the operator of eq. (14) is local. Indeed, it has been shown by Kikukawa  that both the operators in eqs. (14) and (12) are exponentially localized for smooth enough gauge fields, satisfying a plaquette bound .

In realistic simulations of domain-wall fermions, $`N`$ is finite. In this situation, the chiral symmetry is broken by the residual terms $`\delta DD_ND`$. It may be speculated that one could include a small additive quark mass renormalization, in order to get rid of these chirality breaking effects. This is, however, only justified by universality arguments if the subleading corrections in $`N`$ in the action are local. The result of  shows that this is indeed the case since $`\delta D`$ is also local. However it is important to stress that the exponential localization of $`\delta D`$ only sets in at distances of $`\mathrm{O}(N)`$. This can be shown already in the free case. On the other hand, in practical simulations the typical lattice sizes used are not much larger than $`N`$ and consequently $`\delta D`$ is not local at the distances probed. In this situation, a quark mass renormalization cannot cancel the chirality-breaking effects induced by $`\delta D`$.

The convergence rate in $`N`$

For gauge field configurations with a restricted plaquette value, the operator $`A^{}A`$ has a spectral gap :

$$0<uA^{}Av,$$
(16)
ensuring the exponential convergence in $`N`$ of $`D_N`$. The minimum rate of convergence is given by

$$\omega =\mathrm{min}_\mathrm{i}[\omega _i],\omega _i\mathrm{ln}\frac{1+\sqrt{\lambda _i}}{|1\sqrt{\lambda _i}|},$$
(17)
where $`\lambda _i`$ are the eigenvalues of $`A^{}A`$.

However, in realistic simulations the plaquette bound is not satisfied and it is important to study the convergence rate $`\omega `$ for the values of $`\beta `$ and $`m_0`$ at which large-scale numerical simulations are performed nowadays.

The eigenvalues $`\lambda `$, which determine $`\omega `$, can be obtained through the generalized 4D eigenvalue equation

$$a_s^2M^{}M\psi =\lambda (2+a_sM)^{}(2+a_sM)\psi .$$
(18)
It is either the minimum or maximum eigenvalue of $`A^{}A`$ that minimizes $`\omega _i`$. These eigenvalues can be obtained by minimizing (maximizing) the generalized Ritz functional

$$\frac{\psi |a_s^2M^{}M|\psi }{\psi |(2+a_sM)^{}(2+a_sM)|\psi }$$
(19)
using a straightforward generalization of the algorithm described in  777A more detailed description of the algorithm can be obtained from the authors on request..

Eigenvalues above the lowest one can be computed by modifying the operator $`M^{}M`$ in the numerator of eq. (19) in such a way that the already computed eigenvalues are shifted to larger values. For example, this can be achieved by substituting $`M^{}M`$ by $`M^{}M+_i\alpha _iM^{}M|\psi _i\psi _i|(2+a_sM)^{}(2+a_sM)`$ with $`\alpha _i(1\lambda _i)/\lambda _i`$ and $`\lambda _i`$, $`\psi _i`$ being the already computed eigenvalues and vectors. Notice that in this method no inversion of the matrix $`(2+a_sM)^{}(2+a_sM)`$ is needed.

We have studied numerically the convergence rate in the quenched approximation. We find that it is always controlled by the lowest eigenvalue $`\lambda _{\mathrm{min}}`$ of $`A^{}A`$. In Figs. 1 and 2 we show the inverse convergence rates $`\omega _i^1`$ corresponding to the five lowest eigenvalues of $`A^{}A`$ at $`\beta =5.85`$ on an $`8^316`$ lattice, and at $`\beta =6.0`$ on a $`16^332`$ lattice, respectively. In both cases we have set $`a_sm_0=am_0=1.8`$, which is a typical value used in previous simulations.

Figures 1 and 2 give a rather pessimistic view of the convergence of domain-wall fermions to the exact operator: they imply that several hundred or even thousand slices in the extra dimension would be needed to achieve a reasonable approximation. Clearly, this would render domain-wall fermions impracticable. Our results are consistent with the findings in .

It should be noted, however, that a very similar situation was found for Neuberger’s operator . Also in this case very small eigenvalues of the corresponding $`A^{}A=M^{}M`$ occur, turning the computation of its inverse square root extremely costly.

Acceleration of convergence

In the case of Neuberger’s operator, the bad convergence behaviour resulting from the low-lying eigenvalues of $`M^{}M`$ could be cured by treating these modes exactly . It is natural to look for a similar trick also for domain-wall fermions, given the similarity of the two constructions. We have found two ways of achieving this. The first method is equivalent to the one described in , so we will not give any details here. The corresponding improved 5D operator differs from the standard one by boundary terms.

We have tested the inversion of this improved operator, $`𝒟_{\mathrm{impr}}`$, by solving the linear equation $`𝒟_{\mathrm{impr}}X=Y`$ for a given source $`Y`$. As a numerical solver we have used a conjugate gradient method with a stopping criterion $`r/X<ϵ`$, where $`r=𝒟_{\mathrm{impr}}XY`$ is the residual vector and $`ϵ`$ was set to $`ϵ=10^8`$. We found that when using such a relatively low value of $`ϵ`$ the conjugate gradient method behaves very poorly: for a number of configurations at $`\beta =5.85`$ on an $`8^316`$ lattice, the norm of the residual vector developed a very long tail at rather small values of $`r\mathrm{O}(10^5)`$. This resulted in a very large number of iterations in the conjugate gradient algorithm before it converged to the desired accuracy. We suspect that some subtle cancellations occur in the improved operator leading to unexpectedly large rounding error effects.

Since this behaviour of the conjugate gradient algorithm was rather unsatisfactory, we developed an alternative improvement method. The key observation for the new improved 5D operator is that the relations in eq. (15) and eq. (14) hold true for any choice of $`M`$ as long as

$$M^{}=\gamma _5M\gamma _5,\mathrm{det}(2+a_sM)0.$$
(20)
This fact may be used to construct an improved $`M`$ for which the very low eigenvalues of $`A^{}A`$ disappear. A possible form of $`M`$ that achieves this is given by

$$a_s\widehat{M}=a_sM\underset{k,l=1}{\overset{r}{}}X_{kl}w_kw_l^{}\gamma _5,$$
(21)
where

$$\gamma _5Av_k=\alpha _kv_k,k=1,\mathrm{},r,(v_k,v_l)=\delta _{kl}.$$
(22)
and

$$w_k=(2+a_sM)\gamma _5v_k.$$
(23)
Finally

$$(X^1)_{kl}=2\delta _{kl}(\widehat{\alpha }_k\alpha _k)^1+(v_k,(2+a_sM)\gamma _5v_l).$$
(24)
The corresponding 5D operator $`\widehat{𝒟}`$ is given by eq. (1) after substituting $`M`$ by $`\widehat{M}`$. Notice that the improved operator differs from the original one also in the bulk and not just at the boundary.

After some algebra it can be shown that

$$\widehat{A}a_s\widehat{M}(2+a_s\widehat{M})^1=A+\underset{k=1}{\overset{r}{}}(\widehat{\alpha }_k\alpha _k)\gamma _5v_kv_k^{}.$$
(25)
It is now easy to see that $`\gamma _5\widehat{A}`$ has the same eigenvectors as $`\gamma _5A`$, but all eigenvalues $`\alpha _k,k=1,\mathrm{},r,`$ are replaced by $`\widehat{\alpha }_k`$. The limit $`N\mathrm{}`$ of the corresponding improved 4D operator $`\widehat{D}_N`$ is the same as that of the original $`D_N`$ provided $`\mathrm{sign}(\widehat{\alpha }_k)=\mathrm{sign}(\alpha _k)`$. However, the approach to this limit is faster for $`\widehat{D}_N`$ if the lowest eigenvalues of $`\widehat{A}^{}\widehat{A}`$, $`\widehat{\lambda }_k\widehat{\alpha }_k^2`$, are larger than those of $`A^{}A`$, i.e. if $`|\widehat{\alpha }_k|>|\alpha _k|`$.

The concrete choice of $`|\widehat{\alpha }_k|`$ has to be taken with some care to optimize the convergence of the inverter. For example, taking $`|\widehat{\alpha }_k|=1`$ led to a bad convergence behaviour of the conjugate gradient algorithm. It is our experience that choosing $`|\widehat{\alpha }_k|`$ not much higher than $`|\alpha _r|`$, $`r`$ being the index of the largest eigenvalue projected out, see eq. (22), leads to a normal behaviour of the conjugate gradient algorithm.

As an example of the effect of the improvement using $`\widehat{M}`$, eq. (21), we show in Fig. 3 the behaviour of the pion propagator at zero distance $`\mathrm{\Gamma }_\pi (0)`$ at $`\beta =5.85`$ on an $`8^316`$ lattice and for a quark mass, $`am=0.1`$. A similar behaviour is obtained for the pion propagator at other distances.

Already the projection of only three low-lying eigenvalues is sufficient to accelerate the convergence substantially: similar approximations to the $`N\mathrm{}`$ limit are obtained for $`N150`$ in the unimproved case and $`N30`$ in the improved one. It would, of course, be interesting to see the effect also on other physical quantities.

Conclusions

In this note we presented numerical evidence that in practical simulations domain wall fermions need an unacceptably large number of slices in the extra dimension to ensure that the chiral symmetry-breaking terms are suppressed. The reason is that very small eigenvalues of a 4D operator appear, which are directly related to the convergence rate of the 5D domain-wall operator. These results confirm the findings in .

As in the case of Neuberger’s operator, it is possible, however, to separate a number of eigenvalues of the 4D operator and treat them exactly or shift them to larger harmless values. We tested two different implementations of this idea. The first has already appeared in , the second, which is described in detail above, is new. We found that numerically the second implementation works much better.

It is our overall impression, however, that there is no particular advantage in using domain-wall fermions instead of Neuberger’s operator. Theoretical considerations demonstrate that both approaches to realize a chiral symmetry on the lattice are equivalent, to the extent that they satisfy the Ginsparg–Wilson relation. However, in practical applications, it is our present experience that it is easier to control chiral symmetry violations with Neuberger’s operator.

Acknowledgements
We would like to thank Giulia de Divitiis for her help at an early stage of this project and the computer centres at NIC (Jülich) and CIEMAT (Madrid) for providing computer time and technical support.


text:
# The Relative Age of the Thin and Thick Galactic Disks

## 1 Introduction

The dating of the oldest stars in the Milky Way allows us to infer the early history of star formation in our Galaxy, and hence, provides us with important information regarding the formation of the Milky Way. A great deal of attention has been paid to the ages of globular clusters (see for example the reviews by Stetson, VandenBerg & Bolte (1996); Sarajedini, Chaboyer & Demarque (1997)) in the halo and there have been a number of studies which determine the ages of the oldest open clusters in the thin disk (e.g. Carraro, Girardi & Chiosi (1999); Chaboyer, Green & Liebert (1999, hereafter CGL)). Typically, there appears to be a gap of several Gyr between the ages of the globular clusters and the ages determined for the oldest stars in the thin disk. We have elected to investigate this further, in order to quantify the age difference (if any) between the thin disk, the thick disk and the halo of the Milky Way.

In recent years, numerous techniques have been employed to determine the ages of the oldest stars in the Galactic (thin) disk. Jimenez, Flynn & Kotoneva (1998) examined the color magnitude diagram (CMD) of Hipparcos field stars and arrived at a minimum disk age of 8 Gyr; Carraro, Girardi & Chiosi (1999) used isochrone fits to determine the ages of several old open clusters. Bergeron, Ruiz & Leggett (1997) and Knox, Hawkins & Hambly (1999) used observations of local white dwarfs and theoretical white dwarf cooling curves to determine a local disk age of $`6.510`$ Gyr and $`913`$ Gyr respectively. Oswalt et al. (1996) also used observations of white dwarfs to place a minimum age of 9.5 Gyr on the galactic disk. CGL studied the old open cluster NGC 6791 and determined an age of $`8.0\pm 0.5`$ Gyr.

Globular clusters in the thick disk and halo have been the subject of a number of recent investigations. For example, Salaris & Weiss (1998) determined the ages of three thick disk globular clusters, 47 Tuc, M71, and NGC 6352, and arrived at an age of 9.2 Gyr for all three. Harris et al. (1997) examined two halo clusters, NGC 2419 and M92, in order to determine their relative and absolute ages. Their study determined that the clusters had the same relative age to within 1 Gyr, and an absolute age of $`1415`$ Gyr.

In this paper, the relative age of the thin and thick disks is measured, comparing the open clusters NGC 188, NGC 6791 and selected Hipparcos field stars to the thick disk globular cluster 47 Tuc. The results of this analysis are related to the results of Chaboyer, Sarajedini & Armandroff (1999, hereafter CSA), in which the relative age of 47 Tuc to three halo clusters was determined using the same input physics as this study. Compared to other studies which have investigated the relative age of the thin disk and the halo (Carraro, Girardi & Chiosi, 1999, eg. ) this study has the advantage of using the same assumptions and methodologies in dating the thin and thick disks and the halo. This allows us to determine precise relative ages and to provide a reliable estimate of the error in the relative age estimates.

Section 2 of this paper will discuss the stellar model and isochrones used in this study. Section 3 details the age estimates for NGC 188 and the Hipparcos field stars (thin disk). The age of 47 Tuc (thick disk) is determined in §4. The results are summarized in §5 where the relative age between the the thin and thick disks and the halo is discussed.

## 2 Stellar Models and Isochrones

Stellar evolution tracks were constructed using Chaboyer’s stellar evolution code (CGL) for a range of masses of $`0.501.30M_{\mathrm{}}\text{ }`$ in increments of 0.05 $`M_{\mathrm{}}`$ . The input physics, including low and high temperature opacities, nuclear reaction rates, helium diffusion coefficients and the equation of state were identical to those described in CSA. The models were evolved in 6000 time steps from the zero age main sequence through the red giant branch.

The values of $`[\mathrm{Fe}/\mathrm{H}]`$ for each of set of stellar models were chosen to reflect the observed values in each of the clusters. To determine the uncertainty in the age due to the uncertainty in $`[\mathrm{Fe}/\mathrm{H}]`$, stellar models were calculated for a range of $`[\mathrm{Fe}/\mathrm{H}]`$ values, centered upon the preferred value for each of the clusters. Most of the models used a scaled solar heavy element composition (Grevesse & Noels, 1993). The models for 47 Tuc were calculated with $`[\alpha /\mathrm{Fe}]=0.0`$ (scaled solar) and $`[\alpha /\mathrm{Fe}]=+0.40`$.

The calibrated solar model for this study had a helium mass fraction of $`Y_{\mathrm{}}=0.263`$ and a heavy element mass fraction of $`Z_{\mathrm{}}=0.0182`$. The primordial helium abundance was taken to be $`Y_p=0.234`$ (Olive et al., 1997), and the helium abundance for each model was calculated using the relation $`Y=Y_p+(\mathrm{\Delta }Y/\mathrm{\Delta }Z)Z`$ with $`\mathrm{\Delta }Y/\mathrm{\Delta }Z=(Y_{\mathrm{}}Y_p)/Z_{\mathrm{}}`$.

The isochrones were constructed for ages between 5 and 14 Gyr, in 1 Gyr increments for each of the compositions. Color transformations to the observational plane were done in a manner identical to that described in CGL. Figure 1 compares the 10 Gyr isochrone for $`[\mathrm{Fe}/\mathrm{H}]=0.70`$ for both the $`\alpha `$-enhanced and the non $`\alpha `$-enhanced cases. The resulting change in the isochrone corresponds to a change in $`[\mathrm{Fe}/\mathrm{H}]`$ of nearly 0.20 dex. These results are similar to those found in Salaris & Weiss (1998).

Additional isochrones were generated in order to ascertain the effect that changes in the low temperature opacities would have on the isochrones. Figure 2 plots the 10 Gyr isochrone for $`[\mathrm{Fe}/\mathrm{H}]=0.70`$ with Alexander & Ferguson (1994) opacities, as well as the 9, 10 and 11 Gyr isochrones with Kurucz (1991) opacities used in this study. This figure demonstrates that ages determined using isochrones generated with Kurucz opacities will indicate ages about 0.5 Gyr older than isochrones using Alexander & Ferguson opacities. This analysis was done in order to compare the age of NGC 6791 as found in CGL, which used Alexander & Ferguson opacities in its models, with the ages of the clusters in this study and CSA, which utilized Kurucz opacities. Kurucz (1991) opacities were used in this study and by CSA as these opacities are available for both scaled solar and $`\alpha `$-element enhanced compositions. The Alexander & Ferguson opacities were only available to us for scaled solar compositions.

## 3 Thin Disk Ages

### 3.1 NGC 188

Isochrones were fit to the CMD for NGC 188 obtained from observations by Sarajedini et al. (1999, hereafter S99). The CMD was cross-correlated to include only cluster members with a greater than 75 percent membership probability (Dinescu, private communication). Membership criteria was based upon a proper motion study by Dinescu, et al. (1996). NGC 188 has $`[\mathrm{Fe}/\mathrm{H}]=0.05`$ based upon several determinations (Twarog, Ashman & Anthony-Twarog, 1997; Thogersen, Friel & Fallon, 1993; Hobbs, Thorburn & Rodriguez-Bell, 1990). Isochrones with $`[\mathrm{Fe}/\mathrm{H}]=0.05`$ were fit simultaneously in B V and $`VI`$ using the main sequence, allowing the distance modulus to vary around the value established by S99. Figure 3 shows the 5, 6, and 7 Gyr isochrones fit to the observational data. The reddening was taken to be $`\mathrm{E}(\mathrm{B V})=0.09`$ (S99) and transformed to $`\mathrm{E}(VI)`$ by $`\mathrm{E}(VI)=1.25\mathrm{E}(\mathrm{B V})=0.1125`$. The distance modulus which yielded an excellent main sequence fit was $`(mM)_V=11.43`$, well within the error bars for the value determined by S99 of $`(mM)_V=11.44\pm 0.08`$. The age indicated by both B V and $`VI`$ fits was 6.5 Gyr.

Isochrones for a more metal poor composition of $`[\mathrm{Fe}/\mathrm{H}]=0.15`$ were fit to the CMD, again assuming a reddening of $`\mathrm{E}(\mathrm{B V})=0.09`$. The fits produced a distance modulus of $`(mM)_V=13.34`$, slightly less than S99, and an age of 7.5 Gyr. Figure 4 shows the 6, 7, and 8 Gyr isochrones for this lower metallicity fit to the observational data. Isochrones with a higher metallicity, $`[\mathrm{Fe}/\mathrm{H}]=+0.05`$, were fit to the data as well and are shown in Figure 5, indicating $`(mM)_V=11.48`$ and an age of 5.5 Gyr. However, it would appear unlikely that a younger age as a result of higher metallicity is the case, as the CMD does not show a hook at the main sequence turn off which is seen in the isochrones.

In order to determine the uncertainty due to reddening, the $`[\mathrm{Fe}/\mathrm{H}]=0.05`$ isochrones were fit for $`\mathrm{E}(\mathrm{B V})=0.07`$ and $`0.11`$. These fits indicate ages of 7.5 Gyr and 5.5 Gyr, respectively. The three sets of isochrone fits for NGC 188 are summarized in Table 1. The best age estimate was found to be $`6.5\pm 1.0`$ Gyr, allowing for uncertainty in the metallicity and reddening of the cluster. This error estimate does not include the uncertainty in the input physics in the stellar models and isochrones, and so should be viewed as the error on the relative age of NGC 188.

### 3.2 Berkeley 17

Berkeley 17 (Be 17) has been suggested to be the oldest open cluster (Phelps, 1997). It has a metallicity of $`[\mathrm{Fe}/\mathrm{H}]=0.29\pm 0.13`$ from moderate resolution spectroscopy (Friel et al., 1995), while Carraro et al. (1999) estimate $`[\mathrm{Fe}/\mathrm{H}]0.35`$ based upon the slope of the red giant branch in the infrared111We note that the JK photometry presented by Carraro et al. (1999) contains a great deal of scatter on the main sequence and so we are unable to use this data in our isochrone fitting procedure..

Phelps (1997) obtained BVI photometry of this cluster which was used in our isochrone fits. A simultaneous fit to the B V and $`VI`$ photometry was attempted, assuming $`\mathrm{E}(VI)=1.25\mathrm{E}(\mathrm{B V})`$. We attempted to fit isochrones with $`[\mathrm{Fe}/\mathrm{H}]=0.10,0.24,`$ and $`0.40`$ with $`[\alpha /\mathrm{Fe}]=+0.0`$ and $`+0.40`$ (a total of 6 different compositions) without success. In all cases, the Be 17 main sequence was redder than our isochrones in B V and bluer than our isochrones in $`VI`$. Our best attempt at a fit is shown in Figure 6. As this is clearly not an acceptable fit to the data, we are unable to assign an age to Be 17. There are a number of possible explanations for the inability of the isochrones to simultaneously fit the B V and $`VI`$ data. These include: (1) a non-standard extinction law in the direction of Be 17 (ie. $`\mathrm{E}(VI)1.25\mathrm{E}(\mathrm{B V})`$), (2) a helium abundance significantly different than in our models, (3) errors in our isochrones (4) errors in the photometry (of order 0.05 mag in the color).

### 3.3 Hipparcos Field Stars

Open clusters are expected to dissipate/disrupt in the galactic plane. Thus, only the most massive open clusters or those with orbits which keep them far away from the plane for most of their lifetimes are expected to survive for signficant periods of time. For this reason, the age of the oldest field stars in the thin disk may provide a better estimate for the age of the thin disk than open clusters. An analysis of the age of oldest thin disk stars in the solar neighborhood was made by comparing solar metallicity isochrones to a selected sample of field stars from the Hipparcos database (ESA, 1997).

The Hipparcos database was searched for single, non-variable stars with good parallaxes ($`\sigma _\pi /\pi <0.12`$). These stars were cross referenced with the 1996 $`[\mathrm{Fe}/\mathrm{H}]`$ catalogue (Cayrel De Strobel et al., 1997) to find stars with $`[\mathrm{Fe}/\mathrm{H}]`$ within 0.07 dex of solar. Choosing stars with metallicities near the solar value should ensure that these stars are members of the thin disk. In addition, choosing a relatively narrow metallicity range makes the age determination process simple. The selected stars were then plotted on a CMD and compared to solar metallicity isochrones. As many of the stars do not have I band photometry, the comparison to the isochrones was only done in B V.

We are interested in determining the onset of star formation in the thin disk, and so concentrated our attention on the oldest stars in our sample for which it is possible to determine ages. These are the faintest stars in the main sequence turn-off or subgiant branch regions of the CMD. There were 21 stars which appeared to be old in our search of the Hipparcos database. These 21 stars were checked for other published metallicities and space velocities. Some stars were omitted from the sample which were shown in other sources to have metallicities of more than 0.07 dex from solar (Edvardsson et al., 1993; McWilliam, 1990) or shown to have space velocities which are not characteristic of the thin disk (Eggen, 1998; Edvardsson et al., 1993). The resulting sample of stars are shown in Figure 7 along with our zero-age main sequence and 6 – 8 Gyr isochrones (all with $`[\mathrm{Fe}/\mathrm{H}]=0.0`$).

Figure 7b reveals that there are 3 stars which are much brighter/redder than our main sequence isochrones. There a few possible reasons for this discrepency, including that these stars are binaries, or are more metal-rich than indicated by the 1996 $`[\mathrm{Fe}/\mathrm{H}]`$ catalogue. We will not determine the age of these three stars.

An inspection of the oldest stars in the main sequence turn-off and subgiant branch region in Figure 7b reveals that 2 stars lie on the 8 Gyr isochrone, 4 stars lie on or very near the 7 Gyr isochrone and 1 star lies between the 7 and 8 Gyr isochrone. The derived ages for the stars in the turn-off region are quite sensitive to possible errors in the colors of the stars. For example, if the B V color of the 8 Gyr star with $`M_V=4.53`$ was changed from and $`\mathrm{B V}=0.670`$ to $`\mathrm{B V}=0.655`$ then its age would be changed to 7 Gyr. Given the uncertainties in the photometry, reddening, and in our color transformation we are reluctant to base our estimate for the age of the oldest stars in the thin disk on just two stars. For this reason, we estimate that the oldest stars are $`7.5\pm 0.5`$ Gyr old. Isochrones for $`[\mathrm{Fe}/\mathrm{H}]=0.05`$ and +0.05 yielded ages of 8.0 Gyr and 7.0 Gyr respectively and we conclude that the oldest solar metallicity stars in the thin disk have an age of $`7.5\pm 0.7`$ Gyr. In determining this error, we simply added in quadrate the $`\pm 0.5`$ Gyr error due to the uncertainty in the metallicity of the stars with the $`\pm 0.5`$ Gyr uncertainty due to the photometry/reddening and color transformation.

There is a range in metallicities in the thin disk and Edvardsson et al. (1993) found that (in the local solar neighborhood) the thin disk has metallicities ranging down to $`[\mathrm{Fe}/\mathrm{H}]0.20`$. For this reason, we elected to search the Hipparcos database for single, non-variable stars with good parallaxes ($`\sigma _\pi /\pi <0.12`$) with $`0.08[\mathrm{Fe}/\mathrm{H}]0.25`$ ($`[\mathrm{Fe}/\mathrm{H}]`$ values from the 1996 $`[\mathrm{Fe}/\mathrm{H}]`$ catalogue Cayrel De Strobel et al. (1997)). This sample of stars allows us to determine the age of the somewhat metal-poor stars in the thin disk. In order to focus on the oldest stars, color and magnitude cuts were made in order to select stars with $`3.4<\mathrm{M}_\mathrm{V}<5.0`$ and $`0.54<\mathrm{B V}<0.74`$. These values were chosen based upon an inspection of our 7 Gyr and older isochrones with $`[\mathrm{Fe}/\mathrm{H}]=0.10`$ and $`[\mathrm{Fe}/\mathrm{H}]=0.24`$. These cuts were made rather generous to ensure that all potentially older stars are included in the sample. These cuts resulted in a sample of 21 stars which were plotted on a CMD with isochrones of appropriate metallicities. Stars which were less than $`5`$ Gyr old were deleted from further study. This left a sample of 16 stars.

Basic data for each of these 16 stars was retrieved from the SIMBAD database at CDS. Four stars which were classified as variables or spectrospic binaries were removed from list as it would not be possible to determine their ages. The derived ages are quite sensitive to the metallicity of the star. For this reason, the $`[\mathrm{Fe}/\mathrm{H}]`$ references were checked for each star and 3 stars whose recent (post 1995) $`[\mathrm{Fe}/\mathrm{H}]`$ determinations were significantly different than those listed in the 1996 $`[\mathrm{Fe}/\mathrm{H}]`$ catalogue were removed from the list. The basic data for each of these stars are given in Table 2. The parallaxes, proper motions and colors are from the Hippacos catalog. The $`[\mathrm{Fe}/\mathrm{H}]`$ values are from Cayrel De Strobel et al. (1997), and the radial velocities were obtained from the SIMBAD database.

The age of each of these nine stars was determined from our isochrones (interpolating between the $`[\mathrm{Fe}/\mathrm{H}]=0.10`$ and $`[\mathrm{Fe}/\mathrm{H}]=0.24`$ isochrones). The error in the derived age of each star is due to the error in the absolute magnitude (due to the error in the parallax), an assumed color error of $`\pm 0.005`$ in B V and an assumed error in $`[\mathrm{Fe}/\mathrm{H}]`$ of $`\pm 0.05`$ dex. The ages of these stars are given in Table 3, which also includes some derived kinematic data and orbital parameters for each star. The kinematic data and orbital parameters allow us to separate thin and thick disk stars.

The space motion of each star was calculated using a routine kindly provided to us by Dr. Dana Dinescu. The motions were calculated in a cylindrical coordinate system (origin at the Galactic center) by adopting a solar radius of $`R_{\mathrm{}}=8.0`$ kpc and a rotation velocity of the local standard of rest (LSR) of $`\mathrm{\Theta }_0=220`$ km/s. In this cooridinate system, the $`\mathrm{\Pi }`$ component is positive outward from the Galactic center, $`\mathrm{\Theta }`$ is positive in the direction of Galactic rotation and $`W`$ is positive towards the north Galactic pole. Errors in the derived velocities include errors in the proper motions, radial velocities and distances.

Table 3 includes some basic orbital parameters based upon intergation of the orbits in two models of the Galaxy’s potential. These integrations were made by Dr. Dana Dinescu, and full details of the integration routine may be found in Dinescu et al. (1999). The orbits were integrated in Galactic potential models given by Johnston et al. (1995) and Paczyński (1990). The two potentials yield similar orbital parameters and only the orbital parameters from the Paczyński (1990) potential are shown in Table 3. The orbital parameters listed in Table 3 are: $`L_z`$ the $`z`$-component of the angular momentum (a conserved quantity), pericentric ($`R_{per}`$) and apocentric ($`R_{apo}`$) radii, the maximum distance from the plane $`z_{max}`$, the eccentricity of the orbits $`e`$ and the inclination angle with respet to the Galactic plane $`\mathrm{\Psi }`$.

The stars in Table 3 all have $`0.14[\mathrm{Fe}/\mathrm{H}]0.25`$ and in this metallicity range one finds both thin and thick disk stars in the local solar neighborhood. The dispersion in the kinematics of the thin and thick disk make it impossible to definitavily assign a single star to either the thin or thick disk. For example, in the thin disk $`\sigma (W)=20`$ km/s while in the thick disk $`\sigma (W)=40`$ km/s (Edvardsson et al., 1993). Thus, the $`W`$ velocity of HD 41330 ($`W=25.1`$ km/s) implies that it could be member of the thick or thin disk. By considering all of the kinematic and orbital parameters ($`\mathrm{\Pi }`$, $`\mathrm{\Theta }`$, $`W`$, $`L_z`$, $`e`$ and $`\mathrm{\Psi }`$) one can classify a star as either a thin or thick disk star, bearing in mind that these classifications will never be 100% accurate. In general, thin disk stars have small $`\mathrm{\Pi }`$ and $`W`$ velocities, $`\mathrm{\Theta }220`$ km/s, $`L_z1700`$ kpc km/s, $`z_{max}0.3`$ kpc, $`e0.2`$ and small values of $`\mathrm{\Psi }`$.

In Table 3 we have indicated the most propobable classificiation of each star (thin or thick disk) based on the kinematic and orbital parameters alone. HD 52711 and 210918 are prototypical thick disk stars with rotation ($`\mathrm{\Theta }`$) velocities significantly different than the LSR, a low $`z`$-component to their orbital angular momentum ($`L_z`$) and a high eccentricity. In contrast, HD 207129 is a prototypical thin disk star with a rotation velocity and $`L_z`$ similar to the LSR, a low eccentricity, a small $`z_{max}`$ and a small angle to the Galactic plane.

Three stars which are likely thin disk stars (HD 15335, 202628 and 207129) have similar ages of $`6.4`$ Gyr. The oldest stars which appear to be thin disk stars are HD 32923 ($`10.0\pm 0.5`$ Gyr) and HD 41330 ($`9.3\pm 0.6`$ Gyr). Both of these stars have rotation velocities and $`L_z`$ similar to the LSR and low eccentricities. The $`z_{max}`$ of these stars is not too large (0.46 and 0.29 kpc) given that thin disk stars have $`0.325`$ kpc scaleheight exponentials (Majewski, 1993). The angle of their orbits to the Galactic plane ($`\mathrm{\Psi }=3.45^\mathrm{°}`$ and $`2.14^\mathrm{°}`$) are somewhat larger than typical for thin disk stars but not extremely so. As a result, we believe that it is likely that at least one of these two stars is a true member of the thin disk. These stars are located in a region of the CMD where the derived ages are relatively insensitive to the metallicities, colors and absolute magnitudes. Consequently,the derived ages have very small error bars. Averaging the ages of these two stars together we find that the oldest, somewhat metal-poor thin disk stars in the solar neighborhood have an age of $`9.7\pm 0.6`$ Gyr.

The two stars with kinematic and orbital parameters most representive of the thick disk (HD 52711 and 210918) have quite different ages of $`6.7\pm 1.0`$ Gyr and $`11.7\pm 1.1`$ Gyr. The older age is similar to 47 Tuc (see below) while the younger age would suggest that there is a considerable overlap in the ages of the thin and thick disks. However, we are reluctant to reach such a conclusion based only on one star. Two other stars which might be thick disk stars (HD 11007 and 67458) are also fairly young ($`6.6\pm 0.5`$ and $`6.5\pm 1.0`$). However, one could argue that the kinematics and orbital parameters of these two stars are not too different from the thin disk, and so there identification as thick disk stars is debateable. Further age determinations of more metal-poor stars ($`[\mathrm{Fe}/\mathrm{H}]<0.25`$) are needed before one can conclude that their is a significant spread in the age of the thick disk.

## 4 47 Tuc

### 4.1 Isochrone fitting ages

Isochrone fitting ages were determined for the thick disk globular cluster 47 Tuc. The photometric data for the cluster was obtained from Kaluzny et al. (1998). Heavy element abundances in the literature indicate values of $`[\mathrm{Fe}/\mathrm{H}]=0.70\pm 0.07`$ (Carretta & Gratton, 1997), and $`[\mathrm{Fe}/\mathrm{H}]=0.81`$ (Brown & Wallerstein, 1992). Additionally, Brown & Wallerstein (1992) indicate an enhancement in $`\alpha `$-capture elements of $`[\alpha /\mathrm{Fe}]=0.22`$. However, for the stellar models used in this study, opacities were only available for $`[\alpha /\mathrm{Fe}]=0.00`$ and $`0.40`$. Thus isochrones were fit for abundances of $`[\mathrm{Fe}/\mathrm{H}]=0.70`$ and $`0.80`$ both with $`[\alpha /\mathrm{Fe}]=0.40`$.

Isochrones were fit in a manner identical to NGC 188, simutaneously in B-V and V-I. The value of reddening was fixed at $`\mathrm{E}(\mathrm{B V})=0.04`$ (Harris, 1996) and $`\mathrm{E}(VI)=1.25\mathrm{E}(\mathrm{B V})=0.05`$. The distance modulus was varied around $`(mM)_V=11.37`$ (Harris, 1996) in order to obtain a good main sequence fit. Figure 9 shows the isochrone fits for $`[\mathrm{Fe}/\mathrm{H}]=0.70`$ and $`[\alpha /\mathrm{Fe}]=0.40`$. The poor simultaneous fit is likely due to the difference in $`\alpha `$-enhancement by nearly 0.2 dex between the literature value and the models. The fits for $`[\mathrm{Fe}/\mathrm{H}]=0.80`$ and $`[\alpha /\mathrm{Fe}]=0.40`$, shown in Figure 8, likewise indicate a poor simultaneous fit to B-V and V-I. In order to better reflect the total heavy element abundance ($`Z`$) in the cluster, isochrones were generated with $`[\mathrm{Fe}/\mathrm{H}]=0.95`$ and $`[\alpha /\mathrm{Fe}]=0.40`$. This abundance was chosen to match the heavy element mass fraction for the values of $`[\mathrm{Fe}/\mathrm{H}]`$ and $`[\alpha /\mathrm{Fe}]`$ found in Brown & Wallerstein (1992), corresponding to $`Z=0.0040`$. The result, shown in Figure 10, was an acceptable fit which indicated an age of $`12.5\pm 0.5`$ Gyr for a reddening of $`\mathrm{E}(\mathrm{B V})=0.04`$ and distance modulus of $`(mM)_V=11.35`$. Both the reddening and distance modulus were in good agreement with the values quoted in Harris (1996). Uncertainty due to reddening was determined by producing fits with $`\mathrm{E}(\mathrm{B V})=0.03`$ and $`0.05`$. These fits indicated ages of $`14\pm 0.5`$ and $`11.5\pm 0.5`$ Gyr, respectively. Table 4 summarizes the isochrone sets and parameters used in fitting to 47 Tuc. The age used for comparison to the other clusters and field stars was $`12.5\pm 1.5`$ Gyr, as indicated and includes the uncertainty due to reddening.

### 4.2 $`\mathrm{\Delta }\mathrm{V}_{\mathrm{HB}}^{\mathrm{BTO}}`$ ages

Age determinations were also made using the $`\mathrm{\Delta }\mathrm{V}_{\mathrm{HB}}^{\mathrm{BTO}}`$ method described in CSA. The $`\mathrm{\Delta }\mathrm{V}_{\mathrm{HB}}^{\mathrm{BTO}}`$ determination method is theoretically robust, but requires a well defined horizontal branch not present in the CMDs of open clusters. In order to compare the $`\mathrm{\Delta }\mathrm{V}_{\mathrm{HB}}^{\mathrm{BTO}}`$ ages in CSA with the isochrone fitted ages in this paper, the age of 47 Tuc was determined using both methods, and all ages were calculated relative to 47 Tuc. The observed value for $`\mathrm{\Delta }\mathrm{V}_{\mathrm{HB}}^{\mathrm{BTO}}`$ in the CMD for 47 Tuc was 3.18 $`\pm `$0.04 mag. The theoretical value of $`V_{\mathrm{HB}}`$ was calculated using $`\mathrm{M}_\mathrm{v}(\mathrm{RR})=0.23([\mathrm{Fe}/\mathrm{H}]+1.6)+0.56`$ (CSA) and corrected for the fact that the HB was only apparent redward of the RR Lyrae instability strip in 47 Tuc. This correction was determined using the theoretical HB models of Demarque et al. (2000). Using the $`[\mathrm{Fe}/\mathrm{H}]=0.95`$, $`[\alpha /\mathrm{Fe}]=+0.40`$ isochrones yielded an age of $`13.5`$ Gyr. The $`[\mathrm{Fe}/\mathrm{H}]=0.70`$, $`[\alpha /\mathrm{Fe}]=+0.40`$ isochrones yielded an age of 12.2 Gyr. 47 Tuc has $`[\alpha /\mathrm{Fe}]=+0.22`$ (Brown & Wallerstein, 1992) and $`[\mathrm{Fe}/\mathrm{H}]=0.70\pm 0.07`$ (Carretta & Gratton, 1997) which gives a $`Z`$ value intermediate between the above 2 isochrones sets. For this reason, we average the 2 ages and adopt an $`\mathrm{\Delta }\mathrm{V}_{\mathrm{HB}}^{\mathrm{BTO}}`$ age for 47 Tuc of $`12.9\pm 0.7`$ Gyr.

The age determined from the $`\mathrm{\Delta }\mathrm{V}_{\mathrm{HB}}^{\mathrm{BTO}}`$ is $`0.4`$ Gyr older than that determined from the isochrone fits. This small difference in age is well within the estimated error in our isochrone fit age. As the $`\mathrm{\Delta }\mathrm{V}_{\mathrm{HB}}^{\mathrm{BTO}}`$ age estimate is expected to yield more accurate absolute ages than isochrone fitting, we will assume an absolute age of $`12.9`$ Gyr for 47 Tuc. However, when determining the relative age of 47 Tuc to the old thin disk, we will use the isochrone fitting age of 47 Tuc, as the open clusters and Hipparcos field stars had their ages determined using isochrone fitting.

## 5 Summary and Discussion

The age of 47 Tuc ($`12.9\pm 0.7`$) Gyr provides a reliable estimate for the age of the thick disk. The results for the various relative stellar ages determined in this study are summarized in Table 5. This table includes the age for NGC 6791 from CGL and the ages of three inner halo globular clusters (NGC 6652, NGC 1851 and M107) investigated by CSA. The results from CSA indicate that star formation in the inner halo began roughly at the same time as the thick disk. All three of the studies discussed here (CSA, CGL, and this study) used the same input physics and distance scale, allowing for a relative age comparison between all clusters in the studies. The relative age between the thin disk stars and 47 Tuc in this study and CGL were determined using the isochrone fitting ages. Relative ages between 47 Tuc and the inner halo were determined using the $`\mathrm{\Delta }\mathrm{V}_{\mathrm{HB}}^{\mathrm{BTO}}`$ method. The globular cluster 47 Tuc was used as a “bridge” between the studies, as its age was determined using both the isochrone fitting and $`\mathrm{\Delta }\mathrm{V}_{\mathrm{HB}}^{\mathrm{BTO}}`$ methods.

The relative ages determined in this paper indicate that NGC 188 was formed $`6.0\pm 1.8`$ Gyr after 47 Tuc. The oldest solar metallicity Hipparcos field stars were formed $`5.0\pm 1.7`$ Gyr after 47 Tuc. Both of these results suggest that the solar metallicity stars in the thin disk were formed signifiantly later than the thick disk. Simply averaging these results suggest that solar metallicity stars (in the solar neighborhood) started to form $`5.5\pm 1.2`$ Gyr after the thick disk. In contrast, the slightly metal-poor ($`[\mathrm{Fe}/\mathrm{H}]0.22`$) solar neighborhood stars were formed $`2.8\pm 1.6`$ Gyr after 47 Tuc. As this difference in age is less than $`2\sigma `$ it is conceivable that thin disk stars were being born in the solar neighborhood at the same time, or shortly after thick disk objects like 47 Tuc were forming. The small relative ages between the thick disk and inner halo (Table 5) provide evidence that star formation in these regions began at roughly the same time. Figure 11 shows the relative ages versus metallicity of all the clusters discussed in this study.

It is interesting to note that the very metal-rich open cluster NGC 6791 is $`4.0\pm 1.9`$ Gyr younger than 47 Tuc. This old open cluster has an orbit which makes its classification as a thin disk object somewhat problematic. In particular Scott, Friel & Janes (1995) report that this cluster has a rotation velocity $`\mathrm{\Theta }`$ which is 60 km/s lower than the thin disk and it is moving radially outward with $`\mathrm{\Pi }=100`$ km/s resulting in a orbit with a large eccentricity. However, this work assumed a considerably larger distance than was found by CGL (5.3 vs. 4.2 kpc). For this reason, we have calculated the kinematic and orbital parameters of NGC 6791 using the distance found by CGL (4.2 kpc), the radial velocity from Peterson & Green (1998) ($`47\pm 1`$ km/s) and the preliminary proper motions from Cudworth (1994). We find the following: $`\mathrm{\Pi }=88\pm 14`$ km/s; $`\mathrm{\Theta }=165\pm 10`$ km/s; $`W=6\pm 17`$ km/s; $`L_z=1266`$ kpc km/s; $`R_{apo}=9.1`$ kpc; $`R_{per}=4.2`$ kpc; $`z_{max}=0.8`$ kpc; $`e=0.37`$ and $`\mathrm{\Psi }=6.6^\mathrm{°}`$. The orbit is fairly eccentric with a low rotation velocity and a reasonably large $`z_{max}`$ suggesting that NGC 6791 is a thick disk object. However, Scott, Friel & Janes (1995) caution that clusters like NGC 6791 may simply be in the wings of the kinematic distribution of the thin disk, and its peculiar orbit is what has allowed it to survive.

Table 5 includes the median distance of the stars/clusters from the galacic center ($`R_m`$) based upon the orbits of the stars/clusters. We see that all of the thick disk and thin disk objects we have discussed in this paper cover a relatively modest range in their distance from the Galactic center ($`5.9R_m10.1`$). Thus, all of these object were born at reasonably similar Galacto-centric distances at so their ages reflect the star formation in small region of the galaxy. The 3 inner halo clusters discussed by CSA turn out to have a fairly wide range in $`R_m`$ (2 to 18 kpc). The fact that the 47 Tuc has a similar age to these clusters suggests that globular cluster formation was occuring at about the same time over a wide range of Galacto-centric distances

There have been many proposed scenarios of galactic formation. These scenarios include the free-fall collapse of a protogalactic cloud, originally proposed by Eggen, Lynden-Bell & Sandage (1962), and the merger of independently formed stellar systems (Searle, 1977; Searle & Zinn, 1978). Although no definative conclusions about galactic formation can be drawn from the results of this study, some possibilities can be presented. The similarity in age between the inner halo and thick disk indicates that if the thick disk was formed from the gravitational collapse of the halo, this collapse occurred on short timescales ($`<1`$ Gyr). There is a clear age gap of several Gyr between solar metallicity stars in the nearby thin disk and the thick disk (47 Tuc). The age gap between moderately metal-poor thin disk stars and the thick disk is much smaller, and the present uncertainties in the age determinations allow for the possibility that the thin disk formed immediately after ($`<1`$ Gyr) after the thick disk and halo.

We are extremely gratefull to Dana Dinescu for providing us with the routines to calculate space velocities of the field stars, for calculating the orbits of the stars and for cross-referencing her membership list for NGC 188 with the photometry from Sarajedini et al. We would like to think the anonymous referee and the editor Jim Liebert whose numerous suggestions led to a considerable improvement in this paper. This research was supported by NASA through the NASA/New Hampshire Space Grant, an LTSA award NAG5-9225 (BC), by a Burke Research Grant from Dartmouth College (BC) and used the SIMBAD database, operated by operated at CDS, Strasbourg, France. BC acknowledges the hospitality of the Aspen Center of Physics, where he revised the paper.


text:
# Thermal conductivity and specific heat of the linear chain cuprate Sr2CuO3: Evidence for thermal transport via spinons

## Abstract

We report measurements of the specific heat and the thermal conductivity of the model Heisenberg spin-1/2 chain cuprate Sr2CuO3 at low temperatures. In addition to a nearly isotropic phonon heat transport, we find a quasi one-dimensional excess thermal conductivity along the chain direction, most likely associated with spin excitations (spinons). The spinon energy current is limited mainly by scattering on defects and phonons. Analyzing the specific heat data, the intrachain magnetic exchange $`J/k_B`$ is estimated to be $``$2650 K.

There is a considerable theoretical interest in one-dimensional (1D) Heisenberg spin-1/2 systems because they exhibit a number of properties that are entirely dominated by quantum-mechanical behavior and have no analogues in three-dimensional systems. In particular, it has been shown that the Heisenberg $`S`$=1/2 chain represents an integrable system characterized by a macroscopic number of conservation laws. One important conserved quantity is the energy current, implying an ideal (infinite) thermal conductivity along the chains at nonzero temperatures, if perturbations from impurities, phonons, or an interchain coupling, which always lead to non-integrable models, are negligible. It is an open question, to what extent a real material may be regarded as an ideal integrable system. Probably, the most obvious evidence for the predicted anomalous heat transport is the recent observation of an unusually high quasi-1D magnon thermal conductivity in the series (Sr,Ca,La)14Cu24O41. The structure of these materials contains two building blocks with 1D character, namely CuO2 chains and Cu2O3 ladders, both oriented along the same direction. Unfortunately, the dimerisation within the chains and a non-negligible interchain interaction in this system complicate the analysis of the observed features in terms of an integrable model.

In this work, we have searched for anomalies in the thermal transport of Sr2CuO3, which is often considered as the best physical realization of the 1D Heisenberg $`S`$=1/2 model. The crystal structure of Sr2CuO3 contains chains formed by CuO4 squares sharing oxygen corners. The chains run along the $`b`$ axis and, as shown in the inset of Fig. 2, the CuO4 squares lie in the $`ab`$ plane. The intrachain exchange interaction between neighboring Cu2+ ions connected via 180 Cu-O-Cu bonds, measured as $`J/k_B`$, is between 2150 and 3000 K. The ratio $`k_BT_N/J`$, where $`T_N`$ is the 3D Néel temperature, is as small as $`2\times 10^3`$, reflecting an extremely small ratio $`J^{}/J`$, $`J^{}`$ representing the interchain interaction.

Our observations indicate an excess thermal conductivity along the chain direction, provided by quasi-1D spin excitation (spinons). According to our analysis presented below, its magnitude is limited by scattering of spinons on defects and phonons. We find no evidence for a mutual scattering between spin excitations and hence it seems to be absent or at least negligibly small, in agreement with theoretical predictions for integrable models.

The specimens used in these experiments were cut from a single crystal that had been grown by the traveling solvent floating zone method. The details of crystal growth and structural characterization are described elsewhere. For thermal transport measurements, three rectangular-bar-shaped samples of typical dimensions $`2.5\times 1\times 1`$ mm3 with the longest dimension parallel to either the $`a`$, $`b`$ or $`c`$ axis were prepared. Two additional samples, #1 and #2, cut from the same piece, were used for specific heat measurements. The thermal conductivity was measured using a conventional steady-state method as described in Ref. . A standard relaxation technique was employed for the specific heat measurements. The magnetic susceptibility $`\chi `$ was measured with a commercial SQUID-magnetometer.

Small amounts of excess oxygen are known to be present in as-grown crystals of Sr2CuO3, giving rise to a Curie-Weiss term in the temperature dependence of the magnetic susceptibility due to uncompensated Cu $`S`$=1/2 spins. In order to study the influence of excess oxygen, we annealed sample #2 at 870C for 72 h under argon atmosphere, as described in Ref. . From the results of our measurements of $`\chi (T)`$, the ratio of the number of residual spin-1/2 impurities to the total number of Cu ions was estimated to be 1.8$`\times `$10-4 for the unannealed sample #1 and 6$`\times `$10-5 for the annealed sample #2.

The results of the specific heat ($`C_p`$) measurements in the temperature range between 1.5 and 22 K are presented in Fig. 1 as a plot of $`C_p/T`$ versus $`T^2`$.

The solid lines in Fig. 1 are fits to the data above 4 K using the approximation

$$C_p=\gamma T+\beta T^3+\delta T^5.$$
(1)
The parameter values are $`\gamma `$=2.12 $`\times `$10-3 J mole-1 K-2, $`\beta `$=1.359 $`\times `$10-4 J mole-1 K-4, and $`\delta `$=1.310 $`\times `$10-7 J mole-1 K-6 for sample #1 and $`\gamma `$=2.06 $`\times `$10-3 J mole-1 K-2, $`\beta `$=1.531 $`\times `$10-4 J mole-1 K-4, and $`\delta `$=1.310 $`\times `$10-7 J mole-1 K-6 for sample #2. The sum $`\beta T^3+\delta T^5`$ in Eq. (1) is a common low-temperature approximation for the lattice specific heat. The fit values of the parameter $`\beta `$ result in values of the Debye temperature of $`\mathrm{\Theta }_D`$=441$`\pm `$10 K and 424$`\pm `$10 K for samples #1 and #2, respectively. Since Sr2CuO3 is an insulator, the linear in $`T`$ contribution is not due to itinerant electrons but is ascribed to spin degrees of freedom.

The elementary excitations of a 1D Heisenberg spin one-half antiferromagnetic system are not $`S`$=1 magnons but $`S`$=1/2 topological excitations, now commonly called ”spinons”. The applicability of the spinon model has been demonstrated experimentally for several $`S`$=1/2 chain systems, including Sr2CuO3. The corresponding specific heat at $`TJ/k_B`$ is given by

$$C_s=\frac{2Nk_B^2}{3J}T,$$
(2)
where $`N`$ is the number of magnetic ions in the system. The fit values of the parameter $`\gamma `$ give $`J/k_B`$=2620$`\pm `$100 K and 2690$`\pm `$100 K for samples #1 and #2, respectively. This may be compared with $`J/k_B`$2200 K deduced from magnetic susceptibility data and a somewhat larger value of $`J/k_B`$=2850 K, as obtained from the analysis of an optical absorption spectrum. The cited susceptibility measurements covered a wide temperature range between 5 and 800 K, whereas the absorption spectrum presented in Ref.  has been recorded at low temperatures (32 K). Hence the discrepancy between the values of $`J`$ may be ascribed to its possible decrease with increasing temperature. Our rather large low-temperature values of $`J`$ are compatible with this suggestion.

The inset of Fig. 1 reveals anomalies of the specific heat with onsets below approximately 3.5 K for both samples, indicating some sort of phase transition. For the annealed sample the anomaly is shifted to lower temperatures with respect to the peak for the as-grown sample. Both anomalies occur at lower temperatures than the Néel temperatures $`T_N`$ found by $`\mu `$SR ($`4.15<T_N<6`$ K) and neutron scattering ($`T_N=5.4`$ K) measurements, respectively. One possible explanation for this disagreement is that, besides the transition to an antiferromagnetically (AFM) ordered state at $`T_N=5.4`$ K, not reflected in $`C_p(T)`$, there is another transition at $`T_c<3.5`$ K. Recently, two subsequent magnetic phase transitions at $`T_{c1}=5.0`$ K and $`T_{c2}=1.5`$ K were observed for SrCuO2, containing similar spin-1/2 chains but assembled pairwise in arrays of zigzag chains. The more likely possibility that is consistent with our observations is that our anomalies reflect the AFM transition reported in Refs.  and , but now shifted to lower temperatures because of a smaller number of impurities in the samples. It has been shown that nonmagnetic impurities interrupting the spin-1/2 chains enhance staggered spin-spin correlations, a common feature of various low-dimensional Heisenberg spin systems. A convincing manifestation of this feature is the stabilisation of the long-range AFM order by replacing Cu2+ with non-magnetic ions in the spin-ladder system SrCu2O3 (Ref. ) and the spin-Peierls system CuGeO3 (Refs. ). Recently, a field-induced staggered magnetisation near impurities was also observed in Sr2CuO3 by NMR measurements. Therefore, it seems quite likely that the experimentally observed Néel temperatures are always enhanced via the influence of impurities and exceed the value that is given by the ratio $`J^{}/J`$ itself. Indeed, for Sr2CuO3 calculations considering only dipolar interchain coupling yield a value of $`T_N`$ as low as 0.028 K.

The results of the thermal conductivity $`\kappa `$ along the $`a`$, $`b`$, and $`c`$ axes are presented in Fig. 2.

For both directions perpendicular to the chain direction, $`\kappa (T)`$ shows a peak at $`T_{\mathrm{max}}20`$ K and a decrease with increasing temperature, tending to a $`T^1`$ variation above $`T200`$ K. This behavior is typical for phonon thermal transport. The thermal conductivity along the chain direction, $`\kappa _b`$, exhibits the same temperature dependence at $`TT_{\mathrm{max}}`$ but obviously not so at higher temperatures. We suggest that this difference is caused by an additional quasi-1D heat transport along the chain direction, provided by spin excitations. For insulators the phonon-phonon scattering mechanism leads to $`\kappa T^n`$ ($`n1`$) at $`T\mathrm{\Theta }_D`$. For layered structures, such as Sr2CuO3, where the layers are perpendicular to the $`a`$ axis, one may expect that in this temperature region the ratio between an in-plane and the out-of-plane phonon conductivity is larger than the anisotropy of $`\kappa (T)`$ along two different in-plane directions. Since the difference between $`\kappa _c`$ (in-plane) and $`\kappa _a`$ (out-of-plane) is very small, the difference between $`\kappa _b`$ and $`\kappa _c`$ (both in-plane) should even be smaller. This argument is, of course, not valid for the temperature region near and below $`T_{\mathrm{max}}`$, where the influence of sample boundaries and various defects is important and, therefore, the behavior of $`\kappa (T)`$ is, to some extent, sample-dependent. Our fitting of $`\kappa _c(T)`$ and $`\kappa _a(T)`$ of Sr2CuO3, employing the Debye model of phonon thermal conductivity in a similar way as described in Ref. , has shown that the influence of boundary scattering on the thermal conductivity is negligible at $`T50`$ K. The same calculation indicates that an enhanced phonon contribution to $`\kappa _b`$ below 50 K may mostly be traced back to an enhanced phonon mean free path which is limited by boundary and defect scattering. Based on all these arguments presented above, we make the crucial assumption that the phonon thermal conductivity at $`T50`$ K is almost isotropic.

In order to single-out the heat transport due to spin excitations, the phonon thermal conductivity, averaged over the $`a`$\- and $`c`$-directions at $`T50`$ K, was subtracted from the experimental data of $`\kappa _b`$. The resulting spin part $`\kappa _s`$ is also plotted in Fig. 2. In a first approximation, valid for any system of quasiparticles, the 1D thermal conductivity is given by the simple kinetic expression $`\kappa _s=C_sv_sl_s`$ where $`C_s`$ is the specific heat, $`v_s`$ is the velocity, and $`l_s`$ is the mean free path of the spin excitations. The velocity of spinons is $`v_s=Ja\pi /2\mathrm{}`$, where $`a`$ is the distance between the spins along the chain direction. Since $`TJ/k_B`$ still holds, Eq. (2) for the specific heat is valid, and thus the thermal conductivity of spinons is given by the simple equation

$$\kappa _s=N_sa\frac{k_B^2\pi }{3\mathrm{}}l_sT,$$
(3)
where $`N_s`$ is the number of spins per unit volume. Eq. (3) has been shown to also be valid for 1D magnon systems.

We calculated $`l_s`$ using Eq. (3) and taking into account small (maximum 6.5% at 300 K) deviations of $`C_s(T)`$ from linearity. The calculated spinon mean free path is shown in Fig. 3.

It tends to reach a constant value at low temperatures and decreases with increasing temperature. Assuming that the different scattering mechanisms act independently, the inverse total mean free path of spinons $`l_s^1`$ may be written as a sum $`l_{s,i}^1`$ of the inverse mean free paths produced by each scattering mechanism. It turns out that

$$l_s^1=AT\mathrm{exp}(T^{}/T)+L^1,$$
(4)
with the parameters $`A`$=8.2$`\times 10^5`$ m-1K-1, $`T^{}`$=186 K, and $`L`$=7.86$`\times 10^8`$ m, reproduces our results in the whole temperature region between 50 and 300 K with high accuracy, as may be seen in Fig. 3. If the first term on the right-hand side of Eq. (4) is to represent Umklapp processes, spinon-spinon scattering may be ruled out because this case would require $`T^{}`$ to be of the order of $`J/k_B`$2600 K. Since $`T^{}`$ is close to $`\mathrm{\Theta }_D/2`$=220 K, spinon-phonon Umklapp processes are the most likely choice for this contribution.

The second and constant term in Eq. (4) is attributed to scattering of spin excitations on defects interrupting the Cu-O chains. In that case the value of the parameter $`L`$ is a measure for the mean distance between the defects. Surprisingly, this distance is much shorter than the average distance of 2$`\times `$104 $`\mathrm{\AA }`$ between $`S`$=1/2 impurities, estimated from the Curie-Weiss-type term of the magnetic susceptibility. A similar discrepancy has also been encountered in the interpretation of NMR measurements on Sr2CuO3. Recently, in order to explain some peculiar features of NMR spectra of Sr2CuO3, Boucher and Takigawa introduced the concept of mobile “bond-defects”, which are not related to interstitial excess oxygen. The calculated mean distance between two neighboring bond-defects, consistent with the NMR data between 20 and 60 K, is shown in Fig. 3. We note a reasonably good overlap with our data of $`l_s`$ at $`T=60`$ K. The model of Boucher and Takigawa predicts that at lower temperatures the interaction between defects is important and the number of bond-defects decreases with increasing temperature. If our parameter $`L`$ in Eq. (4) indeed represents the distance between bond-defects, it ought to be temperature-dependent below 50 K. Unfortunately, a quantitative check of this conjecture is difficult because of the uncertain subtraction of the phonon background at $`T50`$ K. However, the idea that the bond-defects are the main source of spinon scattering at low temperatures is, at least qualitatively, consistent with the temperature dependence and the anisotropy of the thermal conductivity also below 50 K. To demonstrate this, we calculated the total thermal conductivity ($`\kappa _{\mathrm{ph}}+\kappa _s`$) along the $`b`$ axis at temperatures between 20 and 60 K assuming, first that $`\kappa _{\mathrm{ph}}`$ is isotropic also in this temperature range and equal to the average of $`\kappa _c`$ and $`\kappa _a`$ and, second that $`l_s`$ is equal to the distance between neighboring bond-defects given in Ref.  (the dashed line in Fig. 3). The resulting temperature dependence of $`\kappa _b`$ shown by the solid line in Fig. 2 is in qualitative agreement with the experiment.

Concluding this paper, we return to the question whether the present results are relevant vis à vis of integrable models. We argue that a sizeable quasi-1D thermal transport mediated by spin excitations does exist in Sr2CuO3. Its magnitude is not exceptional but the scatterers, i.e., defects and phonons, limiting the mean free path of spin excitations are extrinsic to the magnetic system. Our analysis indicates the absence or negligibly small influence of spinon-spinon scattering on the thermal conductivity, in agreement with the predictions made for integrable models. Our results imply that a dissipationless energy current, expected for systems that fulfill the assumptions of the integrable Heisenberg 1D $`S`$=1/2 model, is not robust if perturbations like defects and lattice excitations interfere. The very high value of the magnetic exchange interaction within the chains and the low temperature magnetic phase transition, identified via our specific heat measurements, confirm that Sr2CuO3 may, nevertheless, be considered as an excellent realization of a 1D $`S`$=1/2 Heisenberg antiferromagnet.

We acknowledge useful discussions with X. Zotos and F. Naef. This work was financially supported in part by the Schweizerische Nationalfonds zur Förderung der Wissenschaftlichen Forschung.


text:
# PG2131+066: A Test of Pre-White Dwarf Asteroseismology

## 1 Introduction

Asteroseismology of pre-white dwarfs (PWDs) provides fresh insights into their interior structure and bulk properties. However, because of their extreme temperatures and surface gravities, spectroscopic verification of seismic results has been limited. An independent distance determination would provide a comparison between model and observed luminosities, providing a more confident interpretation of seismic and spectroscopic results. Unfortunately, none of the GW Vir type stars are close enough to obtain a distance via parallax. However, PG 2131+066 (hereafter referred to as PG 2131) has a red main sequence companion, which provides an independent test of its distance.

PG 2131 is a member of the PG 1159 spectral class of stars (the pulsating PG 1159 stars are referred to as GW Vir stars). These hot, ($`T_{\mathrm{eff}}100,000K`$) luminous ($`L101000L_{}`$) stars represent a link in evolution between planetary nebula nuclei and white dwarfs. As such, by understanding their interior conditions, we provide constraints on both AGB and white dwarf stars.

Spectroscopy by Wesemael et al. (1985) revealed a red component which they proposed as originating from a companion star (hereafter referred to as the secondary star). Wesemael et al. (1985) deconvolved their spectrum into its constituent components by using a Rayleigh-Jeans approximation to describe the red side of the primary star spectrum. The remaining flux was attributed to the secondary, to which they assigned a spectral type of K5-M0 V, and synthesized a $`V`$ magnitude. From this, they determined a distance to the secondary of $`d=1047_{500}^{+1000}`$ pc.

PG 2131 was discovered to be a variable star by Bond et al. (1984). Pulsations of PG 2131 were investigated in 1992 by the Whole Earth Telescope (WET) collaboration. Seismological analysis by Kawaler et al. (1995) of the WET data, combined with the best available spectroscopic temperature of Teff$`=80,000\pm \mathrm{\hspace{0.17em}10},000`$ K (Werner *et al.* 1991) implied a luminosity of log($`L/L_{}`$)$`=\mathrm{\hspace{0.17em}1.0}\pm \mathrm{\hspace{0.17em}0.2}`$. Adopting a crude bolometric correction, Kawaler et al. (1995) derived a $`V`$ magnitude of 16.6, which provided a “seismic” distance to the primary of $`470_{130}^{+180}`$ pc.

PG 2131 provides an important way to “calibrate” seismic results and confirm (or challenge) the reliability of the procedures by which information is extracted from the stellar pulsations. The asteroseismological distance is the end inference of seismic analysis (depending on the other parameters such as stellar mass, luminosity, temperature etc. determined through earlier steps). Independent testing of this procedure is required to establish the reliability of this distance determination technique. Such a confirmation has been made for a cooler pulsating white dwarf, GD 358 (Bradley & Winget 1994) via trigonometric parallax, but an independent test of the distance determination of a pulsating GW Vir star is vital.

Though the two distances quoted above are within their rather large error estimates we sought to improve both estimates using the best available observations, atmospheric models, and seismic models. Bond and colleagues imaged PG 2131 in 1993 with the first-generation Planetary Camera on HST. Even with the point-spread function problems, they resolved the binary into its two components, and reported concordance between the distance determined using the red companion and the asteroseismic distance from Kawaler et al. (1995) (Bond, private communication). In this paper, we report on our photometric analysis of two HST frames that resolve the binary, compare the photometric colors with models of both stars, and reevaluate the seismological analysis. Section 2 describes the HST photometry and our analysis. Section 3 uses the photometric determinations to compare with stellar atmosphere models and obtain spectroscopic parallax distances. In Section 4, we review the seismological analysis using updated parameters for the pre–white dwarf component. We compare the results and conclude in Section 5.

## 2 Photometry

In 1993, HST obtained images of PG 2131 through the F555W and F785LP filters with WFPC1. These images were originally obtained by Howard Bond as part of a snapshot project to search for binary PN central stars. These images confirm the binary nature of PG 2131. As expected, the images show a visual companion that is as extremely red as the primary star is blue. The separation of the pair is 0.3”, at a position angle of 21 degrees.

We used the IRAF package PHOT to obtain aperture photometry on the stellar cores of both images. To get relative photometry, we chose an aperture that included the stellar core but excluded the first diffraction ring. The small angular separation of 0.3 arc seconds results in the core counts being polluted by their neighbor. This is particularly severe in the F555W image, where nearly half of the secondary core counts can be attributed to the primary. To remove this contamination, we took advantage of the symmetry of the PSF (see Figure 1). We measured the counts on the opposite side of each star from its companion, and subtracted this from the core counts of the other star; the apertures are illustrated in Figure 1c.

With no other stars in the field, we needed to approximate absolute photometry from this pair alone. To do this, we accumulated all the combined flux by using a large aperture ($`2.4\mathrm{"}`$) reaching out to points where the signal was indistinguishable from the sky background. Then we used the relative photometry to distribute the total flux to each star. This provided the “absolute” photometry given in Table 1.

In the analysis that follows, the values we obtained were subject to extensive Monte Carlo simulations. These provided most–likely values for all parameters from instrumental magnitudes through the final distance determinations. We simulated all photometric measurements (i.e. counts and count rates) using Poisson statistics. Uncertainties in color corrections and other expressions were modeled as Gaussian distributions with the standard deviations given in the sources of these relations, as cited. Errors were not propogated analytically, but were determined throughout the sequence by the distribution of values returned by the Monte Carlo trials. The errors quoted are the 1 $`\sigma `$ errors derived from the Monte Carlo probability distributions.

Using filter conversions from Holtzman et al. (1993), we iteratively determined the $`V`$ and $`I`$ magnitudes listed in Table 1 from the HST F555W and F785LP filter bandpasses. It should be noted however that the Holtzman et al. (1993) conversions do not include stars as blue or red as the components of PG 2131, but it does contain the widest range available in the literature - from $`VI`$ of 0 to approximately 1.7. We estimate conversion errors of 0.1 magnitudes for each filter.

If the secondary is a main sequence star, the $`VI`$ color of the secondary may be used to determine the absolute magnitude using an empirical main-sequence $`(M_V,VI)`$ relation. Several are available from the recent literature. Clemens et al. (1998) used nearby stars, combined with Hipparcos data, to determine such a relationship for the lower main sequence. Using their composite fit, we find that $`M_V=10.01\pm 0.68`$, yielding a distance of $`560_{134}^{+200}`$ pc. Error estimates combine the uncertainty for the relationship quoted by Clemens et al. (1998) as well as the photometric errors. It is important to note here that if the companion is not a main sequence star, but is a luminosity class III star of the same color, it would have an absolute magnitude of -0.6, placing it over 78,000 pc. away!

We also obtained an absolute magnitude for the secondary from Figure 6 of Baraffe et al. (1998). Using this figure, we find $`M_V=9.68\pm 0.1`$. With this absolute magnitude we calculate a distance for the secondary of $`688\pm 65`$ pc. Even though this error is smaller, it is based on a region where the Baraffe et al. (1998) sample contains few stars, and therefore is more prone to systematic effects that are not in the error estimate. (We use the value of $`M_V=10.01\pm 0.68`$ from the Clemens et al.  relation for the rest of the paper.)

The only published ground–based photometry of PG 2131 is a magnitude given by Green et al. (1986) of 16.63; this entry in their tabulation is denoted as uncertain. This is a multichannel $`V`$ magnitude, which approximates Johnson $`V`$. For comparison, the combined V magnitude of the two stars on the HST frame is $`16.79\pm 0.10`$. Given the large uncertainty in the Green et al. (1986) value, this represents consistency between the ground-based photometry and the PC photometry. We note that Bond et al. (2000) have obtained newer ground-based photometry for comparison with the HST images.

## 3 Model Fitting

Given the $`VI`$ color and distance for the secondary, we were ready to determine if the pair was an optical binary or a physical binary. To do this, we used the data for the secondary to calibrate models to the HST WFPC1 system. Then we determined the distance for the primary by convolving models with the now calibrated system where the only free parameter is distance.

Starting from the $`VI`$ vs $`T_{\mathrm{eff}}`$ relationship of Monet et al. (1992), we found Teff$`=3600\pm 250`$K for the secondary. Next we processed a grid of Kurucz models111The models included in our grid ranged from models corresponding to K5 V stars to M5 V stars. from the tabulation by Lejeune et al. (1997) for comparison with our derived $`VI`$ color and chose a model with $`T_{\mathrm{eff}}=3500K`$ and $`\mathrm{log}g=4.5`$ as the best match; in good agreement with our Monet et al. (1992) derived temperature.

To calibrate our synthetic HST system, we first adopted stellar radii for the two stars. For the secondary, this was a radius appropriate for an M2 V star, while the radius used for the primary is $`0.0182R_{}`$. We convolved our M2 V model with the WFPC1 system throughputs for the two filters, assuming the distance determined in Section 2. Comparing this flux with the observed flux provided the appropriate zero point so the model magnitudes match the observed magnitudes for the secondary. With these zero points, we ”observed” models for the primary, adjusting the distance until our synthetic magnitudes matched the HST observations. Using this method, we derived a distance of $`681_{137}^{+170}`$ pc for the primary.

Though Paunzen et al. (1998) claim evidence for a possible third component in this system, the HST photometry shows that such a close companion to the white dwarf would produce much redder colors than are observed for it. The spectroscopic signature of H$`\alpha `$ most likely comes from the M2V star visible on the HST images - this is reasonable considering the relative youth of this system (the “newly minted” white dwarf started out as a star of about 2.3 $`M_{}`$ about $`7.5\times 10^8`$ years ago).

## 4 Seismological Analysis

The primary of PG 2131 is a member of the GW Vir class of pulsating pre-white dwarfs. These stars are multiperiodic nonradial $`g`$-mode pulsators. The periods present allow asteroseismological determination of bulk stellar properties such as mass and rotation rate. In some cases, the period spectrum is rich enough to allow determination of the internal structural parameters such as the thickness of various composition layers. Perhaps the best example of this procedure is that shown for PG 1159 by Kawaler & Bradley (1994).

Knowledge of stellar properties from other sources (such as spectroscopy or broadband photometry) frequently provides key input in the seismic analysis. For example, a given set of pulsation periods can be matched by a model of a given mass (within a limited mass range) at discrete values of Teff. Figure 2 shows the least-squares fit between model periods and the observed periods in PG 2131 as a function of $`T_{\mathrm{eff}}`$ for four different stellar masses. Clearly, the fit is temperature dependent within a given mass, but each shows a pronounced minimum at a well-defined effective temperature. To see the quality of the fit in a different way, we present Figure 3. This figure shows the observed periods of PG 2131 as dark horizontal lines; periods of successive overtones of three models (0.60 M, 0.61 M, and 0.62 M) are shown as thinner lines. The intersection of the model periods with the observed are marked as points; a perfect fit would be a precise vertical alignment of points at the temperature of the best model. The departures from a perfect fit could be the result of mode trapping by a subsurface composition transition, as discussed in Kawaler & Bradley (1994). We did not make an effort to match these small departures with the models, but clearly these departures could be reduced with such efforts without affecting the Teff fit significantly.

Considering the range of possible masses, as limited by the spectroscopic temperature determination and its uncertainty, the period–matching constraint results in a relation between bolometric luminosity and effective temperature (as shown by the dashed line in Figure 4). In the case of PG 2131, a specific value of $`T_{\mathrm{eff}}`$, combined with the requirement of matching the observed pulsation period, breaks the degeneracy of the models and isolates a model at a given mass and luminosity.

Kawaler et al. (1995) used this procedure, along with the then–current spectroscopic temperature of 80,000K, to determine the mass and luminosity of PG 2131, and therefore its seismic distance. Later, Dreizler & Heber (1998) revised the temperature upward to 95,000K. We used the revised temperature, and the models shown in Figures 2-4, to recalculate the parameters of PG 2131. Using Dreizler’s atmospheric models, we calculated bolometric corrections to determine the “seismic” distance to PG 2131 of $`668_{83}^{+79}`$ pc. Most of the uncertainty in the seismic distance lies in the bolometric correction. Table 2 shows the properties of the best-fit seismic model of PG 2131.

## 5 Summary

Though the HST images show the 0.3” separation between the two stars, obtaining accurate photometry with these images requires separating the flux into its constituent components. We have used aperture photometry on the stellar cores and the entire field to obtain apparent magnitudes. From the photometry, and a transformation to $`VI`$ colors, we deduced that the secondary is an $`M2\pm 1`$ V star and determined it’s distance to be $`560_{134}^{+200}`$ pc. We used this data to calibrate models of the WFPC1 system for use with stellar models.

By convolving atmosphere models for the primary with a synthetic WFPC1 system, we were able to obtain a reasonable match to the observed magnitudes and colors of the primary star. The best model fit is a GW Vir type star at an effective temperature of 95,000K. For the radius of a pre–white dwarf star, the observed magnitude yields a distance of $`681_{137}^{+170}`$ parsecs, matching the distance determined using the secondary to within the uncertainties.

A re-analysis of seismic data also indicate that the best fit for the primary is a PG 1159 type star at an effective temperature of 95,000K. However, seismic analysis also provides constraints for mass, radius and luminosity. Applying a theoretical bolometric correction to the primary, the luminosity can be used to determine the distance. We find that the seismic distance is $`668_{83}^{+79}`$ pc, entirely consistent with the photometry.

Table 3 summarizes the distance estimates for PG 2131. We conclude that the system is a physical binary. Taking a weighted mean of the seismic distance to the primary and the photometric distance to the secondary provides a distance to the system of $`632_{111}^{+150}`$ pc.

Though the photometric and seismic results are comforting, the large errors are still problematic. The error budget for the spectroscopic parallax is dominated by the $`M_V`$$`(VI)`$ relation. It is not the photometry that is at fault. Our knowledge of the late main sequence is lacking. Even though Clemens et al. (1998) used Hipparcos data to derive the relation between $`M_V`$ and $`VI`$, the absolute magnitude has an uncertainty of 0.35 magnitudes. Add to that an error of 0.1 magnitudes in the filter conversions, and the slim possibility remains that this is not a physical pair. A new image taken with WFPC2 would be useful in two ways. First, if any differential proper motion is present, then that would show that the pair is not a physical pair. The new image would also greatly reduce the photometric errors, but we would still depend on tighter constraints on the late main sequence to refine the distance through spectroscopic parallax.

Still, the determination of asteroseismological distances to white dwarfs represents the “end product” of the analysis pipeline — thus any systematic or internal errors in the analysis will result in an erroneous distance. The remarkable agreement between the asteroseismological distance to PG 2131+066 and the more traditional determinations, coupled with the smaller uncertainties, suggests that the asteroseismological techniques may be the most precise and reliable distances yet obtained for these stars.

The authors wish to thank Howard Bond for numerous discussions during the course of this work. Stefan Dreizler kindly provided details of his models for PG 1159 stars. We gratefully acknowledge the National Science Foundation for support under the NSF Young Investigator Program (Grant AST-9257049) to Iowa State University, and the NASA Astrophysics Theory Program through award NAG5-4060 and NAG5-8352. The HST Archive provided a valuable resource for obtaining and reducing the images.


text:
# Low-Lying Excitations from the Yrast Line of Weakly Interacting Trapped Bosons

## Abstract

Through an extensive numerical study, we find that the low-lying, quasi-degenerate eigenenergies of weakly-interacting trapped $`N`$ bosons with total angular momentum $`L`$ are given in case of small $`L/N`$ and sufficiently small $`L`$ by $`E=L\mathrm{}\omega +g[N(NL/21)+1.59n(n1)/2]`$, where $`\omega `$ is the frequency of the trapping potential and $`g`$ is the strength of the repulsive contact interaction. The last term arises from the pairwise repulsive interaction among $`n`$ octupole excitations and describes the lowest-lying excitation spectra from the Yrast line. In this case, the quadrupole modes do not interact with themselves and, together with the octupole modes, exhaust the low-lying spectra that are separated from others by $`N`$-linear energy gaps.

The Yrast state, which is a subject of active research in nuclear physics , is the lowest-energy state of a system of particles for a given total angular momentum (AM). This state is important in that the whole of the excitation energy is used up for the rotation of the system and hence the system is at zero temperature, opening up the possibility of performing precise spectroscopic measurements of energy levels close to the Yrast state. Recently, Mottelson  has pointed out that similar problems arise in rotating Bose-Einstein condensates (BECs) of trapped atomic vapor . Two recent observations of vortices in trapped BECs  have lent a great impetus to theoretical studies of this subject .

The Yrast state is obtained by distributing a given AM $`L`$ over $`N`$ bosons so as to minimize the total energy. When the trapping potential is harmonic, there would be a huge number of such partitions for $`L1`$ and hence a huge degeneracy, were it not for the interactions . The problem thus reduces to finding how a weak interaction lifts the degeneracy, selects the Yrast state, and determines the low-lying excitations from this state. Wilkin et al.  have found that when the interaction is attractive, all the AM in the Yrast state resides in the center-of-mass motion of the system. Mottelson  discussed that low-lying excitations are well described by collective modes whose excitation operators are given by

$`\widehat{Q}_\lambda ={\displaystyle \frac{1}{\sqrt{N\lambda !}}}{\displaystyle \underset{p=1}{\overset{N}{}}}(x_p+iy_p)^\lambda ,`$ (1)
where $`\lambda `$ is an integer, and $`x_p`$ and $`y_p`$ are the position operators of the $`p`$-th particle. The corresponding excitation energy is given by

$$ϵ_\lambda =2gN\left(1\frac{1}{2^{\lambda 1}}\right),$$
(2)
where $`g`$ denotes the strength of the contact interaction. Equation (2) suggests that for the case of attractive interactions (i.e., $`g<0`$), the AM of the Yrast state is carried by the dipolar ($`\lambda =1`$) mode, in agreement with Ref. , while for the case of repulsive interactions it is carried by the quadrupole ($`\lambda =2`$) or octupole ($`\lambda =3`$) modes , as they have the greatest energy gain per unit of AM. A more elaborate study taking the mode-mode interaction into account has shown that the quadrupole modes have slightly larger energy gains than the octupole modes .

Bertsch and Papenbrock  have performed numerical diagonalization for small systems ($`N=25`$ and $`50`$), finding that the energy of the Yrast line for $`LN`$ can be given by $`E=L\mathrm{}\omega +gN(NL/21)`$, where $`\omega `$ is the frequency of the trapping potential. The corresponding eigenstate was very recently shown to exist  (see also Refs. ). This state partitions the AM equally among the particles  and smoothly crosses over to the many-body single vortex state  as $`L`$ approaches $`N`$.

In this paper, we report on the results of our extensive numerical study of this system with up to 640000 particles and for a total AM up to 30. Our primary finding is that of the lowest-lying excitation spectra from the Yrast line that arise from the pairwise repulsive interaction between octupole modes; the corresponding interaction energy is given by $`1.59gn(n1)/2`$, where $`n`$ is the number of excited octupole modes . These energy levels are separated from the other excitations, including the one discussed in Ref. , by $`N`$-linear energy gaps. We have also found that for small $`L/N`$ and sufficiently small $`L`$ all the many-body eigenstates have surprisingly large overlap ($`0.99`$) with trial wavefunctions defined by

$`{\displaystyle \underset{\lambda =1}{\overset{L}{}}}(\widehat{Q}_\lambda )^{n_\lambda }|0\mathrm{with}{\displaystyle \underset{\lambda =1}{\overset{L}{}}}\lambda n_\lambda =L,`$ (3)
where $`n_\lambda `$ is a non-negative integer and $`|0`$ describes the exact many-body ground state of the system with $`L=0`$.

The model we study is the same as that of Refs. . We consider a system of weakly-interacting bosons subject to a given AM and trapped in a parabolic confining potential. The system is isotropic in the $`xy`$ plane and strongly confined in the $`z`$ direction so that all the particles occupy the lowest-energy state in this direction. The problem is thus essentially two-dimensional. The many-body Hamiltonian for this system can be written as $`\widehat{H}=_i\widehat{h}_i+\widehat{V}`$. Here $`\widehat{h}_i`$ denotes the single-particle Hamiltonian for the $`i`$-th particle

$$\widehat{h}_i=\frac{\mathrm{}^2}{2M}\left(\frac{^2}{x_i^2}+\frac{^2}{y_i^2}\right)+\frac{M\omega ^2}{2}(x_i^2+y_i^2),$$
(4)
where $`M`$ is the mass of the boson, $`\omega `$ is the frequency of the trap, and $`\widehat{V}=(4\pi \mathrm{}^2\stackrel{~}{g}/M)_{i<j}\delta ^{(2)}(𝐫_i𝐫_j)`$ describes the interaction between particles ($`\stackrel{~}{g}`$ is the dimensionless ‘scattering length’). In this paper we consider the case of $`\stackrel{~}{g}>0`$, i.e., repulsive interactions.

The single-particle spectrum is given by $`E_{n,m}=\mathrm{}\omega (2n+|m|+1)`$, where $`n`$ and $`m`$ denote the radial and AM quantum numbers, respectively. When only low-lying states of the many-boson system are concerned, we may assume that all the particles are in states with $`n=0`$. The single-particle states are then labeled only by $`m`$ and described by $`\varphi _m(z)=(z^m/\sqrt{\pi m!})\mathrm{exp}(|z|^2/2)`$, where $`zx+iy`$, and the lengths are measured in units of $`(\mathrm{}/M\omega )^{1/2}`$. We use these single-particle states as a basis set to rewrite the many-body Hamiltonian as

$`\widehat{H}`$ $`=`$ $`\mathrm{}\omega {\displaystyle \underset{m0}{}}mb_m^{}b_m`$ (6)
$`+g{\displaystyle \underset{m_1m_4}{}}V_{m_1m_2m_3m_4}b_{m_1}^{}b_{m_2}^{}b_{m_3}b_{m_4},`$
where $`g\stackrel{~}{g}\mathrm{}\omega `$, the operators $`b_m`$ and $`b_m^{}`$ annihilate and create one boson in the single-particle state $`\varphi _m`$, respectively, and $`V_{m_1m_2m_3m_4}`$ is the two-body matrix element given as

$$V_{m_1m_2m_3m_4}=\frac{\delta _{m_1+m_2,m_3+m_4}(m_1+m_2)!}{2^{m_1+m_2}\sqrt{m_1!m_2!m_3!m_4!}}.$$
(7)
Given a total AM $`L`$, the state may be spanned by the Fock states $`|n_0,n_1,\mathrm{},n_L`$ with $`_{i=0,L}n_i=N`$ and $`_{j=0,L}jn_j=L`$, where $`n_i`$ denotes the occupation number of the $`i`$-th single particle state $`\varphi _i`$, i.e., $`b_i^{}b_i|n_0,n_1,\mathrm{},n_L=n_i|n_0,n_1,\mathrm{},n_L`$. Numerical diagonalization of the two-body interaction $`\widehat{V}`$ is performed within the Hilbert subspace subject to these constraints.

The energy eigenvalue corresponding to the $`i`$-th many-body eigenstate with total AM $`L`$ may be written as $`L\mathrm{}\omega +ϵ_{L,i}`$, where $`ϵ_{L,i}`$ denotes the interaction energy. Because the lowest-energy eigenvalue with total AM $`L`$ ($`2LN`$) is given by $`ϵ_L=gN(NL/21)`$, the unique state with $`L=0`$, which is composed of $`N`$ bosons in a single-particle state $`\varphi _0`$, is the ground state for $`gN\mathrm{}\omega `$, and it will therefore be denoted as $`|0`$. Our trial wavefunctions are constructed by acting on the ground state $`|0`$ with the collective operators, $`\{\widehat{Q}_\lambda \}`$, whose second quantized forms are given by $`\widehat{Q}_\lambda =(1/\sqrt{N})_m\sqrt{(m+\lambda )!/m!\lambda !}b_{m+\lambda }^{}b_m`$.

Figure 1 shows the energy spectra for $`N=10000`$ and $`L8`$, where the energy is measured in units of $`gN`$. Our numerical results show that all the energy eigenstates for $`L8`$ have remarkably large overlaps ($`0.994`$) with the corresponding trial wavefunctions given in Eq. (3). Among the eigenstates, only “intrinsic” ones, i.e., those that do not include dipole excitations, are shown in Fig. 1 (the only exception being the state $`\widehat{Q}_1|0`$ for $`L=1`$), as the dipole mode corresponds to the center-of-mass motion that does not affect the interaction in a harmonic trap. For example, for $`L=6`$ only the “intrinsic” eigenstates, $`(\widehat{Q}_2)^3|0`$, $`(\widehat{Q}_3)^2|0`$, $`\widehat{Q}_2\widehat{Q}_4|0`$, and $`\widehat{Q}_6|0`$, are shown among all the eigenstates for $`L=6`$, and they are identified in the ascending order of energy as $`(\widehat{Q}_2)^3|0`$, $`(\widehat{Q}_3)^2|0`$, $`\widehat{Q}_2\widehat{Q}_4|0`$, $`\widehat{Q}_1\widehat{Q}_2\widehat{Q}_3|0`$, $`(\widehat{Q}_1)^2(\widehat{Q}_2)^2|0`$, $`\widehat{Q}_6|0`$, $`\widehat{Q}_1\widehat{Q}_5|0`$, $`(\widehat{Q}_1)^2\widehat{Q}_4|0`$, $`(\widehat{Q}_1)^3\widehat{Q}_3|0`$, $`(\widehat{Q}_1)^4\widehat{Q}_2|0`$, $`(\widehat{Q}_1)^6|0`$. We should note here that two states, $`(\widehat{Q}_2)^3|0`$ and $`(\widehat{Q}_3)^2|0`$, are seen as ‘a single circle’ in Fig. 1 due to their quasi-degeneracy. The energy difference between them is caused by mode-mode interactions and is on the order of $`g`$ . Two states $`(\widehat{Q}_2)^4|0`$ and $`\widehat{Q}_2(\widehat{Q}_3)^2|0`$ for $`L=8`$ are also quasi-degenerate. The Yrast line is shown by a solid line in Fig. 1, and the other states are separated from this line by energy gaps roughly equal to or greater than $`gN/4`$. The energy difference $`gN/4`$ corresponds to the replacement of the factor $`(\widehat{Q}_2)^2`$ included in the Yrast states with $`\widehat{Q}_4`$.

In Fig.1, the collective excitations discussed by Mottelson, $`\widehat{Q}_\lambda |0`$ ($`\lambda =1,2,\mathrm{},8`$), are linked by a dashed curve. The energies of these states are given in Table I. These values are very well described by Eq. (2); in particular, for $`\lambda 3`$ the agreement is perfect within the machine accuracy .

For higher (but not very large) AM, the low-lying states are again well described by the collective excitations. In particular, the quadrupole ($`\lambda =2`$) and octupole ($`\lambda =3`$) modes well describe the low-lying, quasi-degenerate states as shown in Table II. These quasi-degenerate states for $`2L30`$ are shown in Fig. 2, where the energies are measured from the Yrast line and are on the order of $`g`$. We note that the other excitations are separated from these in Fig. 2 by $`N`$-linear energy gaps ($`gN/4`$).

However, the overlap integrals between the collective trial wavefunctions and the numerically-obtained low-lying eigenstates become small even for small $`L/N`$ as $`L`$ increases. In fact, in Table III, the deviations of these overlaps from unity seem to include not only a factor of $`L/N`$ but also another factor like $`L^s`$ ($`s2`$) for small $`L`$. Thus these trial wavefunctions describe the low-lying states well for small $`L/N`$ and sufficiently small $`L`$.

Moreover, for small $`L/N`$ and sufficiently small $`L`$, the energies of these quasi-degenerate states are found to be described by $`1.59g\times n(n1)/2`$, where $`n`$ is the number of excited octupole modes in each state . This energy ‘quantization’ can also be seen in the data for $`N=640000`$, $`L=18`$ in Table IV as $`1.5901.59\times 2(21)/2`$, $`9.5301.59\times 4(41)/2`$ and $`23.821.59\times 6(61)/2`$. Here, $`2`$, $`4`$ and $`6`$ are the number of excited octupole modes in the first, second, and third excited states for $`L=18`$, as seen in Table II.

These results show that the quadrupole modes do not interact with themselves or with the octupole mode in the limit of $`L/N0`$, while the octupole modes undergo a pairwise repulsive interaction. Although the mode-mode interactions were discussed in Refs. , such concrete quasiparticle features of these modes were not pointed out previously. Because of such quasiparticle features, the Yrast states are well described by $`(\widehat{Q}_2)^n|0`$ for even $`L`$ and by $`(\widehat{Q}_2)^n^{}\widehat{Q}_3|0`$ for odd $`L`$ in that limit ($`n`$, $`n^{}`$: integers). The lowest energy for sufficiently small AM $`L`$ is therefore simply given by the sum of the excitation energies of the quadrupole modes (and the energy of an octupole mode for odd $`L`$).

As the total AM increases, however, the interconversion between these two modes becomes significant. This is seen in Table V and is also suggested by the level repulsion between the lowest and first-excited states for larger $`L`$ in Fig. 2. Although contributions by the modes, $`\widehat{Q}_4`$, $`\widehat{Q}_5`$ and $`\widehat{Q}_6`$ to the Yrast states have been claimed , the interconversion between modes $`\widehat{Q}_2`$ and $`\widehat{Q}_3`$ makes a greater contribution than these three modes, as seen in Table V.

Whereas the Yrast state is analytically given for $`LN`$ , its relations to collective excitations for small $`L/N`$ and sufficiently small $`L`$ (especially the interconversion effects between the quadrupole and octupole modes) remain unclear. An understanding of these relations may help us clarify the integrability in a two-dimensional bose system with a harmonic potential that includes weak repulsive delta-function interactions or a hidden symmetry of the Hamiltonian .

T.N. was supported for the research in the USA by Japan Society for the Promotion of Science. M.U. acknowledges support by a Grant-in-Aid for Scientific Research (Grant No. 11216204) by the Ministry of Education, Science, Sports, and Culture of Japan, and by the Toray Science Foundation.


text: # 1 Mira Variables