text:
# The Origin of Soft X-rays in DQ Herculis

## 1 Introduction

A classical nova eruption (see Shara 1989 for a review) is a thermonuclear runaway in a cataclysmic variable (CV; a semi-detached binary in which a white dwarf primary accretes from a late-type, usually M dwarf, secondary). The runaway is caused when $`10^4`$ M of hydrogen-rich gas is accreted onto the white dwarf surface; this material becomes degenerate, until the pressure and temperature at the base of this envelope starts the runaway reaction. A nova eruption ejects most of the accreted envelope. After a poorly known recurrence time (103 – 105 years), the underlying binary accumulates enough fresh fuel to erupt again. Novae are of interest to nuclear astrophysicists; they are useful as extragalactic distance indicators; they play a role in Galactic chemical evolution; and they are likely sources of $`\gamma `$-rays. To understand a nova, it is essential to study the underlying system away from its eruption. Conversely, the study of nova eruptions is essential to understand the secular evolution of CVs.

The subject of this paper, DQ Herculis (Nova Herculis 1934), has played a key historical role in the studies of nova eruptions and of post-nova binaries (hereafter old novae) due to its brightness and fortuitous viewing geometry. The 1934 classical nova eruption of DQ Her reached $`m_v1.3`$ and enabled extensive observations (see Martin 1989 for a review). The discovery of eclipses in DQ Her (Walker, 1956) was pivotal in establishing the standard picture of CVs. Walker (1956) also discovered the coherent 71 sec photometric oscillation: we now understand that DQ Her contains a magnetic white dwarf, accreting onto the magnetic polar regions which are offset from the rotational poles, i.e., an oblique rotator. The ejecta of the 1934 eruption are visible in H$`\alpha `$ and other emission lines as an elliptical shell: as of 1978, 43 years after the eruption, the equatorial and polar radii (expansion velocities) of the DQ Her shell were $`6.3\pm 0.2`$ arcsec ($`384\pm 21`$ km s-1) and $`8.6\pm 0.2`$ arcsec ($`528\pm 25`$ km s-1), respectively (Herbig & Smak, 1992). They found no evidence for a deceleration of the expansion, and derived a revised distance estimate of 561$`\pm `$19 pc (we use this value throughout this paper). This makes DQ Her among the nearest old novae known. By 1995 September (61 years after eruption), the equatorial and polar radii have grown to $`9.1\pm 0.8`$ arcsec and $`12.5\pm 0.8`$ arcsec, respectively, still consistent with no deceleration, according to our analysis of an archival Hubble Space Telescope WFPC2 image.

The eclipsing nature of DQ Her in principle allows us to place tight constraints on its system parameters. The optical eclipse width, measured from mid-ingress to mid-egress as is customary, is 0.110$`\pm `$0.003 cycles (as compiled from literature by Horne et al. 1993). Assuming that the source of eclipsed light is centered on the white dwarf, this can be used to constrain the relationship between the mass ratio $`q=M_2/M_1`$ (where $`M_1`$ and $`M_2`$ denote the masses of the white dwarf and the mass donor, respectively) and the inclination angle $`i`$. Horne et al. (1993) measured the radial velocity amplitude of the secondary $`K_2`$ through the detection of the Na I infrared doublet absorption lines. This, combined with the eclipse width and the $`K_1`$ from literature, have led to their estimates of $`q=0.66\pm 0.04`$, $`i=86^{}.5\pm 1^{}.6`$, $`M_1=0.60\pm 0.07`$ M, and $`M_2=0.40\pm 0.05`$ M.

About 20 X-ray bright CVs with coherent spin modulations have been discovered since the late 1970s, and are classified as intermediate polars (IPs) or DQ Her stars (Patterson, 1994). They are among the brightest CVs in the 2–10 keV band, with luminosities generally in excess of $`10^{32}`$ ergs s-1. However, the prototype DQ Her was not detected in an Einstein Imaging Proportional Counter (IPC) observation with a 2$`\sigma `$ upper limit of 0.0046 ct s-1 (Córdova et al., 1981). One of their hypotheses, that high inclination may be the cause of this non-detection, found support in a model of the phase shifts of the 71-s oscillation during eclipse (Petterson, 1980). In this model, a beam of high energy photons originating on the white dwarf is reprocessed by the accretion disk and results in the optical oscillation. Petterson (1980) finds a good fit at $`i89^{}`$, although the actual value of $`i`$ can be lower, as long as the entire front half of the disk including the white dwarf is obscured by the outer edge of the disk. Such permanent obscuration of the white dwarf would greatly reduce the X-rays reaching Earth.

Despite the high inclination angle, X-rays were detected from DQ Her with the ROSAT Position Sensitive Proportional Counter (PSPC; Silber et al. 1996)111Patterson (1994) was already aware of this detection, citing “Patterson and Eracleous 1994, in preparation.”. The inferred luminosity (0.1–2.0 keV) of $`4.5\times 10^{30}`$ ergs s-1 is modest for an IP. We note that, using a thermal plasma model with $`kT=0.3`$ keV with NH=1.0 $`\times 10^{20}`$ cm-2, consistent with the PSPC spectrum, and the revised distance of 561 pc (both are different from the values assumed by Córdova et al. 1981), the Einstein IPC upper limit corresponds to 2.8 $`\times 10^{30}`$ ergs s-1 in the IPC band (0.4–4 keV), or 3.3 $`\times 10^{30}`$ ergs s-1 in the 0.1–2.0 keV band. These numbers may reflect source variability or the limit of accuracy of cross-instrument comparisons; neither would be surprising. A deep eclipse is not observed, consistent with the permanent obscuration of the white dwarf, suggesting an alternative origin of the observed X-rays. Although Silber et al. (1996) considered a likely origin to be the secondary, it has to be emitting X-rays at significantly above the saturation limit ($`10^3`$ Lbol, or a few times 1029 erg s-1 in this case) seen in rapidly rotating late type dwarfs (see, e.g., Singh et al. 1999). If not the secondary, what is the origin of the soft X-rays observed from DQ Her?

Only $``$150 photons were detected in the ROSAT observation, insufficient for further progress. We have therefore secured a Chandra X-ray Observatory Advanced CCD Imaging Spectrometer (ACIS) observation of DQ Her. We describe the observation and data reduction in §2, present the results in §3, and discuss the implications in §4.

## 2 Observations

We observed DQ Her with Chandra (Weisskopf et al., 1996) ACIS without a grating, placing the object on the back-illuminated S3 chip222For details of the ACIS instrument, see the web site at http://www.astro.psu.edu/xray/acis.. The observations were performed as two pointings, one from 2001 July 26 13:00 UT to July 27 02:31 UT and the other from 2001 July 29 17:09 UT to July 30 02:27 UT, for a total integration time of 68.8 ksec. These observations were performed at near identical pointing positions and roll angles. We took the outputs from the processing pipeline at the Chandra X-ray Center, and extracted images, spectra, and light curves from these observations individually, before combining the results.

We detect a source at RA = 18h07m30.24s, Dec=+45d51m32.6s using the Ciao routine celldetect. Although this is offset by 1.7$`\mathrm{"}`$ from the position of DQ Her given in Downes et al. (2001), this is likely to reflect the limitation of comparison of positions measured using different astrometric reference frames (Downes 2002, private communication). We identify this X-ray source with DQ Her based on this positional coincidence and based on the orbital variability (see §3.2 below).

We extracted source spectra and light curves from a 5 arcsec radius circular region, and background from an annulus of 40 arcsec outer radius and 20 arcsec inner radius centered on the source. The size of the background annulus is chosen to exclude the nova shell emission (see §3.3). The source count rate is 0.025 cts s-1 (3 $`\times 10^{30}`$ ergs s-1 in the 0.2–5 keV band, based on the spectral model discussed below), while the background in the source region (scaled from the counts in the background region by the detector areas) is estimated to be 8.2$`\times 10^4`$ cts s-1. DQ Her did not show a significant difference in count rate between the two observations.

## 3 Results

### 3.1 Spectral Characteristics

There is a bump in the observed ACIS-S spectrum of DQ Her around 1 keV (left panel, Figure 1) which is not in the instrument response. Because of this, the spectrum cannot be fitted with smooth continuum models (power-law, bremsstrahlung or black body) alone. Using models of optically thin thermal plasma emission, such as mekal (Mewe et al., 1985, 1986; Liedahl et al., 1995), the 1 keV bump is interpreted as a complex of lines in the 0.6–1.2 keV range, notably lines of O VIII (0.64 keV), Fe XVII (0.73 & 0.83 keV), and Ne IX (0.90 keV). These lines, however, are weaker relative to the continuum than the predictions of single-temperature plasma models. Moreover, the continuum above 2 keV is much stronger than expected for plasma temperatures that give rise to these lines (kT$``$0.6 keV). As a result, a single component mekal model fit also fails to reproduce the observed spectrum.

The simplest successful representation ($`\chi _\nu ^2`$=1.4; see Figure 1) is a two-component model333We have also included the acisabs model in the fit, to account for instrumental absorption due to contaminant, at the appropriate level. consisting of kT=0.59$`\pm `$0.04 keV mekal plasma model plus a power law of photon index 2.7$`\pm `$0.2, with NH = 2.9 $`\pm 1.2\times 10^{20}`$ cm-2. The power-law component in our model provides the added continuum necessary to achieve an acceptable fit, both near the lines (0.8–1.2 keV) and at higher energies (2–5 keV). Inferred intrinsic flux in the 0.2–5 keV band is 8.2 $`\times 10^{14}`$ ergs cm-2s-1, or a luminosity of 3.1 $`\times 10^{30}`$ ergs s-1. DQ Her was not significantly detected above 5.0 keV (we obtain a count rate and an estimated 1$`\sigma `$ error in the 5–8 keV band of $`0.3\pm 2.0\times 10^4`$ cts s-1), although this is in part due to the lower effective area of the Chandra ACIS instrument at high energies. We take $`2.0\times 10^4`$ cts s-1 as our upper limit, which corresponds to $`4\times 10^{29}`$ ergs s-1 assuming a power law spectral shape with a photon index of 2 and a source distance of 561 pc.

Note, however, that our adopted model for the spectrum is not unique. Given the modest number of detected photons ($``$1700), and instrumental limitations (both in bandpass and in spectral resolution), there can be a wide range of models that fit the data. We are only confident that DQ Her shows a spectral feature around 1 keV which mimics a collisionally excited, optically-thin plasma of kT$``$0.6 keV. The presence of such a plasma would also lead to a Si line at 1.85 keV, which may also be present in the data. Consequently, spectral modeling alone of the Chandra data is insufficient to uncover the origin of the soft X-rays in DQ Her.

### 3.2 Timing Characteristics

We show, in Figure 2, the ACIS light curve of DQ Her folded on the orbital ephemeris of Zhang et al. (1995) into 16 bins per cycle. The ephemeris has an estimated accuracy of $`\pm `$11 s for this epoch (see §4.1 for further details). A partial eclipse, lasting for about 3 bins and about 30% deep, is clearly seen. In this representation, the X-ray eclipse center appears to be offset from the ephemeris prediction by one half of the bin size, i.e., 0.03 cycle or 500 s. Moreover, the X-ray eclipse appears to be wider (3 bins, or 0.1875 cycles) than the optical eclipse (0.110$`\pm `$0.003, according to a compilation in Horne et al. 1993), which should be 2 bins in Figure 2.

