text:
# The asteroseismological potential of the pulsating DB white dwarf stars CBS 114 and PG 1456+103

## 1 Introduction

The helium-atmosphere DB white dwarf stars are very interesting from the standpoint of stellar structure and evolution. The chemical evolution of their atmospheres cannot be satisfactorily explained to date (see Shipman 1997 for a review). In particular, the presence of the so-called DB gap, the absence of DB white dwarfs between temperatures of $``$30000 to 45000 K (Liebert 1986), is poorly understood. It is also not clear whether DBs are mostly produced by single-star evolution or whether a significant fraction originate from binary progenitors.

The examination of these problems can be aided by asteroseismology – the study of the interiors of pulsating stars via the analysis of their normal mode spectra. Fortunately, a class of pulsating DB white dwarf stars (hereinafter DBVs) exists, and their prototype, GD 358, is one of the classical examples for the successful application of asteroseismological methods (Winget et al. 1994, Vuille et al. 2000).

Compared to other classes of pulsating star, the DBVs seem quite promising candidates for asteroseismology: their mode spectra are shown or believed to be rich, and pulsation theory for these stars is quite advanced and well-tested. For instance, the central oxygen abundance of the models has been shown to have a measurable effect on the pulsation frequencies (Metcalfe, Winget & Charbonneau 2001), and the possible presence of a 3He/4He transition zone due to chemical diffusion has also been shown to produce a measurable effect (Montgomery, Metcalfe & Winget 2001), although recent observations have ruled out this possibility (Wolff et al. 2002). In addition, the recently developed application of genetic-algorithm-based model-fitting to DBV mode spectra (Metcalfe, Nather & Winget 2000) allows an objective and effective exploration of parameter space to find an optimal model for the observations.

The application of these methods to DBVs other than GD 358 is however hampered by a lack of suitable observational data. Most of these stars are rather faint (around 16th magnitude), which requires the use of 2-m class telescopes to acquire suitable data with photomultiplier detectors. However, time on these telescopes is generally oversubscribed and long runs are difficult to obtain. With the advent of CCD detectors and with the development of systems capable of acquiring measurements with high time resolution, the DBVs are now in reach of smaller telescopes. We have therefore started to acquire measurements aiming at making as many DBVs as possible accessible to asteroseismology.

The pulsations of the DB white dwarf star CBS 114 were discovered by Winget & Claver (1988, 1989), who reported multiperiodic oscillations with peak-to-peak light variations up to 0.3 mag and with a time scale of about 650 s. PG 1456+103 (hereinafter referred to as PG 1456) was reported as a DBV by Grauer et al. (1988). These authors reported multiperiodic pulsations with dominant periodicities between 420 and 850 seconds and maximum peak-to-peak amplitudes around 0.15 mag. As no further observations of either of these stars have been published, we deemed it worthwhile to obtain larger data sets to study the pulsations of CBS 114 and PG 1456 in more detail.

## 2 Observations

Most of our measurements were acquired as differential CCD photometry with the 0.75-m telescope at the Sutherland station of the South African Astronomical Observatory (SAAO), at which we had been allocated three weeks of observing time.

We used the University of Cape Town CCD camera (O’Donoghue 1995). Its CCD was operated in full-frame mode to acquire a maximum number of local comparison stars. We used 2$`\times `$2 or 3$`\times `$3 prebinning depending on seeing, which results in readout times of about 2–3 seconds. For the faint CBS 114 ($`B17.2`$ as estimated from our CCD frames compared to PG 1456+103), we used integration times of 27 or 28 seconds to obtain one frame per 30 seconds, and we used 20-second integrations for PG 1456. No filter was used to maximize the number of photons detected. On every clear night, sky flatfields were taken during twilight.

In addition, we acquired measurements of PG 1456 with the 0.9-m SARA telescope at Kitt Peak National Observatory (Arizona), to reduce the aliasing problem in the frequency analysis. The CCD used is an Apogee AP7p, with a back-illuminated SITe SIA-502AB 512$`\times `$512 detector. The pixels are 24$`\mu `$m square, or 0.75 arcsec/pixel. Read noise for the camera is 12.2 electrons r.m.s., and the gain is 6.1 electrons per ADU. We used integration times of typically 30 seconds, and our readout time was about 7 seconds. Again, we took sky flats during each twilight. Our measurements are summarized in Table 1.

We started data reduction with the corrections for bias, dark counts (SARA measurements only) and flatfield. Mean weekly flatfields were used for the SAAO data sets and combined nightly flats for the SARA frames. Photometric measurements on these reduced frames were made with the program MOMF (Kjeldsen & Frandsen 1992), which applies a combination of PSF and aperture photometry relative to a user-specified ensemble of comparison stars. No variability of any star other than the targets in the different CCD fields was found, and the comparison star ensemble resulting in the lowest scatter in the target star light curves was chosen.

The resulting differential light curves were corrected for differential colour extinction and for correlations with seeing. For CBS 114, it was sometimes also found necessary to de-trend the data with (x, y) position of the star on the chip because of flatfield errors. Residual low-frequency trends in the data (not found to be coherent over the time span of the observations, hence judged not to be intrinsic to the stars) were removed by means of low-order polynomials.

Finally, the times of measurement were transformed to a homogeneous time base. We chose Terrestrial Time (TT) as our reference for measurements on the Earth’s surface and applied a correction to account for the Earth’s motion around the solar system’s barycentre. As this barycentric correction varied up to $`\pm `$1 s throughout a run, we applied it point by point. Our final time base therefore is Barycentric Julian Ephemeris Date (BJED). The reduced time series were subjected to frequency analyses; we show example light curves in Figs. 1 and 2. The typical rms scatter per single data point is around 40 – 50 mmag for CBS 114, and for PG 1456 about 15 mmag (SAAO data) and 11 mmag (SARA data).

## 3 Frequency analysis

Our frequency analysis was performed with the programme PERIOD 98 (Sperl 1998). This package applies single-frequency power spectrum analysis and simultaneous multi-frequency sine-wave fitting, but also has some advanced options. In particular, it can be used to calculate optimal solutions for multiperiodic signals including harmonic, combination, and equally spaced frequencies, which are often found in the analysis of the light curves of pulsating white dwarf stars.

