text:
# Timing analysis of the isolated neutron star RX J0720.4-3125Based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA).

## 1 Introduction

RXJ0720.4-3125 is a nearby, isolated neutron star (NS) that was originally discovered with ROSAT during a systematic survey of the galactic plane (Haberl et al., 1997). The source exhibits all the common characteristics of the other six ROSAT NS candidates (hereafter dim NSs, see e.g. Treves et al., 2000 for a review): a blackbody-like, soft spectrum with $`kT80`$ eV; an exceedingly large X-ray to optical flux ratio; a low X-ray luminosity, $`L_X10^{30}10^{31}\mathrm{erg}\mathrm{s}^1`$; a low column density and no evidence for a binary companion. An optical counterpart has been identified (Motch & Haberl, 1998; Kulkarni & van Kerkwijk, 1998).

Pulsation associated with the spin period has been observed in four dim NSs: RXJ0720.4-3125 ($`P8.39`$ s), RXJ0420.0-5022 ($`P22.69`$ s), RXJ0806.4-4123 ($`P11.37`$ s) and RBS1223 ($`P5.16`$ s). The pulse shape of RXJ0720.4-3125 is almost sinusoidal (Haberl et al., 1997; Cropper et al., 2001).

The true nature of the mechanism powering the X-ray emission from dim NSs is still unclear. Until a few years ago, these sources were thought to constitute a stand alone class and two mechanisms were proposed for their emission: either accretion from the interstellar medium onto an old neutron star or release of thermal radiation from a younger, cooling object (see e.g. Pavlov et al., 1996; Treves et al., 2000; Motch, 2001 and references therein). More recently, it has been noted that the period of RXJ0720.4-3125 is somewhat unusual and falls in the very narrow range in which the anomalous X-ray pulsars (AXPs) periods cluster ($`P612`$ s, see e.g. Mereghetti, 2000). This hints at a possible evolutionary link between dim NSs, AXPs, and soft gamma-ray repeaters (SGRs, see e.g. Heyl & Kulkarni, 1998; Heyl & Hernquist, 1998; Colpi et al., 2000; Alpar, 1999, 2000, 2001). Two kind of “unified” scenarios have been then proposed. In the first one, the common mechanism powering the three classes of objects is dissipation of a decaying, superstrong magnetic field ($`B>\text{ }10^{14}10^{15}`$ G). In this case dim NSs represent the descendants of SGRs and AXPs, and RXJ0720.4-3125 may be one of the closest old magnetars. Alternatively, all three classes may contain neutron stars with lower (canonical) magnetic field ($`B10^{12}G`$) endowed by a fossil disk (Alpar, Ankay & Yazgan, 2001; Alpar, 2001). In this interpretation, dim NSs would be in the propellor phase and would be the progenitors of AXPs and SGRs, the latter having entered an accretion phase.

In any model involving magnetars, proton cyclotron features are expected to lie in the X-ray range, while electron cyclotron lines are expected in the same band for canonical magnetic field. According to Zane et al. (2001), for surface magnetic fields strengths of $`10^{14}10^{15}`$ G spectra exhibit a distinctive absorption feature at the proton cyclotron energy $`0.63(B/10^{14}\mathrm{G})`$ keV. The required resolving power is $`100`$, therefore the detection of this feature is well within the capabilities of XMM-Newton grating spectrometers. Recently, Paerels et al. (2001) have presented spectra of RXJ0720.4-3125 using XMM-Newton, and found that there is no evidence for absorption lines and edges in the X-ray spectrum. Unless different atmospheric effects may concur in suppressing the feature, when taken straightforward the absence of electron or proton cyclotron resonances in the RGS band excludes a range of average magnetic field strengths, $`B(0.32)\times 10^{11}`$ G and $`(0.52)\times 10^{14}`$ G.

Based on the same XMM-Newton observation, Cropper et al. (2001) presented the pulse-shape analysis of RXJ0720.4-3125, modelling spin pulse intensity and hardness ratio profiles. By assuming that the source of the flux variation is the changing visibility of the heated magnetic polar caps, they derived an upper limit on the cap size and showed that a polar cap larger than $`60^{}65^{}`$ can be rejected at a confidence level of 90%. Whatever the mechanism, the X-ray emitting region is therefore confined to a relatively small fraction of the star surface. It is worth noticing that such a small hot region can not be explained only in terms of non-uniformity in the distribution of the surface cooling temperature induced by an high magnetic field (Greenstein & Hartke 1983; Page 1995). The temperature gradient which is induced by the $`B`$-dependence of the thermal conducivity of the neutron star crust is relatively smooth and, even in the (unrealistic) limit case $`B\mathrm{}`$ the associated blackbody luminosity drops only by a factor 1/2 at $`60^{}`$ and by an order of magnitude at $`77^{}`$. By performing pulse phase spectroscopy, Cropper et al. (2001) also found that the observed hardness ratio is softest slightly before flux maximum. The same characteristic was later discovered by Perna et al. (2001) in the spectra of AXPs. Cropper et al. (2001) suggested two possible explanations for this effect: their best-fitting model is based on radiation beaming, while an alternative one assumes a spatially variable absorbing matter, co-rotating in the magnetosphere. The latter may be indeed the case if the star is propelling matter outward (Alpar, Ankay & Yazgan 2001).

Further information about the nature of this puzzling source can be obtained by the spin history. In magnetars magneto-dipolar radiation will cause rapid spin-down at a rate $`\dot{P}10^{11}(B/10^{14}\mathrm{G})^2/P`$ ss-1, and it has been the positive detection of a secular spin-down of this order that has suggested the association of AXPs and SRGs with ultra-magnetized NSs (Kouveliotou et al., 1998, 1999; Thompson et al., 2001). The preliminary measure of $`\dot{P}`$ published by Haberl et al. (1997) for RXJ0720.4-3125 is uncertain to a considerably large value, and does not allow spin-up and spin-down to be discriminated. Very recently, Hambaryan et al. (2002) presented the first evidence of large spin-down rate in a dim NS, RBS1223. Their measure is partially based on ROSAT data where the detection of the period was not highly significant, and also the value ($`\dot{P}=1.35_{0.67}^{+0.69}\times 10^{11}`$ s/s) is still compatible with being due to torque from a fossil disk. However, it is certainly worth noticing that, when taken face value and interpreted in terms of dipolar rotational losses, this implies the presence of a superstrong magnetic field, $`B(3.56.5)\times 10^{14}`$ G.

A similar measure of $`\dot{P}`$ for the closest pulsating candidate, RXJ0720.4-3125, is therefore of extreme importance, as well as the accurate tracking of its spin history.

This is the first opportunity to investigate the pulse timing in the XMM-Newton data, since at the time of the previous papers (Paerels et al. 2001, Cropper et al., 2001) the timing correlation files were not sufficiently well determined. Here we present a combined analysis of XMM-Newton, Chandra and Rosat data, spanning a period of $`7`$ years.

## 2 Observations

The log of observations used for this study is given in Table 1.

RXJ0720.4-3125 has been observed twice by XMM-Newton: first on 2000 May 13, during the Calibration/Performance Verification phase (Paerels et al. 2001, Cropper et al. 2001) and later on 21 November 2001, again for calibrations. For the analysis presented in this paper we have used data from all three EPIC cameras on both epochs. The data were reduced using the XMM-Newton Scientific Analysis System (SAS) version 5.1. Source counts in the range 0.12 to 1.2 keV (PN) and 0.1 to 2.0 keV (MOS) were extracted in an aperture of 40 arcsec for the PN X00a, 30 arcsec for PN X00b, and 30 arcsec for the MOS cameras. Events with patterns 0-4 (single and double events) were selected for MOS, while patterns 0-12 (all the valid PN patterns) were retained for PN. Data arriving during episodes of high background which were experienced in the 2000 May 13 observation were excluded. Also in this observation, some PN event times were noted not to be integer multiple of the frame times, but delayed by 1 second: such events (in total about 34 ks from the 50 ks) were corrected in the event list. Finally, times were converted to Barycentric Dynamical Time (TDB) at solar system barycentre.

The Chandra observations were made over three days starting on 2000 Feb 1 using the HRC-S with LETG. The events were extracted from zero-order images of the source (namely, from circles of 2 arcsec radius centered at the star’s position). The dispersed data (1st and higher orders) are strongly contaminated by background (about 50-60%) and are not useful for this task. On the other hand, the background contamination in the zero-order data is negligible, about 2-3%. Also, the zero-order data are not piled-up: this effect can appear only at larger count rates ($`10^5`$ ct/s).

The Rosat observations were taken over a period of five years in 1993, 1996 and 1998. An earlier report of the analysis of the 1993 and 1996 data is in Haberl et al. (1997). The data were re-extracted and re-analysed using the EXSAS software system to make use of the most recent knowledge of spacecraft clock corrections. The 1998 data were taken shortly before a spacecraft emergency safe-hold, and do not benefit from a clock calibration: timing corrections are therefore extrapolated from the last available calibration. Timings may therefore be incorrect up to several seconds (W. Becker, private communication): this means that the reference of this observation to the rest of the dataset is problematic. However, no significant drifts are expected over the duration of the observation itself, so that the data are still useful in a standalone analysis.

Although the countrate is too low for our main analysis, we have also extracted the SAX LECS observation of RX J0720.4–3125 from the SAX data archive (the source was not visible in the other SAX instruments) for use as an a posteriori check on our derived timing parameters. Here we used the event list from the pipeline, correcting the times to the solar system barycentre using the FTOOL earth2sun. This neglects the light travel time from Earth centre to the satellite in low Earth orbit, $`20`$ ms, which is negligible for the purposes of our a postiori check.