To quantify this impression, we have fit the 32 bin per cycle version of the light curve using a piecewise linear function with symmetric ingress and egress; we show our best-fit model in Figure 2. The eclipse center is found to be at phase 0.033$`{}_{0.027}{}^{}{}_{}{}^{+0.030}`$ (90% confidence), confirming the offset. However, the eclipse width (mid-ingress to mid-egress) is measured to be 0.18$`{}_{0.12}{}^{}{}_{}{}^{+0.05}`$ (90%) cycle, thus we cannot claim that the X-ray eclipse is definitely wider than in the optical. In addition to the formal uncertainties, there is additional uncertainty arising from our choice of specific model, thus both these numbers need confirmation. Our Chandra data are of insufficient quality for us to be able to derive further details of the eclipse profile, such as ingress/egress durations.

### 3.3 Imaging Characteristics

At first glance, the Chandra ACIS image of DQ Her looks like that of a point source. Indeed, over 70% of the flux is contained within 1 arcsec ($``$560 AU, or 8.4$`\times 10^{15}`$ cm) of the image centroid, as expected for a soft point source. However, a closer inspection of the image shows that there is a region about 2 arcsec by 2 arcsec in area and about 10 arcsec NE of the point source with more photons (28) than are expected by chance ($``$8). These excess counts appear equally distributed over this area, and does not display the sharply peaked core of the point spread function (PSF) expected of a point source.

We have investigated this further by constructing radial profiles of the ACIS image of DQ Her in several energy bands. We then fit a model consisting of a flat background plus the PSF available in the Chandra calibration database (interpolated in energy and off-axis angle as appropriate). We find notable excesses in the azimuthally-averaged radial profiles in the 0.2–0.5 and 0.5–0.8 keV bands (Figure 3; we find no significant features at energies above 0.8 keV, although this may in part due to poor counting statistics). The statistical significance of this is unfortunately difficult to establish: our experiences with other bright sources is that PSF fit is usually less than perfect because of effects of pile-up, background non-uniformity, uncertainties in the point source location, and imperfect calibration of the PSF. However, the structures (clusters of several points above the best fit models) seen in the radial profiles of DQ Her is unusual.

The two features marked with horizontal bars (one at 5–9 arcsec in the 0.2–0.5 keV band, and another at 8–10 arcsec in the 0.5–0.8 keV band) contain 18.0$`\pm `$4.6 and 12.2$`\pm `$4.1 counts, respectively. Other possible features either have too few photons (the 6–7 arcsec bin in the 0.5–0.8 keV profile has 5) or outside the optically detected nova shell (12–20 arcsec excess in the 0.2–0.5 keV profile) and will not be discussed further. For the two features, we also construct an azimuthal histogram in 30 degree bins of photons falling within the appropriate radial and energy bins. We plot this in the insets of Figure 3 as a histogram in polar coordinates (the position angle on the graph corresponds to the position angle on the sky, and radial distance from the center is proportional to the number of photons in that bins). We find that the 0.2–0.5 keV feature is about 8 arcsec away S – SW of DQ Her, while the 0.5–0.8 keV feature is about 9 arcsec away NE of the binary.

## 4 Discussion

### 4.1 Phase Offset

The folded X-ray light curve appears to be offset from the ephemeris prediction by about 0.03 cycles, or 500 s. The phase offset is greater than 0.006 cycles, or 100 s, at 90% confidence level. Our check of the data reduction processes (such as the barycentric correction) has not revealed any potential systematics that could account for such a large offset. The linear ephemeris, taken from Zhang et al. (1995), has a nominal error of 11 s for this epoch. Moreover, their linear and quadratic ephemerides have diverged by less than 15 s. Even though a long-term sinusoidal modulation in O$``$C eclipse timings has been found (Patterson et al., 1978), the magnitude of this effect is 70 s. Thus, we believe that the phase offset likely does not result from an inaccurate ephemeris.

Instead, we consider it likely that the phase shift of the X-ray light curve is caused by departure from axial symmetry of the DQ Her system. We see a considerable scatter in the times of the optical eclipse, up to 170 s in the O$``$C timings, although some of it can be modeled by the sinusoidal term (see Figure 4 of Zhang et al. 1995). Moreover, the eclipse light curves (see their Figure 1) are obviously skewed, with ingress beginning around phase $``$0.08 and egress completing around phase +0.12. This clearly demonstrates the departure from axisymmetry of the eclipsed object. We might still argue that the deepest point in the eclipse light curve marks the position of the unseen white dwarf. However, this is not necessarily the case, particularly if the eclipsed object is the back side of an asymmetric accretion disk (Petterson, 1980). Interestingly, the radial velocity study of the secondary, using the Na I infrared doublet, shows that the inferior conjunction of the secondary is at phase 0.013$`\pm `$0.009 on the eclipse ephemeris (Horne et al., 1993). This possible phase offset between the radial velocity curve of the secondary and the optical eclipse should be studied further. In any case, a non-axisymmetric disk can result in a phase offset between the optical and X-ray eclipses.

Our best-fit X-ray eclipse width is also larger than the optical eclipse width, although in this case they can be equal within 90% confidence limits. We do not necessarily expect the two to be equal, since the model of Petterson (1980) for the phase shift of the 71-s oscillations suggests that the emission-weighted center of eclipsed light is somewhere on the back half of the accretion disk, not centered on the white dwarf. If true, the inclination angle $`i`$ derived with the latter assumption (Horne et al., 1993) would be an underestimate for a given $`q`$. We cannot solve this difficulty with the data in hand, so we treat this as an extra source of uncertainty in the analysis to follow.

### 4.2 Wind-Scattered X-rays

At this inclination, a partial eclipse is a signature of an extended emission region. Any X-ray emission or scattering from the surface of the accretion disk should be as deeply eclipsed as the optical light from the disk. Since the X-ray eclipse is much shallower than in the optical, we conclude that the X-rays must originate in a vertically extended region. Assuming the best-fit system parameters of Horne et al. (1993), we estimate that materials that are $`>0.28a`$ (where $`a`$ is the binary separation) above the orbital plane should remain uneclipsed. The limit could be somewhat larger if the inclination is underestimated by Horne et al. (1993).

Is this a single, extended X-ray source, or can there be an eclipsed component and an uneclipsed component? From the data in hand, we cannot rule out the latter possibility. We could test this if we had a higher quality X-ray light curve of the eclipse. A two-component model would predict a sharp ingress and egress, while a single extended source would result in a more gradual eclipse transitions. Here we proceed assuming a single extended source, since we do not have a physical model corresponding to the two-component model.

The partial eclipse of UV emission lines in DQ Her also demonstrates the presence of an extended emission region (Córdova & Mason, 1985; Silber et al., 1996; Eracleous et al., 1998). The last authors, in particular, have analyzed the time-resolved profile of the C IV $`\lambda `$1550 line through the eclipse. At mid-eclipse, this line is clearly redshifted, which they argue is a strong indication that the residual line flux originates in an accretion disk wind (see also Figure 7. of Córdova & Mason 1985), even though an accretion disk wind is usually seen only in non-magnetic CVs with high accretion rate (see, e.g., Drew & Proga 2000). These systems have an accretion disk that extends all the way down to the white dwarf surface and an optically thick boundary layer. The presence of wind in DQ Her, a system with no boundary layer and a luminous but truncated disk, suggests that the accretion geometry does not dictate whether a CV has a wind or not (Eracleous et al., 1998). DQ Her appears to be sufficient luminous to drive a wind, as the fluxes of the He II recombination lines in DQ Her requires an ionizing luminosity of order $`10^{34}`$ ergs s-1 (Silber et al., 1996; Eracleous et al., 1998).

The presence of the accretion disk wind in DQ Her might make this an unusual IP, since we are not aware of any reports of accretion disk winds in other IPs. Instead, it suggests a natural parallel with the dwarf nova OY Car in superoutburst and the nova like system UX UMa. Both are bright non-magnetic CV with accretion disk wind; both are deeply eclipsing in the optical yet no X-ray eclipses were observed in OY Car in supreoutburst (Naylor et al., 1988; Pratt et al., 1999) or in UX UMa (Wood et al., 1995). This is interpreted to imply that the central X-ray source in these systems cannot be seen directly, but instead is seen indirectly after scattering in an extended corona or wind. The EUVE observation of OY Car in superoutburst (Mauche & Raymond, 2001) confirms this basic picture, since the observed spectrum can be modeled as resonant scattering of continuum photons in an accretion disk wind. A wind, rather than a static corona, is strongly favored since the EUV emission lines are seen to have FWHM of $``$2300 km s-1. We believe the same resonant scattering model can be applied to the observed soft X-rays from DQ Her, except that the nature of the unseen central source must be different (see §4.3).

A static corona can also be the scattering site for the X-ray photons observed in DQ Her. However, that would make it different from OY Car, where the EUV line widths prove the wind to be responsible for scattering. Secondly, the UV line profile during the eclipse is indicative of a wind, while we do not know if a static corona exists in this system. Finally, the depth of the partial eclipse demands a large vertical extent ($`>0.28a`$, see above) for the scattering region. This is well away from the disk and probably is an unlikely location for a static corona.

We therefore interpret the observed soft X-rays as due to resonant scattering of an unseen central source. Furthermore, we prefer the accretion disk wind as the scattering site, pending final confirmation in the form of X-ray line widths. It is interesting to note that a partial X-ray eclipse is detected in DQ Her, while none has been found in OY Car or UX UMa. It is not clear, however, whether this indicates some difference in the wind geometry compared with between DQ Her and non-magnetic systems, or if this is simply due to the more extreme inclination angle in DQ Her.

### 4.3 Non-Detection of Reflected X-rays

The observed X-ray luminosity of 3–5 $`\times 10^{30}`$ ergs s-1 is likely to be a small fraction of the intrinsic luminosity of the central source, although this depends strongly on the unknown efficiency of scattering in the wind. The flux of the He II $`\lambda 1640`$ line has been used to estimate an ionizing luminosity of 1–2 $`\times 10^{34}`$ ergs s-1 (Silber et al., 1996; Eracleous et al., 1998).

Can this be in the form of a hard component (kT$``$20 keV bremsstrahlung) usually seen in IPs? To investigate this question, we have modeled the reflection of such a component off the surface of the Roche-lobe filling secondary, applying the method of Still et al. (2001) developed for Hercules X-1. X-rays hitting the cold material on the surface of the secondary can be scattered or absorbed; since softer photons are more likely to be absorbed, the scattered component is harder than the incident component. We have evaluated the 5–8 keV flux expected from reflection from the secondary as a function of the orbital phase, taking into account the shadowing by the accretion disk and the Roche-lobe filling shape of the secondary, assuming the system parameters of Horne et al. (1993), as listed in §1. For the computational details, see Still et al. (2001).