### 3.1 CBS 114

A frequency analysis of single-site data of a multiperiodic variable in the presence of noise is difficult because of aliasing; extreme caution is required. However, the structure and extent of our data set aids us in this effort.

We started by calculating nightly amplitude spectra of our data. The amplitude spectra of the individual nights were quite similar, with a number of dominant peaks always occurring at similar frequencies. However, some variations in the individual amplitudes were noted.

As the next step, we combined the three weekly data sets (which are very similar in length and extent) and we computed their amplitude spectra as well as that of all three weeks combined. They are shown in Fig. 3, which shows six dominant structures between 1.4 and 2.6 mHz and some further signals which could be combination frequencies.

Some trial prewhitening within the strongest features convinced us that all of them are dominated by one signal each. Consequently, we determined their frequencies by examining the combined, the weekly and the average of the weekly amplitude spectra and by least-squares fitting to the light curves. In this way, we could rule out a number of alias frequencies because they gave unreasonable results if adopted. However, in some cases, no unambiguous decision could be made.

Having determined the best set of frequencies, we refined their values by calculating an optimal solution for the whole data set, prewhitened it from the data, and calculated the amplitude spectrum of the residuals (upper panel of Fig. 4). Because of the temporal changes in the weekly amplitude spectra, it is not surprising that some residual mounds of amplitude remained near the dominant frequencies.

We examined the frequencies of these features and their aliases and performed trial prewhitening to check for multiplet structure possibly present. No such evidence was found. We then fit our best frequencies to the individual nights longer than 3 h (assuring the dominant structures are resolved) and determined their amplitudes, phases, and corresponding error estimates (following Montgomery & O’Donoghue 1999). We found that significant amplitude and phase changes occurred over time scales of 3 – 4 days, and that either amplitude or phase variability alone is not sufficient to explain the observations; both are present.

As we cannot determine the cause of these variations, we treated them phenomenologically and calculated our final multifrequency solution with constant frequencies as determined before, but with amplitudes and phases variable over half-week periods. In this way, we can determine ranges of the amplitudes of the dominant signals, and most of the residual mounds in the amplitude spectrum are thus removed (lower panel of Fig. 4). In Table 2, we list the best multifrequency solution determined with this method.

#### 3.1.1 Re-analysis of Winget & Claver’s data

Although the measurements by Winget & Claver (1988, 1989) only comprise 5.2 h of data obtained during two nights, we can still use them for a comparison with our results. The amplitude spectrum of these measurements (upper panel of Fig. 5) shows dominant features at frequencies agreeing very well with the ones determined from our data. Consequently, we fitted the first six frequencies from Table 2 to these measurements and noticed that they were consistent with the older data. We therefore determined corresponding frequencies, amplitudes and phases, subtracted this fit from the data and calculated a residual amplitude spectrum, which is shown in the lower panel of Fig. 5, revealing the presence of yet another signal.

A variation at this frequency is also present in our new observations, although its amplitude is below our detection threshold. Consequently, we include it in our final frequency solution for the combined data sets (Table 2). We note some more suspected signals in the lower panel of Fig. 5, but due to the small amount of data we do not push the analysis further. We have therefore detected altogether seven independent pulsation frequencies and two combination signals in the available light curves of CBS 114.

### 3.2 PG 1456

Due to our two-site coverage of the light variations of PG 1456, the frequency analysis should be comparatively easy. However, initial trials showed that this star’s pulsational amplitudes and frequencies are also somewhat variable, and as our light curves from the first two weeks are rather short, they complicate the analysis if the total data set is used. Consequently, we used the JD 2451987 - JD 2452002 two-site measurements as our primary data set for the frequency analysis; the remaining data were used for comparison purposes and consistency checks. Fig. 6 shows some prewhitening steps in these two main subsets of data.

A signal ($`f_B`$) not present in the early SAAO observations reached considerable amplitude in the combined two-site data obtained about 7 weeks later; we note that this is not due to an aliasing problem. As a matter of fact, the amplitude of this signal was still increasing during the latter measurements. It had about 12 mmag in the SARA data, but 18 mmag in the last set of SAAO measurements, in which it also showed a significant phase variation. A combination frequency ($`f_A+f_B`$) not present initially also appeared in the two-site data. On the other hand, the signal $`f_C`$ was found to be constant in amplitude and phase over the whole data set, and $`f_A`$ was constant in amplitude but showed a slight frequency shift of about 0.4 $`\mu `$Hz between our two subsets of data.

The residual amplitude spectra in the lowest panel of Fig. 6 still show a number of interesting peaks, several of which may be real. To estimate the richness of the frequency spectrum of PG 1456, we performed further prewhitening in the two-site data (still checking with the early SAAO measurements), taking the temporal variations of $`f_B`$ into account. The resulting amplitude spectra are shown in Fig. 7.

It becomes clear that the light variations of PG 1456 are quite complicated and that our data set is not sufficient for finding all the frequencies present. A preliminary multifrequency solution is given in Table 3.

A few comments are still necessary. Signals in the region around 1970 $`\mu `$Hz are present in both data sets with similar amplitude (e.g., the 1973.4 $`\mu `$Hz signal in Table 3), but their frequencies are different by 5 $`\mu `$Hz, even taking aliasing into account. The dominant frequency region in the residual single site-data is around 1304 $`\mu `$Hz, but neither this frequency nor its aliases are present in the two-site data. It seems that only intensive multisite observations will make a good understanding of the frequency spectrum of PG 1456 possible.

The periods we determined are in very good agreement with those detected by Grauer et al. (1988). Regrettably, a detailed comparison of these periods is impossible as the discovery data are no longer available (Grauer, private communication).

## 4 Interpretation of the observational results

### 4.1 CBS 114

The independent pulsation periods of CBS 114 as listed in Table 2 seem to have similar spacing. We therefore searched for and determined their mean spacing by calculating the Fourier power spectrum of the period values with unit amplitude (Handler et al. 1997), which we show in Fig. 8. A mean period spacing (significant at the 98% level) is indicated, amounting to 37.1$`\pm `$0.7 s. This result is corroborated by a Kolmogorov-Smirnov test (not shown).