The major datasets in Table 1 are from the two XMM Newton observations which have high count-rates (particularly the EPIC-PN), and the long 1996 Nov 3 Rosat data. The 1998 Rosat and the Chandra observations are valuable by nature of their several day durations.

## 3 Time Series Analysis

Our data originate from instrumentation with widely differing sensitivities: typical count rates vary from 1 count every $`3`$ s for Rosat HRI to $`6`$ counts/s for XMM-Newton PN. However, none of these count rates is sufficiently high for a normal distribution of counts to be expected in a time bin of adequate time resolution. The body of literature addressing discrete Fourier Transforms (for example Deeming 1975, Scargle 1982, Schwarzenberg-Czerny 1998) is not directly applicable for our analysis. The treatment of sparse data and event list data is generally accomplished using the Rayleigh Transform (for example de Jager 1991, Mardia 1972). It should be noted that this transform is numerically the same as the discrete Fourier Transform.

It will be important in our analysis to define confidence intervals to the derived quantities, in particular the period. In the case where sufficient counts per sample are available, this is generally achieved by least squares fits to the data (for example Kurtz 1985), but this is not possible for sparse or event list data. In this case it is necessary to use the more general maximum likelihood technique, which makes no assumptions on the data distribution (Cash 1979, Lampton, Margon & Bowyer 1976). Using this, Cash (1979) defined the $`C`$\- and $`\mathrm{\Delta }C`$-statistics. The $`\mathrm{\Delta }C`$-statistic (Cash 1979) has the considerable advantage of being distributed as $`\chi ^2`$, so that the required confidence intervals can be identified directly from the plot.

The maximum likelihood technique depends on a model probability distribution, which must be specified a priori. Haberl et al. (1997) and Cropper et al. (2001) find the flux variation of RXJ0720.4-3125 to be almost sinusoidal, so in our particular case we can take a sinusoidal variation superimposed on a mean value as an appropriate model. The single frequency and absence of harmonics in the Rayleigh transform of the X00a PN dataset shown in Figure 1 provides the justification explicitly for this approach. The maximum likelihood is then computed simply as equation (35) of Bai (1992). The normalisation of the model is unimportant, as it results in a constant term in the C-statistic, which then is eliminated in the $`\mathrm{\Delta }C`$-statistic. In our case we have two free parameters: normalised amplitude and phase. In general, when performing maximum likelihood fits, the best fit free parameters have to be determined by whatever method may be most appropriate. However, in the case of sinusoids, as noted by Bai (1992), their orthogonality allows them to be determined directly from the Rayleigh Transform (or Lomb-Scargle periodogramme).

The computation of the $`\mathrm{\Delta }C`$-statistic is therefore only slightly more complicated than the periodogramme. It starts in the same way, by a scan through period, calculating the Rayleigh Transform to provide the amplitude and phase at each period. Then at each period the maximum likelihood is calculated using equation (35) of Bai (1992). Finally from equations (4) and (7) of Cash (1979), $`\mathrm{\Delta }C`$ is simply twice the difference between the maximum likelihood at a particular period, and the maximum value over the whole frequency range. This results in an inverted “periodogramme”, with minimum value of zero. The permitted period range is that between the confidence levels appropriate for $`\chi ^2`$ with the appropriate number of degrees of freedom (1 or 2 depending on whether $`\dot{P}`$ is included as a parameter – see section 4).

Details of the procedure are given in Appendix A.

With a period of 8.391 s (Haberl et al. 1997), more than $`10^7`$ cycles have elapsed over the time span of our data set. It is therefore essential to incorporate a $`\dot{P}`$ term in our analysis techniques when combining datasets at different epochs. The search then becomes a $`(P_0,\dot{P})`$ search, where $`P_0`$ is the period at the start of the first dataset (the 1993 Rosat observations), and the instantaneous period $`P=P_0+\dot{P}(tt_0)`$ is computed at each point in the grid. Here $`t_0`$ is the time at the start of the R93 run.

As a more general comment, we recommend this maximum likelihood periodogramme (hereafter MLP) for time series analysis. The uncertainty in a period determination is not evident directly from the width of a peak in a power spectrum computed by classical discrete Fourier transform, Lomb-Scargle periodogramme or Rayleigh transform, as this is set by the window function (c.f. Scargle 1982 appendix D). However, it is directly available from the MLP $`\mathrm{\Delta }C`$ statistic: the inverted peak from a sinusoidal signal in data with high signal-to-noise ratio will be narrower than the peak from low signal-to-noise data and the $`\chi ^2`$ can be read directly from the vertical axis (see e.g. Figure 2). Furthermore, as noted above, the MLP technique requires an appropriate model, which may be unknown a priori. In this case the signal can be reconstructed from the harmonic content of the Rayleigh transform or Lomb-Scargle periodogramme, or from the phase folded data, and an appropriate model generated. This model can then be used in the calculation of the MLP, incorporating all of the harmonic information explicitly.

## 4 Measure of The Spin Period

The only measure of $`P`$ and $`\dot{P}`$ presented in the literature for RXJ0720.4-3125 was given by Haberl et al. (1997) using Rosat data. The period was determined to be $`8.39115\pm 0.00002`$ seconds from the 1996 Nov 3 dataset. Based on a linear fit of the data spanning 3 years, these authors derived a 90 % confidence range of $`6.0\times 10^{12}`$ to $`0.8\times 10^{12}`$ for $`\dot{P}`$, compatible with no period change.

It is computationally unfeasible to explore this range in the $`P_0\dot{P}`$ plane in a single analysis of the entire dataset to determine refined values for $`P_0`$ and $`\dot{P}`$. A reduction in the range is required. This can be accomplished initially by extending the Haberl et al. (1997) analysis to include the XMM-Newton data. This is sufficiently accurate to identify the correct peak in the alias patterns in the extended duration R98 and Ch00 data, leading to a further limiting of the permitted range in the $`P\dot{P}`$ plane. The next stage is a grid search in $`P`$ and $`\dot{P}`$ over this restricted range for the 1993 and 1996 combined Rosat data only. The final stage is a grid search using all data in the restricted zone for which there are acceptable solutions from this Rosat analysis.

### 4.1 Individual datasets

We begin by performing an MLP assuming $`\dot{P}=0`$ on each of the longer pointing (R93, R96d, X00a, X00b) in the period range (8.375-8.400) s, incorporating the value of 8.39115 s found by Haberl et al. (1997). The results are shown in Figure 2 (left). As we can see, there is no ambiguity in the period determinations from the R96d, X00a and X00b datasets; the best fit frequencies and their uncertainties are reported in Table 2. A linear least square fit using the 68% formal errors from the MLP results in a value of $`P_0=8.39113\pm 0.00011`$ s, $`\dot{P}=0.0\pm 5.5\times 10^{13}`$ s/s, where, as throughout, $`P_0`$ is referenced to the start of the R93 run.

The proximity in time between the Ch00 and X00a observations, together with the above upper limit in $`\dot{P}`$, permits an unambiguous determination of the second highest peak in the Ch00 power spectrum (Figure 2, right). The same is true for the R98 observations, in which the highest peak is selected. Adding these to the linear least squares fit results in a value of $`P_0=8.39107\pm 0.00005`$ and $`\dot{P}=2.7\times 10^{13}\pm 2.5\times 10^{13}`$. As may be expected, the $`P_0`$ and $`\dot{P}`$ are highly correlated. The 68, 90 and 99% confidence intervals are shown in Figure 3.

The duration and high count-rate of the X00a PN run alone provides a strong constraint in the $`(P_0,\dot{P})`$ plane. The MLP 68% and 90% confidence intervals are also shown in Figure 3. This indicates that acceptable values for $`P_0`$ and $`\dot{P}`$ lie within the range 8.39106 to 8.39115 s and $`1`$ to $`+2\times 10^{13}`$ s/s, respectively. We checked the effect of a $`10^7`$ error in the XMM-Newton clock calibration, as reported recently (Kuster, M., 2001, reported at “New Visions of the X-ray Universe in the XMM-Newton and Chandra Era”). The change in the confidence region is negligible.

### 4.2 Combined Rosat dataset

The $`P_0`$ $`\dot{P}`$ range identified in Figure 3 is sufficiently restricted for an MLP to be performed on the combined 1993 and 1996 Rosat datasets. The use of combined datasets requires coherent phasing to be maintained over the whole dataset (coherent phasing was required only within each of the individual datasets used until this point). The $`(P_0,\dot{P})`$ plane was searched over between $`8.39095P_08.39120`$ s, and $`4\times 10^{13}\dot{P}8\times 10^{13}`$ s/s, exceeding the parameter range for the 90% confidence limits from the linear least squares fit. The resulting 68% and 90% confidence contours break up into small regions (aliases) distributed along lines in the $`(P_0,\dot{P})`$ plane. These are also shown in Figure 3

As is evident, there is overlap between the Rosat 90% confidence interval contours and the 99% contours of the linear least squares fit to the individual data (but not quite at the 90% level). This overlap region is also consistent with the X00a PN confidence intervals.

### 4.3 Full dataset