The reflection component should be most prominent at the highest energy band of Chandra, 5–8 keV. The viewing geometry for the reflection component is most favorable at orbital phase 0.5; it would result in a strong, sinusoidal orbital modulation. However, no such component is detected in our Chandra data in the 5–8 keV data. This can be converted to the upper limit of the bolometric luminosity of the hard component as a function of the disk thickness (Figure 4). Taking our upper limit in the phase bin around phase 0.5, 1/16th cycle wide, to be 7.0$`\times 10^4`$ cts s-1, we estimate an upper limit to the kT=20 keV bremsstrahlung luminosity of 1033 ergs s-1 for a disk $`H/R`$ (height over radius) of less than 0.23, or half-opening angle $`<13^{}.0`$ ($`<14^{}.6`$ for a kT=50 keV bremsstrahlung due to the different bolometric correction). If the intrinsic hard component luminosity is 1034 ergs s-1, the reflected luminosity can be reduced to an undetectable level if the disk half opening angle is $`>17^{}.7`$ for kT=20 keV ($`>18^{}.0`$ for kT=50 keV), almost completely shadowing the secondary.

In comparison, the eclipsing and dipping low-mass X-ray binary EXO 0748$``$676 has an inclination of $`i75^{}`$ (Parmar et al., 1986): The presence of X-ray eclipses require the opening angle of the disk to be less than 15 toward the secondary in this system. Although azimuthal structures at such heights are commonly inferred (including one in U Gem 25 above the disk; Szkody et al. 1996), we are not aware of claims for disks thicker than 15 in the direction toward the secondary. We therefore conclude that the central hard X-ray component is less luminous than inferred from the He II $`\lambda 1640`$ line, unless the disk in DQ Her has an unusually large thickness. Further studies of the disk thickness, however, are required to tighten the upper limit of the hard component luminosity.

In contrast, the lack of reflection component does not allow us to put a useful constraint on the central luminosity, if it has a soft spectrum. Accretion luminosity in DQ Her would be dominated by a soft, blackbody-like component, if the shock is buried within the white dwarf atmosphere. In addition, the entire photosphere of the primary may still be hot from the thermonuclear runaway of 1934. The soft component in AM Her type magnetic CVs have a typical temperature of 20–30 eV (Mauche, 1999); even using a kT=100 eV blackbody, the spectrum is too soft to result in appreciable reflection signature. Therefore, our results are consistent with the ionizing luminosity of 1–2 $`\times 10^{34}`$ ergs s-1 (Silber et al., 1996; Eracleous et al., 1998) if this is in the form of a soft component, regardless of the disk height.

### 4.4 The Nova Shell around DQ Her

We appear to have detected extended emission features, one about 8 arcsec S to SW of, and another about 9 NE of DQ Her. The latter location corresponds to a bright \[N II\] knot in the DQ Her shell, while the former has no obvious optical feature associated with it (Slavin et al., 1995). We consider it likely that these small excesses are real and represent X-ray knots in the DQ Her shell. A similar, but much brighter, X-ray shell has been seen around GK Per (Nova Persei 1901; Balman & Ögelman 1999). This shell has a soft thermal spectrum, with a prominent Ne IX line at 0.9 keV (Balman, 2002). We do not have enough photons for a spectral analysis for the DQ Her shell, but it appears much softer than the GK Per shell. This could result from a different speed class (DQ Her was a slow nova, while GK Per was a fast nova), as well as from a different environment around DQ Her than around GK Per. The luminosity in the DQ Her knot appears to be a few times 1028 ergs s-1, subject to a large uncertainty in count rate to flux conversion, due to the unknown spectral shape.

Given the probable discovery of X-rays from the shell, and the apparent lack of strong hard X-rays from the central source, it may be useful to revisit the models of UV and optical spectra of the shell. For example, Ferland & Truran (1981) assumed a hard X-ray emitting magnetic CV with a luminosity in excess of 1034 ergs s-1 as the source of ionizing photons. Our work casts some doubt on such a picture. Petitjean et al. (1990), on the other hand, assumed a 105 K blackbody source, which is consistent with our X-ray data of the central binary. However, how do the shell X-rays fit in?

Balman (2002) interprets the X-rays from the GK Per shell as a miniature supernova remnant. In this model, the nova ejecta are plowing into the surrounding interstellar medium and shock heated into X-ray emitting temperatures. There certainly is ample kinetic energy in the nova ejecta to power shell X-ray emission, and the ejecta velocities are fast enough to result in X-ray emitting temperatures. The photoionization model of UV and optical lines, on the other hand, suggest a cold shell. We regard this not as an outright contradiction, but as an additional evidence for the co-existence of hot and cold matter in nova shells. We have long known of spectra that show both hot and cold components in nova shells (Payne-Gaposchkin, 1957; Williams, 1992). More recent imaging results (Slavin et al., 1995) and theoretical work (Lloyd et al., 1997) show the importance of Rayleigh-Taylor instability. Nova shells are complex, multi-phased entities, including, as it now appears, X-ray emitting knots.

## 5 Conclusions

Our Chandra observation of the old nova DQ Her have revealed a partial X-ray eclipse. We believe that the most likely explanation is scattering in the accretion disk wind, with a central X-ray source which is hidden from our view by the accretion disk at all orbital phases. In this model, DQ Her is similar to OY Car in superoutburst and to UX UMa. With a higher signal-to-noise observation, such as we can hope to achieve with XMM-Newton, we should be able to measure the X-ray eclipse profile in detail, and check the reality of the phase offset between X-ray and optical eclipses. We also detect an apparent weak extended feature, which we interpret as an X-ray emitting knot in the nova shell. This would make DQ Her the second classical nova, after GK Per, with an X-ray emitting shell.


format: text

text:
# New Pseudo-Phase Structure for 𝛼-Pu

## Abstract

In this paper we propose a new pseudo-phase crystal structure, based on an orthorhombic distortion of the diamond structure, for the ground-state $`\alpha `$-phase of plutonium. Electronic-structure calculations in the generalized-gradient approximation give approximately the same total energy for the two structures. Interestingly, our new pseudo-phase structure is the same as the Pu $`\gamma `$-phase structure except with very different b/a and c/a ratios. We show how the contraction relative to the $`\gamma `$ phase, principally in the $`z`$ direction, leads to an $`\alpha `$-like structure in the  plane. This is an important link between two complex structures of plutonium and opens new possibilities for exploring the very rich phase diagram of Pu through theoretical calculations.

Pure plutonium (which also includes Pu with small amounts of other solute atoms) is an exotic metal that is poorly understood. It appears to be going through a Mott-insulator transition in its $`5f`$ electronic orbitals as it transforms into different crystal structures as a function of temperature. Actinide elements lighter than Pu have itinerant $`5f`$ electrons; elements heavier than Pu form a second rare-earth series with localized $`5f`$ electrons. The $`\alpha `$-Pu phase (zero-temperature ground state) is believed to have mainly itinerant $`5f`$ electrons. However, $`\delta `$-phase (fcc) Pu is believed to be in neither limit (neither completely itinerant nor localized), but is instead in some kind of exotic electronic state that we may describe as partially localized, for want of any better terminology to apply to this systemAlbers (2001).

The understanding of this phase transformation, with its associated large volume per atom expansion, has been an exciting topic of discussion within band-structure theory and related correlated-electron extensions to it for several decades now. But almost all of the more sophisticated recent calculations have been done for the fcc $`\delta `$ phaseSavrasov and Kotliar (2000); Savrasov et al. (2001); Bouchet et al. (2000); Söderlind et al. (2002). This crystal structure has one atom per unit cell and is amenable to better, more computationally intensive techniques for handling the correlated electron physics, which is believed to be the fundamental issue controlling the $`f`$-shell Mott-insulator transformation. Because the physics is so difficult, we believe that studying the strongly correlated $`\delta `$ phase in isolation is less useful than also including the more weakly correlated $`\alpha `$ phase in the study. A correct theory should get both limits right, and is a strong constraint on the theoretical approaches that are used. The chief problem with including the $`\alpha `$ structure is its complicated crystal structure, which is a low symmetry monoclinic structure with space group P21/m and involves 16 atoms per unit cellZachariasen (1963); Espinosa (2001). Because of this complexity, for example, a fully relaxed structure (volume and internal parameters) has not yet been obtained with DFT.

We believe that we have now found a solution to this problem. In this paper we propose a new pseudo structure for $`\alpha `$ Pu that involves only two atoms per unit cell and is a relatively symmetrical orthorhombic crystal structure. We show that the real $`\alpha `$ crystal structure only involves small distortions of various atoms around this pseudo structure and that the electronic-structure energy of both systems (at least at the band-structure level) is nearly identical.

By a pseudo structure we mean a simpler crystal structure where the atoms are quite close to the positions of the more complex real structure. Pseudo-structures are interesting for two reasons: (1) they provide insight into the geometrical structure of the atoms in the original structure and possible phase-transformation pathways, and (2) they provide an alternative less expensive crystal structure for complicated electronic and microstructural calculations. From a physics point of view, one can view the pseudo structure as an approximation to the real structure that captures the dominant total energy of the electronic-structure (heat of formation) and local bonding effects; the real structure can then be viewed as just a minor distortion of the pseudo structure, with only a small energy difference between the two structures.

In the case of Pu, insight into the crystal structure is extremely important. Unlike, for example, the cubic diamond structure of Si where it is easy to understand the tetrahedral nearest-neighbor structure as a consequence of directional s-p bonds, or the close packed fcc (face-centered-cubic) crystal structures favored by Madelung energies in metals, it is extremely difficult to understand why the lowest temperature Pu phases are so extremely anisotropic and complex. In $`\alpha `$ Pu, for example, the nearest-neighbor structure involves a complicated arrangement of short and long bondsDonohue (1974). How can we understand this, and why do the crystal structures become more isotropic and metallic like at high temperatures? Good pseudo structures give us the opportunity to look at these crystal structures in new ways that may eventually stimulate insights that will resolve some of the questions.

Understanding the complex phase diagram of Pu is likely to involve consideration of strong directional bonding and high density of states effects of the $`f`$ electrons and how that bonding weakens under strong electron-electron correlations. Sophisticated electronic-structure calculations are necessary to make progress in this area. Unfortunately, the $`\alpha `$ and $`\beta `$ phases, with 16 and 34 atoms in the unit cell, respectively, are just too computationally expensive for such calculations to be practical. To have a simple pseudo-phase crystal structure that can mimic the actual bonding arrangement of the extremely complicated crystal structure of the real phase makes it possible to do these very sophisticated electronic-structure calculations with little loss of accuracy or insight. If accurate pseudo phases exist, they also make possible complicated microstructural models of the phase transformations that captures the dominant components of free energy differences without requiring huge numbers of parameters. For example, Ginzburg-Landau models become much more accessible through the use of pseudo phases. An attempt is currently being made to do this for the very important $`\alpha `$ to $`\delta `$ phase transformation of PuSaxena .

What is even more interesting is that our proposed pseudo structure is identical with the $`\gamma `$ phase of Pu, except that the unit cell dimensions are greatly changed. This crystal structure is an orthorhombic distortion of the cubic diamond lattice. It has two atoms in the smallest unit cell, one at the origin, and the second 1/4 of the way along the a, b, and c axes (which are now all different in the orthorhomic unit cell). Because both atoms have a face-centered configuration, there are 8 atoms in the primitive orthorhombic unit cell (just like cubic diamond). As we will show, our pseudo structure is far better than the only other previously proposed pseudo structure for $`\alpha `$ Pu, the Crocker pseudo structureCrocker (1971). This structure was based on the insight that the the $`\alpha `$ structure can be viewed as two repeating planes of a distorted hexagonal structure. Although this structure also has two atoms in its unit cell (each of the distorted hexagons are replaced by symmetric hexagons), it is also monoclinic, and has not proved useful for electronic structure calculations. We will also show that its near-neighbor shells are a poor approximation to the real $`\alpha `$ phase.