We interpret this preferred period separation as a sign of the excitation of a number of pulsational radial overtones with the same spherical degree $`\mathrm{}`$. Some small deviations around its mean value (see Fig. 9 in Sect. 5) are an indication of mode trapping. The presence of a preferred period separation is predicted by asymptotic theory (Tassoul 1980), and was already observed within e.g. the $`\mathrm{}=1,m=0`$ modes of GD 358 (Winget et al. 1994, Vuille et al. 2000), for which a mean period spacing of 39.2 s was found. In addition, Bradley & Winget (1994) showed that the systematic deviations from this mean spacing were due to mode trapping.

We cannot directly identify the $`\mathrm{},m`$ values of the modes we observed in CBS 114 because of the lack of any rotational $`m`$-mode splitting in our frequency spectra. However, it seems very likely that the independent modes of CBS 114 are all $`\mathrm{}=1`$: if we saw a mixture of $`\mathrm{}`$ values, we would not expect to find such a significant mean period spacing. Therefore the modes we observed must originate exclusively (or at least predominantly) from the same $`\mathrm{}`$. They should be $`\mathrm{}=1`$ or 2 as the effects of geometrical cancellation (Dziembowski 1977) are expected to render modes of higher $`\mathrm{}`$ photometrically undetectable.

The mean period spacing between consecutive overtones of a pulsating white dwarf star is a measure of its mass (Kawaler 1987). If the modes of CBS 114 were $`\mathrm{}=1`$, their mean period spacing would be consistent with the star being a bit more massive than GD 358. However, if the modes were all $`\mathrm{}=2`$, CBS 114 would need to have a mass below 0.3$`M_{\mathrm{}}`$ (see Bradley, Winget & Wood 1993). Such a low mass is inconsistent with the spectroscopic gravity of the star (Beauchamp et al. 1999) and would be quite unusual considering the mass distribution of the DB white dwarf stars (Beauchamp et al. 1996).

Based on these arguments, we suggest that CBS 114 is a DBV star pulsating predominantly in nonradial g-modes of spherical degree $`\mathrm{}=1`$. We note that its individual pulsation periods, especially the longer period modes, are very similar to those of GD 358.

### 4.2 PG 1456

The interpretation of the frequency spectrum of PG 1456 is more difficult. A search for a mean period spacing such as in the previous section did not give a significant result. This suggests that the modes we detected are not of the same $`\mathrm{}`$ or that the star is a fast rotator. The modes of $`m0`$ would then mask the possible patterns of equally spaced periods.

About half of the indpendent periods we detected or suspect are similar to periods of GD 358, but the others are not. We are therefore unable to interpret the mode spectrum of PG 1456 at this stage. Extensive multisite observations, e.g. with the Whole Earth Telescope (Nather et al. 1990), are required to understand this star’s pulsations.

## 5 Asteroseismology of CBS 114

Using the optimization method developed by Metcalfe, Nather, & Winget (2000) and Metcalfe, Winget, & Charbonneau (2001), we performed a global search for the optimal model parameters to fit the 7 independent pulsation periods of CBS 114 listed in Table 2. The method uses a parallel genetic algorithm to minimize the root-mean-square (rms) differences between the observed and calculated periods for models with effective temperatures ($`T_{\mathrm{eff}}`$) between 20,000 and 30,000 K, total stellar masses ($`M_{}`$) between 0.45 and 0.95 $`M_{}`$, and helium layer masses with $`\mathrm{log}(M_{\mathrm{He}}/M_{})`$ between 2.0 and $``$7.0. This technique has been shown to find the globally optimal set of parameters consistently among the many possible combinations in the search space, but it requires between $``$10 and 4000 times fewer model evaluations than an exhaustive search of parameter-space to accomplish this (depending on the number of free parameters), with a failure rate $`<10^5`$.

We assumed that all of the observed modes had spherical degree $`\mathrm{}=1`$ (as suggested in Sect. 4.1) and azimuthal order $`m=0`$. The latter assumption can bias our determination of the particular set of model parameters that produces the optimal fit to the data, but if the rotation period is $``$1 day, any set of periods drawn from $`m=(1,0,+1)`$ will produce essentially the same overall picture of the parameter-space. We demonstrated this general behavior by generating C-core fits to 100 data sets for GD 358 that used randomly-selected $`m`$-components from those identified in Winget et al. (1994), and the $`m=0`$ values for $`k=12,18`$ (for which no triplet structure was found). In every case, the optimal set of model parameters fell within the same families of good solutions that were identified using the $`m=0`$ modes. Thus, when the spherical degree of the modes is known, the model-fitting procedure correctly identifies the families of possible solutions even when the values of $`m`$ are unknown. Furthermore, we found that the optimal solutions fell into the various families within parameter-space in proportion to the relative fitness of that family when only $`m=0`$ modes were used for the fit. For example, if a family has a peak fitness in the $`m=0`$ case that is 3$`\sigma `$ better than any other family, then the use of $`m0`$ modes will identify an optimal solution in this same family with $``$50% probability. This means that if one of the families produces much better solutions than the others, we can be reasonably confident that we have at least identified the correct family of solutions.

### 5.1 Carbon-core fit

Assuming a pure C core extending to a fractional mass of $`0.95m/M_{}`$, the optimal set of model parameters found by the genetic algorithm for CBS 114 were

$$T_{\mathrm{eff}}=24,600\mathrm{K},M_{}=0.655M_{},\mathrm{log}(M_{\mathrm{He}}/M_{})=3.96,$$

with rms period residuals $`\sigma _P=1.02s`$. The observed and calculated periods are shown in the top panel of Fig. 9 plotted against the deviations from the mean period spacing (dP), which were calculated using the same set of periods in both cases. Each point in this representation of the data is independent of the others, unlike a period spacing diagram using $`\mathrm{\Delta }P`$ ($`P_{k+1}P_k`$; cf. Bradley & Winget 1994). Note that the genetic algorithm only fits the periods of the pulsation modes, and the agreement between the deviations from the mean period spacing is simply a reflection of the overall quality of the match.