With the further restriction available from the fit to the Rosat 1993 and 1996 data subset in Figure 3, a coherent MLP search was made to all of the data (excluding the uncertain R98 run). The search was made for $`8.391050P_08.391150`$ s and $`3\times 10^{13}\dot{P}3\times 10^{13}`$ s/s within the parallel lines $`P_0=a\dot{P}+b`$ defined by intercepts $`b=8.391097`$ s and 8.391105 s and slope $`a=1.75\times 10^8`$ s, encompassing the Rosat 1993/1996 contours. Care was taken to ensure adequate sampling of the $`(P_0,\dot{P})`$ plane. This identified two pairs of values, $`(P_0,\dot{P})`$ which cannot be discriminated between on statistical grounds. These are $`(P_0,\dot{P})=8.39109273`$ s, $`5.409\times 10^{14}`$ s/s and $`8.39109148`$ s, $`3.749\times 10^{14}`$ s/s (Table 3). All other periods can be excluded at the 95.4% level, but two other possibilites cannot be excluded at the 99% level: these are also given in Table 3 for completeness. The contours for the two best fits are given in Figure 4.

We have folded the data on both $`P_0\dot{P}`$ solutions (1) and (2) in Table 3 to check that the relative phasing of all individual runs is correct. We show the folded data on these periods in Figure 5. The data from runs R96a and R96b have been combined to increase the signal-to-noise ratio, and we have not shown the data from R98 because of the uncertainty in the time reference for this run (section 2).

Now performing the a postiori check with the BeppoSax data, we find that these data phase correctly with the main datasets for $`P_0\dot{P}`$ (1), but not (2). This is shown in the lowest panel of Figure 5 in each case. This allows us to select $`P_0\dot{P}`$ (1) as the most likely timing parameters for RX J0720.4–3125. In any case, both acceptable fits to the data have $`3\times 10^{14}<\dot{P}<6\times 10^{14}`$. For the purposes of our further discussion, the difference between the $`\dot{P}`$ in (1) and (2) is not important.

## 5 Discussion

The refined value of $`\dot{P}`$ reported in this paper is consistent with, but two orders of magnitude lower than the extrema of the range reported by Haberl et al. (1997). The first implication is that RXJ0720.4-3125 is unlikely to be spinning down under propeller torques. If the propeller is due to the torque exerted by the interstellar medium, it would be (e.g. Colpi et al., 1998)

$$\dot{P}_{prop}10^8n^{9/13}v_{10}^{27/13}B_{12}^{8/13}P^{21/13}\frac{\mathrm{s}}{\mathrm{yr}},$$
(1)
where $`n`$ is the external density in cm-3 ($`n1`$ for the interstellar medium), $`B_{12}=B/(10^{12}\mathrm{G})`$ and $`v_{10}`$ is the star’s velocity normalized to 10 km/s. By using the values of $`P`$ and $`\dot{P}`$ of RXJ0720.4-3125, we obtain $`B_{12}16n^{9/8}v_{10}^{27/8}`$. Propeller spin-down dominates over dipolar losses if $`B_{12}<25\sqrt{n}v_{10}^{3/2}`$ and this, combined with the previous equation, constrains the star’s velocity to extremely small values: $`v_{10}n^{1/3}`$. Moreover, this scenario makes the matching of thermal and dynamical time scales more difficult. On the other hand, the fact that the propeller contribution is negligeable supports larger values of the star’s velocity.

In the “unified” model suggested by Alpar (2001), the observed spin-down is associated with propeller effect by a fossil disk. In this case, under the assumption that the X-ray luminosity of the source is supplied by energy dissipation (frictional torque), upper and lower bounds on $`\dot{\mathrm{\Omega }}`$ can be derived by eq.(2) of Alpar (2001). For RXJ0720.4-3125 the observed luminosity is $`L2\times 10^{31}(d/100\mathrm{pc})^2`$ erg cm-2 s-1, where $`d_{100}=d/100\mathrm{pc}`$ and $`d`$ is the distance (Haberl et al. 1997). This gives:

$$2\times 10^{12}d_{100}^2<\text{ }\dot{\mathrm{\Omega }}<\text{ }2\times 10^{10}d_{100}^2\frac{\mathrm{rad}}{\mathrm{s}^2},$$
(2)
that in turn translates in $`\dot{P}=(P^2/2\pi )\dot{\mathrm{\Omega }}`$

$$2\times 10^{11}d_{100}^2<\text{ }\dot{P}<\text{ }2\times 10^9d_{100}^2\frac{\mathrm{s}}{\mathrm{s}}.$$
(3)
While this scenario is still consistent with the spin-down measured for RBS 1223 (Hambarayan et al., 2002), the value of $`\dot{P}`$ reported here for RXJ0720.4-3125 is well below this range. Thus, for RXJ0720.4-3125, the associated energy dissipation cannot, alone, account for the source luminosity (unless the source is at $`5`$ pc, which is unrealistic). That also make less plausible an interpretation of the hardness ratio profile in terms of a spatially variable absorbing matter, co-rotating in the magnetosphere (Cropper et al. 2001). The observed behaviour is more probably explained by the angle-dependent properties of the emitted radiation.

On the other hand, the slow spin-down rate of RXJ0720.4-3125 is still considerable. The other plausible mechanism which may account for such large and stable values of $`\dot{P}`$ is rotational loss by emission of magnetic radiation. For a dipolar magnetic field,

$$\dot{P}10^{15}B_{12}^2/P\frac{\mathrm{s}}{\mathrm{s}},$$
(4)
which gives a present value of the magnetic field of $`B=2.13\times 10^{13}`$ G 111Note that a twisted magnetosphere may lead to a reduction up to an order of magnitude in the inferred polar value of the magnetic field. The considerations about the age, however, remain unchanged (Thompson, Lyutikov & Kulkarni, 2001).222Here and in the following we specify the discussion to solution (1) of Table 3. Despite this value being below those typically quoted for magnetars ($`B10^{14}10^{15}`$ G), it is still extreme and close to the critical value $`B_Q=4.41\times 10^{13}`$ G at which quantizing effects start to be important in shaping the atmospheric emission. Furthermore, a scenario in which this source is powered by accretion from the interstellar medium must be ruled out: for the present values of $`P`$ and $`B`$ the corotating magnetosphere will prevent the incoming material to penetrate below the Alfven radius.

The corresponding spin-down age is

$$t_{sd}=\dot{P}/\left(2P\right)2.48\times 10^6\mathrm{yr},$$
(5)
which, given the numerous uncertainties, is marginally compatible with that inferred by the cooling curves: a few $`10^5`$ yrs for a surface temperature of $`80`$ eV (see e.g. Kaminker, Haensel & Yakovlev, 2001, Kaminker, Yakovlev & Gnedin, 2001; Schaab et al. 1997, 1999). The discrepancy is less significant if we notice that: a) the cooling curves are strongly dependent on the neutron star mass (Kaminker, Haensel & Yakovlev, 2001), and b) what we are probably observing in X-rays is a region of limited size which is kept hotter than the average star surface, as inferred by the analysis of the pulse-shape (Cropper et al., 2001).

On the other hand, $`t_{sd}`$ is representative of the true age of the source only in the case in which the magnetic field remained almost constant during the star evolution. The same condition applies for the validity of the cooling curves mentioned above, which do not include the extra input of energy released in the neutron star in case of field decay. It is therefore of fundamental importance to address the field evolution: the strength of $`B`$ and its variation during the neutron star history determine, in fact, the actual luminosity, the lifetime and even the nature of the energy loss from the star. A related question is where the field is anchored: in the core, penetrating the whole star or confined in the crust. There are three meachanisms which are typically proposed for inducing field-decay: ambipolar diffusion in the solenoidal or irrotational mode and Hall cascade (Goldreich & Reisenegger 1992, Heyl & Kulkarni 1998, Colpi, Geppert & Page 2000). In reality all the three processes co-exist with different timescales, and each of them may dominate the evolution depending on $`B`$, $`L`$ at any given time. In absence of more detailed computations, Heyl & Kulkarni (1998) and Colpi, Geppert & Page (2000) tentatively isolated the three processes and proposed some simple, phenomenological rules to mimic the evolution in the three regimes. We stress that these descriptions are far from being comprehensive of all the effects that can modify substantially the results: they are used in a first approximation and have the main advantage of having a simple analytical form. By using the expressions of Colpi et al. (2000), we have estimated the source age and the value of the magnetic field at the birth of the neutron star, $`B_0`$. Results are shown in Table 4. In all cases the star is assumed to be born with a period of 1 ms: results are not strongly dependent on this exact value, provided it is far less than the present period.

As we can see, allowing for a mechanism involving very fast decay, such as the Hall cascade, we find that the source is now $`4\times 10^4`$ yr old, and is born with a superstrong field $`B_010^{15}`$ G. Such a young age is only marginally compatible with the absence of a remnant and, more important, is not compatible with the low observed X-ray luminosity (as we compared using not only the standard cooling curves mentioned above, but also some cooling curves kindly provided by Geppert & Colpi, private communication; these latter are computed allowing for the extra-heating due to $`B`$-decay from Hall cascade and predict larger luminosities than the former, making the discrepancy even higher). Underluminous models have been recently presented by Kaminker, Haensel & Yakovlev (2001), who accounted for the enhanced neutrino cooling in presence of strong neutron superfluidity. These solutions may match an age of $`10^4`$ yrs for RXJ0720.4-3125, but, as discussed by the same authors, they must probably be rejected since they fail in the comparison with observational data of a sample of other neutron stars.

On the other hand, both mechanisms involving ambipolar diffusion predict a magnetic field quite stable over the source lifetime and close to the actual value. Accordingly, the predicted age is $`10^6`$ years in all cases, close to $`t_{sd}`$. Below $`10^{14}`$ G the cooling curves are not significantly influenced by decay through ambipolar diffusion (Heil & Kulkarni, 1998), thus, as in the case of constant $`B`$ discussed above, the scenario is compatible with the observed luminosity. The larger age is also compatible with the absence of a remnant.