What makes Pu fascinating from an electronic-structure point of view is the instability in the $`5f`$ shell, which have itinerant electrons for $`\alpha `$ Pu, and localize into the core for Am and higher-Z actinides. This Mott-like transition for the $`5f`$ electrons involves very strong electron-electron correlations, and requires a much more complicated theory for the electrons than the usual local density approximation (LDA) for density functional theories (DFT), which are like a mean-field treatment of the correlationsLundqvist and March (1983). If the parabolic decrease of the volume of the light actinides as a function of atomic number can be explained with a simple Friedel model, plutonium deviates from this theoretical estimation, since Pu has a larger volume than Np. Until now there is no definitive explanation of this upturn. In addition, it is well known that DFT in its local density or generalized gradient approximation (GGA), due to the presence of strong correlation effects, has been unable to describe the $`\delta `$-phase of plutonium, underestimating its atomic volume by about 20%, for example. But recently, several attemps to go beyond LDA, have given some understanding of the $`\alpha `$-$`\delta `$ transition. The LDA+$`U`$ methodBouchet et al. (2000); Savrasov and Kotliar (2000) and the LDA+DMFT (dynamical mean-field theory)Savrasov et al. (2001), by taking into account the localization of the $`f`$-states, have shown the role of spin and orbital polarization in describing the $`\delta `$-phase.

Söderlind et alSöderlind et al. (1997) have performed ab initio relativistic total energy calculations and offered an explanation for the crystal structure of $`\alpha `$-Pu in terms of a Peierls-Jahn-Teller-like distortion. Comparing the total energy of the $`\alpha `$ phase with total energies of symmetric (fcc, bcc, hcp) structures and other actinides-like structures (such as $`\alpha `$-U and $`\alpha `$-Np) they showed that the $`\alpha `$-Pu structure minimizes the total energy for plutonium, but they couldn’t explain the upturn in equilibrium volume for Pu. More recently Jones et alJones et al. (2000) have performed calculations with two different methods, FLAPW and LCGTO-FF for all the light actinides. Although their theoretical volumes agree well (within about 2%) with the experimental volumes for Th to Np, they show a larger discrepancy (7%) for $`\alpha `$-Pu, and no upturn between Np and Pu. No attempt was made to relax the structural parameters of the $`\alpha `$ phase. They also speculated that some localization effects could already be present in $`\alpha `$-Pu accounting for the discrepancy between the computed and measured atomic volume and causing the upturn. To our knowledge, no fully relaxed electronic structure calculations have been performed on the $`\alpha `$, $`\beta `$ and the $`\gamma `$ phase and a complete and comprehensible picture of the plutonium phase diagram still does not exist.

In our calculations, we have mainly used the projector augmented-wave (PAW) method, which combines ideas from pseudopotentials and the linearized augmented-plane-wave (LAPW) method Blaha et al. (1999). Forces and the full stress tensor can be easily calculated and used to relax atoms into their instantaneous ground state. This is implemented in the plane-wave basis VASP codeKresse and Furthmuller (1996) with spin-polarized DFT and the GGA of Perdew and WangPerdew and Wang (1992). Ultrasoft Vanderbilt pseudopotentials were used to represent the core electrons.

This allowed us to perform for the first time a fully relaxed calculation of the $`\alpha `$-phase. Results are presented in Table 1. As in previous calculationsSöderlind et al. (1997); Jones et al. (2000), the atomic volume is slightly lower compared to the experimental value. This discrepancy is partially due to the absence of spin-orbit coupling in our calculationsKollar et al. (1997); Jones et al. (2000). The relaxed structure is very closed to the experimental one (the energy difference is about 4 mRy), and the positions of the atoms relaxed by only 2 or 3 %, showing that DFT-GGA reproduces the ground state of plutonium with a good precision.

We have compared the total energies of the $`\alpha `$ structure, Crocker’s structure, and the $`\gamma `$ structure. As expected, the $`\alpha `$ phase has the lowest energy, while Crocker’s structure is about 20 mRy larger, and the $`\gamma `$ phase 36 mRy larger. The theoretical volume of the $`\gamma `$ phase is very far from the experimental volume (23.14 Å3), but very close to the experimental volume of the $`\alpha `$ phase, see Table 1. If we minimize the total energy of this structure as a function of $`b/a`$ and $`c/a`$ we obtain a large contraction in the $`z`$ direction and a smaller one in the $`y`$ direction to obtain a structure very different from $`\gamma `$ (we will call this new structure “pseudo-$`\alpha `$” in the rest of the paper). The volume increases slightly and is roughly equal to the theoretical volume for the $`\alpha `$ structure, see Table 1. Surprisingly the total energy is also very close to the total energy obtained for the $`\alpha `$ structure, about 3 mRy higher, almost within the accuracy of calculations. Using a full-potential linear augmented-plane wave method (FLAPW) we have checked the influence of the spin-orbit coupling on our results (see Table 1). The total energy of the pseudo-structure is still very close to the $`\alpha `$ one and the spin-orbit coupling doesn’t help to differentiate the two structures.

It is interesting to speculate on why we find nearly the same total energy for two structures that appear to be so completely different. To understand this it is useful to compare the interatomic distances in the $`\gamma `$, the pseudo-$`\alpha `$, and the $`\alpha `$ structure, which are shown in Table 2. Although the atomic environment around each atom in the $`\alpha `$ structure is very complex, we can divided it into roughly two categories: short and long bonds, well separated. Each atom has four first-nearest neighbors and ten second-nearest neighbors, except atom I which has 3 and 13 and atom VIII which has 5 and 7. During the transformation from $`\gamma `$ to pseudo-$`\alpha `$ the first shell of neighbors contracts around the atom while the second and the third shell expand to be very close to the fourth shell. This gives a local environment of 4 nearest neighbors and 10 second-nearest neighbors for the pseudo-$`\alpha `$ , exactly what is found in average in the $`\alpha `$ structure. In addition, the interatomic distances are very close. So the relaxation of the $`\gamma `$ phase gives us a new structure, with nearly exactly the same distribution of interatomic distances than in the $`\alpha `$ phase. Crocker’s structure, however, is quite different with only two first-nearest neighbors, and an intermediate shell with four atoms. So, the local atomic environment is more important than the monoclinic symmetry in determining the energy of the $`\alpha `$ structure.

In Fig. 1, we show the density of states (DOS) obtained for the three structures at their theoretical volumes. They are very similar, with the Fermi level in a pseudo gap and two broad features, one above and one below the Fermi level. The DOS of the $`\gamma `$-Pu with the experimental parameters has been compared to the pseudo-$`\alpha `$ structure and we can see the presence of a peak at the Fermi energy, so during the contraction of the cell along the $`b`$ and $`c`$ axis there is the opening of a large pseudo gap with the reduction of the density of states at the Fermi level.

Looking at the  plane of the pseudo-$`\alpha `$ structure, see Fig. 2, one can almost visualize the $`\alpha `$ phase. The atoms in dark and grey are the two distinct atoms in the pseudo-$`\alpha `$ structure, and we have superimposed in blue the $`\alpha `$-structure atoms. They are very similar despite the complexity of the $`\alpha `$ pattern and the relative simplicity of the distorted diamond structure. We have calculated the distance between related atoms in the two different crystal structures, which has an average of about 0.25 Å. These distortions are responsible for the small total-energy difference between the two structures. In addition, the  planes have the same interlayer distance as in the $`\alpha `$ phase (2.31 Å and 2.26 Å for pseudo-$`\alpha `$ and $`\alpha `$ respectively). As we indicated earlier the new structure respects the average surrounding of an atom in the $`\alpha `$ phase, with four nearest neighbours, but exactly as in the $`\alpha `$ phase, there is two in the same plane, one in the plane above and one in the plane below. The angle between one atom and its nearest neighbours from the two other planes is also very similar (around 137 and $`130^{}`$ for pseudo-$`\alpha `$ and $`\alpha `$ respectively, in diamond it is 109). In addition, the connection of our pseudo phase with the gamma phase suggests a deep structural relationship between the low-temperature phases of Pu that needs further examination. It is interesting to observe that the b/a and c/a ratios of the real $`\gamma `$ phase of Pu have a tendency to move toward the values for our new pseudo phase as the temperature is lowered. Zachariasen and Ellinger measured the lattice constants in function of temperature. The lengths of $`b`$ and $`c`$ both decrease with temperature, while the length of $`a`$ increasesDonohue (1974). If the $`\beta `$ phase did not interrupt this trend, it is likely that $`\gamma `$ Pu would slowly continue in this way until it would reach the phase transformation into the $`\alpha `$ phase that would then require just a tiny distortion.

In conclusion we have discovered a very simple pseudo structure reproducing the main features of the complex ground state of plutonium. Our results show that the local atomic environment in the $`\alpha `$ structure with four first-nearest neighbors and 10 second-nearest neighbors is its crucial characteristic. We have found a very simple path (following the experimental variation of the lattice constants of $`\gamma `$ Pu with temperature) between $`\gamma `$ and $`\alpha `$ plutonium. These results show the predictive power of DFT calculations and our new pseudo structure of the $`\alpha `$ phase opens new possibilities for exploring the behavior of plutonium. It will allow approximate calculations of the $`\alpha `$ phase with far more sophisticated techniques than have been previously possible.


text:
# The core/envelope symmetry in pulsating stars

## 1 Astrophysical context

Asteroseismology is the study of the internal structure of stars through their pulsation frequencies. The distribution of the observed pulsation periods is theoretically determined entirely by the run of two fundamental frequencies in stellar models: the buoyancy, or Brunt-Väisälä frequency and the acoustic, or Lamb, frequency. For nonradial $`g`$-mode pulsations in white dwarfs – the only modes so far observed in these stars – the buoyancy frequency is the dominant physical quantity that determines the pulsation periods, reflecting the internal thermal and mechanical structure of the star. A considerable body of work by numerous investigators has focussed on the signature that core properties, such as the crystallized mass fraction and the C/O abundance distributions, and envelope properties, such as diffusion profiles and surface layer masses, have on the pulsation frequencies. Unambiguous interpretation of the distribution of periods is hindered by the difficulty of disentangling the signatures of structures in the deep core from those in the envelope.

It is well known both analytically and numerically that the separation between the periods of consecutive radial overtones approaches a constant value for high-overtone modes, so one might think that these periods do not actually contain information about the internal structure of the star. The resolution of this apparent paradox is that in general there are sharp features in the Brunt-Väisälä frequency for which the modes may not be considered to be in the asymptotic limit, and the presence of these features perturbs the periods so that the period spacings are no longer uniform, an effect which is commonly referred to as ‘mode trapping’. Thus, these deviations from uniform period spacing contain information about sharp features in the Brunt-Väisälä frequency, and can be used to discern internal structure in the star, such as the locations of composition transition zones. To date, almost all analyses of observed mode trapping in white dwarf pulsators have focussed on determining the thicknesses of the various chemically pure envelope layers.