As noted in Table 2, the daily aliases of several of the identified modes in CBS 114 cannot be ruled out entirely. We can assess the impact of this uncertainty on the final set of optimal model parameters by repeating the fitting procedure using one or more of the alias periods. However, since the genetic algorithm fitting method with 3 free parameters is only about ten times more efficient than calculating the entire grid of models, it is better to calculate the pulsation periods of all $`10^6`$ models if we intend to repeat the procedure more than a few times. As there are four periods with a total of five possible aliases, we chose to calculate the entire grid of models. This allows us to check the answer that resulted from the genetic algorithm fit, and will enable us to generate C-core fits very quickly in the future for any additional data sets on DBV white dwarf stars.

All C-core models with rms period residuals smaller than 2.22 seconds are shown in the top half of Fig.10, which includes front and side views of the entire search space. Each point in the left panel corresponds to a point in the right panel, and the darkness of a point indicates the relative quality of the match with the observations. Black points are within $`\mathrm{\Delta }\sigma _P=0.03s`$ of the optimal model, and the four progressively lighter shades of grey correspond to models within 3, 10, 25, and 40 times this difference. The parameter-correlations explained by Metcalfe, Nather, & Winget (2000) are clearly visible, causing the good models to fall along lines in the plot rather than on a single point.

As expected, the optimal parameters from comparison of the observations to the complete grid of models were identical to those found using the genetic algorithm method. In fact, the identified periods led to lower rms residuals than any set of periods that included one of the possible aliases. Of the combinations with two periods replaced by their aliases, fits with lower residuals were achieved in three cases: (1) $`f_1f_1^{}`$ and $`f_7f_7^+`$, (2) $`f_1f_1^{}`$ and $`f_7f_7^{}`$, and (3) $`f_6f_6^{}`$ and $`f_7f_7^{}`$. When either three or four periods were replaced with aliases, only one combination led to lower residuals: $`f_1f_1^{}`$, $`f_5f_5^+`$, and $`f_7f_7^+`$. Note that aliases of $`f_1`$ and $`f_7`$ appear in most of these alternative period lists, implying that they are relatively important to the outcome of the fit. Although lower residuals were possible using these various combinations of two or more alias periods, the optimal models in every such case are too massive and too cool to be reconciled with the spectral line fits for CBS 114 by Beauchamp et al. (1999).

### 5.2 Oxygen-core fit

The optimal set of model parameters found by the genetic algorithm for a pure O core extending to $`0.95m/M_{}`$ had a mass that was inconsistent with the measurements of Beauchamp et al. (1999) and a helium layer mass at the edge of our search range, near the theoretical limit where nuclear burning will occur at the base of the envelope. The second-best model had residuals only $``$0.03 seconds higher and was similar to the optimal C-core model. The parameters of this model were

$$T_{\mathrm{eff}}=25,800\mathrm{K},M_{}=0.640M_{},\mathrm{log}(M_{\mathrm{He}}/M_{})=3.96,$$

with rms period residuals $`\sigma _P=0.91s`$. The middle panel of Fig. 9 shows the calculated periods of this model along with the observations, plotted against the deviations from the mean period spacing. The fit is only slightly better than the C-core model, and is qualitatively similar.

The O-core models that were calculated by the genetic algorithm during the optimization process with rms period residuals smaller than 2.08 seconds are shown in the bottom half of Fig. 10. Although this is not a complete sampling of the parameter-space, it is heavily sampled in regions where models produce better than average residuals. Again, the shade of each point indicates the relative quality of the match—black points are within $`\mathrm{\Delta }\sigma _P=0.03s`$ of the optimal model, and the shades of grey are at 3, 10, 25, and 40 times this difference. Note that because the optimal O-core model has residuals $``$0.1 seconds lower than the optimal C-core model, the darkest grey points in the bottom half of Fig. 10 are actually better than the the black points in the top half.

### 5.3 C/O-core fit

Since it may seem dubious to fit the 7 observed periods of CBS 114 using a model with five free parameters, we proceed cautiously. To allow a systematic exploration of models with various internal C/O profiles, Metcalfe, Winget, & Charbonneau (2001) used a simple parameterization that explored a general class of profiles similar to those used by Bradley, Winget, & Wood (1993). The parameterization fixes the oxygen mass fraction to its central value ($`X_\mathrm{O}`$) out to some fractional mass ($`q`$) where it then decreases linearly in mass to zero oxygen at $`0.95m/M_{}`$.

Using this method, the optimal model parameters for the observed pulsation periods of CBS 114 were

$$T_{\mathrm{eff}}=21,000\mathrm{K},M_{}=0.730M_{},$$

$$\mathrm{log}(M_{\mathrm{He}}/M_{})=6.66,X_\mathrm{O}=0.61,q=0.51,$$

with rms period residuals $`\sigma _P=0.43s`$. All models within Delta sigma = 0.03 s of this fit had the same mass, temperature, and helium layer mass at our sampling resolution of 0.005 $`M_{\mathrm{}}`$, 100 K, and 0.05 dex respectively. Only models with a central oxygen mass fraction within Delta $`X_O=\pm 0.01`$ and the optimal value of q had rms residuals within this range.

The calculated periods of this model are shown in the bottom panel of Fig. 9 along with the observed periods, plotted against the deviations from the mean period spacing. To evaluate whether or not this fit is better by an amount that justifies the addition of two free parameters, we can use the Bayes Information Criterion (BIC, following Koen & Laney 2000):

$$\mathrm{BIC}=N_p\left(\frac{\mathrm{log}N_{\mathrm{obs}}}{N_{\mathrm{obs}}}\right)+\mathrm{log}\sigma ^2$$

where $`N_p`$ is the number of free parameters, $`N_{\mathrm{obs}}`$ is the number of observed periods, and $`\sigma `$ is the rms period residual of the optimal fit. The value of BIC must be lower for a decrease in $`\sigma `$ to be considered significant. Our best 3-parameter fit was the O-core model, which leads to a value of $`\mathrm{BIC}=0.28`$. This leads us to expect the residuals of a 5-parameter fit to decrease to 0.69 seconds without being considered significant. The rms residuals of our 5-parameter fit are substantially lower than 0.69 seconds, so the addition of the extra parameters seems to be justified, and the C/O fit is significantly better than the O-core model.