If our conclusions are valid, the connection between dim INS and AXPs is not so obvious. RXJ0720.4-3125 has a strong, but not superstrong, field which is compatible with those of the canonical radio-pulsars which have passed the death line. On the other hand, having excluded accretion, what mechanism causes an X-ray emission concentrated in a fraction of $`60`$% of the star surface remains a mystery, as well as the related question about the validity of using the observed blackbody temperature to locate the source in the cooling diagram. The variation of the surface cooling temperature with the latitude predicted so far for strong fields is not large enough and more complicated explanations are required.

## 6 Acknowledgements

We are grateful to Mat Page for his advice on the maximum likelihood methods, and to Darragh O’Donoghue for pointing us to the Raleigh transform, and for the use of his Eagle Fourier transform code used in the preliminary stages of this work. We are grateful to Monica Colpi for lots of useful discussion, to Ulrich Geppert for providing the cooling curves and to Sandro Mereghetti for reading the manuscript.

## Appendix A Details of the Maximum Likelihood Periodogramme

The $`C`$-statistic is derived from the maximum likelihood ratio and is

$$C=2(E\underset{i=1}{\overset{N}{}}n_i\mathrm{ln}I_i)$$
(6)
(Cash 1979), where $`E`$ is the sum of the expected number of counts according to the model distribution, $`n_i`$ is the observed number of counts in an interval and $`I_i`$ is the expected number of counts at the time of event $`t_i`$ according to the model distribution. $`N`$ is the total number of counts. For event list data, the sampling duration tends to zero and $`n_i=1`$. The $`\mathrm{\Delta }C`$-statistic is

$$\mathrm{\Delta }C=(C_{min})_{pq}^T(C_{min})_p$$
(7)
where $`p`$ is the total number of parameters, and $`q`$ is the number of free parameters. This is distributed as $`\chi ^2`$ with $`q`$ degrees of freedom (Cash 1979). $`(C_{min})_p`$ is the minimum $`C`$-statistic when all parameters are varied – the minimum over all points in the search grid in the free parameter space. $`(C_{min})_{pq}^T`$ is the minimum $`C`$-statistic at any particular point in the search grid. Substituting,

$$\mathrm{\Delta }C=2\left[\left(\underset{i=1}{\overset{N}{}}n_i\mathrm{ln}I_i\right)_{best}\left(\underset{i=1}{\overset{N}{}}n_i\mathrm{ln}I_i\right)_{particular}\right].$$
(8)
The expected number of counts $`I_i`$ in equation (A1) is calculated from the model distribution. In our case of a pure sinusoidal variation, because the expected arrival of an event is directly proportional to the model prediction for time $`t_i`$, we have

$$I_i=a_0(1+A\mathrm{cos}(2\pi t_i/P+\theta _0)$$
(9)
where $`A`$ is the normalised amplitude, $`P`$ is the period and $`\theta _0`$ is the phase. The scaling of this model $`a_0`$ is unimportant, as constant factors are eliminated when substituting in equation (A3). This model has parameters $`(A,P,\theta _0)`$ so $`p=3`$. The periodogramme will scan in $`P`$, with best fit values for $`A`$ and $`\theta _0`$, so $`q=1`$. In the case where we allow a period change,

$$P=P_0+\dot{P}t_i$$
(10)
where $`P_0`$ is the period at $`t_i=0`$ and $`\dot{P}`$ is the period change, the periodogramme will consist of a scan in $`P_0`$ and $`\dot{P}`$. In this case $`p=4`$ and $`q=2`$.

In the general case of a maximum likelihood statistic, an optimisation search is required to obtain the best fit for the $`pq`$ fitted parameters. Here, in the case of a sinusoid, these can be derived directly from the Rayleigh transform, as pointed out by Bai (1992),

$$z=\frac{1}{N}\left[\left(\underset{i=1}{\overset{N}{}}n_i\mathrm{cos}2\pi t_i/P\right)^2+\left(\underset{i=1}{\overset{N}{}}n_i\mathrm{sin}2\pi t_i/P\right)^2\right]$$
(11)
where $`z`$ is the Rayleigh power, and again, $`n_i=1`$ for event list data. The amplitude is then

$$A=2\sqrt{z/N}$$
(12)
and the phase is

$$\theta _0=\mathrm{arctan}\left[\left(\underset{i=1}{\overset{N}{}}n_i\mathrm{sin}2\pi t_i/P\right)/\left(\underset{i=1}{\overset{N}{}}n_i\mathrm{cos}2\pi t_i/P\right)\right].$$
(13)
This is equivalent numerically to the discrete Fourier transform (DFT) (Deeming (1975), as can be ascertained by reference to Kurtz (1985), equation (A2).

The procedure is then for each $`P`$ (which may be calculated from equation (A5) in a $`P_0,\dot{P}`$ search) to calculate the amplitude and phase through equations (A7) and (A8) (which amounts to calculating the Rayleigh or discrete Fourier transform), then calculate $`\mathrm{\Delta }C`$ through equations (A4) and (A3). As the $`P`$ grid is searched, the maximum $`n_i\mathrm{ln}I_i`$ is stored: $`\mathrm{\Delta }C`$ is then simply twice the difference between the calculated value at each period, and this recorded maximum. This difference is distributed as $`\chi ^2`$ for one or two degrees of freedom, depending on whether $`P`$ is used directly, or calculated in a $`(P_0,\dot{P})`$ search through equation (A5).

An illustration of the power of the MLP is given in Figure 6. Here we have taken the X00a PN and X00b PN data, and computed a discrete Fourier transform in the narrow period range 8.390 to 8.392 s (a much expanded scale by comparison with Figure 2), assuming no $`\dot{P}`$. This (upper plot) shows a multiplicity of aliases created by the long gap in the data, superimposed on the broad peak corresponding to the duration of the longest observation (X00a). These are barely resolved in the plot. The lower plot is the MLP with the 68% and 90% confidence levels indicated. This eliminates most of the aliases, as the narrowness of the distribution is related to the precision with which the best fit period can be identified within the broad peak in the upper plot. In the case of the high signal-to-noise ratio X00a and X00b data, this is significantly narrower than the width of the peak, which is set by the window function (Scargle 1982).


format: text

text:
# Heavy quark action on the anisotropic lattice

## I Introduction

Recent experimental progress in flavor physics to look for the effect of new physics requires precise theoretical predictions from the standard model. However, model independent calculation of hadronic matrix elements is difficult because of nonperturbative nature of QCD. The lattice QCD simulation is one of most promising approaches in which the systematic uncertainties can be reduced systematically Lattice . The ultimate goal of this work is to construct a framework for lattice calculations of hadronic matrix elements in a few percent systematic accuracy, as required by the experiments in progress. More practically, the level of accuracy requested from the CLEO-c experiment CLEOC  as well as from the B Factories KEKB ; SLACB  are about 2 %. This paper intends to propose a systematic program to achieve this accuracy making use of the anisotropic lattice QCD.

In a lattice calculation of heavy quark systems such as the charm and bottom, one needs to control the large discretization error of $`O((am_Q)^n)`$. Extensive studies in various approaches have achieved much progress in understanding the heavy quark systems. However, it still seems difficult to achieve aforementioned systematic accuracy with the current techniques for several reasons. Here we categorize them into three types and summarize their advantages and disadvantages from the viewpoint of precision study of matrix elements in flavor physics.

1) Effective theories. One approach is to use descriptions of the heavy quark based on the heavy quark effective theory using nonrelativistic QCD TL91 . As an advantage, the action incorporates higher order correction in spatial momenta and can easily remove the all the mass dependent errors at the tree or one-loop level. However, since the theory does not have the continuum limit, one cannot remove the discretization errors by extrapolating the results on lattices with finite spacings. Another disadvantage is that the nonperturbative renormalization is difficult due to strong mass dependence. Precise measurements of the weak matrix elements for the heavy-light mesons below 10% are therefore difficult.

2) Relativistic framework. The most straightforward approach is to use the $`O(a)`$ improved Wilson action in natural way for relatively lighter heavy quark, and to extrapolate the results in $`1/m`$ according to the heavy quark effective theory. Since the theory has the continuum limit, the discretization errors can be removed by an extrapolation. The perturbative errors can also be avoided by employing the nonperturbative renormalization. In practice, however, lattice artifacts of $`O((am_Q)^2)`$ give dominant errors KS02  so that the precise calculation will remain difficult for the next few years even in quenched approximation. Brute force improvement by decreasing lattice spacing and employing large lattices quickly increase the simulation cost, and therefore not a realistic solution for the calculation of matrix elements.

3) Fermilab approach. The Fermilab approach links the above two approaches EKM97 ; Sro00 . In the heavy quark mass region, the $`O(a)`$ improved Wilson action with asymmetric parameters is reinterpreted as an effective theoretical description just like the nonrelativistic QCD action. Since the action reduces to the conventional $`O(a)`$ improved Wilson action for small masses $`am_Q1`$, it has in principle smooth continuum limit. The disadvantage is that it is not known how to extrapolate the results on lattices with $`am_Q>1`$, which are currently unavoidable in particular for systems with $`b`$-quark. For such a quark mass region, a mass dependent tuning of parameters in the action is required for proper improvement, although a systematic tuning beyond the perturbation theory is still a challenge AKT01 . Therefore precise calculation below 10% level are nontrivial for the weak matrix elements of the heavy-light mesons.

Therefore, it is desirable to develop a new framework for the heavy quark with the following features. (i) The continuum limit can be taken; (ii) A systematic improvement program, such as the nonperturbative renormalization technique NRimp , can be applied not only to the parameters in the action but also to the operators; (iii) A modest size of computational cost is required for systematic computation of matrix elements.