In particular, initial estimates of the He layer mass of the DBV GD 358 yielded $`10^{5.7}M_{}`$ (Bradley & Winget, 1994), while more recent analyses (Metcalfe et al., 2000; Metcalfe et al., 2001) indicated a globally optimal fit with $`M_{\mathrm{He}}=10^{2.7}M_{}`$ and a changing C/O profile at $`M_r0.5`$$`0.9M_{}`$. However, the Metcalfe et al. analyses still found a local minimum in parameter space near $`M_{\mathrm{He}}10^6M_{}`$. Most recently, Fontaine & Brassard (2002) calculated a grid of carbon core models that included a double-layered envelope structure, a result similar to the earlier time-dependent diffusion calculations of Dehner & Kawaler (1995). They were able to fit the observed periods of GD 358 down to a level of precision comparable to the fit of Metcalfe et al. (2001), without including any structure in the core. Their model had composition gradients at two locations in the envelope, near $`10^3`$ and $`10^6M_{}`$, and a chemically uniform core.

In this paper, we show that there is an inherent symmetry in the way in which high-overtone pulsations sample the cores and envelopes of the models, with the result that it may be possible to reconcile the differences between the fits mentioned above. Specifically, by expanding on the work of Montgomery (2003), we show that a sharp feature in the Brunt-Väisälä frequency in the core of the models can produce qualitatively the same mode trapping patterns as a corresponding feature in the Brunt-Väisälä frequency in the envelope, leading to a potential ambiguity in the interpretation of the observed periods. In the following sections, we demonstrate how this core/envelope mapping works for both DBV and DAV models, and we use it to suggest a re-interpretation of the previous and current fits for GD 358. As a further example, we show the results of cross-fitting between two structurally dissimilar models. Finally, we show what considerations allow this ‘symmetry’ to be broken.

## 2 Simple analogy: the vibrating string

In order to better illustrate the problem of $`g`$-mode oscillations in a white dwarf, we first imagine doing asteroseismology of a far simpler object, the vibrating string. As it turns out, all of the major results can be carried over to the astrophysical case.

### 2.1 The uniform string

If we take a uniform string with a constant ‘sound speed’ $`c`$, and assume a sinusoidal time dependence, then we obtain the following eigenvalue equation:

$$\frac{^2\psi }{x^2}+\frac{\omega ^2}{c^2}\psi =0,\psi (0)=0=\psi (L),$$
(1)
where $`\psi (x)`$ is the spatial part of the eigenfunction, $`\omega `$ is the eigenfrequency, and $`L`$ is the length of the string. The solution is given by

$$\psi (x)=A\mathrm{sin}(k_nx),k_n=\frac{n\pi }{L},\omega _n=k_nc,n=1,2,3,\mathrm{}$$
(2)
Thus, the spectrum of eigenfrequencies $`\{\omega _n\}`$ is equally spaced. If we now make a small position-dependent perturbation to the sound speed $`c`$, call it $`\delta c(x)`$, then the shift in frequencies may be calculated from a variational principle, yielding

$$\frac{\delta \omega _n}{\omega _n}=\frac{2}{L}_0^L𝑑x\left(\frac{\delta c(x)}{c}\right)\mathrm{sin}^2(k_nx).$$
(3)
Finally, we may ask what the perturbation to the eigenfrequencies is for a new $`\delta \stackrel{~}{c}(x)`$ which is the reflection of $`\delta c(x)`$ about the midpoint of the string, i.e., $`\delta \stackrel{~}{c}(x)=\delta c(\stackrel{~}{x})`$, where $`\stackrel{~}{x}Lx`$:

$`{\displaystyle \frac{\delta \omega _n}{\omega _n}}`$ $`=`$ $`{\displaystyle \frac{2}{L}}{\displaystyle _0^L}𝑑x\left({\displaystyle \frac{\delta c(\stackrel{~}{x})}{c}}\right)\mathrm{sin}^2(k_nx)`$ (4)
$`=`$ $`{\displaystyle \frac{2}{L}}{\displaystyle _0^L}𝑑\stackrel{~}{x}\left({\displaystyle \frac{\delta c(\stackrel{~}{x})}{c}}\right)\mathrm{sin}^2(k_nLk_n\stackrel{~}{x})`$
$`=`$ $`{\displaystyle \frac{2}{L}}{\displaystyle _0^L}𝑑\stackrel{~}{x}\left({\displaystyle \frac{\delta c(\stackrel{~}{x})}{c}}\right)\mathrm{sin}^2(k_n\stackrel{~}{x}).`$
From comparison of Eqs. 3 and 4, we see that $`\delta c(x)`$ and $`\delta c(\stackrel{~}{x})`$ produce identical sets of perturbed frequencies. Thus, from an asteroseismological standpoint, the two perturbations are degenerate and cannot be distinguished from one another. This is hardly surprising since we know that reflection about the midpoint of the string is a symmetry of the problem, i.e., it doesn’t matter which end of the string we call $`x=0`$, and how we choose it had better have no effect on the observed frequencies. Thus, the above symmetry must be present even if the perturbation $`\delta c(x)`$ is not small.

### 2.2 The nonuniform string

For the case when $`c=c(x)`$ is a function of $`x`$, Eq. 1 is still separable, but the $`x`$-dependence is no longer sinusoidal. However, if we consider only modes with relatively high overtone number ($`n1`$), then we may use the JWKB approximation to obtain

$$\psi _n(x)=A\frac{1}{\sqrt{k_n(x)}}\mathrm{sin}[\varphi _n(x)],k_n(x)=\frac{\omega _n}{c(x)},$$
(5)
where

$$\omega _n=\frac{n\pi }{_0^L𝑑xc^1(x)},\varphi _n(x)n\pi \frac{_0^x𝑑x^{}c^1(x^{})}{_0^L𝑑x^{}c^1(x^{})},n=1,2,3,\mathrm{}$$
(6)
Just as before, a small perturbation to the sound speed, $`\delta c`$, may be related to a change in the eigenfrequencies by

$$\frac{\delta \omega _n}{\omega _n}=\frac{2}{_0^L𝑑xc^1(x)}_0^L𝑑xc^1(x)\left(\frac{\delta c(x)}{c(x)}\right)\mathrm{sin}^2[\varphi _n(x)].$$
(7)
Since the string is not uniform, it no longer possesses reflection symmetry about its midpoint. However, in analogy with the uniform case, we consider a reflection coordinate $`\stackrel{~}{x}`$ defined by

$`\varphi _n(\stackrel{~}{x})`$ $``$ $`\varphi _n(L)\varphi _n(x)`$ (8)
$`=`$ $`n\pi \varphi _n(x).`$
As before, if we consider a ‘reflected’ perturbation $`\delta \stackrel{~}{c}(x)/\stackrel{~}{c}(x)=\delta c(\stackrel{~}{x})/c(\stackrel{~}{x})`$, then from Eq. 8 we have $`dxc^1(x)=d\stackrel{~}{x}c^1(\stackrel{~}{x})`$, so we find that

$`{\displaystyle \frac{\delta \omega _n}{\omega _n}}`$ $`=`$ $`{\displaystyle \frac{2}{_0^L𝑑xc^1(x)}}{\displaystyle _0^L}𝑑xc^1(x)\left({\displaystyle \frac{\delta c(\stackrel{~}{x})}{c(\stackrel{~}{x})}}\right)\mathrm{sin}^2[\varphi _n(x)]`$ (9)
$`=`$ $`{\displaystyle \frac{2}{_0^L𝑑xc^1(x)}}{\displaystyle _0^L}𝑑\stackrel{~}{x}c^1(\stackrel{~}{x})\left({\displaystyle \frac{\delta c(\stackrel{~}{x})}{c(\stackrel{~}{x})}}\right)\mathrm{sin}^2[n\pi \varphi _n(\stackrel{~}{x})]`$
$`=`$ $`{\displaystyle \frac{2}{_0^L𝑑xc^1(x)}}{\displaystyle _0^L}𝑑\stackrel{~}{x}c^1(\stackrel{~}{x})\left({\displaystyle \frac{\delta c(\stackrel{~}{x})}{c(\stackrel{~}{x})}}\right)\mathrm{sin}^2[\varphi _n(\stackrel{~}{x})],`$
which is identical to Eq. 7 if $`\stackrel{~}{x}`$ is replaced by $`x`$. Thus, the reflected perturbation produces the same set of frequency perturbations as did the original perturbation, so there again exists an ambiguity in inferring the structure of the string from its eigenfrequencies; this is illustrated in Fig. 1.

Finally, we note that this symmetry bears a close resemblance to the phenomenon of frequency aliasing in time-series analysis. If we consider the frequency perturbation to be a continuous function of the overtone number, such as the solid line in the lower panel of Fig. 1, then the fact that the spectrum of eigenfrequencies is discrete means that this curve is sampled only at an evenly spaced set of points. Now if the curve in the lower panel of Fig. 1 represents the frequency changes due to a perturbation in $`\delta c/c`$ at $`x/L=0.30`$, then the curve corresponding to the symmetric perturbation at $`x/L0.87`$ would be a highly oscillatory curve passing through the same set of integer values of $`n`$. However, since the eigenmodes sample this oscillatory curve only at integer values of $`n`$, this high frequency signal is ‘aliased’ back to lower frequencies, i.e., the solid curve in the lower panel of Fig. 1. In this way, a perturbation $`\delta c/c`$ at the point $`\stackrel{~}{x}`$ can, through aliasing, appear as if it originates at the point $`x`$.

## 3 White dwarf core/envelope symmetry

### 3.1 Analytical approach

We would now like to apply the ideas of the previous section to the pulsating white dwarf stars. We believe this should be possible since the adiabatic equations of oscillation of a spherically symmetric stellar model (in the Cowling approximation) can be reduced to a second-order differential equation, i.e., the $`g`$-mode pulsations can be reduced to a form which mimics that of the vibrating string. In particular, the oscillation equation may be written as (Deubner & Gough, 1984; Gough, 1993)

$$\frac{d^2}{dr^2}\psi (r)+K^2\psi (r)=0,K^2\frac{\omega ^2\omega _c^2}{c^2}\frac{L^2}{r^2}\left(1\frac{N^2}{\omega ^2}\right),$$
(10)
where $`N`$ is the Brunt-Väisälä frequency, $`L^2\mathrm{}(\mathrm{}+1)`$, $`c`$ is the sound speed, and $`\omega _c`$ is the acoustic cutoff frequency, which is usually negligible except near the stellar surface. For high overtone $`g`$-modes, we have $`KLN/\omega r`$, and applying the JWKB approximation we find (Gough, 1993)

$$\psi _k(r)=A\frac{1}{\sqrt{K_k(r)}}\mathrm{sin}\left[\varphi _k(r)+\frac{\pi }{4}\right],\omega _k=\frac{L}{\left(k\frac{1}{2}\right)\pi }_{r1}^{r2}𝑑r\frac{N}{r}$$
(11)
and

$$\varphi _k(r)\left(k\frac{1}{2}\right)\pi \frac{_{r1}^r𝑑r^{}|N|/r^{}}{_{r1}^{r2}𝑑r^{}|N|/r^{}},k=1,2,3,\mathrm{}$$
(12)
where $`r_1`$ and $`r_2`$ are the inner and outer turning points of the mode, respectively, and where we have switched from $`n`$ to $`k`$ to denote the radial overtone number. The major difference between this case and that of the string is that for the string the turning points were the same for every mode (i.e., the fixed endpoints), whereas for $`g`$-modes, the inner and outer turning points of the modes are weak functions of the mode frequency, so they are different for every mode. Thus, the reflection symmetry may not be as exact as it was for the string problem.

