text:
# Untitled Document

agbppr3.tex, version of 11 March 2002 - revised by JHT

FUEL-SUPPLY-LIMITED STELLAR RELAXATION OSCILLATIONS: APPLICATION TO MULTIPLE RINGS AROUND AGB STARS AND PLANETARY NEBULAE

Hugh M. Van Horn1 Current address: Department of Terrestrial Magnetism, Carnegie Institution of Washington, 5241 Broad Branch Road, Washington, DC 20015; on leave from the National Science Foundation, John H. Thomas2 Also Department of Mechanical Engineering, University of Rochester, Adam Frank, and Eric G. Blackman

Department of Physics and Astronomy

University of Rochester

Rochester, NY 14627-0171

Abstract

We describe a new mechanism for pulsations in evolved stars: relaxation oscillations driven by a coupling between the luminosity-dependent mass-loss rate and the H fuel abundance in a nuclear-burning shell. When mass loss is

included, the outward flow of matter can modulate the flow of fuel into the shell when the stellar luminosity is close to the Eddington luminosity $`L_{\mathrm{Edd}}`$. When the luminosity drops below $`L_{\mathrm{Edd}}`$, the mass outflow declines and the shell is re-supplied with fuel. This process can be repetitive. We demonstrate the existence of such oscillations and discuss the dependence of the results on the stellar parameters. In particular, we show that the oscillation period scales specifically with the mass of the H-burning relaxation shell (HBRS), defined as the part of the H-burning shell above the minimum radius at which the luminosity from below first exceeds the Eddington threshold at the onset of the mass loss phase. For a stellar mass $`M_{}0.7M_{}`$, luminosity $`L_{}10^4L_{}`$, and mass loss rate $`|\dot{M}|10^5M_{}`$ yr-1, the oscillations have a recurrence time $`1400`$ years $`57\tau _{\mathrm{fsm}}`$, where $`\tau _{\mathrm{fsm}}`$ is the timescale for modulation of the fuel supply in the HBRS by the varying mass-loss rate. This period agrees very well with the $``$ 1400-year period inferred for the spacings between the shells surrounding some planetary nebulae. We also find the half-width of the luminosity peak to be $`0.39`$ times the oscillation period; for comparison, the observational shell thickness of $``$ 1000 AU corresponds to $`0.36`$ of the spacing between pulses. We find oscillations only for models in which the luminosity of the relaxation shell is $``$ 10–15% of the total stellar luminosity and for which energy generation occurs through the $`pp`$ chain. We suggest this mechanism as a natural explanation for the circumnebular shells surrounding some planetary nebulae, which appear only at the end of the AGB phase.

1. Introduction

Asymptotic giant branch (AGB) stars exhibit pulsations on a variety of timescales, ranging from long periods associated with thermal pulses of the He-burning shells ($`\tau 10^5`$ years; cf., Iben and Renzini 1983, Schönberner 1983, Mazzitelli and D’Antona 1986, Iben 1991, Dorman et al., 1993, Vassiliadis and Wood 1993) to far shorter oscillations associated with acoustic modes of the star ($`\tau 1`$ year). Recently, a new class of pulsational behavior has been inferred for such stars. The evidence takes the form of nested circumnebular shells around a number of planetary nebulae (cf., Mauron and Huggins 2000; Terzian and Hajian 2000; Balick, Wilson, and Hajian 2001). These appear to have been ejected in pulses during the AGB phase or succeeding post-AGB but pre-planetary-nebula (pre-PN) phases of stellar evolution. For example, surrounding the Cat’s Eye Nebula (NGC 6543), Balick et al. (2001) find rings that appear to have been ejected quasi-periodically, with the estimated time between successive pulses ranging from $``$ 1000 to $``$ 1900 years; the mean spacing for eight successive shells is approximately 1400 years. The authors take the rings to be $``$ 1000 AU thick, indicating that the ejection mechanism operates during 0.36 of the time between successive pulses. Each shell is estimated to contain $`0.01M_{}`$ of material, and a total mass $`0.1M_{}`$ appears to have been ejected in the process that produced the shells.

To date, four models have been proposed to explain these enigmatic shells. Soker and Harpez (1988) suggested a binary model, in which the secondary follows a highly elliptical orbit. As the secondary approaches periastron, it interrupts the mass loss from the primary, leading to the formation of a shell. The intershell timescale is then set by the orbital period $`\tau _{\mathrm{orb}}`$ of the binary. Mastrodemos and Morris (1998, 1999) proposed a different binary model, in which the shells are the result of spiral shocks in the AGB wind. In these models the intershell timescale is also $`\tau _{\mathrm{orb}}`$. Simis, Icke, and Dominik (2001) proposed that the shells are produced by quasi-periodic variations in dust formation rates in dust-driven AGB winds. They argue that non-linear coupling between pulsation of the star, dust formation in the atmosphere, and radiation-driving in the wind produces shells with a quasi-regular period $`10^3`$ years. Finally, Garcia-Segura et al. (2002) have proposed an MHD model in which a dynamo operating in the star sets the period of the oscillations, and magnetic pressure in the wind creates the shells. It is difficult at present to evaluate the validity of these models. The eccentric binary model seems unlikely to account for the several PN discovered with multiple shells; it is not clear why a $`10^3`$ year period should be favored; and it is hard to understand how the observed variation from $``$ 1000 to 2000 years in shell spacings could be produced by a perfectly periodic binary model. The MHD model may be based on sound physics, but existing AGB dynamo models (Blackman et al. 2001) predict far shorter dynamo periods.

A number of indications seem to us to point instead to fuel-supply-limited relaxation oscillations as the explanation for the observed circumnebular mass shells:

(1) At high luminosities, stars lose mass at exceptionally high rates. Indeed, the brightest precursors of the planetary nebulae may have $`L_{}>L_{\mathrm{Edd}}`$, driving the so-called “superwind” that terminates evolution up the AGB. However, the effect of the stellar wind upon the internal structure of a star has not, to our knowledge, been included self-consistently in previous calculations. Models of the wind assume that the stellar interior is dynamically unaffected by the mass flux through the surface layers. Conversely, models of stellar evolution that include wind mass loss simply strip mass away from the stellar surface. Because the rate of mass loss increases dramatically as a star climbs the AGB (cf., Vassiliadis and Wood 1993), any effect due to the wind will be accentuated near the end of the AGB phase. As the observed circumnebular shells were evidently produced during the final stages of evolution just prior to the ejection of a planetary nebula, some effect of mass loss upon the internal stellar structure thus seems a natural candidate to explain the pulsed ejection of the circumnebular shells.

(2) Comparison of the rate at which H burning consumes mass with the rate of mass loss from the stellar surface also is suggestive. Let $`Q^{}`$ be the amount of energy released per gram of H consumed. This quantity is (cf., Clayton 1968, p. 287 ff)

$$Q^{}=\frac{(4m_\mathrm{H}m_{\mathrm{He}})c^2}{4m_\mathrm{H}}=6.398\times 10^{18}\mathrm{erg}\mathrm{g}^1.$$
$`(1)`$
If $`X`$ is the mass fraction of H in the unburned material, $`\dot{M}_s`$ the rate at which an H-burning shell consumes mass, and $`L_s`$ the resulting luminosity produced by this shell, then

$$\dot{M}_sX=\frac{L_s}{Q^{}}.$$
$`(2)`$
With $`L_s10^4L_{}`$ and $`X0.7`$, the rate at which mass must be consumed by a H-burning shell in order to provide the stellar luminosity is thus $`|\dot{M}_s|10^7M_{}`$ yr-1. For comparison, the observed rates of mass loss range from $`5\times 10^{11}<(|\dot{M}|/M_{}\mathrm{yr}^1)<2\times 10^6`$ for the central stars of the planetary nebulae (Perinotto 1989) to $`10^4M_{}`$ yr-1 for long-period variables (cf., Cox et al. 2000, p. 415). We note that the OH/IR stars which have the highest mass loss rates may have $`L_{}>L_{\mathrm{Edd}}`$. Thus, at most stages in an AGB star’s evolution, mass loss from the surface is comparable to or higher than mass consumption at the H-burning shell. Consequently, there is a competition between the need to supply fuel downward into the H-burning shell and upward into the stellar wind.