The anisotropic lattice, on which the temporal lattice spacing $`a_\tau `$ is finer than the spatial one $`a_\sigma `$, is a candidate of such a framework Kla98a ; Aniso01a . The approach is fundamentally along the line of Fermilab approach on the anisotropic lattice. However, large temporal lattice cutoff is expected to improve the above problems drastically. The most crucial one is that the mass dependence of the parameters in the action may become so mild in the region of practical interest that one can adopt the $`O(a)`$ improving clover coefficient determined in the nonperturbative renormalization technique. On the other hand, standard size of spatial lattice spacing keeps a total computational cost modest. Therefore, extrapolation of the result to the continuum limit may be possible with keeping systematic uncertainties under control with sufficient accuracy. Whether these observations practically hold should be examined numerically, as well as in the perturbation theory.

Our form of quark action on the anisotropic lattice is founded in Refs. Aniso01a ; Aniso01b , where the perturbative results for light and heavy quarks and the simulation results of the light quark were studied. In this paper, we focus on the heavy quark and study the mass dependence of the the Lorentz symmetry breaking effect in order to understand the mass range for which a consistent description is possible. This signals the applicability of anisotropic lattice to the heavy quark systems. We investigate the heavy-heavy and heavy-light mesons on a quenched lattice of anisotropy $`\xi =4`$ and spatial cutoff $`a_\sigma ^11.6`$ GeV. The heavy quark mass is varied from the charm quark mass to about 6 GeV to examine an applicability of the action to these systems.

This paper is organized as follows. In the next section, we first give our quark action, which is discussed in detail in Ref. Aniso01a . We then give our conjecture on the lattice spacing dependence of the anisotropy parameter toward the continuum limit, and discuss the advantage of the anisotropic lattice compared to the isotropic one. In Section III, we observe the tree-level expectation of the mass dependence of the $`O(a^2)`$ terms in the quark dispersion relation and study how the anisotropic parameters or the breaking effect of relativity behave as a functions of the heavy quark mass. Section IV describes the numerical simulation. We compute the heavy-light and heavy-heavy meson spectra and dispersion relations with two sets of heavy quark parameters. In one set (Set-I), the bare anisotropy is set to the value for the massless quark. The other set (Set-II) adopt the result of the mass dependent tuning using the heavy-heavy meson dispersion relation as obtained in Ref. Aniso01b . We observe the mass dependence of the renormalized anisotropy in order to probe the breaking of relativity for heavy quark. We also observe how the inconsistency among the binding energies of heavy-heavy, heavy-light and light-light mesons SAnomaly1 ; SAnomaly2  grows as a function of mass. The last part of this section discusses the hyperfine splitting of the heavy-light meson. In Section V, we summarize the result of simulation, and draw our perspective on further development of this framework, as our conclusion in this paper.

## II Formulation

### II.1 Quark action

We adopt the following quark action in the hopping parameter form Aniso01a ; Ume01 :

$$S_F=\underset{x,y}{}\overline{\psi }(x)K(x,y)\psi (y),$$
(1)
$`S_F`$ $`=`$ $`\delta _{x,y}\kappa _\tau [(1\gamma _4)U_4(x)\delta _{x+\widehat{4},y}`$ (2)
$`+(1+\gamma _4)U_4^{}(x\widehat{4})\delta _{x\widehat{4},y}]`$
$`\kappa _\sigma {\displaystyle \underset{i}{}}\left[(r\gamma _i)U_i(x)\delta _{x+\widehat{i},y}+(r+\gamma _i)U_i^{}(x\widehat{i})\delta _{x\widehat{i},y}\right]`$
$`\kappa _\sigma c_E{\displaystyle \underset{i}{}}\sigma _{4i}F_{4i}(x)\delta _{x,y}r\kappa _\sigma c_B{\displaystyle \underset{i>j}{}}\sigma _{ij}F_{ij}(x)\delta _{x,y},`$
where $`\kappa _\sigma `$ and $`\kappa _\tau `$ are the spatial and temporal hopping parameters, and are related to the bare quark mass and bare anisotropy as given below. The parameter $`r`$ is the spatial Wilson coefficient, and the parameters $`c_E`$ and $`c_B`$ are the clover coefficients for the $`O(a)`$ improvement. Although the explicit Lorentz symmetry is not manifest due to the anisotropy in lattice units, it can be restored in principle for physical observables in physical units at long distances up to errors of $`O(a^2)`$ by properly tuning $`\kappa _\sigma /\kappa _\tau `$, $`r`$, $`c_E`$ and $`c_B`$ for a given $`\kappa _\sigma `$. The action is constructed in accord with the Fermilab approach EKM97  and hence applicable to an arbitrary quark mass, although a mass dependent tuning of parameters is difficult beyond the perturbation theory. This may be circumvented by taking $`a_\tau ^1m_Q`$, with which the mass dependence of parameters are expected to be small so that $`O(a)`$ Symanzik improvement program for the heavy quark can be applied. To study whether this expectation really holds or not is the main subject of this paper.

In present study, we vary only two parameters $`\kappa _\sigma `$ and $`\kappa _\tau `$ with fixed other parameters. We set the Wilson parameter as $`r=1/\xi `$ and the clover coefficients as the tadpole-improved tree-level values, $`c_E=1/u_\sigma u_\tau ^2`$, and $`c_B=1/u_\sigma ^3`$. The tadpole improvement LM93  is achieved by rescaling the link variable as $`U_i(x)U_i(x)/u_\sigma `$ and $`U_4(x)U_4(x)/u_\tau `$, with the mean-field values of the spatial and temporal link variables, $`u_\sigma `$ and $`u_\tau `$, respectively. Instead of $`\kappa _\sigma `$ and $`\kappa _\tau `$, we introduce $`\kappa `$ and $`\gamma _F(=1/\zeta )`$ as

$`{\displaystyle \frac{1}{\kappa }}`$ $``$ $`{\displaystyle \frac{1}{\kappa _\sigma u_\sigma }}2(\gamma _F+3r4)=2(m_0\gamma _Fa+4),`$
$`\gamma _F`$ $``$ $`1/\zeta \kappa _\tau u_\tau /\kappa _\sigma u_\sigma .`$ (3)
The former controls the bare quark mass and the latter corresponds to the bare anisotropy.

### II.2 Mass dependence of anisotropy parameter

On an anisotropic lattice, one must tune the parameters so that the hadronic states satisfy the relativistic dispersion relations. In general, the lattice dispersion relation for arbitrary values of $`\zeta `$ is described in lattice units as

$$E^2(𝒑)=m^2+𝒑^2/\xi _F^2+O(p^4).$$
(4)
In the above expression the energy $`E`$ and the rest mass $`m`$ are in temporal lattice units while the momentum $`𝒑`$ is in spatial units. The parameter $`\xi _F`$ introduced in this equation characterizes the anisotropy of the quark fields. The difference between the quark field anisotropy $`\xi _F`$ and the gauge field anisotropy $`\xi `$ probes the breaking of relativity. Therefore, the calibration is nothing but the tuning of $`\zeta (=1/\gamma _F)`$ for a given $`\kappa `$ so that $`\xi _F`$ equals $`\xi `$ and hence the relativistic dispersion relations are satisfied. Let us call the tuned parameter as $`\zeta ^{}(=1/\gamma _F^{})`$. $`\zeta ^{}`$ depends on the quark mass and can in general be different from unity even for the case of isotropic lattice. The requirement of the relativity relation automatically enforces that the rest and kinetic masses are equal to each other. In this sense, our calibration procedure of anisotropic lattice action is a natural generalization of the Fermilab approach.

Now we consider the tuning of anisotropic parameter either on isotropic or anisotropic lattice, for fixed physical quark mass. The anisotropy parameter $`\zeta =1/\gamma _F`$ can be tuned using the dispersion relation of either heavy-heavy or heavy-light meson in mass dependent way. Alternatively, one can also adopt the value tuned for the light-light mesons neglecting the mass dependence. These procedures of calibration can give different results due to the discretization errors. We would like to know: 1) the value of lattice spacing $`a_\sigma ^{hh\text{-}hl}`$ above which the calibrated parameters $`\zeta ^{}`$ using the heavy-heavy and heavy-light mesons differ from each other by more than a certain accuracy $`ϵ_{acc}`$ (say, 2%), and 2) the value of lattice spacing $`a_\sigma ^{hl\text{-}ll}`$ above which the calibrated parameters $`\zeta ^{}`$ using the heavy-light and light-light mesons are different by more than a certain accuracy $`ϵ_{acc}^{}`$.

Naively speaking, the difference between $`\zeta ^{}`$’s with heavy-heavy and heavy-light mesons originates from the $`O((a_\sigma p)^2)`$ effects, where $`p`$ is the typical quark momenta in heavy-heavy meson, $`p\alpha m_Q`$. On the other hand, the deviation of $`\zeta ^{}`$ with heavy-light mesons from the $`\zeta ^{}`$ with light-light mesons is due to the $`O((a_\tau m_Q)^2)=O((a_\sigma m_Q/\xi )^2)`$ effects. With the criteria that these errors stay within the required accuracies, we obtain