If we consider a perturbation to the Brunt-Väisälä frequency $`\delta N/N`$, and we assume that the inner and outer turning points are fixed, then the problem is very similar to that of the nonuniform string, with the perturbations to the eigenfrequencies given by

$$\frac{\delta \omega _k}{\omega _k}=\frac{2}{_{r1}^{r2}𝑑rN/r}_{r1}^{r2}𝑑r\left(\frac{\delta N}{N}\right)\frac{N}{r}\mathrm{sin}^2\left[\varphi _k(r)+\frac{\pi }{4}\right].$$
(13)
Similarly, the approximate reflection mapping between the points $`r`$ and $`\stackrel{~}{r}`$ is given by

$`\varphi _k(r)`$ $`=`$ $`\varphi _k(r_2)\varphi _k(\stackrel{~}{r})`$ (14)
$`=`$ $`\left(k{\displaystyle \frac{1}{2}}\right)\pi \varphi _k(\stackrel{~}{r}),`$
which may be explicitly written as

$$_{r1}^r𝑑r\frac{|N|}{r}=_{\stackrel{~}{r}}^{r2}𝑑r\frac{|N|}{r}.$$
(15)
Essentially, the above equation says that points which are the same number of ‘wavelengths’ (nodes in the radial eigenfunction) from the upper and lower turning points are reflections of one another.

We remark that if we take $`r_1=0`$ and $`r_2=R_{}`$ in Eq. 12, then $`\varphi _k(r)`$ becomes a monotonic increasing function between the centre and the surface, and could itself be used as a radial coordinate. In particular, if we define

$$\mathrm{\Phi }(r)\frac{_0^r𝑑r\frac{|N|}{r}}{_0^R_{}𝑑r\frac{|N|}{r}},$$
(16)
then $`\mathrm{\Phi }(0)=0`$ and $`\mathrm{\Phi }(R_{})=1`$; since $`\mathrm{\Phi }`$ is directly linked to the Brunt-Väisälä frequency, we will call it the ‘normalised buoyancy radius’. The radial coordinate $`\mathrm{\Phi }`$ has the advantage that it makes the problem look quite similar to that of the uniform string: the ‘reflection mapping’ between the points $`r`$ and $`\stackrel{~}{r}`$ is given by $`\mathrm{\Phi }(r)=1\mathrm{\Phi }(\stackrel{~}{r})`$, with the reflection point (‘centre of the string’) having $`\mathrm{\Phi }=\frac{1}{2}`$. In addition, it can be shown that in the JWKB approximation the kinetic energy density per unit $`\mathrm{\Phi }`$ is constant. Since this approximation is valid everywhere except at sharp features such as composition transition zones, a change in the kinetic energy density (as well as the ‘weight function’) plotted as a function of $`\mathrm{\Phi }`$ gives us a clear indication of the amount of mode trapping which a mode experiences. We show examples of this in section 5.

For completeness, we mention that the analysis of the preceding paragraphs may also be carried out for high-order $`p`$-modes. In this case, we have $`K\omega /c`$, with $`\mathrm{\Phi }`$ given by

$$\mathrm{\Phi }(r)\frac{_{r_1}^r𝑑rc^1}{_{r_1}^{r_2}𝑑rc^1}.$$
(17)
Since $`\mathrm{\Phi }`$ is a function of the sound speed, $`c`$, we refer to it as the ‘normalised acoustic radius’. The reflection mapping is unchanged from the previous result, namely $`\mathrm{\Phi }(r)=1\mathrm{\Phi }(\stackrel{~}{r})`$. If a range of $`\mathrm{}`$ values are observed, as is the case for the Sun, then this symmetry can be broken due to the fact that the lower turning point, $`r_1`$, is a function of $`\mathrm{}`$.

### 3.2 Full numerical treatment

In order to investigate this core/envelope symmetry in detail, we use a fiducial DB white dwarf model with the parameters $`M_{}=0.6\mathrm{M}_{}`$, $`T_{\mathrm{eff}}=24,000\mathrm{K}`$, $`M_{\mathrm{He}}=10^6M_{}`$, and we consider pulsation modes having $`\mathrm{}=1`$ with periods between 400 and 900 s. We choose the diffusion exponents for the C/He transition zone so as to make this transition zone as smooth as possible. The reason for this is that we wish the background model to be smooth so that the only bumps in the Brunt-Väisälä frequency are the ones which we put in by hand.

In Fig. 2, we show the change in mode periods caused by placing a bump in the Brunt-Väisälä frequency at a point in the core as well as those produced by a bump in the envelope: the filled circles in the upper panel correspond to a core bump at $`M_r/M_{}=0.5622`$, while the open circles are for a bump in the envelope at $`\mathrm{log}(1M_r/M_{})=5.331`$. We see that the shapes of the perturbations are qualitatively the same, and that both sets of perturbations are well fit by the same asymptotic formula (solid curve). In the lower panel, we show how the residuals change as we vary the position of the (induced) envelope bump in order to try to reproduce the period changes due to the (intrinsic) core bump, i.e., we move the envelope bump through a range of radii in order to examine whether we have found a global best fit. Clearly, not only is there one unambiguous global minimum which is an order of magnitude smaller than the other local minima, the residuals are a smooth function of the bump position, so we are justified in using a nonlinear fitting algorithm in order to find the local (and global) minimum at $`\mathrm{log}(1M_r/M_{})=5.331`$.

Using the above procedure, we can map out the corresponding pairs of core/envelope points in our equilibrium model, obtaining the solid curve in Fig. 3 (the dashed line is the result of using the analytical relation given by Eq. (15) with $`r_1=0`$ and $`r_2=R_{}`$). Thus, we see that a feature in the core at $`M_r/M_{}0.5`$ can mimic a feature in the envelope at $`\mathrm{log}(1M_r/M_{})5.5`$. This is a very significant result because we a priori expect there to be both envelope bumps due to chemical diffusion as well as core bumps due to the prior nuclear burning history of the white dwarf progenitor.

For instance, in Fig. 4 we show evolutionary C/O profiles (Salaris et al., 1997) which result from assuming either standard semiconvective mixing (solid line) or complete mixing in the overshoot region (dashed line). Since the Brunt-Väisälä frequency depends on the radial derivative of these profiles, the regions of large slope at $`M_r/M_{}0.5`$ and $`M_r/M_{}0.98`$ will produce bumps in the Brunt-Väisälä frequency at these points. For the shaded regions in this figure, we have taken the quoted ranges of the He transition zones in the envelope of the model of Fontaine & Brassard (2002) and used the reflection mapping of Fig. 3 to indicate which regions in the core correspond to these envelope ranges. Almost eerily, these ranges correspond quite closely to those in which we would expect to see structure in the core. Thus, the possibility exists that Fontaine & Brassard are fitting ‘real’ structure in the core with assumed structure in the envelope. The reverse is also possible, of course, and at the very least the potential exists for the two signatures to become entangled with one another.

### 3.3 DA models

The above analysis can be directly applied to DAV models. For our fiducial model, we take $`M_{}=0.60\mathrm{M}_{}`$, $`T_{\mathrm{eff}}=12,000\mathrm{K}`$, $`M_{\mathrm{He}}=10^3M_{}`$, $`M_\mathrm{H}=10^6M_{}`$. In order to make our background model as smooth as possible, we ignore the Ledoux term in the computation of the Brunt-Väisälä frequency; this reduces but does not eliminate the bumps due to the chemical transition zones. Finally, since these stars are cooler and therefore more degenerate, we consider longer period $`\mathrm{}=1`$ modes, of between 800 and 1300 s.

For the DAVs, we find a couple of surprises. First, as shown in Fig. 5, the mapping is no longer single-valued: there is more than one point in the envelope which is approximately mapped to a point in the core, and vice versa. In fact, there appear to be at least two and possibly more families of mappings which have approximately the same residuals. This is due to the fact that our background model already has two ‘bumps’ in it corresponding to the C/He and He/H transition zones, which our test bump is interacting with.

Second, as shown in Fig. 6, the core/envelope mapping is weighted much more toward the envelope than it is for the DBV models, i.e., a given point in the core is mapped to a point farther out in the envelope than it is in the DBV models. Physically, this is because the DAV model is cooler and has a more degenerate core, and therefore $`N^2`$ is smaller in its core. Thus, one has to move farther out from the centre in order to accumulate a given amount of phase in the sense of Eq. (15) or (16). The implications of this could be important, since a changing C/O profile in the range 0.5–0.9 $`M_{}`$ maps to an envelope point of $`10^{10}`$$`10^8`$. Thus, such a C/O profile cannot mimic a hydrogen layer mass greater than $`10^6M_{}`$, so the ambiguity in interpreting the origin of perturbations to the periods may be alleviated in the DAVs.

Finally, we wish to mention that at least one of the DAVs, BPM 37093, is theoretically predicted to have a partially crystallized core. In terms of seismology, the main effect of such a core would be to exclude the pulsations from the crystallized region (Montgomery & Winget, 1999a, b), effectively making the lower turning point, $`r_1`$, the upper boundary of the crystallized core. Thus, instead of running from 0 to 1, the horizontal axis in Fig. 6 would run from $`M_{\mathrm{Cr}}/M_{}`$, to 1, where $`M_{\mathrm{Cr}}/M_{}`$ is the crystallized mass fraction. This would tend to push the symmetry point for a given location in the core farther out into the envelope. Also, the higher mass of BPM means that its core will be more degenerate than that of the other DAVs, which will also push its symmetry mapping yet farther out into the envelope.

## 4 Cross-fitting two structural models

Empirical evidence that this symmetry in the white dwarf models may lead to some ambiguity in the interpretation of model-fits from various descriptions of the stellar interior was recently published in Fontaine & Brassard (2002) and Metcalfe (2003). To investigate the effects of this symmetry more directly, we performed several cross-fitting experiments with two different structural models – attempting to match the calculated pulsation periods from one model by using a structurally distinct model to do the fitting. In particular, we attempted to fit the carbon core double-layered envelope model periods from Table 1 of Fontaine & Brassard (2002) using the 5-parameter model of Metcalfe et al. (2001), which includes a single-layered envelope and an adjustable C/O core; we were able to find a match with root-mean-square period residuals of only $`\sigma _P=1.26`$ seconds (the optimal model parameters were identical to those found by Metcalfe (2003) for GD 358).

We may also use Fontaine & Brassard’s published results to turn the question around and see what goodness of fit they would obtain by attempting to fit the periods of our model. An upper limit on the root-mean-square residuals that may arise from such a comparison can be obtained by comparing the periods of the two optimal models for GD 358 from Metcalfe (2003) and Fontaine & Brassard (2002), leading to $`\sigma _P1.26`$ seconds. It is important to note that the double-layered envelope models may be able to match the periods from the adjustable C/O models even better than this with slightly different model parameters – a possibility that would even more clearly demonstrate a core/envelope symmetry.

Given the discussion of the previous sections, the results from the cross-fitting are not too surprising, since we are in essence fitting structure in the envelope of the Fontaine & Brassard model using structure in the core of our model located at the appropriate reflection point, and vice versa. The level of the residuals from our model-to-model comparison tells us directly what we have suspected for some time: the model-to-observation residuals for GD 358 ($`\sigma _P1`$ second) are dominated by structural uncertainties in the current generation of models, regardless of which type of model is used to do the fitting. This does not necessarily mean that the conclusions based on fitting from either of these models should be thrown out; it simply means that neither model is a complete description of the actual white dwarf stars, a statement that can hardly be considered controversial.