We are thus led to suggest that an interaction between radiation-driven mass loss from the stellar surface and mass consumption by H-burning in a highly evolved AGB star may produce a relaxation oscillation that ejects shells of matter quasi-periodically. Near the tip of the AGB, a star has a luminosity $`L_{}10^4L_{}`$, and it contains a dense, degenerate C/O core with mass $`0.6M_{}`$ and radius $`0.02R_{}10^9`$ cm, surrounded by a tenuous H/He envelope extending outward to radius $`R_{}4\times 10^{13}\mathrm{cm}500R_{}`$ (cf., Iben 1987). Suppose that such a star is initially in a steady state. We define the H-burning relaxation shell (HBRS) as a thin outer part of the H-burning region above the smallest radius for which the luminosity from below can exceed the Eddington luminosity. This region can be quite a small fraction of the total H-burning shell, but is the region which is important for relaxation oscillations that we now discuss.

We assume that the HBRS consumes mass at a rate $`|\dot{M}_s|`$. Suppose also that the star undergoes radiation-driven mass loss from the stellar surface at a steady rate $`|\dot{M}|`$. A perturbation that increases the temperature of the HBRS, leading to a thermonuclear runaway, greatly increases the luminosity. Because the surface mass-loss rate is proportional to the luminosity (cf., §2), the increased luminosity increases the rate of mass loss from the surface. Conservation of mass, however, requires that the mass flowing out through the surface must be replaced from below; that is, mass loss is not simply a stripping away of the outer layers of the star. Some form of gradual readjustment must occur along the AGB and into the post-AGB and pre-PN phases. The increased mass-loss rate thus at least impedes the supply of fuel to the HBRS, and it may even temporarily interrupt that supply. In any event, the increased surface mass loss has the effect of reducing the luminosity produced by the HBRS, which in turn reduces the rate of radiation-driven mass loss from the stellar surface, allowing the cycle to repeat itself. This is a previously unexplored coupling between the HBRS luminosity, the mass-loss rate, and the fuel supply available in the HBRS. Such fuel-supply-limited relaxation oscillations are familiar in other contexts, such as the pulsation of a fuel-starved flame from a guttering candle.

In the present paper, we consider in more detail the hypothesis that the source of the circumnebular shells is just such a fuel-supply-limited relaxation oscillation in a star near the end of its AGB phase. In §2, we first review briefly the relation between the rate of mass loss $`\dot{M}`$ and the stellar luminosity $`L_{}`$ for the optically thick, radiation-driven winds appropriate to high-luminosity stars. In §3, we next collect together the fundamental Eulerian equations that describe the hydrodynamic interaction between matter and radiation in a spherically symmetric star with a spherically symmetric mass outflow. In §4 we develop a simplified model that accounts qualitatively for the main features of the relaxation oscillations described above. We emphasize that this is a simplified model, intended to illustrate the essential physics of these oscillations. A more complete calculation, beyond the scope of this paper, would require a detailed stellar evolution model including self-consistent mass loss. We describe the results of numerical calculations using our simplified model in §5, and we conclude in §6 with a summary.

2. Mass Loss on the AGB

From his observations of mass loss from luminous K and M giants and supergiants, Reimers (1975a, b) proposed an empirical formula to express the dependence of the rate of mass loss upon the stellar parameters:

$$|\dot{M}|=4\times 10^{13}M_{}\mathrm{yr}^1\eta \frac{(L/L_{})(R/R_{})}{(M/M_{})},$$
$`(3)`$
where $`\eta `$ is a parameter of order unity. This expression has been used extensively, both in analyzing observations and in theoretical studies. For example, for an AGB star with $`M_{}M_{}`$, $`L_{}10^4L_{}`$, and $`R_{}500R_{}`$, this expression yields the result $`|\dot{M}|2\times 10^6M_{}`$ yr-1. However, as originally pointed out by Renzini (1981; see also Iben and Renzini 1983), this equation appears to underestimate the rate of mass loss in the very high luminosity, thermally pulsing AGB phase, requiring a postulated “superwind” to provide the necessary mass-loss rates.

Bowen (1988) and Bowen and Willson (1991) have described hydrodynamic calculations of the atmospheres of long-period, Mira variables – pulsating AGB stars of low to intermediate masses – that lead to significant rates of mass loss. In models that include dust formation, radiation pressure acting on the dust grains produces a mass loss extending up to $`\dot{M}10^6M_{}`$ yr-1 for a star with $`M_{}=1.2M_{}`$, $`R_{}=270R_{}`$, and $`L_{}=5315L_{}`$. Bowen suggests this as the “superwind” that appears to terminate AGB evolution. Willson (2000) has recently re-interpreted the empirical Reimers mass-loss rate (3) in terms of these results.

In their computations of the evolution of intermediate-mass stars, Vassiliadis and Wood (1993) included an empirical mass-loss rate that reproduces the very large values of $`\dot{M}`$ observed during the AGB phase of evolution. It reaches an essentially constant “superwind” rate that is within a factor of two of the value

$$|\dot{M}|=\frac{L_{}}{cv_\mathrm{w}}10^5M_{}\mathrm{yr}^1\times \frac{L_{}}{10^4L_{}}$$
$`(4)`$
appropriate for a radiation-driven wind, where $`v_\mathrm{w}15`$ km s-1 is the asymptotic wind speed far from the central star. For stars with $`M<2.5M_{}`$, the “superwind” mass-loss rate is achieved only during the last few thermal pulses, and the vast majority of the mass loss occurs in these episodes.

Koesterke, Dreizler, and Rauch (1998) and Koesterke and Werner (1998) measured the mass-loss rates for four PG 1159 stars – post-AGB stars that are among the immediate precursors of the white dwarfs – finding $`8<\mathrm{log}(|\dot{M}|/M_{}\mathrm{yr}^1)<7`$, consistent with the theory of radiation-driven winds.

3. Fundamental Equations

The fundamental Eulerian equations that describe the structure and time variation of a spherically symmetric star with spherically symmetric mass loss are (cf., Landau and Lifshitz 1959; Hansen and Kawaler 1994) the continuity equation,

$$\frac{\rho }{t}+\frac{1}{r^2}\frac{(r^2\rho v)}{r}=0,$$
$`(5a)`$
where $`\rho `$ is the density and $`v`$ is the radial velocity; the definition of the mass $`M_r`$ interior to radius $`r`$,

$$\frac{M_r}{r}=4\pi r^2\rho ;$$
$`(5b)`$
the equation of conservation of radial momentum,

$$\frac{P}{r}+\rho \left(\frac{v}{t}+v\frac{v}{r}\right)=\rho \frac{GM_r}{r^2},$$
$`(5c)`$
where $`P=P_g+P_r`$ is the sum of the gas pressure $`P_g=\rho k_BT/\mu H`$ and radiation pressure $`P_r=\frac{1}{3}a_rT^4`$, and where $`T`$ is the temperature; and the heat-flow equation