$`a_\sigma ^{hh\text{-}hl}`$ $``$ $`\sqrt{ϵ_{acc}}/(\alpha m_Q),`$
$`a_\sigma ^{hl\text{-}ll}`$ $``$ $`\sqrt{ϵ_{acc}^{}}\xi /m_Q.`$ (5)
Figure 1 conjectures the expected behavior of the tuned anisotropy parameter $`\zeta ^{}`$ for the isotropic ($`\xi =1`$) and anisotropic ($`\xi =4`$) lattices. The tuned anisotropy normalized by the value at massless limit is displayed for fixed physical quark mass, around charm quark mass, as an example. In the case of isotropic lattice, as shown in the left panel of Fig. 1, $`a_\sigma ^{hl\text{-}ll}<a_\sigma ^{hh\text{-}hl}`$ is expected to hold. In this case, for a sufficiently small lattice spacing $`a_\sigma <a_\sigma ^{hl\text{-}ll}`$, both the heavy-heavy and heavy-light mesons can be successfully described by setting the anisotropy parameter to the value in the massless limit. For a coarser lattice spacing $`a_\sigma ^{hl\text{-}ll}<a_\sigma <a_\sigma ^{hh\text{-}hl}`$ , mass dependent tuning, namely the Fermilab approach, is necessary, but the description of both the heavy-heavy and heavy-light mesons still works. For an even coarser lattice spacing $`a_\sigma ^{hh\text{-}hl}<a_\sigma `$, simultaneously consistent description of heavy-heavy and heavy-light mesons no longer works, because of severe discretization effect in the heavy quarkonia.

Now let us consider the anisotropic lattice case with $`\xi =4`$. The expected mass dependence of $`\zeta ^{}`$ is schematically represented in the right panel of Fig. 1. As is obvious in Eq. (5) and will be studied in the following sections, the anisotropy does not improve the situation for heavy-heavy meson. Therefore $`a_\sigma ^{hh\text{-}hl}`$ remains roughly of the same size as in the isotropic lattice case. On the other hand, Eq. (5) and the studies in the following sections show that for $`\xi =4`$, $`a_\sigma ^{hl\text{-}ll}`$ is expected to increase and $`a_\sigma ^{hh\text{-}hl}<a_\sigma ^{hl\text{-}ll}`$ occasionally holds. This latter situation is particularly probable for quark with not very large mass, and expected to hold up to the bottom quark mass region. In this case, for sufficiently small lattice spacing $`a_\sigma <a_\sigma ^{hh\text{-}hl}`$, both the heavy-heavy and heavy-light mesons can be successfully described by setting the anisotropy parameter to the value in the massless limit. For a coarser lattice spacing $`a_\sigma ^{hh\text{-}hl}<a_\sigma <a_\sigma ^{hl\text{-}ll}`$ consistent description of heavy-heavy and heavy-light mesons no longer works. Nevertheless the heavy-light meson is successfully described with the anisotropy parameter tuned for the massless limit. Finally for an even coarser lattice spacing $`a_\sigma ^{hl\text{-}ll}<a_\sigma `$, the genuine Fermilab approach, namely mass dependent tuning of $`\zeta `$, is required to describe the heavy-light mesons.

In the range of lattice spacing $`a_\sigma <a_\sigma ^{hl\text{-}ll}`$, if we give up the description of heavy quarkonia and only consider the heavy-light systems, one can successfully adopt the tuned value of $`\zeta `$ at the massless limit. This is the most important region for our approach, since there the mass dependence of parameters in the action other than $`\zeta =1/\gamma _F`$, such as $`c_E`$ and $`c_B`$, are also expected to be small, and hence the result of tuning at the massless limit is also considered to be valid. For the tuning of clover coefficient for the massless quark, it may be possible to apply the nonperturbative renormalization technique NRimp , which is currently most efficient procedure to remove the $`O(a)`$ effect beyond perturbation theory. The most advantage of the anisotropic lattice is that this region $`a<a_\sigma ^{hl\text{-}ll}`$ is expected to be extended by a factor $`\xi `$ as is found in Eq. (5), so as to cover the heavy quark mass region of practical interests. Therefore the anisotropic lattice is a candidate of framework for heavy-light systems, which satisfies the required characters (i), (ii), and (iii) mentioned in the Introduction. Also for heavy-heavy systems, it would be important to study quantitatively the heavy quark mass at which the action starts to fail in correct description.

This paper is dedicated to a numerical study whether these expectations actually holds. Using Eq. (5), similar expectation can be deduced for the heavy quark mass dependence with a fixed lattice spacing. In this paper, since we work with only one lattice spacing, instead of studying the lattice spacing dependence with a fixed heavy quark mass, we study the mass dependence with a fixed lattice spacing. This means to specify in which region considered above the present $`a_\sigma `$ is located for each value of quark mass. This subject is discussed in more detail with the tree level analysis in the next section, and with the numerical simulation in the successive section.

## III Expectation from tree level analysis

In this section, we review the tree level analysis of the heavy quark action in the Fermilab formulation. In the following, we express the energy and the momenta in physical units so that the dependence on the lattice spacings $`a_\tau `$ and $`a_\sigma `$ is explicitly shown. The dispersion relation of the heavy quark on anisotropic lattice is given as EKM97 ; Aniso01a

$`\mathrm{cosh}(a_\tau E)`$ $`=`$ $`1+`$ (6)
$`{\displaystyle \frac{[\overline{m}_0+2r\zeta _i\mathrm{sin}^2(a_\sigma p_i/2)]^2+\zeta ^2_i\mathrm{sin}^2(a_\sigma p_i)}{2[1+\overline{m}_0+2r\zeta _i\mathrm{sin}^2(a_\sigma p_i/2)]}},`$
where $`\overline{m}_0=a_\tau m_0`$ . This leads to the following dispersion relation.

$`E^2`$ $`=`$ $`M_1^2+\left({\displaystyle \frac{\xi ^{\mathrm{tree}}}{\xi _F^{\mathrm{tree}}}}\right)^2𝒑^2`$ (7)
$`+A_1a_\sigma ^2(𝒑^2)^2+A_2a_\sigma ^2{\displaystyle \underset{i}{}}p_i^4+\mathrm{},`$
where $`M_1`$ is the pole mass,

$`M_1`$ $`=`$ $`\mathrm{log}(1+\overline{m}_0)a_\tau ^1,`$ (8)
and the anisotropy of the quark field at the tree level,

$$\left(\frac{\xi ^{\mathrm{tree}}}{\xi _F^{\mathrm{tree}}}\right)^2=\mathrm{log}(1+\overline{m}_0)\xi ^2\left[\frac{2\zeta ^2}{\overline{m}_0(2+\overline{m}_0)}+\frac{r\zeta }{1+\overline{m}_0}\right].$$
(9)
The third and fourth terms in Eq (7) represent the $`O((ap)^2)`$ errors, where $`A_1`$ and $`A_2`$ are given as

$`A_1`$ $`=`$ $`{\displaystyle \frac{\xi ^2}{4}}[({\displaystyle \frac{2\zeta ^2}{\overline{m}_0(2+\overline{m}_0)}}+{\displaystyle \frac{r\zeta }{1+\overline{m}_0}})^2`$ (10)
$`\mathrm{log}(1+\overline{m}_0)({\displaystyle \frac{8\zeta ^4}{\overline{m}_0^3(2+\overline{m}_0)^3}}`$
$`+{\displaystyle \frac{4\zeta ^3[\zeta +2r(1+\overline{m}_0)]}{\overline{m}_0^2(2+\overline{m}_0)^2}}+{\displaystyle \frac{r^2\zeta ^2}{(1+\overline{m}_0)^2}})],`$
$$A_2=\frac{\xi ^2}{3}\mathrm{log}(1+\overline{m}_0)\left[\frac{2\zeta ^2}{\overline{m}_0(2+\overline{m}_0)}+\frac{r\zeta }{4(1+\overline{m}_0)}\right].$$
(11)
The above coefficients are derived using the $`M_1`$, $`M_2`$, $`M_4`$ and $`w_4`$ in Ref. EKM97  extended to anisotropic lattice by replacing certain parameters as given in Ref. Aniso01a .

Fig. 2 shows the mass dependence of $`\zeta `$ tuned by requiring $`\xi _F=\xi `$, both on the isotropic lattice and on the anisotropic lattice with $`\xi =4`$. In the latter case, it is clear that the mass dependence is drastically reduced so that taking the value of $`\zeta `$ with massless limit is a good approximation over a wide range of quark mass, in contrast to the case of the isotropic lattice.

In the following, let us consider either the case with mass dependent tuning of $`\zeta `$ (denoted by “Fermilab” in the figures) or with the $`\zeta `$ at the massless limit (“massless”). Figure 3 shows the breaking effect of relativity relation, $`\xi _F/\xi `$, according to Eq. (9). The kinetic mass $`M_2`$ of quark is related to $`\xi _F/\xi `$ as $`M_2=M_1\xi _F^2/\xi ^2`$. For the mass dependently tuned $`\zeta `$, $`\xi _F`$ trivially equals $`\xi `$. Without tuning, the mass dependence of $`\xi _F`$ is drastically reduced in the case of $`\xi =4`$, as is expected from that the deviation of $`\xi _F/\xi `$ from unity is the $`O(a_\tau m_Q)`$ error. In Figure 4, one find that the coefficient $`A_1`$, whose deviation from zero signals the $`O(a^2)`$ effect, is also drastically reduced on the anisotropic lattice, either with or without mass dependent tuning of $`\zeta `$. On the other hand, $`O(a^2)`$ error from $`A_2`$ is not reduced but somewhat worse on the anisotropic lattice as is seen in Figure 5. However, the coefficient $`A_2`$ is not larger than the value at the massless limit. As a general feature for all of the $`\xi _F/\xi `$, $`A_1`$ and $`A_2`$, the mass dependences are drastically reduced on the anisotropic lattice. On the other hand, mass dependent tuning of $`\zeta `$ on isotropic lattice does not reduce the coefficients $`A_1`$ and $`A_2`$, while it completely removes the discrepancy between $`\xi _F`$ and $`\xi `$. In order to reduce these errors one has to introduce higher order terms as is pointed out in Ref. EKM97 .