Even so, we might worry that our optimal model parameters may be less accurate than the values obtained for GD 358, due to the smaller number of observed pulsation periods. To determine the magnitude of any systematic uncertainties due to the smaller number of observed modes, we performed a new fit to the pulsation periods of GD 358, but using only the 7 modes corresponding to those observed in CBS 114: $`k=8,9`$ and 11-15. When we compared the result of this fit with the one using all 11 modes observed in GD 358, we found only small shifts to the optimal model parameters:

$$\mathrm{\Delta }T_{\mathrm{eff}}=500\mathrm{K},\mathrm{\Delta }M_{}=+0.015M_{},$$

$$\mathrm{\Delta }\mathrm{log}(M_{\mathrm{He}}/M_{})=0.05,\mathrm{\Delta }X_\mathrm{O}=0.02,\mathrm{\Delta }q=0.01.$$

By analogy, we might expect our 5-parameter fit to CBS 114 to represent a slight underestimate of the temperature and overestimate of the mass. Both potential biases help to explain part of the discrepancy between the mass and temperature of our optimal model for CBS 114 and the values inferred from the spectroscopic analysis of Beauchamp et al. (1999) \[$`T_{\mathrm{eff}}=23,300`$; $`\mathrm{log}g=7.98`$\]. Additional differences are expected since we have used different mixing-length parameters (ML3) than Beauchamp et al. (ML2/$`\alpha =1.25`$). Metcalfe, Salaris, & Winget (2002) quantified the offsets between fits to GD 358 using ML3 and ML2, and found them to be approximately the same size as the shifts due to the smaller number of data points. Unfortunately, CBS 114 does not have a published parallax, so an independent constraint on the mass and temperature from the luminosity is not presently available.

## 6 Conclusions

We have presented time-resolved CCD photometry of the pulsating DB white dwarf stars CBS 114 and PG 1456+103. Our data, obtained with telescopes of only 0.75 m and 0.9 m apertures, are sufficient for the extraction of useful asteroseismological information even for a star as faint as CBS 114 ($`B17.2`$). We have also shown that it is possible to understand the mode spectra of multiperiodic pulsating white dwarf stars from single-site observations in suitable cases (see Handler 2001 for another example). In other cases, like PG 1456, this is not possible, and worldwide observing campaigns, e.g. with the Whole Earth Telescope, are necessary.

The frequency analysis of our measurements of CBS 114 resulted in the discovery of a mean period spacing of $`37.1\pm 0.7`$ s in the independent modes of the star. Although no convincing evidence for rotational splitting within the modes was detected, we argued that the star pulsates in nonradial $`\mathrm{}=1`$ g-modes.

We used our data set to check the sensitivity of the genetic-algorithm-based model-fitting method (Metcalfe et al. 2000) to artifacts of data analysis, like aliasing problems and our lack of assignments of the azimuthal order $`m`$ to the observed modes. Encouragingly, our results are hardly affected by these shortcomings.

Most significantly, the genetic algorithm found a C/O fit to the pulsation periods of CBS 114 that proved to be significantly better than the pure-core fits, with an accuracy comparable to a similar fit to GD 358. The optimal mass and central oxygen mass fraction from this fit can provide an independent measurement of the astrophysically-important rate for the $`{}_{}{}^{12}\mathrm{C}(\alpha ,\gamma )^{16}\mathrm{O}`$ reaction, as was recently done for GD 358 (Metcalfe, Salaris, & Winget 2002). In general, higher mass white dwarfs are expected to have a lower central oxygen mass fraction, because the 3$`\alpha `$ rate rises faster with increasing density than the $`{}_{}{}^{12}\mathrm{C}(\alpha ,\gamma )^{16}\mathrm{O}`$ rate. Our optimal model for CBS 114 has a higher mass and a lower central oxygen mass fraction than the optimal model for GD 358 (Metcalfe, Winget, & Charbonneau 2001), so it is consistent with this general trend. A model of the internal chemical profile with the same mass as our fit to CBS 114 requires a rate for the $`{}_{}{}^{12}\mathrm{C}(\alpha ,\gamma )^{16}\mathrm{O}`$ reaction near $`S_{300}=180`$ keV b to produce a central oxygen mass fraction of 0.61 (M. Salaris, private communication). This value is close to the rate derived from recent high-energy laboratory measurements ($`S_{300}=165\pm 50`$ keV b; Kunz et al. 2002). By contrast, the rate derived from the optimal model of GD 358 by Metcalfe, Salaris, & Winget (2002) was significantly higher ($`S_{300}=370\pm 40`$ keV b). This suggests either that presently unknown sources of systematic uncertainty in our models must affect the analysis of GD 358 and CBS 114 in different ways, or that the two stars have different evolutionary origins, or both. An asteroseismological determination of the central oxygen mass fraction for additional DBV white dwarfs will help us to decide which of these scenarios is most likely.

## ACKNOWLEDGEMENTS

We thank Don Winget and Chuck Claver for permission to use their data of CBS 114, Maurizio Salaris for providing a new internal chemical profile to match our optimal model and Mike Montgomery for carefully reading a draft version of this paper.


format: text

text:
# High-energy magnon dispersion in the half-filled Hubbard model: A comparison with La2CuO4

\[

## Abstract

We use quantum Monte Carlo methods and single-mode approximation to study the magnon dispersion in the 2D half-filled Hubbard and phonon-coupled Heisenberg models. We find that in the Hubbard model with $`U/t<8`$, high-energy magnon dispersion is similar to those observed in inelastic neutron scattering experiments in $`\text{La}_2\text{CuO}_4`$. On the other hand, our studies of a 2D Heisenberg model coupled to dynamic optical bond phonons, fails to reproduce the experimental dispersion. These results can be interpreted as evidence for intermediate $`U/t`$ and charge fluctuations in the cuprate materials.

\]