We can attempt to determine which of the two structural descriptions is closer to reality (or better describes the interior structure as sampled by the pulsations) by comparing the absolute level of the residuals for GD 358 from each model, corrected for the number of free parameters. The 4-parameter model of Fontaine & Brassard (2002) leads to $`\sigma _P=1.30`$ seconds when compared to the periods observed in GD 358. The Bayes Information Criterion (Koen & Laney, 2000) would lead us to expect the residuals of a 5-parameter fit to be reduced to 1.17 seconds just from the addition of an extra parameter. The 5-parameter model of Metcalfe (2003) does much better than this, with residuals of $`\sigma _P=1.05`$ seconds (an improvement equivalent to $`4\sigma _{\mathrm{obs}}`$) – suggesting that the internal C/O profile may be the more important of the two possible structures. Based on the results of our own experiments with double-layered envelope models (Metcalfe et al., in preparation) this possibility may be even more likely.

## 5 ‘Breaking’ the symmetry

In the previous sections, we have demonstrated that a core/envelope symmetry exists for high-overtone modes. Fortunately, this symmetry is approximate, and there are many ways in which it may be lifted or broken. For example, many of the DBVs appear to be higher overtone pulsators ($`k10`$), although this is not necessarily true for all of the observed modes in a given star, nor is it true of every member of the class. The DAVs as a class are not high-overtone pulsators, so this symmetry will be less of an issue for them. In addition, as shown in Fig. 6, the core/envelope mapping for DAV models is weighted more towards the envelope, so that the reflection point for many points in the core is in a region of the envelope where it may not overlap with the expected position of the chemical transition zones.

In addition, the core/envelope symmetry can be broken if we make additional assumptions. For instance, due to the different physical processes which produce them, the generic shape expected for the C/O profile in the core will be different from the expected shape of the C/He profile in the envelope. Using this information (parametrized in some form, for instance), we should be able to discern a core feature from an envelope feature. We show an example of this in Fig. 7, in which we see that the reflected bumps (lower panel) do indeed have different amplitudes and shapes from the actual bumps (middle panel) in the Brunt-Väisälä frequency \[note that we have used the buoyancy radius, $`\mathrm{\Phi }`$, as defined in Eq. 16, as our radial coordinate; along the top axis we use the more familiar $`\mathrm{log}(1M_r/M_{})`$ \].

The physical inputs used to generate Fig. 7 include a Salaris-like C/O profile (e.g., Salaris et al., 1997) and a two-tiered C/He Dehner-like envelope profile (Dehner & Kawaler, 1995). Furthermore, we have defined a ‘bump’, $`\delta N/N`$, as the fractional difference between the Brunt-Väisälä frequency calculated both with and without the Ledoux term (this term explicitly takes account of the effect which composition changes have on the value of $`N^2`$). While this yields reasonable results for the inner transition zones, for the outer C/He transition zone near $`\mathrm{log}(1M_r/M_{})6`$ this prescription would greatly overestimate the importance of this bump. Instead, for this case we have taken $`\delta N/N`$ to be the fractional difference of the actual Brunt-Väisälä frequency and a third-order polynomial which smoothly joins onto the Brunt-Väisälä frequency on each side of the transition zone. With these definitions, we see that the inner C/O transition zone should produce the most dominant mode trapping feature, and that the inner C/He transition zone should produce the next most important feature, since these bumps are both the highest and the narrowest features present. In contrast, the outer C/He transition zone should have a relatively minor effect on the mode trapping.

We illustrate the mode trapping ability of this model in Fig. 8, in which we plot, again as a function of $`\mathrm{\Phi }`$, the weight functions for $`N^2`$ (see Kawaler et al., 1985, their Eq. 8c) of the full numerical problem: the middle panel is the weight function for an $`\mathrm{}=1`$, $`k=11`$ mode, and the lower panel is that of an $`\mathrm{}=1`$, $`k=21`$ mode. In the upper panel we show the bumps corresponding to the composition transition zones, and in all panels we have shaded the regions containing the two largest bumps, since these should produce the strongest mode trapping. As we can see, this is certainly the case for the two modes we have selected: the modes show a marked change in amplitude in these shaded transition regions, but otherwise they propagate with essentially ‘constant amplitude’, i.e., their amplitudes evolve according to the JWKB formula.

So how does mode trapping help break the core/envelope symmetry? First, imagine that we wish to calculate the shift in periods due to all the bumps, not just to the two large ones we have highlighted. The $`k=21`$ mode is trapped between the two transition zones, so it will sample the bump at $`\mathrm{\Phi }0.37`$ more strongly than modes such as the $`k=11`$ mode, which is not trapped in this region. Conversely, the $`k=11`$ mode has an enhanced amplitude in the outer layers ($`\mathrm{\Phi }0.6`$), so it will sample the bump at $`\mathrm{\Phi }0.77`$ more strongly than modes such as the $`k=21`$ mode, which is not trapped in this outer region. Thus, the effect which the smaller bumps have on the frequencies is strongly influenced by the mode trapping which the larger bumps produce. As a result, it should in principle be possible to discern whether there are in fact two bumps rather than one large one, what their relative locations are (e.g., both in the core, both in the envelope, or one each in the core and envelope), and what their relative strengths are.

Additionally, the presence of modes of different $`\mathrm{}`$ may help in resolving this core/envelope symmetry, since the outer turning point is a function of $`\mathrm{}`$ for many of the modes. Also, constraints on the rotational splitting kernels (e.g., Kawaler et al., 1999) derived from rotationally split multiplets may be of help. In future calculations, we will attempt to address quantitatively the conditions under which it is possible to resolve the core/envelope symmetry, both in terms of the number of modes required as well as the physical assumptions needed regarding the shapes of the features to be resolved in the Brunt-Väisälä frequency.

## 6 Discussion & Conclusions

We have demonstrated, both analytically and numerically, that a symmetry exists which connects points in the envelope of our models with points in the core of our models. Specifically, we find that a sharp feature (‘bump’) in the Brunt-Väisälä frequency in the envelope of our models can produce the same period changes as a bump placed in the core, and we have numerically calculated this core/envelope mapping. While we have restricted ourself to the case of high-overtone white dwarfs which pulsate in $`g`$-modes, such a symmetry should be a generic feature of all high-overtone pulsators, whether they pulsate with $`g`$\- or $`p`$-modes.

The specific motivation for much of our analysis has been the well-studied DBV, GD 358. Given its mass, stellar evolution theory leads us to expect that it has a C/O core, and that in some region of the core there should be a transition from a C/O mixture to a nearly pure C composition. If this transition begins at $`0.5M_{}`$, then it will produce a bump in the Brunt-Väisälä frequency at this point. As discussed in the previous sections, such a bump could be mimicked by an envelope transition zone with a depth of $`10^{6.0}M_{}`$. Therefore, any incompleteness in the modelling of the core bump could be ‘corrected’ by a transition zone in the envelope placed at a depth of $`10^{6.0}M_{}`$. This could explain the initial He layer determination of $`10^{6.0}M_{}`$ by Bradley & Winget (1994) as well as the continuing presence of a local minimum near $`10^{6.0}M_{}`$ in the current generation of models (Metcalfe et al., 2000; Metcalfe et al., 2001).

In addition, there is reason to believe that the actual C/He profile may be two-tiered, with a transition from pure C to a C/He mixture at $`10^2M_{}`$ and another transition to pure He at $`10^6M_{}`$ (Dehner & Kawaler, 1995; Córsico et al., 2002; Fontaine & Brassard, 2002); if this is the case, it is quite possible that the mode trapping effects of a C/O transition zone in the core could become entangled with those of the outer He transition zone. While additional considerations may break this core/envelope symmetry, we need to be aware of this aspect of the problem in order to make progress in modelling these stars.

In summary, we have shown that, for moderate to high overtone pulsators, there exists a symmetry in the mode trapping produced by features in the core and those in the envelope which can lead to ambiguity in determining the location of features such as composition transition zones. This may explain the present and previous fits for GD 358, and at the very least is something which must be taken into account in future asteroseismological fits.

## ACKNOWLEDGMENTS

The authors would like to thank D. O. Gough for useful discussions, and the referee for his help comments. This research was supported by the UK Particle Physics and Astronomy Research Council, by the Smithsonian Institution through a CfA Postdoctoral Fellowship, by the National Science Foundation through grant AST-9876730, and through grant NAG5-9321 from NASA’s Applied Information Systems Research Program.


text:
# Ratcheting of granular materials

## Abstract

We investigate the quasi-static mechanical response of soils under cyclic loading using a discrete model of randomly generated convex polygons. This response exhibits a sequence of regimes, each one characterized by a linear accumulation of plastic deformation with the number of cycles. At the grain level, a quasi-periodic ratchet-like behavior is observed at the contacts, which excludes the existence of an elastic regime. The study of this slow dynamics allows to explore the role of friction in the permanent deformation of unbound granular materials supporting railroads and streets.

preprint: HEP/123-qed

A particularly intriguing phenomenon in driven systems is the so-called ratchet effect. There is already an extensive body of work on this subject, driven by the need to understand biophysical systems such as molecular motors  and certain mechanical and electrical rectifiers . The classic ratchet is a mechanical device consisting of a pawl that engages the sloping teeth of a wheel, permitting motion in one direction only. Ratchet-like motion have been proposed as a mechanism to explain the convective motion and size segregation in vibrated granular materials . The understanding of this phenomenon is crucial in the investigation on the permanent deformation in structures subjected to cyclic loading. In soils, this loading can be induced by earthquakes, sea waves, road traffic, etc. .

The classical theory of elasto-plasticity describes the cyclic loading response by postulating an elastic region in the stress space, which changes during the deformation . This elastic region, however, is not easy to identify because the onset of the plastic deformation is gradual and not sharply defined. A great variety of modifications have been proposed in order to provide a more appropriate description which, however, make the theory too complicated, and require too many material parameters which are difficult to calibrate .

This research was motivated by experiments of cyclic loading tests on unbound granular material used to support railroads and streets. A slow deformation is observed during the service life of these structures, where the breakage, corrosion and the friction between the grains play an important role . We emphasise some recent experiments performed in Darmstadt . These experiments show that when the samples consist of very wear resistant grains, the long time cyclic loading behavior is given by a linear accumulation of plastic deformation. This surprising result suggests that the grains attain a periodic irreversible motion at the sliding contacts, which could, in principle, be detected using numerical simulations. Here we report on the first micromechanical observation using molecular dynamic simulations. Ratcheting motion was detected in the sliding contacts on a polygonal packing subjected to quasi-static cyclic loading.

The polygons representing the particles in this model are generated by using a simple version of Voronoi tessellation: First, we set a random point in each cell of a regular square lattice of side $`\mathrm{}`$, then each polygon is constructed assigning to each point that part of the plane that is nearer to it than to any other point. This method gives a diversity of areas of polygons following a Gaussian distribution with mean value $`\mathrm{}^2`$ and variance of $`0.36\mathrm{}^2`$. The number of edges of the polygons is distributed between $`4`$ and $`8`$ for $`98.7\%`$ of the polygons, with a mean value of $`6`$.