From these results, one can expect the following about the error of the heavy quark systems. For the heavy-light systems, since heavy quark momenta typically takes the size of the hadronic scale, $`\mathrm{\Lambda }_{QCD}`$ 200–500 MeV, the $`O((ap)^2)`$ errors are well under control even for large $`a_\sigma M_2`$, if one keep $`a_\sigma ^1\mathrm{\Lambda }_{QCD}`$. While the largest error originates from $`A_2`$, this is of the same order as that in the light quark systems. For the heavy-heavy system, the situation is quite different. Since the heavy quark momenta are typically $`p\alpha m_Q`$, The $`O((ap)^2)`$ errors are expected to be as large as $`A_1(\alpha a_\sigma m_Q)^2`$, $`A_2(\alpha a_\sigma m_Q)^2`$. Again the $`A_2`$ gives the largest contribution, but the size of the error is uncontrolled when $`a_\sigma m_Q`$ is large.

To summarize, the anisotropic lattice largely reduces the discretization effects represented by $`\xi _F/\xi 1`$ and $`A_1`$, while it does not improve the $`A_2`$. For the heavy-light systems, this suffices for a computation with discretization effects under control. On the other hand, when $`a_\sigma M_2`$ is of order of unity, the anisotropic lattice does not improve the situation for heavy quarkonia, because of severe effect of $`A_2`$ in these systems. Although mass dependent tuning of $`\zeta `$ further removes the deviation of $`\xi _F`$ from $`\xi `$, the $`\zeta `$ tuned for massless quark also stay a good approximation if $`a_\tau m_01`$ holds. This is the most advantage of anisotropic lattice compared with the isotropic Fermilab approach, in which the mass dependence of the parameter is much stronger. These observations are in accord with our conjecture in the last section for developing the anisotropic Fermilab formulation.

## IV Numerical simulation

This section numerically examines the idea described in previous sections. For this purpose, we perform simulations with two series of heavy quark parameters. For Set-I the anisotropy parameter is set to the value at the chiral limit, $`\gamma _F^{}(m_q=0)`$, just as same as for the light quark. The Set-II adopts the fully tuned anisotropy parameter $`\gamma _F^{}`$ using heavy-heavy mesons. This calibration is done in the first subsection. For the heavy-heavy and heavy-light mesons, the rest and kinetic masses are obtained with two sets of parameters. Two quantities are used to probe the breaking effect of the relativity: the fermionic anisotropy $`\xi _F`$ for heavy-heavy and heavy-light mesons, and the inconsistency among the binding energies of heavy-heavy, heavy-light, and light-light mesons. These behaviors give us a hint to judge the regions in which our framework can be applied. The last subsection treats the heavy-light meson spectrum, and observe how the hyperfine splitting behaves with meson kinetic mass.

### IV.1 Calibration in heavy quark region

The simulation parameters used in this paper is the second set in Ref. Aniso01b : the quenched anisotropic lattice of the size $`16^3\times 128`$ generated with the standard plaquette gauge action with $`(\beta ,\gamma _G)=(5.95,3.1586)`$, which correspond to the renormalized anisotropy $`\xi =4`$ Kla98b  and spatial lattice cutoff $`a_\sigma ^1=1.623(9)`$ GeV set by hadronic radius $`r_0`$ Som94 . The mean-field values in the quark action are set to the mean values of link variables in the Landau gauge: $`u_\sigma =0.7917`$ and $`u_\tau =0.9891`$.

In Ref. Aniso01b , the optimum bare anisotropy $`\gamma _F^{}`$ is determined using the dispersion relation of mesons with degenerate quark masses, and the resultant values of $`\gamma _F^{}`$ are well represented by a linear function

$$\frac{1}{\gamma _F^{}}=\zeta _0+\zeta _2m_q^2,m_q=\frac{1}{2\xi }\left(\frac{1}{\kappa }\frac{1}{\kappa _c}\right),$$
(12)
where for the present lattice $`\zeta _0=0.2490(8)`$, $`\zeta _2=0.189(15)`$, and $`\kappa _c=0.12592(6)`$.

In this paper, we use seven values of $`\kappa `$ for heavy quark ($`\kappa _h`$), which roughly cover 1–6 GeV. Three of them already appeared in Ref. Aniso01b . We start with the calibration of remaining four values of $`\kappa `$ in the heavy quark region in the same manner as in Ref. Aniso01b . The values of $`\kappa `$ used are listed in Table 1 together with the result of calibration. The second column is the naive estimate of bare quark mass according to Eq. (12). For the heaviest case, $`m_0`$ is almost unity in temporal lattice units, and therefore the breaking effect of relativity may be visible. Here we note that for heavier quark masses the difference of $`\gamma _F^{}`$ for pseudoscalar and vector mesons,

$$\delta \gamma _F^{}=\gamma _F^{(V)}\gamma _F^{(PS)},$$
(13)
is sizable beyond the statistical fluctuations. This signals that the quarkonium system is not properly described within the present framework at this lattice spacing. This problem will again be discussed later in terms of the fermionic anisotropies for heavy-heavy and heavy-light mesons.

Figure 6 shows the result of calibration. The result is well fitted to a linear form in $`m_q^2`$ or a quadratic form in $`m_q`$, in spite of large quark mass. The fits including previous data in Ref. Aniso01b  result in

Linear: $`\zeta _0=0.2510(6),\zeta _2=0.1437(26).`$ (14)
Quadratic: $`\zeta _0=0.2473(12),\zeta _1=0.0267(72),`$ (15)
$`\zeta _2=0.1151(81).`$
In both cases, the values of $`\zeta _0`$ is close to the value at the mean-field tree level, $`\xi ^1=0.25`$. Since there is no reason that the small quark mass dependence persist up to $`m_qa_\tau ^1`$, the difference between $`\zeta _2`$’s in linear fits with present data and with only previous data in lighter quark mass region is irrelevant, and above result just shows that the quark mass dependence of $`1/\gamma _F^{}`$ is still tractable at this mass region.

### IV.2 Relativity breaking effect

In this subsection we compute the heavy-heavy and heavy-light meson spectra and the dispersion relations using two sets of heavy quark parameters. For the first one, Set-I, the bare anisotropy is set to the value of $`\gamma _F^{}`$ in the chiral limit obtained in Ref. Aniso01b , namely $`\gamma _F=4.016`$ at $`\beta =5.95`$, for all quark masses. The second one, Set-II, adopts the above result of mass dependent calibration. We use the same hopping parameter values $`\kappa _h`$ for the heavy quark as given in the previous subsection. For the light quark, we use single value, $`\kappa _l=0.1235`$. The value of $`\gamma _F`$ at this $`\kappa _l`$ is set to the value in the chiral limit, as same as in the light hadron spectroscopy in Ref. Aniso01b . We regard that at $`\kappa =\kappa _l`$ the quark mass is sufficiently light for present purposes, and do not extrapolate to the chiral limit. The numbers of configurations used are 200 and 500 for heavy-heavy and heavy-light meson masses, respectively.

The meson spectra are listed in Table 2 and 3 for heavy-heavy and heavy-light mesons, respectively. In these tables, we also list the results of $`\xi _F`$ for each meson channel, and difference between them, $`\delta \xi _F=\xi _F^{(V)}\xi _F^{(PS)}`$. If the anisotropic lattice action does not describe the quarks inside mesons in respecting the relativity relation, the breaking effect appears in the dispersion relations of mesons. Therefore the deviation of fermionic anisotropy $`\xi _F`$ from $`\xi `$ signals the breaking effect of relativity. In the following, we first discuss on the result of $`\xi _F`$ for Set-I parameters, namely with $`\gamma _F`$ tuned for massless limit, and then briefly summarize the result for Set-II parameters.

Figure 7 displays the heavy quark mass dependence of $`\xi _F`$ for Set-I. The horizontal axis is the bare quark mass $`m_q`$ in temporal lattice units. The behaviors of $`\xi _F`$’s are well in accord with the expectation in Sec. II. For quantitative discussion, let us consider the case that the required accuracies to define the $`a_\sigma ^{hh\text{-}hl}`$ and $`a_\sigma ^{hl\text{-}ll}`$ are 2 %, namely $`ϵ_{acc}=ϵ_{acc}=0.02`$. $`\xi _F`$’s from the heavy-heavy and heavy-light mesons disagree beyond this accuracy at $`m_q>0.2`$ ($`a_\sigma m_Q>0.8`$). Therefore one must keep $`m_q<0.2`$ to avoid large systematic uncertainty in the heavy quarkonia. On the other hand, the $`\xi _F`$ from the heavy-light mesons are rather close to $`\xi `$, and $`\gamma _F^{}(m_q=0)`$ can be applied up to $`m_q0.3`$ within presently required accuracy. For quark mass larger than this value, the discrepancy between $`\xi _F`$’s from the pseudoscalar and vector mesons gradually grows beyond the statistical error. This signals the growth of systematic error. In the region of $`m_q<0.3`$ such effect is sufficiently small.

For the charm quark mass, the present lattice spacing $`a_\sigma `$ is already well less than $`a_\sigma ^{hl\text{-}ll}`$ so that the $`\gamma _F`$ tuned for massless quark is applicable to the charmed hadron systems. On the other hand, bottom quark mass is not the case, and one need finer lattice spacing or larger anisotropy $`\xi `$. The present lattice would be also sufficient for the charmonium states, since the region $`m_q<0.2`$ also covers the charm quark mass. Another striking feature is that the observed $`\xi _F`$’s for heavy-light mesons are close to the tree level expectation. This implies that the deviation of $`\xi _F`$ from $`\xi `$ may be largely removed by the tree level tuning of $`\gamma _F`$. Such an approach would work for the spectroscopy of hadrons containing single bottom quark. This is a good alternative procedure to the mass dependent calibration using heavy-light meson, since the statistical error of $`\xi _F`$ from heavy-light mesons rapidly grows as heavy quark mass.