While there is still no consensus on the microscopic origin of superconductivity in high-temperature copper oxide superconductors, it is widely believed that magnetic fluctuations play an important role. This has led to extensive studies of magnetic fluctuations in the parent compounds of these materials, such as $`\text{La}_2\text{CuO}_4`$. The low-energy long-wavelength properties of these layered compounds are well described by the non-linear sigma model and the renormalized classical theory . However, short wavelength or high energy spin-fluctuations may also be relevant to superconductivity, which only sets in after long-range antiferromagnetic order is lost. The description of these excitations has raised many questions. Are there well-defined magnons at short wavelengths? Is there significant spectral weight in multimagnon excitations? Do the spinons present a better description of short-wavelength excitations? Is there a coexistence of spinons and magnons? The most quantitative theoretical studies have been done for the nearest-neighbor Heisenberg model, where one finds magnons throughout the Brillouin Zone, with significantly reduced spectral weight along the antiferromagnetic zone boundary .

On the experimental side, Raman scattering has long provided evidence for anomalous short-wavelength spin fluctuations . The line-shape of the two-magnon Raman spectra has not been adequately explained. Similar results are inferred from optical absorption spectra . Understanding these spectra is complicated by the fact that light couples directly to the charge degrees of freedom and thus their understanding within a spin-only picture requires many assumptions . In contrast, neutron scattering provides a direct probe of the spin excitations . Recently, the high energy magnon dispersion has been measured in La2CuO4 by inelastic neutron scattering. One finds magnon dispersion at short wavelengths, with ($`\pi /2,\pi /2`$) magnon energy about $`13\%`$ lower than that at ($`\pi ,0`$) . This is in contrast to results for the nearest-neighbor Heisenberg model, where one finds the magnon energy at ($`\pi ,0`$) to be $`710\%`$ lower than at ($`\pi /2,\pi /2`$) . The latter results have also been seen experimentally in the materials $`\text{Cu(DCOO)}_2.4\text{D}_2\text{O}`$  and $`\text{Sr}_2\text{Cu}_3\text{O}_4\text{Cl}_2`$ .

Explanation for the anomalous dispersion in $`\text{La}_2\text{CuO}_4`$ has evoked considerable interest. Adding a second neighbor antiferromagnetic exchange takes the spectra in the opposite direction, ruling out the simplest possibility. One mechanism that has been suggested is ring exchange of four electrons in a plaquette . Such a term can be directly obtained from the Hubbard model and was first proposed for the high-$`T_c`$ materials to explain the Raman spectra more than a decade ago . One interesting aspect of the ring-exchange term is that, when treated fully, it does not change the spin-wave dispersion within linear spin-wave theory. Thus its effect on the dispersion is a purely quantum effect. Recent numerical studies have focused on whether such a term can explain the experimental dispersion without destabilizing Néel order .

In a recent paper, Peres and Araújo have studied the 2D Hubbard model using mean-field theory and obtained a dispersion relation similar to that observed in $`\text{La}_2\text{CuO}_4`$ for an intermediate value of the on-site interaction parameter, $`U/t=6`$, indicating that charge fluctuations – ignored in the Heisenberg model – need to be taken into account to explain the properties of this material. Given the inherently approximate nature of the mean-field theory, it is important to confirm and extend this result using independent approaches. We have used the determinant quantum Monte Carlo (det QMC) method and single mode approximation (SMA) to study the magnon dispersion in the 2D Hubbard model at half-filling. The det QMC method has been extensively used to study the ground state properties of the 2D Hubbard model at half-filling where the fermion sign problem can be avoided. We have studied the magnon dispersion along the magnetic zone boundary as a function of the on-site interaction parameter, $`U/t`$. We find that with decreasing $`U`$, the zone-boundary dispersion changes sign and for $`U/t<8`$ it becomes similar to those observed in $`\text{La}_2\text{CuO}_4`$.

Coupling to phonons is another potential source of anomalous dispersion in the cuprates, because of the rather high magnon energies. We have used the stochastic series expansion QMC to study a Heisenberg model coupled to optical bond phonons. We find that such a model fails to reproduce the experimental spectra. While suggestive, this, however, does not rule out the possibility that a more realistic treatment of the spin-phonon couplings can mimic the experimental results.

The Hubbard model in two dimensions is given by the Hamiltonian

$`H`$ $`=`$ $`t{\displaystyle \underset{𝐢,𝐣\sigma }{}}(c_{𝐢,\sigma }^{}c_{𝐣,\sigma }+c_{𝐣,\sigma }^{}c_{𝐢,\sigma })`$ (2)
$`+U{\displaystyle \underset{𝐢}{}}(n_{𝐢,}{\displaystyle \frac{1}{2}})(n_{𝐢,}{\displaystyle \frac{1}{2}})+\mu {\displaystyle \underset{𝐢}{}}n_𝐢,`$
where $`c_{𝐢,\sigma }^{}(c_{𝐢,\sigma })`$ creates (annihilates) an electron with spin $`\sigma `$ at lattice site $`𝐢`$. The kinetic energy term includes a sum over nearest neighbors $`𝐢,𝐣`$ and $`t`$ is the hopping integral between adjacent sites. $`U`$ is the on-site interaction, and $`\mu `$ is the chemical potential. We shall be dealing solely with the half-filled band, ie., $`n_{𝐢,}+n_{𝐢,}=1`$. With the interaction term written in a particle-hole symmetric from as above, the half-filled band corresponds to setting $`\mu =0`$. Henceforth we set $`t=1`$ and express the interaction parameter $`U`$ in units of $`t`$.

In the limit of large $`U`$, the Hubbard model at half-filling maps on to the Heisenberg model with exchange parameter $`J=4t^2/U`$. The Heisenberg model in 2D is given by the Hamiltonian

$$H=J\underset{𝐢,𝐣}{}𝐒_𝐢.𝐒_𝐣$$
(3)
The spin operators are related to the electron creation and annihilation operators by the relation

$`𝐒_𝐢={\displaystyle \frac{1}{2}}c_{𝐢,\sigma }^{}\stackrel{}{\tau }_{\sigma ,\sigma ^{}}c_{𝐢,\sigma }`$

where $`\stackrel{}{\tau }`$ are the Pauli matrices.