The interaction between the polygons is modelled as follows: when two polygons overlap, two points can be defined by the intersection of their edges. The segment connecting these two intersection points defines the contact line. Then, the elastic part of the contact force is calculated as $`\stackrel{}{f^e}=k_n\mathrm{\Delta }x_n^e\widehat{n}^c+k_t\mathrm{\Delta }x_t^e\widehat{t}^c`$, where $`\widehat{n}^c`$ and $`\widehat{t}^c`$ denote the normal and tangential unitary vectors with respect to the contact line, and $`k_n`$ and $`k_t`$ are the stiffnesses in the respective directions. $`\mathrm{\Delta }x_n^e`$ is the overlapping area divided by a characteristic length of the interacting polygon pair . $`\mathrm{\Delta }x_t^e`$ defines the elastic tangential displacement of the contact, that is given by the time integral starting at the beginning of the contact

$$\mathrm{\Delta }x_t^e=_0^tv_t^c(t^{})\mathrm{\Theta }(\mu f_n^e|f_t^e|)𝑑t^{},$$
(1)
where $`\mathrm{\Theta }`$ is the Heaviside step function and $`\stackrel{}{v}_t^c`$ denotes the tangential component of the relative velocity $`\stackrel{}{v}^c`$ at the contact . This equation defines a limit of elasticity in the contact force. When the contact force satisfies $`f_t^e=\mu f_n^e`$ the contact slides. This produces an irreversible motion between the polygons which contributes to a plastic deformation in the assembly.

In each contact a viscous force $`\stackrel{}{f}^v=m\nu \stackrel{}{v}^c`$ is included. Here $`m`$ is the relative mass of the polygons in contact and $`\nu `$ is the coefficient of viscosity. This force reduces the acoustic waves produced during the loading. Inertia effects were excluded by using a very long time of loading. indeed, doubling it affects the strain response by less than $`5\%`$. This correspond to the quasi-static approximation.

The parameters of our model have been reduced to a minimum set of dimensionless parameters: The period of cyclic loading $`t_0`$ is chosen $`4000`$ times the characteristic period of oscillation $`t_s=\sqrt{k_n/\rho \mathrm{}^2}`$, where $`k_n`$ is the normal stiffness of the contacts and $`\rho `$ the density of the grains. The relaxation time $`t_r=1/\nu `$ is chosen $`10`$ times $`t_s`$. The tangential stiffness is $`k_t=0.33k_n`$; the initial pressure $`p_0=0.001k_n`$ and the friction coefficient $`\mu =0.25`$. The normal stiffness is chossen $`k_n=160MPa`$.

The simulations are performed on six different samples of $`400`$ polygons. First, the polygons are placed into a rectangular box such that they do not overlap each other. Next, a gravitational field is applied and the sample is allowed to consolidate. An external load is imposed by applying a force $`\sigma _1W`$ and $`\sigma _2H`$ on the horizontal and vertical walls, respectively. Here $`\sigma _1`$ and $`\sigma _2`$ are the horizontal and vertical stresses. $`W`$ and $`H`$ are the width and the height of the sample. When the velocity of the polygons vanishes gravity is switched off. A fifth-order predictor-corrector algorithm is used to solve the equations of motion.

In order to calculate the strain we select the polygons whose centers of mass are less than $`10\mathrm{}`$ from the center of the sample. Then, the strain tensor is calculated as the displacement gradient tensor averaged over the area enclosed by the initial configuration of these polygons. From the eigenvalues $`ϵ_1`$ and $`ϵ_2`$ of the symmetric part of this tensor we obtain the shear strain as $`\gamma =ϵ_1ϵ_2`$. The volume fraction is calculated as $`\mathrm{\Phi }=(V_pV_0)/V_b`$, where $`V_p`$ is the sum of the areas of the polygons, $`V_0`$ the sum of the overlapping area between them, and $`V_b`$ the area of the rectangular box. The initial volume fraction is $`\mathrm{\Phi }_0=84.49\pm 0.05\%`$.

Initially, the sample is isotropically compressed until the pressure $`p_0`$ is reached. Then, the sample is subjected to vertical load-unload cycles as $`\sigma _1=p_0+\frac{\sqrt{2}}{2}\mathrm{\Delta }\sigma [1+\mathrm{cos}(\pi t/t_0)]`$ whereas $`\sigma _2=p_0`$ is kept constant. Part (a) of Fig. 1 shows the relation between the shear stress $`q=(\sigma _1+\sigma _2)/2`$ and the shear strain $`\gamma `$ in the case of a loading amplitude $`\mathrm{\Delta }\sigma =0.6p_0`$. This relation consists of open hysteresis loops which narrows as consecutive load-unload cycles are applied. This hysteresis produces an accumulation of strain with the number of cycles which is represented by $`\gamma _N`$ in the part (b) of Fig 1. We observe that the strain response consists of short time regimes, with rapid accumulation of plastic strain, and long time ratcheting regimes, with a constant accumulation rate of plastic strain of around $`2.4\times 10^6`$ per cycle.

Part (c) of Fig. 1 shows the relation between the shear stress and the volume fraction. This consists of asymmetric compaction-dilation cycles which makes the sample to compact during the cyclic loading. This compaction is shown in part (d) of the Fig. 1. We observe a slow variation of the volume fraction during the ratcheting regime, and a rapid compaction during the the transition between two ratcheting regimes. The evolution of the volume ratio seems to be rather sensitive to the initial random structure of the polygons. Even so we found that after $`8\times 10^3`$ cycles the volume fraction still slowly increasing in all the samples, without any evidence of a saturation level.

This extremely slow dynamics in the evolution of the granular packing shows an astonishing analogy with the behavior of glassy systems . Even more surprising is that no elastic regime is detected by decreasing the amplitude of the loading cycles. Fig. 2 shows the accumulation rate of strain $`\mathrm{\Delta }\gamma /\mathrm{\Delta }N`$ for different loading amplitudes $`\mathrm{\Delta }\sigma `$. A constant accumulation of strain is observed during the cyclic loading, even when the amplitude is as slow as $`10^3`$ times the applied pressure.

The existence of these ratcheting regimes proves to be a consequence of the particular disorder of the distribution of contact forces. It is well know that the stress applied on the boundary is transmitted through an heterogeneous network of contact forces . Both tangential and normal force distributions turn out to be very broad. The part (a) of Fig. 3 shows the distribution of contact forces in a polygonal packing that has been isotropically compressed. The components of the force are limited by the condition $`|f_t|\mu f_n`$, where $`\mu `$ is the friction coefficient. Most of the contacts satisfy the elastic condition $`|f_t|<\mu f_n`$, but due to the heterogenieties, some contacts are able to reach the sliding condition $`|f_t|=\mu f_n`$ during the compression. Under arbitrarily small load stress cycles these contacts slide, given irreversible motion between the grains. Those sliding contacts are the most relevant micromechanic rearrangements during the cyclic loading. Opening and closure of contacts are quite rare events, and the coordination number of the packing keeps approximately it initial value $`4.43\pm 0.08`$ in all the simulations.

After some loading cycles a quasi-periodic irreversible motion is observed at these sliding contacts. Two different sliding modes are represented in Fig. 3: (b) a contact sliding forward in the load phase and backward in the the unload one and (c) a contact sliding in the load phase and sticking in the unload one. A measure of the plastic deformation of the sliding contact is given by $`\xi =(\mathrm{\Delta }x_t^c\mathrm{\Delta }x_t^e)/\mathrm{}`$, where $`\mathrm{\Delta }x_t^c`$ and $`\mathrm{\Delta }x_t^e`$ are the total and the elastic part of tangential displacement at the contact, the last one being given by Eq. 1. The part (d) of Fig. 3 shows the plastic deformation $`\xi `$ of the sliding contact shown in (c). Due to the load-unload asymmetry of the contact force loop, a net accumulation of plastic deformation is observed in each cycle. This is given by a slip-stick mechanism which resembles mechanical ratchets.

It is interesting to see the spatial correlation of these ratchets. Fig. 4 shows a snapshot of the field of plastic displacement per cycle at the contacts inside of the specimen. We see that correlated displacements coexist with a strongly inhomogeneous distribution of amplitudes. Localized slip zones appear periodically during each ratcheting regime. Some slip zones are destroyed and new ones are created during the transition between two ratcheting regimes. Moreover, we notice that these ratchets are found as well at the boundaries as in bulk material, without the layering effects observed in vibrated granular materials .

Let us next discuss the correlation between the dynamics of the sliding contacts and the evolution of the stiffness of the material. The last one is given by the slope of the stress strain curve in part (a) of the Fig. 1. The evolution of the fraction $`n_s=N_s/N_c`$ of sliding contacts with the number of loading cycles is shown in Fig. 5. Here $`N_s`$ is the number of sliding contacts and $`N_c`$ is the total number of contacts. During each loading phase, the number of sliding contacts increases, giving rise to a continuous decrease of the stiffness as shown in part (a) of Fig. 1. An abrupt reduction in the number of sliding contacts is observed at the transition from load to unload, producing a discontinuity in the stiffness and hence, a plastic deformation. We can see also that some contacts reach almost periodically the sliding state even for extremely small loading cycles. The ratchet-like behavior of these contacts produces a net displacement of the hysteretic stress-strain loop in each cycle, ruling out an elastic regime.

The fact that there is no evidence of an elastic regime in the cyclic loading response suggests that the hysteresis of the granular materials can not be described by using the classical theory of elasto-plasticity. Actually, the deficiencies of this theory have been addressed in the experimental investigation of soil deformation . Ten years ago, the hypoplastic model was formulated in order to mend these deficiencies. This continuum theory has been supported on the experimental evidence that any load involves plastic deformation .

The hypoplastic approach of granular materials requires to introduce new ingredients in the current micro-mechanical description. The traditional fabric tensor, measuring the distribution of the orientation of the contacts, can not fulfill this description, because it does not make a distinction between elastic and sliding contacts . New structure tensors taking into account the statistics of the sliding contacts, must be introduced in order to give a micro-mechanical basis to the granular ratcheting. The identification of these internal variables and the determination of their evolution equations and their connection with the macroscopic variables would be an important key to develop an appropriate continuous description of granular soils.

To conclude, we have performed a grain scale investigation of the cyclic loading response of a polygonal packing. We have shown the existence of long time regimes with a constant accumulation of plastic deformation per cycle, due to ratcheting motion at the sliding contacts. This phenomenon may have deep implications in the study of the permanent deformation of soils subjected to cyclic loading. More precisely, it may be necessary to introduce internal variables in the constitutive relations, connecting the dynamics of the sliding contacts with the evolution of the continuous variables during cyclic loading.

At this time, a comparison of the dynamic simulations with realistic situations is limited by the computer time of the simulations. Using a computer of $`2.4GHz`$ we are able to simulate only $`3\times 10^3`$ cycles per week. Such amount of cycles corresponds only to some hours in realistic cases such as railroads and highways . However, further investigations of the role of the stiffness, the initial stresses, and the inertial effects on this ratcheting behavior are currently feasible.

We thank G. Gudehus, P. Cundall, D. Potyondy, A. Schuenemann, J. Gallas and S. McNamara for helpful discussions and acknowledge the support of the DFG project Modellierung kohäsiver Reibungsmaterialen and the European DIGA project HPRN-CT-2002-00220.