$$\frac{T}{r}=\{\begin{array}{cc}\frac{3}{4ac}\frac{\kappa \rho }{T^3}\frac{L_r}{4\pi r^2},\hfill & \text{in radiative equilibrium;}\hfill \\ _{\mathrm{ad}}\frac{T}{P}\frac{P}{r},\hfill & \text{in convective equilibrium.}\hfill \end{array}$$
$`(5d)`$
Here $`\kappa `$ is the opacity of the stellar matter, and $`L_r`$ is the luminosity flowing out through a sphere of radius $`r`$. Note that the first form of equation (5d) assumes that energy transport is dominated by radiation flow rather than convection, while the second form is that appropriate for convective energy transport. We also require the equation of energy conservation,

$$\frac{L_r}{r}=4\pi r^2\rho \left(ϵT\frac{S}{t}Tv\frac{S}{r}\right),$$
$`(5e)`$
where $`ϵ`$ is the thermonuclear energy generation rate, and $`S`$ is the entropy. In addition, we need the equation for the time rate of change of the H mass fraction $`X`$ due to nuclear burning in the HBRS:

$$\frac{X}{t}+v\frac{X}{r}=\frac{ϵ}{Q^{}},$$
$`(5f)`$
where $`Q^{}`$ is given by equation (1). Note that if we multiply equation (5f) by $`dM_r`$ and integrate over the HBRS, we recover equation (2), with $`\dot{M}_s4\pi r_s^2\rho _sv_s`$.

4. A Simplified Model

Integrating equation (5b) from the center of the star out to some fixed radius $`r`$ and computing the time rate of change of this quantity gives

$$\dot{M}_r\frac{dM_r}{dt}=_0^r4\pi r^2\frac{\rho }{t}𝑑r=4\pi r^2\rho vF(r,t),$$
$`(6)`$
where we use equation (5a) to eliminate $`\rho /t`$, and the last equivalence defines the local mass-flow rate $`F`$. At the stellar surface, $`\dot{M}_r\dot{M}=F(R_{}[t],t)`$, a quantity that depends only upon time. If the star experiences steady-state mass loss, then equation (5a) shows that the mass flux $`F(r,t)`$ is a constant, independent of position as well as time. In general, $`F`$ may depend upon both position and time, but if the changes in the star occur quasi-statically (as is the case in stellar evolution calculations), then we can expect the mass flux through the star to remain close to its (position-independent) value throughout the entire stellar envelope.

From a scale analysis of the momentum equation (5c), we find that the ratio of the velocity-dependent terms to the remaining terms is of order $`(v/c_s)^2`$, where $`c_s`$ is the sound speed. For a mass-loss rate of $`10^5M_{}`$ yr-1, and assuming the mass flow rate $`F4\pi r^2\rho v`$ to be independent of position in the stellar envelope, the flow velocity at the HBRS is $`v0.5`$ cm s-1. In comparison, the sound speed at the shell is $`c_s4\times 10^7`$ cm -1. Thus, at the HBRS, the velocity-dependent terms in equation (5c) are only $`10^{16}`$ times the other terms in this equation. These terms also appear to be negligible near the photosphere, and for this reason, we neglect them in the momentum equation. This approximation is equivalent to assuming that the variations of interest to us here all take place on timescales long in comparison with the hydrostatic readjustment timescale. For this reason, the hydrodynamic timescale $`\tau _{\mathrm{hyd}}R_{}/v500R_{}/15`$ km s$`{}_{}{}^{1}1`$ year plays no further role in our development of the equations that describe this model.

Dropping these $`v`$-dependent terms reduces the momentum equation to the usual equation of hydrostatic equilibrium. Separating out the radiation pressure $`P_r`$, we find that the equation of hydrostatic equilibrium for the gas pressure $`P_g`$ can be written as

$$\frac{P_g}{r}=\frac{\rho GM_r}{r^2}\left[1\frac{\kappa L_r}{4\pi cGM_r}\right].$$
$`(7)`$
In the envelope of the star we can reasonably assume $`L_rL_{}=\mathrm{constant}`$ and $`M_rM_{}=\mathrm{constant}`$. With these approximations, equations (7) and (5d) can be written in the forms

$$\frac{P_g}{r}=\frac{\rho GM_{}}{r^2}\left[1\frac{\kappa L_{}}{4\pi cGM_{}}\right]\frac{\rho GM_{}}{r^2}(1\lambda ),$$
$`(8a)`$
where the final equivalence defines the parameter $`\lambda `$, and we have assumed radiative energy transport, at least near the HBRS:

$$\frac{T}{r}=\frac{3}{4a_rc}\frac{\kappa \rho }{T^3}\frac{L_{}}{4\pi r^2}.$$
$`(8b)`$
The quantity $`\lambda `$ in equation (8a) has a simple physical interpretation – it is a dimensionless form of the stellar luminosity:

$$\lambda \frac{\kappa L_{}}{4\pi cGM_{}}\frac{L_{}}{L_{\mathrm{Edd}}}.$$
$`(9)`$
where $`L_{\mathrm{Edd}}4\pi cGM_{}/\kappa `$ is the Eddington luminosity.

The assumption that radiative transport dominates in the stellar envelope is not clearly justified a priori, as the highly distended surface layers of AGB stars are well known to contain deep convection zones. We nevertheless believe that equation (8b) is appropriate here for at least two reasons: (1) The surface convection zone in an AGB star does not extend all the way down to the HBRS. We are most concerned with conditions near this shell, and it is well known that the temperature profile below the convection zone in the surface layers of a star rapidly approaches the “radiative zero” approximation we use below. (2) It is not clear whether the relaxation oscillations we wish to study actually occur on the AGB or whether they occur at some later high-luminosity phase, where the surface convection zone may be shallower and equation (8b) correspondingly better justified.

Dividing equation (8a) by equation (8b), assuming a Kramers’ law opacity of the form The more familiar form of this equation is $`\kappa =\kappa _0\rho T^{3.5}`$. The form given in equation (10) is more convenient for our present purposes, however, and the two forms are equivalent if $`a=1`$ and $`b=1.125`$.

$$\kappa =\kappa _{}P_g^aP_r^b,$$
$`(10)`$
where $`\kappa _{}`$, $`a`$, and $`b`$ are constants, and employing the conventional definition

$$P_g/P_r=\beta /(1\beta )$$
$`(11)`$
enables us to write the equation of hydrostatic equilibrium in a form that can be integrated directly, if we assume $`\beta `$ to be independent of position in the star. Using the “radiative zero” surface boundary condition, $`P_g=0`$ at $`r=R_{}`$, we obtain

$$\lambda =\frac{(1\beta )}{(ba+1)}.$$
$`(12)`$
Note that equation (12) includes the full effect of the Kramers’ law opacity, subject only to the approximation that $`\beta `$ is independent of position.

Using equation (10) for $`\kappa `$ in equation (8b), expressing $`\rho `$ in terms of $`P_g`$ and $`T`$ with the ideal gas law, and again using equation (11) with $`\beta `$ assumed independent of position, we can also integrate equation (8b), obtaining for the temperature distribution in the stellar envelope

$$T=\frac{(ba+\frac{1}{4})}{(ba+1)}\frac{GM_{}\mu H}{k}\left(\frac{1}{r}\frac{1}{R_{}}\right),$$
$`(13)`$
where we assume $`T=0`$ at $`r=R_{}`$, consistent with the approximation made in deriving equation (12). Equation (13) expresses the temperature at position $`r`$ in terms of constants and stellar properties, including the stellar radius $`R_{}`$.

Equations (12) and (13) express the conditions of hydrostatic and thermal equilibrium in the stellar envelope. Using these equations, we can also obtain an expression for the mass contained in the stellar envelope above the HBRS at radius $`r_s`$:

$$\mathrm{\Delta }M_{\mathrm{env}}_{r_s}^R_{}4\pi \rho r^2𝑑r=\frac{\mu H}{k}\frac{\beta }{(1\beta )}\frac{4\pi }{3}R_{}^3a_r\left[\left(\frac{ba+\frac{1}{4}}{ba+1}\right)\frac{GM_{}\mu H}{kR_{}}\right]^3_{r_s/R_{}}^1\left(\frac{1}{s}1\right)^3s^2𝑑s,$$
$`(14)`$
where $`sr/R_{}`$. The integral in equation (14) can be evaluated analytically and yields

$$_{r_s/R_{}}^1\left(\frac{1}{s}1\right)^3s^2𝑑s=\mathrm{ln}\frac{r_s}{R_{}}3\left(1\frac{r_s}{R_{}}\right)+\frac{3}{2}\left[1\left(\frac{r_s}{R_{}}\right)^2\right]\frac{1}{3}\left[1\left(\frac{r_s}{R_{}}\right)^3\right].$$
$`(15)`$
Note that the factors of $`R_{}`$ cancel in the constants in equation (14), so that the envelope mass depends only upon $`\beta `$ (and thus, from equation  upon $`\lambda `$), $`M_{}`$, and the fractional radius $`r_s/R_{}`$. Consequently, as the stellar luminosity increases ($`\lambda `$ increases and $`\beta `$ accordingly decreases), the parameter $`r_s/R_{}`$ must vary in order to keep $`\mathrm{\Delta }M_{\mathrm{env}}`$ fixed. From the ratio of the quantity $`\mathrm{\Delta }M_{\mathrm{env}}`$ in the perturbed state to the quantity $`\mathrm{\Delta }M_{\mathrm{env}}^{(0)}\mathrm{\Delta }M_{\mathrm{env}}`$ in the unperturbed state, together with equation (12), we obtain

$$\frac{\lambda }{\lambda ^{(0)}}=\frac{\mathrm{ln}s3(1s)+\frac{3}{2}(1s^2)\frac{1}{3}(1s^3)}{\mathrm{ln}s_03(1s_0)+\frac{3}{2}(1s_0^2)\frac{1}{3}(1s_0^3)},$$
$`(16)`$
where now $`sr_s/R_{}`$, and $`s_0`$ is the value of $`s`$ in the unperturbed state. If $`r_s10^{10}`$ cm and $`R_{}10^{12}`$ cm – e.g., corresponding to a star leaving the AGB – then $`s_010^2`$. If we denote the temperature of the HBRS by $`T_s`$, we can use equation (13) to express it in the form of a ratio to the temperature $`T_s^{(0)}`$ in the unperturbed state:

$$y\frac{T_s}{T_s^{(0)}}=\frac{\left({\displaystyle \frac{1}{s}}1\right)}{\left({\displaystyle \frac{1}{s_0}}1\right)}.$$
$`(17)`$
Note that equations (16) and (17) provide a parametric relation for the quantity $`\lambda /\lambda ^{(0)}`$ in terms of the dimensionless temperature $`y`$ and the constant $`s_0`$.

Next, consider the energy equation, using the form of $`TdS`$ that includes both radiation and gas pressure. The general expression can be written in the form (cf., Cox and Giuli 1968, p. 206 ff)

$$TdS=C_V\left[dT\frac{T}{\rho }(\mathrm{\Gamma }_31)d\rho \right].$$
$`(18)`$
Here the quantity $`(\mathrm{\Gamma }_31)`$ is the thermodynamic derivative

$$(\mathrm{\Gamma }_31)\left(\frac{\mathrm{ln}T}{\mathrm{ln}\rho }\right)_S=\frac{2}{3}\frac{(43\beta )}{(87\beta )},$$

and the last equality gives the explicit expression for a mixture of gas plus radiation. Using equation (18), the energy equation (5e) becomes

$$\frac{L_r}{r}=4\pi r^2\rho \left\{ϵC_V\left[\frac{T}{t}\frac{T}{\rho }(\mathrm{\Gamma }_31)\frac{\rho }{t}\right]vC_V\left[\frac{T}{r}\frac{T}{\rho }(\mathrm{\Gamma }_31)\frac{\rho }{r}\right]\right\}.$$
$`(19)`$
We now wish to integrate this equation over the stellar envelope. If we assume

$$ϵ=ϵ_0X\rho T^\nu ,$$
$`(20)`$
we can integrate the first term on the right-hand side of equation (19) to obtain the luminosity generated by the HBRS,

$$L_s_{r_s}^R_{}ϵ𝑑M_rL_s^{(0)}\left(\frac{X_s}{X_s^{(0)}}\right)\left(\frac{\rho _s}{\rho _s^{(0)}}\right)\left(\frac{T_s}{T_s^{(0)}}\right)^\nu .$$
$`(21a)`$
Here the subscript $`s`$ denotes conditions at the HBRS at some time $`t`$, and the superscript $`(0)`$ denotes the conditions in the unperturbed state. The temperature exponent $`\nu `$ in the H-burning reaction rate has the value $`\nu 26`$ for the $`pp`$ chain, which is appropriate for the relatively low shell-burning temperatures we consider here, while it has the value $`\nu 1020`$ for the CNO cycle (cf., Hansen and Kawaler 1994, p. 238, 241). From equations (11) and (12) we can also obtain the relative density distribution in the envelope,

$$\frac{\rho _s}{\rho _s^{(0)}}=\frac{\lambda ^{(0)}}{\lambda }\left(\frac{T_s}{T_s^{(0)}}\right)^3,$$
$`(22)`$
and this can be used to replace $`\rho _s/\rho _s^{(0)}`$ in equation (21a).

In the term involving time derivatives in equation (19), we use the continuity equation to eliminate $`\rho /t`$. The resulting integrand involves the spatial derivative of the mass-flow rate $`F`$. We assume that $`F`$ is independent of position as argued previously and accordingly neglect this term. The remaining integral involving a time derivative can be written in the approximate form

$$_{r_s}^R_{}C_V\frac{T}{t}𝑑M_rC_V^{}\mathrm{\Delta }M_{\mathrm{env}}\frac{dT_s}{dt},$$
$`(21b)`$
where $`C_V^{}`$ is some average heat capacity for the envelope, $`\mathrm{\Delta }M_{\mathrm{env}}`$ is the envelope mass defined in equation (14), and $`T_s`$ is again the temperature of the HBRS. The remaining terms in equation (19) represent heat advection. We estimate them to be only $`2\times 10^3`$ as large as other terms in this equation, and we have neglected them.

Collecting equations (21) together, we can thus write the integral of equation (19) in the form

$$L_{}=L_1+L_s^{(0)}\left(\frac{X_s}{X_s^{(0)}}\right)\left(\frac{\lambda ^{(0)}}{\lambda }\right)\left(\frac{T_s}{T_s^{(0)}}\right)^{\nu +3}L_{}^{(0)}\tau _{\mathrm{th}}\frac{d}{dt}\left(\frac{T_s}{T_s^{(0)}}\right),$$
$`(23)`$
where $`L_1`$ denotes the luminosity originating interior to the HBRS, and

$$\tau _{\mathrm{th}}\frac{C_V^{}\mathrm{\Delta }M_{\mathrm{env}}T_s^0}{L_{}^0}20\mathrm{years}$$
$`(24)`$
is the thermal relaxation timescale. In evaluating $`\tau _{\mathrm{th}}`$ we have taken $`\mathrm{\Delta }M_{\mathrm{env}}0.1M_{}`$; $`T10^6`$ K, which may be typical of a fuel-starved HBRS; and $`L_{}10^4L_{}`$.

In the final equation (5f), we first carry out a Galilean transformation to a frame moving with the HBRS. For any quantity $`f(r,t)`$, if we define a new position coordinate to be $`r^{}rv_st`$, where $`v_s`$ is the propagation speed of the HBRS, the general relations among the partial derivatives in the old and new frames are

$$\frac{f}{r^{}}|_t=\frac{f}{r}|_t\mathrm{and}\frac{f}{t}|_r^{}=\frac{f}{t}|_r+v_s\frac{f}{r}|_t.$$
$`(25)`$
If we apply the relations (25) we can express equation (5f) in the frame moving with the HBRS:

$$\frac{X}{t}|_r^{}+(vv_s)\frac{X}{r^{}}|_t=\frac{ϵ(r^{},t)}{Q^{}},$$
$`(26)`$
where $`Q^{}`$ is given by equation (1). If we define the position of the HBRS by the condition

$$X(r_s,t)\frac{1}{2}X_0,$$
$`(27)`$
where $`X_00.7`$ is the primordial H abundance, then at the location of this shell, we have $`X/t|_{r_s}0`$. Evaluating equation (26) at this location thus gives the shell speed:

$$[v(r_s)v_s]\frac{X}{r_s}|_t=\frac{ϵ(r_s,t)}{Q^{}},$$
$`(28)`$
where $`v(r_s)`$ is the speed of the flow, as measured in the rest frame of the star, at the position $`r_s`$ of the HBRS, and $`v_s`$ is the desired propagation speed of the shell. If we use equation (28) to describe the presumed unperturbed steady-flow condition and subtract this from equation (26), we obtain for the perturbed flow the result

$$\frac{X}{t}|_r^{}+(vv_s)\frac{X}{r^{}}|_t(v^{(0)}v_s^{(0)})\frac{X^{(0)}}{r^{}}|_t=\frac{[ϵ(r^{},t)ϵ^{(0)}(r^{},t)]}{Q^{}},$$
$`(29)`$
We now need to work out the integrals of the various terms in this equation. The first term can be written in the form

$$_{M_s}^M_{}\frac{X}{t}𝑑M_r^{}\mathrm{\Delta }M_{\mathrm{nuc}}\overline{\frac{X}{t}}\mathrm{\Delta }M_{\mathrm{nuc}}\frac{dX_s}{dt}.$$
$`(30a)`$
This equation defines the spatial average of the time derivative, which we approximate as the time derivative of the H abundance in the HBRS. Note that $`X`$ is constant in time and space except in the very thin HBRS. Accordingly the mass $`\mathrm{\Delta }M_{\mathrm{nuc}}`$ in HBRS regions is relatively small; we use the approximate value $`\mathrm{\Delta }M_{\mathrm{nuc}}5\times 10^4M_{}`$ for this mass. (This value is actually consistent with the final H-burning shell mass found by Iben (1984), Iben and MacDonald (1986), and Herwig et al. (1999) in late AGB or immediate post-AGB phases of stellar evolution. However since our model also requires a significant luminosity internal to the HBRS, $`L_1`$, we do not require that the HBRS is the entire H-burning shell. Some arbitrary ratio of H and He could be burning below the HBRS to provide $`L_1`$.)

In the terms in equation (29) that involve spatial gradients of $`X`$, we write $`dM_r^{}=4\pi (r^{})^2\rho (r^{},t)dr^{}`$ and use equation (6) to replace $`v`$ with the mass-flow rate $`F`$, which we assume to be independent of position. The resulting integral has the form

$$_{M_s}^M_{}\left[(vv_s)\frac{X}{r^{}}(v^{(0)}v_s^{(0)})\frac{X^{(0)}}{r^{}}\right]𝑑M_r^{}=[(FF^{(0)})(F_sF_s^{(0)})]X_0,$$
$`(30b)`$
where $`F_s4\pi r^2\rho v_s`$, and $`X_0`$ is again the primordial H abundance. The integral of the terms on the right-hand side of equation (29) just gives the difference in the HBRS luminosities in the perturbed and unperturbed states. Collecting these results together, and for simplicity neglecting $`F_sF_s^{(0)}`$, we write the integral of the $`X`$ equation in the approximate form

$$\mathrm{\Delta }M_{\mathrm{nuc}}\frac{dX_s}{dt}+(FF^{(0)})X_0=\frac{(L_sL_s^{(0)})}{Q^{}}.$$
$`(31)`$
The system of equations describing the presumed relaxation oscillations now consists of (23) and (31), with the dimensionless luminosity $`\lambda `$ expressed in terms of the dimensionless shell temperature $`y`$ by the parametric equations (16) and (17). For given values of $`y`$ and $`s_0`$, we solve equation (17) to obtain $`s`$, and we use these values in equation (16) to obtain $`\lambda /\lambda ^{(0)}`$. For convenience, we collect the resulting equations together here, defining $`xX_s/X_s^{(0)}`$:

$$L_{}=L_1+L_s^{(0)}\left(\frac{\lambda ^{(0)}}{\lambda }\right)xy^{\nu +3}L_{}^{(0)}\tau _{\mathrm{th}}\frac{dy}{dt},$$
$`(32a)`$
$$X_s^{(0)}\mathrm{\Delta }M_{\mathrm{nuc}}\frac{dx}{dt}+(FF^{(0)})X_0=\frac{L_s^{(0)}}{Q^{}}\left[\left(\frac{\lambda ^{(0)}}{\lambda }\right)xy^{\nu +3}1\right],$$
$`(32b)`$
$$\left(\frac{1}{s}1\right)=y\left(\frac{1}{s_0}1\right),$$
$`(32c)`$
and

$$\frac{\lambda }{\lambda ^{(0)}}=\frac{\mathrm{ln}s3(1s)+\frac{3}{2}(1s^2)\frac{1}{3}(1s^3)}{\mathrm{ln}s_03(1s_0)+\frac{3}{2}(1s_0^2)\frac{1}{3}(1s_0^3)}.$$
$`(32d)`$
For radiation-driven mass loss, we assume the simple form given by equation (4), which we write as

$$\dot{M}=\alpha L_{}.$$
$`(33)`$
Combining (6), (9), and (33), we obtain

$$\frac{\dot{M}}{\dot{M}^{(0)}}=\frac{F}{F^{(0)}}=\frac{L_{}}{L_{}^{(0)}}=\frac{\lambda }{\lambda ^{(0)}}.$$
$`(34)`$
Equation (34) eliminates equation (6) and expresses the mass flux $`F`$ in terms of $`\lambda `$. We can now put the equations in fully dimensionless form. Dividing equation (32a) by $`L_{}^{(0)}`$, we obtain

$$\frac{\lambda }{\lambda ^{(0)}}=A+(1A)\left(\frac{\lambda ^{(0)}}{\lambda }\right)xy^{\nu +3}\tau _{\mathrm{th}}\frac{dy}{dt},$$
$`(35a)`$
where $`AL_1/L_{}^{(0)}`$. Note that, in a truly steady state with $`dy/dt=0`$, we have $`x=1`$, $`y=1`$, and $`\lambda =\lambda ^{(0)}`$, obtaining from equation (23) $`L_{}^{(0)}L_s^{(0)}+L_1`$. Similarly, we divide equation (32b) by the quantity $`X_s^{(0)}\mathrm{\Delta }M_{\mathrm{nuc}}`$ and use equation (34) to obtain

$$\frac{dx}{dt}+\left(\frac{\lambda }{\lambda ^{(0)}}1\right)\frac{1}{\tau _{\mathrm{fsm}}}=\frac{1}{\tau _{\mathrm{ex}}}\left[\left(\frac{\lambda ^{(0)}}{\lambda }\right)xy^{\nu +3}1\right].$$
$`(35b)`$
In equations (35), $`\tau _{\mathrm{th}}`$ is given by equation (24). We define

$$\tau _{\mathrm{fsm}}\frac{X_s^{(0)}\mathrm{\Delta }M_{\mathrm{nuc}}}{X_0F^{(0)}}\frac{X_s^{(0)}\mathrm{\Delta }M_{\mathrm{nuc}}}{X_0|\dot{M}|}25\mathrm{years},$$
$`(36a)`$
as the timescale on which the fuel supply is modulated by mass loss; i.e., the time required for mass loss to strip material out of the HBRS. In evaluating $`\tau _{\mathrm{fsm}}`$ we assume the mass in the HBRS to be $`\mathrm{\Delta }M_{\mathrm{nuc}}5\times 10^4M_{}`$, as noted above, and we assume $`\dot{M}10^5M_{}`$ yr-1 to be the mass-loss rate. We define the fuel-exhaustion timescale similarly:

$$\tau _{\mathrm{ex}}\frac{Q^{}X_s^0\mathrm{\Delta }M_{\mathrm{nuc}}}{L_s^0}1.75\times 10^4\mathrm{years}.$$
$`(36b)`$
An important constraint is provided by the fact that the H abundance can never exceed the primordial value $`X_0`$ nor become negative. From equation (27), this correspondingly constrains our dimensionless parameter $`xX_s/X_s^{(0)}`$ to the range

$$0x(X_0/\frac{1}{2}X_0)=2.$$
$`(37)`$
Equations (35a) and (35b), together with the constraint equation (37), the parametric relations (32c) and (32d), and the definitions of the characteristic timescales (24), (36a), and (36b), constitute the system of equations that characterize this relaxation-oscillator model. In these equations, the quantity $`\lambda ^{(0)}`$ is a dimensionless number presumably just slightly less than unity (i.e., the luminosity is just slightly less than the Eddington luminosity).

The mathematical nature of the governing system of nonlinear differential equations (35a,b) is better seen by writing them in the fully nondimensional form

$$\frac{dy}{dt}=A+(1A)xh(y)g(y),$$
$`(38a)`$
$$\frac{dx}{dt}=\gamma _1[g(y)1]\gamma _2[xh(y)1],$$
$`(38b)`$
where

$$g(y)\frac{\lambda }{\lambda ^{(0)}},h(y)\frac{y^{(\nu +3)}}{g(y)},\gamma _1\frac{\tau _{\mathrm{th}}}{\tau _{\mathrm{fsm}}},\gamma _2\frac{\tau _{\mathrm{th}}}{\tau _{ex}},$$
$`(39)`$
and where $`t`$ is now a nondimensional time, scaled by the thermal time scale $`\tau _{\mathrm{th}}`$. Note that $`g(y)\lambda /\lambda ^{(0)}`$ is a nonlinear function of $`y`$ defined by equations (32c) and (32d). Note also that $`g(1)=h(1)=1`$ and that equations (38) have the equilibrium solution $`x=y=1`$. A linear stability analysis, presented in the Appendix, shows that this equilibrium is unstable to growing oscillations for a range of values of the stellar parameters.

5. Numerical Solutions of the Model Equations

We have carried out a number of numerical solutions of the nonlinear system (38) in order to explore the behavior of the relaxation oscillations and to determine their dependence on the various parameters in the model. As a specific example, consider first a benchmark model with $`M_{}=0.7M_{}`$, $`L_{}=10^4L_{}`$, $`R_{}=10^{13}`$ cm, $`\mathrm{\Delta }M_{\mathrm{env}}=0.1M_{}`$, $`\mathrm{\Delta }M_{\mathrm{nuc}}=5\times 10^4M_{}`$, and $`T_s^{(0)}=10^6`$ K. For these parameters, equation (13) and the definition $`s_0r_s^{(0)}/R_{}`$ yield $`s_0=0.03`$. We adopt the values $`\nu =3`$ and $`AL_1/L_{}^{(0)}=0.9`$ in this model. These values yield $`\tau _{\mathrm{th}}=20`$ years, $`\tau _{\mathrm{fsm}}=25`$ years, and $`\tau _{\mathrm{ex}}=1.75\times 10^4`$ years. For these parameter choices, the numerical solution of equations (38) yields sustained relaxation oscillations with a period $`P72\tau _{\mathrm{th}}`$ 1400 years, as shown in Fig. 1. During the initial rise in temperature and luminosity, the H abundance $`X`$ drops rapidly, reaching $`X=0`$ in $`39\tau _{\mathrm{th}}`$ 780 years. The resulting relaxation oscillations reach a peak temperature $`T_s6\times 10^7`$ K and a peak luminosity $`L_{}3.28L_{}^{(0)}`$. The half-width of a pulse, defined as the time during which the luminosity exceeds half its peak value, is $`\mathrm{\Delta }t28\tau _{\mathrm{th}}`$ 560 years $`=0.39P`$. We define the “duty cycle” $`W`$ of a pulse as the time required for the shell temperature to relax from its peak value to $`T_s^{(0)}`$ (i.e., for $`y`$ to return to the value $`y=1`$). This is approximately the time required for the star to radiate away the heat generated in the thermonuclear shell flash. A simple analytic argument based on our dimensionless model equations gives the duty cycle as $`W26\tau _{\mathrm{th}}`$, while our explicit numerical calculation yields $`W=39.4\tau _{\mathrm{th}}=0.55P`$.

The H abundance in the shell remains vanishingly small until the luminosity – and the accompanying mass-loss rate – has declined sufficiently for the high-temperature shell again to overtake the H fuel supply. The recovery between successive pulses is governed by the time required for the H abundance – and the associated energy-generation rate – to recover to, and then exceed, the equilibrium rate. The time between pulses is $`P72\tau _{\mathrm{th}}=1400`$ years, closely comparable to the timescale between the observed circumnebular shells in NGC 6543. We obtain similar sustained relaxation oscillations starting from shell temperatures $`T_s`$ either slightly less than or slightly greater than $`T_s^{(0)}`$, where $`T_s^{(0)}`$ is the presumed steady-state burning temperature; in this sense, these results are robust.

We have investigated the numerical accuracy of these results in two ways. Calculations with timesteps selected to be a factor of two larger and smaller than the values used in the calculations we report here give values of $`P`$ within $`\pm 0.2`$%, values of $`W`$ within $`\pm 1.8`$%, and values of the peak temperature within $`\pm 2.0`$%. Calculations for this same model obtained using two different algorithms yield values of $`P`$ that agree to within 20%, values of $`W`$ within 20%, and values of the peak temperature within about 3%.

We have explored the variations in the properties of these relaxation oscillations for the range of model parameters listed in Table 1 below. Model F is the illustrative benchmark case discussed above. The numerical results are summarized in Table 2 below. Depending upon the specific choice of parameters, we obtain oscillations with a range of periods, duty cycles, and peak temperatures and luminosities. Note that the peak temperatures for all except model F exceed $`10^8`$ K, violating our simplifying model assumption that energy production in the shell is due solely to the H-burning $`pp`$ reactions.

| Table 1. Model Parameters | | | | | | | | | |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | Model | $`\dot{M}`$ | $`T_s^{(0)}`$ | $`s_0`$ | $`\nu `$ | $`A`$ | $`\tau _{\mathrm{th}}`$ | $`\tau _{\mathrm{fsm}}`$ | $`\tau _{\mathrm{ex}}`$ |
|  | A | $`10^6`$ | $`10^7`$ | 0.003 | 3 | 0.9 | 200 | 500 | $`3.50\times 10^4`$ |
|  | B | $`10^5`$ | $`10^7`$ | 0.003 | 3 | 0.9 | 200 | 50 | $`3.50\times 10^4`$ |
|  | C | $`10^6`$ | $`10^6`$ | 0.03 | 3 | 0.9 | 20 | 500 | $`3.50\times 10^4`$ |
|  | D | $`10^5`$ | $`10^6`$ | 0.03 | 3 | 0.9 | 20 | 50 | $`3.50\times 10^4`$ |
|  | E | $`10^4`$ | $`10^6`$ | 0.03 | 3 | 0.9 | 2 | 5 | $`3.50\times 10^4`$ |
|  | F | $`10^5`$ | $`10^6`$ | 0.03 | 3 | 0.9 | 20 | 25 | $`1.75\times 10^4`$ |

The mass-loss rate $`\dot{M}`$ is in $`M_{}`$ yr-1, $`T_s^{(0)}`$ is in K, and all timescales are in years.

| Table 2. Model Results | | | | | |
| --- | --- | --- | --- | --- | --- |
|  | Model | $`P`$ | $`T_s^{\mathrm{max}}/T_s^{(0)}`$ | $`L_s^{\mathrm{max}}/L_s^{(0)}`$ | $`W`$ |
|  | A | $`2.28\times 10^4`$ | 31.7 | 1.86 | 0.45 |
|  | B | $`1.38\times 10^4`$ | 30.3 | 1.85 | 0.67 |
|  | C | $`9.12\times 10^3`$ | 171.7 | 3.85 | 0.16 |
|  | D | $`2.47\times 10^3`$ | 140.9 | 3.73 | 0.55 |
|  | E | $`9.29\times 10^2`$ | 1,409. | 5.05 | 0.92 |
|  | F | $`1.44\times 10^3`$ | 62.2 | 3.28 | 0.55 |

The oscillation period $`P`$ is in years, and the “duty cycle” $`W`$ of the pulse is defined as the fraction of the pulse period required for the shell temperature $`T_s`$ to decline from its maximum value $`T_s^{\mathrm{max}}`$ to the mean steady-state temperature $`T_s^{(0)}`$.

In order to investigate the dependence of the oscillations on the paramter $`AL_1/L_{}^{(0)}`$ we did a series of calculations using $`\nu =3`$ and $`\mathrm{\Delta }M_{\mathrm{nuc}}=0.001`$. In this case we find that oscillations occur only when the luminosity ratio $`AL_1/L_{}^{(0)}`$ lies in the range $`0.84<A<0.92`$. We can also estimate the range of values of $`A`$ that give oscillations from a linear stability analysis of the governing nonlinear equations (see the Appendix). In this case, the linear analysis predicts growing oscillations in the range $`0.854<A<0.907`$. We have not found relaxation oscillations for $`A0.8`$, even though we have carried out a number of calculations in this range. For example, with $`A=0.3`$ and $`\nu =10`$, the system goes through a single flash and then stabilizes at an HBRS luminosity lower than the initial assumed steady-state luminosity. Evidently this is not a self-consistent solution of the problem. Further, neglect of the luminosity $`L_1`$ emerging from below the HBRS (i.e., setting $`A=0`$) prevents matter in the envelope from re-heating when the luminosity drops and mass loss shuts down. Conversely, if the (presumed constant) interior luminosity $`L_1`$ is too large, the system goes through a single thermonuclear pulse and then settles down into a state of steady burning. For example, with $`L_1/L_{}^{(0)}=0.95`$, $`\nu =3`$, $`\tau _{\mathrm{fsm}}/\tau _{\mathrm{th}}=2.5`$, and $`\tau _{\mathrm{ex}}/\tau _{\mathrm{th}}=3000`$, there is no oscillation, and the shell temperature and luminosity stabilize at the assumed steady-state values. Thus, there appears to be a rather small region of parameter space in which the star exhibits such fuel-supply-limited relaxation oscillations.

We have studied the dependence of the relaxation oscillations upon the assumed value of $`\mathrm{\Delta }M_{\mathrm{nuc}}`$, taking $`A=0.9`$, $`\nu =3`$, and $`s_0=0.01`$ for all models. As shown in Table 3 below, we found that the oscillation period $`P`$ decreases and the peak values of the temperature and luminosity decline as $`\mathrm{\Delta }M_{\mathrm{nuc}}`$ decreases. The period appears to approach an asymptotic value $`(1517)\times \tau _{\mathrm{th}}`$ for small HBRS masses. For the model in Table 3 with the lowest listed value of $`\mathrm{\Delta }M_{\mathrm{nuc}}`$, the amount of energy available from the shell flash is so small that the oscillations actually decay slowly in time, rather than attaining a roughly constant amplitude.

| Table 3. Variations with $`\mathrm{\Delta }M_{\mathrm{nuc}}`$ | | | | |
| --- | --- | --- | --- | --- |
|  | $`\mathrm{\Delta }M_{\mathrm{nuc}}/M_{}`$ | $`P/\tau _{\mathrm{th}}`$ | $`T_s^{\mathrm{max}}/T_s^{(0)}`$ | $`L_s^{\mathrm{max}}/L_s^{(0)}`$ |
|  | 0.001 | 33 | 2.8 | 1.36 |
|  | 0.0006 | 26 | 2.1 | 1.26 |
|  | 0.0005 | 22 | 1.8 | 1.20 |
|  | 0.0004 | 20 | 1.6 | 1.16 |
|  | 0.0003333 | 19 | 1.5 | 1.14 |
|  | 0.0003077 | 18 | 1.25 | 1.08 |
|  | 0.0002857 | 17–15 | 1.1 | 1.03 |

We have also studied the dependence of the relaxation oscillations upon the assumed value of $`\nu `$, taking $`\mathrm{\Delta }M_{\mathrm{nuc}}=0.001M_{}`$, $`A=0.9`$, $`T_s^{(0)}=10^7`$ K, and $`s_0=0.01`$ for all models. As shown in Table 4 below, the system does not depend very much upon $`\nu `$, apparently because the temperature spikes so rapidly. We found no oscillations for $`\nu >7`$. Within the limitations of our simplified model, a stronger dependence of the thermonuclear rate upon $`T`$ apparently liberates energy so rapidly that the burning simply adjusts to a different stable-burning state. For calculations with $`\nu =6`$, the peak temperatures reach $`5\times 10^8`$ K, and our approximation for the nuclear burning rate is invalid.

| Table 4. Variations with $`\nu `$ | | | | |
| --- | --- | --- | --- | --- |
|  | $`\nu `$ | $`P/\tau _{\mathrm{th}}`$ | $`T_s^{\mathrm{max}}/T_s^{(0)}`$ | $`L_s^{\mathrm{max}}/L_s^{(0)}`$ |
|  | 3 | 33 | 2.8 | 1.36 |
|  | 4 | 35 | 3.5 | 1.44 |
|  | 5 | 38 | 3.9 | 1.47 |
|  | 6 | 51 | 3.9 | 1.48 |

One aspect of the results surprised us: the conditions under which we found oscillations correspond to very low values of the nuclear-energy-generation exponent $`\nu `$. The values for which we find oscillations correspond to energy generation by the $`pp`$ chain, whereas we had expected that the CNO cycle would dominate. Conceivably, this might also be a consequence of the terminal stage of shell burning.

Finally, we note that in order to justify our suggested link between these relaxation oscillations and ejection of mass shells, the time scale for acceleration from rest to the escape velocity must be less than the oscillation period. This condition is easily satisfied; the acceleration time for material illuminated from below by luminosity $`L_{}>L_{\mathrm{Edd}}`$ is $`0.13\mathrm{year}\times (M_{}/M_{_{\mathrm{}}})^{1/2}(R_{}/10^{13}\mathrm{cm})^{3/2}(L_{}/L_{\mathrm{Edd}}1)^1`$, assuming the Thomson cross section as a conservative bound. This acceleration time is always much less than the $`>10^3`$ year oscillation period we find, so that we are well justified in applying our model to PN systems with evidence for periodic ejecta.

6. Conclusions

We have shown that fuel-supply-limited relaxation oscillations may develop in some stars near the ends of their nuclear-active lifetimes just prior to the ejection of a planetary nebula. In the benchmark model we discuss in this paper, with $`M_{}<M_{}`$, $`L_{}10^4L_{}`$, and $`|\dot{M}|10^5M_{}`$ yr-1, such oscillations occur when the luminosity generated by the H-burning relaxation shell (HBRS) falls to $``$ 10% of the total luminosity of the star. The pulsation timescale we find to be $`72\tau _{\mathrm{th}}57\tau _{\mathrm{fsm}}1400`$ years, closely comparable to the timescale inferred for the observed circumnebular shells (Balick et al. 2001; Zijlstra et al. 2002), for the present choice of parameters. While our choice of parameters does produce results consistent with the observed circumnebular shells, it is not clear whether such parameters are actually attained during the evolution of a real star. To determine whether such oscillations occur in more realistic stellar models, or whether they do indeed provide the correct explanation for the observed shells surrounding some planetary nebulae, will require much more detailed calculations of stellar evolution with mass loss than have been done to date.

Acknowledgments

HMVH acknowledges support from the National Science Foundation, which made his participation in this work possible. He also wishes to express his gratitude to the University of Rochester and to the Department of Terrestrial Magnetism at the Carnegie Institution of Washington, both of which graciously served as host institutions during the course of this work. In particular, he thanks S. Keiser for her help with the calculations and display of the results. JHT acknowledges support from NSF grant AST95-28398 and thanks Clare Hall and the Department of Applied Mathematics and Theoretical Physics, University of Cambridge, for their hospitality during the final preparation of this manuscript. EB acknowledges support from DOE grant DE-FG02-00ER54600, and the Laboratory for Laser Energetics. Support for AF was provided at the University of Rochester by NSF grant AST-9702484, AST-0098442, NASA grant NAG5-8428 and the Laboratory for Laser Energetics.

Figure Caption

Fig. 1. Fuel-supply-limited relaxation oscillations of the model described in §4 (model F in Tables 1 and 2). The dotted curve represents the dimensionless temperature $`yT_s/T_s^{(0)}`$, the dashed curve gives the dimensionless HBRS luminosity $`L_s/L_s^{(0)}`$, and the solid curve is the dimensionless H abundance $`xX_s/X_s^{(0)}`$ in the HBRS. The abscissa gives the time elapsed from the start of the calculation in units of the thermal timescale defined by equation (24), $`t/\tau _{\mathrm{th}}`$. At the onset of a pulse, the rapid increase in luminosity drives the H abundance in the HBRS to zero in less than $`39\tau _{\mathrm{th}}780`$ years. During this initial phase, the shell temperature rises rapidly to a peak value $`62`$ times the intial value. The energy generated during this pulse is gradually radiated away during the next $`39\tau _{\mathrm{th}}780`$ years, and the luminosity declines toward the initial luminosity on the same timescale. When the surface luminosity and the associated mass-loss rate have declined sufficiently, the H abundance in the HBRS begins to build up again. In $`72\tau _{\mathrm{th}}1400`$ years $`57\tau _{\mathrm{fsm}}`$, the energy generation rate in the HBRS builds up to a level exceeding the level of steady-state burning, and another pulse occurs.

Appendix: Linear Stability Analysis

To investigate the stability of the equilibrium solution $`x=y=1`$ of equations (38), let $`x=1+\xi `$ and $`y=1+\eta `$ and expand the functions $`g(y)`$ and $`h(y)`$ in Taylor series around $`y=1`$ in the form

$$g(1+\eta )=1+b\eta +\mathrm{O}(\eta ^2),h(1+\eta )=1+c\eta +\mathrm{O}(\eta ^2),$$
$`(A1)`$
where the coefficient $`b`$ depends on the parameter $`s_0`$ and the coefficient $`c`$ depends on the parameters $`s_0`$ and $`\nu `$. Substituting these expressions into equations (38) and linearizing yields the equations

$$\frac{d\xi }{dt}=\gamma _2\xi (\gamma _1b+\gamma _2c)\eta ,$$
$`(A2a)`$
$$\frac{d\eta }{dt}=(1A)\xi +[c(1A)b]\eta .$$
$`(A2b)`$
If we assume a nonzero solution of these equations in the form $`\xi =\xi _0\mathrm{exp}(\alpha t)`$, $`\eta =\eta _0\mathrm{exp}(\alpha t)`$, then we find that $`\alpha `$ satisfies a quadratic equation with roots

$$\alpha =[B\pm (B^24C)^{1/2}]/2,$$
$`(A3)`$
where

$$B=(1A)cb\gamma _2,C=[(1A)\gamma _1+\gamma _2]b.$$
$`(A4)`$
The solution is an oscillation when $`B^24C<0`$ and these oscillations are growing when $`B>0`$.

Using this analysis, we can determine the regions of parameter space in which we expect to find growing oscillations. For example, with $`s_0=0.03`$, $`\nu =3`$, $`\gamma _1=0.8`$, and $`\gamma _2=0.001143`$ (as in model F in Table 1), we find that the solutions of the linearized equations are oscillatory for values of the parameter $`A`$ in the range $`0.807<A<0.957`$, and these oscillations are growing for $`A`$ in the range $`0.807<A<0.908`$.

References

Balick, B., Wilson, J., and Hajian, A. R. 2001, A. J., 121, 354.

Blackman, E.G., Frank, A., Markiel, J.A., Thomas, J.H., and Van Horn, H. M. 2001, Nature 409, 485

Bowen, G. H. 1988, Ap. J., 329, 299.

Bowen, G. H., and Willson, L. A. 1991, Ap. J., 375, L53.

Clayton, D. D. 1968, Principles of Stellar Evolution and Nucleosynthesis (McGraw-Hill: New York).

Cox, A. 2000, ed. Allen’s Astrophysical Quantities, 4th ed. (Springer-Verlag: New York).

Cox, J. P., and Giuli, R. T. 1968, Principles of Stellar Structure (Gordon and Breach: New York).

Dorman, B., Rood, R. T., and O’Connell, R. W. 1993, Ap. J., 419, 596.

Garcia-Segura et al. 2002 \*** Need refc. \***

Hansen, C. J., and Kawaler, S. D. 1994, Stellar Interiors: Physical Principles, Structure, and Evolution, (Springer-Verlag: New York).

Herwig, F., Blöcker, T., Langer, N., and Driebe, T. 1999, A. & A., 349, L5.

Iben, I., Jr. 1984, Ap. J., 277, 333.

Iben, Jr., I. 1987, in Late Stage of Stellar Evolution, ed. S. Kwok and S. R. Pottasch (Reidel: Dordrecht), p. 175.

Iben, Jr., I. 1991, Ap. J. Suppl., bf 76, 55.

Iben, I., Jr., and MacDonald, J. 1986, Ap. J., 301, 164.

Iben, Jr., I., and Renzini, A. 1983, Ann. Rev. Astr. Ap., 21, 271.

Koesterke, L., Dreizler, S., and Rauch, T. 1998, A. & A., 330, 1041.

Koesterke, L., and Werner, K. 1998, Ap. J., 500, L59.

Landau, L. D., and Lifshitz, E. M. 1959, Fluid Mechanics, (Pergamon Press: Oxford).

Mastrodemos, N., and Morris, M. 1998, Ap. J., 497, 303.

Mastrodemos, N., and Morris, M. 1999, Ap. J., 523, 357.

Mauron, N., and Huggins, P. J. 2000, Ap. J., 359, 707.

Mazzitelli, I., and D’Antona, F. 1986, Ap. J., 308, 706.

Zijlstra, A.A., Bedding, T.R., & Mattei, J.A. 2002, MNRAS, accepted for publication, astro-ph/0203328

Perinotto, M. 1989, in Planetary Nebulae, ed. S. Torres-Peimbert (Kluwer?: Dordrecht), p. 293.

Reimers, D. 1975a, Mem. Soc. R. Sci. Liege, 6 Ser., 8, 369.

Reimers, D. 1975b, in Problems in Stellar Atmospheres and Envelopes, ed. B. Baschek, W. H. Kegel, and G. Traving (Springer: Berlin?), p. 229.

Renzini, A. 1981, in Physical Processes in Red Giants, ed. I. Iben, Jr., and A. Renzini (Reidel: Dordrecht), p. 431.

Schönberner, D. 1983, Ap. J., 272, 708.

Simis, Y., Icke, V., and Dominik, C. 2001, http://xxx.lanl.gov/abs/astro-ph/0103347.

Soker and Harpaz 1988 \*** Need refc. \***

Terzian, Y., and Hajian, Y. 2000, in Asymmetric Planetary Nebulae. II. From Origins to Microstructure, ed. J. H. Kastner, N. Soker, and S. A. Rappaport, CS-199, (Astron. Soc. Pacific: San Francisco), p. 33.

Vassiliadis, E., and Wood, P. R. 1993, Ap. J., 413, 641.

Willson, L. A. 2000, Ann. Revs. Astron. Ap., 38, 573.