The det QMC used to study the 2D Hubbard model is a finite temperature method that is based on discretizing the imaginary time $`\beta =L\mathrm{\Delta }\tau `$ and employing Trotter approximation to decompose the full imaginary-time evolution. This approach treats the electron-electron interactions exactly and at half-filling, where we focus in this work, is able to produce results with small statistical fluctuations. Ground state expectation values are obtained using sufficiently large values of $`\beta `$. Like most other numerical approaches, the technique is limited to finite-size lattices and we have been able to study primarily two different lattice sizes, $`N`$=6x6 and $`N`$=8x8. While this is not sufficient to do a complete finite-size scaling analysis to obtain thermodynamic expectation values, previous studies have shown that a lattice of size $`N`$=8x8 is large enough to give reasonably good estimates of the thermodynamic limit, especially when measuring quantities away from the antiferromagnetic wave vector.

The Heisenberg model has been studied using the Stochastic Series Expansion (SSE) QMC method. The SSE method is also a finite-temperature QMC method based on the importance sampling of the diagonal matrix elements of the Taylor expansion of $`e^{\beta H}`$. Ground state expectation values can be obtained by using sufficiently large values of $`\beta `$, and there are no approximations beyond statistical errors. With the recently developed “operator loop update”, the method has proven to be very efficient tool for studying several different models.

In both approaches, we have measured the equal-time spin-spin correlation function $`S_𝐢S_𝐣`$ and the associated static spin structure factor in the momentum space given by

$$S(𝐪)=\frac{1}{N}\underset{𝐢,𝐣}{}e^{𝐪.(𝐢𝐣)}(n_𝐢n_𝐢)(n_𝐣n_𝐣)$$
(4)
We also evaluate the static spin susceptibility in the momentum space given by the Kubo integral

$$\chi (𝐪)=_0^\beta 𝑑\tau S(𝐪,\tau )S(𝐪,0)$$
(5)
The magnon energy is calculated using the relation

$$\omega (𝐪)=2S(𝐪)/\chi (𝐪).$$
(6)
We start the discussion of the numerical results by studying the convergence of the calculated quantities as a function of the inverse temperature, $`\beta `$. Figure 1 shows a plot of the magnon energies at the edge, $`(0,\pi )`$, and the center, $`(\pi /2,\pi /2)`$, of the zone boundary as a function of $`\beta `$ for $`U=6.0`$. The magnon energies are seen to converge to their $`T=0`$ values fairly rapidly with $`\beta `$. This convergence at relatively small values of $`\beta `$ is not surprising given that we only deal with high-energy quantities. Henceforth, all the results presented are for $`\beta =3.0`$.

Next we present the results for the magnon energies along the entire magnetic zone boundary for $`U=6.0`$. This value of $`U`$ is consistent with estimates of the effective on-site energy in $`\text{La}_2\text{CuO}_4`$ obtained from photoemission and optical spectroscopy data. Figure 2 shows the variation of the magnon energy for a path along the Brillouin zone similar to the one considered in Figure 3a of ref. for two different lattice sizes. Note that the data near ($`\pi ,\pi `$) and ($`0,0`$) are strongly temperature dependent, the latter being exactly proportional to temperature. However, as shown in Figure 1, the data along the zone boundary from ($`\pi ,0`$) to ($`\pi /2,\pi /2`$) have reached close to their zero temperature values. The qualitative nature of the variation of the magnon energy along the path is in good agreement with the experimental results and mean-field predictions. In particular, the magnon energy at the center of the zone boundary, $`(\pi /2,\pi /2)`$, is found to be lower than that at the edge, $`(0,\pi )`$, as seen in the experimental data and opposite to the dispersion found for the Heisenberg model. The magnitude of the deviation is estimated to be $`12\%`$ of the $`N`$=8x8 system, matching closely the experimental observation. While this quantitative agreement with experimental data is not conclusive– the magnon energies were evaluated in our calculation using the SMA, and the finite-size effects are fairly large – it is safe to conclude that the 2D Hubbard model, with a value of the on-site interaction energy in the range estimated from photoemission and optical spectroscopy data, can reproduce the correct magnon dispersion along the magnetic zone boundary.

We now discuss the variation of the dispersion along the zone boundary in the 2D Hubbard model as a function of the on-site interaction parameter, $`U`$. Figure 3 shows the results of our simulation for a $`N`$=8x8 lattice at $`\beta =3.0`$. We have focused primarily on the deviation of the magnon energy at the center of the zone boundary relative to the edges. To that end, we have scaled the magnon energies by their value at $`(0,\pi )`$. The deviation is seen to be maximum for the smallest value of $`U`$ studied and decreases with increasing $`U`$ – becoming essentially flat within statistical error at $`U=8.0`$. Unfortunately, for $`U>8`$, the results of our simulation for the 2D Hubbard model become too noisy. However, in this limit of large $`U`$, the Hubbard model maps on to the Heisenberg model. Hence the dispersion for larger $`U`$ values is expected to be be qualitatively similar to that for the Heisenberg model. We have also shown in the figure the dispersion obtained for a Heisenberg model with exchange parameter $`J`$ corresponding to $`U=10.0`$. The point $`U=8.0`$ marks a transition to “Heisenberg-like” behavior. This value of the transition point is consistent with that found from specific heat measurements of the 2D Hubbard model .

To investigate other possible sources for the nature of magnon energy dispersion along the magnetic zone boundary, we have studied the effects of spin-phonon coupling within the framework of the 2D Heisenberg model. The model involves coupling of the spins to dynamic optical bond phonons and is given by the Hamiltonian

$$H=J\underset{𝐢,𝐣}{}(1+\alpha (a_b^{}+a_b))𝐒_𝐢.𝐒_𝐣+\omega _0\underset{b}{}a_b^{}a_b$$
(7)
where $`b`$ denotes the bond between the lattice sites i and j and the operator $`a_b^{}(a_b)`$ creates(annihilates) a phonon on bond $`b`$. $`\omega _0`$ represents the bare phonon frequency and $`\alpha `$ measures the strength of the spin-phonon interaction.