Now we summarize the result for the Set-II. The heavy-heavy meson satisfies the relativity relation by definition. However, the heavy-light meson dispersion relation violates the relativity relation so that $`\xi _F/\xi `$ deviates from unity toward smaller value as heavy quark mass increases. This implies that the mass dependent tuning cannot absorb the discrepancy between the $`\xi _F`$’s from the heavy-heavy and heavy-light meson systems. For the quark mass of $`m_q<0.2`$, the result for Set-II are consistent with the result for Set-I. In that region the heavy-heavy meson can be successfully described, and therefore the result of calibration with heavy-heavy mesons is also valid.

We now observe the inconsistency among the binding energies of heavy-heavy, heavy-light, and light-light mesons discussed in Refs. SAnomaly1 ; SAnomaly2 . The inconsistency is probed by a quantity

$$I\frac{2\delta M_{hl}(\delta M_{hh}+\delta M_{ll})}{2M_{2hl}},$$
(16)
where $`\delta M=M_2M_1`$, $`M_1`$ and $`M_2`$ are the rest and kinetic masses of the meson, respectively. The subscripts $`hh`$, $`hl`$, and $`ll`$ represent the quark contents of each meson ($`h`$ for heavy and $`l`$ light quarks). We neglect the last term in the numerator, since the calibration of light quark mass region requires that $`\delta M_{ll}`$ vanishes. Since the rest and kinetic quark masses cancel in each kind of mass, nonvanishing $`I`$ represents the inconsistency in the binding energy, namely the dynamical effect. The anomalous behavior of $`I`$ in large kinetic mass region was first reported in Ref. SAnomaly1  for the $`O(a)`$ improved quark action on isotropic lattice. Ref. SAnomaly2  explained that this behavior originates from the $`O((ap)^2)`$ discretization effect in the heavy quarkonium system and estimate the size of $`I`$ with the help of potential model analysis. It results in $`I0.5`$ at $`aM_{2hl}3.2`$, which is in good agreement with the result in SAnomaly1 .

Figure 8 displays the results of $`I`$ for Set-I and II for the pseudoscalar channel. The behavior of $`I`$ is quite similar to that in Refs. SAnomaly1 ; SAnomaly2 , as is expected from that it originates from $`O((a_\sigma p)^2)`$ error and therefore is not improved by the anisotropy. The behavior of $`I`$ of two sets is very similar to each other. It means that the inconsistency cannot be eliminated by tuning anisotropy parameters with either heavy-heavy or heavy-light system, and a kind of universal quantity to signal the failure in consistent description of heavy quarkonium systems. The $`I`$ rapidly deviate from zero for the heavy quark mass larger than our lightest one. This is consistent with the above observation of $`\xi _F`$.

These results indicate that for $`a_\tau m_Q>0.2`$ the present action fails to describe the dynamics of heavy quarkonium, and cannot be reinterpreted by the subtraction of quark mass effect. Thus we need to either give up to apply the present framework for such a quark mass region, or to improve the action by incorporating higher order correction terms. The remaining part of this section therefore focus on the heavy-light meson spectrum.

### IV.3 Heavy-light meson spectrum

Now we turn our attention to the heavy-light meson, which is of our main interest. According to the heavy quark expansion, the spin flipping interaction of heavy quark in the heavy-light systems is of $`O(1/m_Q)`$, and hence the hyperfine splitting of mesons is proportional to the meson mass at the leading order. In the heavy-light systems, the large mass of heavy quark is not important for the dynamics. Employing the Fermilab formulation, this character is taken into account and the the correct heavy quark expansion is in terms of kinetic mass EKM97 . Therefore, the hyperfine splitting, which is measured as the difference of rest masses, are expected to be reciprocally proportional to the kinetic meson mass. If this is the case, Set-I and Set-II should show the similar behaviors, up to the $`O(a^2)`$ systematic uncertainty.

The hyperfine splitting $`m_Vm_{PS}`$ is displayed in Figure 9. The horizontal axis is spin averaged kinetic meson mass inverse in physical units defined through the hadronic radius $`r_0`$. Since the most serious uncertainty in the kinetic mass is from the systematic uncertainty in $`\xi _F`$, we estimate this error as

$$\delta M_2=2M_2|\delta \xi _F|/\xi _F,$$
(17)
where $`\delta \xi _F=\xi _F^{(V)}\xi _F^{(PS)}`$, and hence it does not include the statistical errors of $`\xi _F`$ and $`M_2`$. On the other hand, the error associated to the hyperfine splitting is the statistical one. Although the data shows a linear dependence in the heavy quark mass region, there appears small negative intercept. We regard this small discrepancy with the heavy quark expansion as the Lorentz symmetry breaking effect brought in the meson dispersion relation and of $`O(a^2)`$. For fixed physical quark mass, this effect is expected to disappear linearly in $`a^2`$ toward the continuum limit. The results of Set-I and Set-II clearly show similar behavior, and therefore the above interpretation of Fermilab formulation works well up to the relativity breaking effect represented by $`\delta \xi _F`$. For more quantitative analysis, one must quantify the size of $`O(a^2)`$ systematic errors and observe how they disappear toward the continuum limit. This is out of the scope of this paper.

## V Conclusion

In this paper, we investigated an applicability of the anisotropic lattice quark action in heavy quark mass region, on the quenched lattice with $`a_\sigma ^11.6`$ GeV and renormalized anisotropy $`\xi =4`$. The anisotropy is expected to extend the region in which the parameters in the action tuned for massless quark are applicable to precision computations of heavy-light matrix elements. For this purpose, we measured the heavy-heavy and heavy-light meson masses and dispersion relations to observe the relativity breaking effect. The calculation was carried out for two sets of parameters, Set-I and Set-II. Set-I adopts the values at the massless limit, while in Set-II the bare anisotropy is tuned using the heavy-heavy mesons. Our main results are as follows.

(a) In the quark mass region $`a_\tau m_Q<0.2`$, the observed fermionic anisotropy $`\xi _F`$’s are consistent for heavy-heavy and heavy-light mesons as well as for Set-I and Set-II parameters within 2 % accuracy. This implies that present framework is applicable to both these systems even without tuning of anisotropy parameter. Beyond this region, the action fails to describe the heavy quarkonium state correctly, as expected from the tree level analysis of the quark dispersion relation.

(b) The mass dependence of the renormalized anisotropy for the heavy-light mesons with $`\gamma _F^{}`$ at the massless limit (Set-I) is so small that one can exploit the massless tuning in the region of $`a_\tau m_Q<0.3`$ with less than 2% errors. This result is essentially important for our strategy, since it implies that the parameters tuned at the massless limit is directly available for this mass region, which already contains charm quark mass with the present lattice and can be extended to the bottom quark with development of computational resources in successive decade.

(c) For $`a_\tau m_Q>0.3`$, relativity breaking effect in the heavy-light mesons seems to grow as a function of the heavy quark mass. This is signaled by the discrepancy of $`\xi _F`$’s from the pseudoscalar and vector mesons, although the present statistics is not enough for a quantitative estimate of this effect. In the scaling of hyperfine splitting, we also found small discrepancy with the expectation from the heavy quark expansion, which is considered as $`O(a^2)`$ systematic effect. It is important to quantify these effects and to observe how they disappear toward the continuum limit, for future high precision computations with this approach.

We therefore conclude that the anisotropic lattice quark action satisfies the required characters (i)–(iii) for high precision calculations of heavy-light matrix elements mentioned in Introduction. Among above results, the small mass dependence of anisotropy parameter is in particular quite encouraging for further development of the framework in this direction.

One of the promising strategies is to calibrate the parameters in the action at the massless limit and use them for all masses. This should work perfectly for $`a_\tau m_Q<0.3`$ which is already sufficient to describe the $`D`$ meson systems. If one wants to treat $`B`$ meson systems, one has to work with heavier quark, $`a_\tau m_Q>0.3`$, in which case the mass dependent errors cannot be neglected. However, since the mass dependence is small, it can be interpreted as an $`O(a^2)`$ error and can be removed by taking the continuum limit. Alternatively, one can also apply the genuine Fermilab approach for the bottom quark. In this case, the mass dependences of the renormalization coefficients are the source of systematic errors. Nevertheless, as long as one obtains such coefficients nonperturbatively in the massless limit first, and use the one-loop perturbation theory only to compute the mass dependent corrections, the perturbative error can be much better controlled compared to the Fermilab approach on the isotropic lattice. The nice agreement between the observed $`\xi _F`$ and the tree level expectation suggest that this idea is also promising in the bottom quark mass region.

On the other hand, present level of systematic errors as well as the statistical one is not sufficient for these ends. As a most significant example, one need to eliminate the $`O(\alpha a)`$ error using nonperturbative renormalization technique in the massless limit Kla98a , and confirm that this suffices also for heavy quark mass region treated in this paper.

## Acknowledgments

We thank T. Umeda for useful discussions. The simulation was done on NEC SX-5 at Research Center for Nuclear Physics, Osaka University and Hitachi SR8000 at KEK (High Energy Accelerator Research Organization). H. M. is supported by Japan Society for the Promotion of Science for Young Scientists. T. O. is supported by the Grant-in-Aid of the Ministry of Education No. 12640279. A. S. is supported by the center-of-excellence (COE) program at Research Center for Nuclear Physics, Osaka University.