The spin-phonon model was studied using the SSE QMC method. Our results indicate that over any reasonable range of spin-phonon coupling strength and bare phonon frequency, the magnon dispersion along the magnetic zone boundary is qualitatively the same as the pure Heisenberg model. Figure 4 shows the variation of the magnon energy along the same path in momentum space as considered in Figure 2 for the pure Heisenberg model and the spin-phonon model with $`\alpha =0.1`$ and $`\omega _0=0.25`$ for a lattice of size $`N`$=8x8 at $`\beta `$=5.0. The value of $`\beta `$ is chosen such that $`\beta J`$ is the same as $`\beta t`$ for the data presented in Figure 2 with $`J=4t^2/U`$. This allows for a direct comparison between the two sets of data. Once again we found that while the data near $`(\pi ,\pi )`$ and (0,0) are strongly temperature dependent, the magnon energies along the magnetic zone boundary from $`(\pi ,0)`$ to $`(\pi /2,\pi /2)`$ have converged close to their ground state values. The data suggests that the deviation in the magnon dispersion observed for $`\text{La}_2\text{CuO}_4`$ from that found in the Heisenberg model cannot be explained by the type of spin-phonon coupling considered here. It remains to be seen whether a more realistic treatment of the electron-phonon coupling in the cuprates can account for this behavior.

To summarize, we have used the det QMC to show that the experimentally observed magnon dispersion along the magnetic zone boundary in $`\text{La}_2\text{CuO}_4`$ can be reproduced by the 2D Hubbard model, using reasonable $`U/t`$ values. The deviation of the magnon energy at the center of the zone boundary is maximum for the smallest value of $`U`$ considered, where the charge fluctuations are the largest, and decreases with increasing $`U`$. For $`U8t`$, the dispersion is flat within statistical errors. The dispersion for larger values of $`U`$ is qualitatively similar to that of the pure Heisenberg model. Thus, at $`U8t`$, there is a transition to “Heisenberg-like” behavior where the effects of the charge fluctuations become negligible. An interesting question is whether this intermediate $`U/t`$ regime of the Hubbard model is equivalent to a Heisenberg plus ring exchange term. At really small $`U/t`$ the electrons will be strongly delocalized, so that a description in terms of a spin Hamiltonian would break down. However, at the intermediate $`U/t`$, found in our calculations, a higher order perturbative treatment maybe justified leading to a Heisenberg plus ring exchange Hamiltonian.

We would like to thank Anders Sandvik for many useful discussions. This work was supported in part by NSF grant number 9986948. The numerical simulations were carried out in part at the Condor Flocks at University of Wisconsin, Madison and NCSA, Urbana, Illinois.


text:
# On an Application of Relative Entropy

## Abstract

We show that in problems of authorship attribution and other linguistic applications, a Markov Chains approach is a more attractive technique than Lempel-Ziv based compression.

We wish to point out a number of inaccurate and misleading statements that Benedetto et al. make in their paper titled “Language Trees and Zipping”Benedetto et al. (2002). First, they claim the technique they used for construction of a language tree does not make use of any a-priori information about the alphabet, but it does, both in the alphabet chosen (Unicode) and in the set of languages they chose to experiment with; second, they propound Lempel-Ziv (LZ, gzip) compression as being applicable to DNA analysis, where the usefulness of LZ is quite doubtful; third, in practice their definition of relative entropy and distance can yield negative values; fourth, the classification performance of the method they use is significantly worse than other entropy-based methods as has been noted in prior work; and fifth, the classification speed is significantly worse as well, which shows that its “potentiality” is questionable. We elaborate on each of these points in more detail in the subsequent paragraphs.

Notice that the “Language Tree”(LT) diagram Benedetto et al. (2002) does not include the Russian language (Slavic family of Indo-European family of languages; 288 million speakers). Our computations show that once Russian is included, it does not cluster with the other members of the Slavic group. Obviously, certain Cyrillic alphabet based languages were left out of the study Benedetto et al. (2002), which “improves” results significantly and shows that a-priori information about the alphabet is being taken advantage of to achieve the results outlined in paper Benedetto et al. (2002).

The LZ compressor makes few assumptions about the input string, but in practice, we do have a-priori information that we can take advantage of. Biologists widely use an amino acid substitution matrix (PAM250 or BLOSUM62) in search for similar biological sequences Gusfield (1997). It is not at all clear how a substitution matrix could be implemented with the LZ algorithm. That is why compression is not widely used for DNA analysis, although first trials for its application go back to 1990 Gusfield (1997).

The quantity $`S_𝒜`$ Benedetto et al. (2002) defined as “relative entropy” in (1) and redefined as “distance” in (2) can take negative values. Negative values indeed appeared in our study which showed that the “LT”Benedetto et al. (2002) reflects significantly the structure of Unicode or vice versa, and its relevance to language classification should be supported additionally.

A traditional definition and estimates for (relative) entropy via $`n`$th order Markov Chain on letters Shannon (1948); Yaglom and Yaglom (1983); Kukushkina et al. (2001) always lead to a proper positive number. Markov Chains are also traditional in text entropy analysis Shannon (1948); Yaglom and Yaglom (1983), compression Cleary and Witten (1984), authorship and subject attribution Teahan (2000); Khmelev (2000). In Kukushkina et al. (2001), the classification performance of compression programs was compared with the Markov Chain approach Khmelev (2000). 82 authors of large enough texts ($`10^5`$ characters) were chosen. Afterwards 82 one-per-author texts were held out and used for control purposes. The classification algorithm Kukushkina et al. (2001) had to determine the author of each control text among 82 alternatives. The corresponding numbers of exact guesses for 15 compression programs and Markov Chains are presented in the following list Kukushkina et al. (2001):

Program(number of guesses): 7zip(39), arj(46), bsa(44), compress(12), dmc(36), gzip(50), ha(47), huff(10), lzari(17), ppmd5(46), rar(58), rarw(71), rk(52); Markov Chain approach (see Khmelev (2000)) 69 guesses.

Clearly, gzip is significantly outperformed by other compression algorithms and the first order Markov chain model Khmelev (2000). Notice also that in practical implementations, the gzip-based approach Benedetto et al. (2002) is significantly slower than the first order Markov chains method Khmelev (2000).

To sum up, in natural language processing (and, perhaps, in other fields) the $`n`$th order Markov chain models Khmelev (2000); Teahan (2000) are more appropriate than an LZ-approach Benedetto et al. (2002).