Title: Thermal effects in stellar pulsations
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00097423/00001
 Material Information
Title: Thermal effects in stellar pulsations
Physical Description: iv, 129 leaves : ill. ; 28 cm.
Language: English
Creator: Pesnell, William Dean, 1957-
Publication Date: 1983
Copyright Date: 1983
Subject: Pulsating stars   ( lcsh )
Physics thesis Ph. D
Dissertations, Academic -- Physics -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
Thesis: Thesis (Ph. D.)--University of Florida, 1983.
Bibliography: Bibliography: 125-128.
Additional Physical Form: Also available on World Wide Web
General Note: Typescript.
General Note: Vita.
Statement of Responsibility: by William Dean Pesnell.
 Record Information
Bibliographic ID: UF00097423
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 000444009
oclc - 11353743
notis - ACK4888


This item has the following downloads:

PDF ( 4 MBs ) ( PDF )

Full Text


B '


FOF THE [EFGR.EE Of( L)CTO'n Of Pti lLOSu-'piu



I would like to thank Dr. J. Robert Buchler for his assistance

during the past four years as well as his great patience. In addition, I

would like to thank A. N. Cox for insight into the physical nature of

stellar pulsations and Dr. Oded Regev for useful suggestions and


Several other people have assisted in parts of the work. They

include Dr. Robert Coldwell, who assisted in the writing of the programs

that are inchludic and, Sharoni Eullivant, who demonstrated the use of the

urd pro.ceSsor. I oulid also iike to thank Marie Jose Goupil for reading

the r..-gh draft .f t-.c z.an'iscript.

.:,rue *f the co.niputations in this work were performed at the Los

ij.-mo) NaLrional L borat:,r' Los Alamos, New Mexico, while I was a Summer

'7raduatE Fesearch .As.istsnt there in the summer of 1982. The bulk of the

comput~tion- w-ere d.,e at the University of Florida, primarily under the

IUSIC s-Is ,'teu, for which I ad grateful.


ACKNOWLEDGMENTS ... ............ ........ ................... i

ABSTRACT ... ........ ...... ....... .... ..... .............. iv


I INTRODUCTION.....................................

II BACKGROUND .......................................3

III THEORY...........................................7

IV THE THERMAL OPERATOR ............................10

Asymptotic Expansions...........................16

Asymptotic Thermal Eigenfunctions...............22

V NUMERIC RESULTS.................................33

BW Vulpeculae...................................39

Cepheid model...................................41

VI DISCUSSION ......................................74

VII CONCLUSIONS....................................91



B THE INITIAL MODEL INTEGRATOR ...................103

BIBLIOGRAPHY ...... .......................................... 125

BIOGRAPHICAL SKETCH ...................................... 129


Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy





Chairman: Dr. Jean Robert Buchler
Major Department: Physics

A Sriurm-Li,':u. lie operator, embedded in the linearized radial

pulSation equations, is ujed to analyze the thermal effects of stellar

pul aioni. ampctoric exF.panions are utilized to demonstrate the nature

of these eicenf.,nicoioni .and, to define the thermal response time of a

scell3r moiel. The eiten.,alues and eigenvectors for models of a Cepheid

and et-5 Crcphel variable are given showing the differences in the two

ilas-se of varlt3rle'. TrIe tie scales introduced agree with earlier ones

in L eir ci m-1i:n u-e of predicting the location of possible driving

zor.e:. in additiolr, uwa., of displaying the temperature variations in

the hydrogen loniiatlion ione without the customary "spike" is

demont rated, ani a. new class of pulsations, called "sudden" modes, is




The use of time scales to interpret phenomena found in stars has

been a favorite tool in astrophysics, (see Cox and Ciuli 1968; Saio and

Wheeler 1982; and Shibahashi and Osaki 1981). Primarily, they are used

to simplify the governing equations by ignoring time dependence that

are much longer or averaging those that are much shorter than the

subject under study. Examples of this include the suppression of

hydrodynamic and possibly thermal phenomena when examining nuclear

evolution (Iben 1965), nuclear evolution effects in stellar pulsations

(Christy 1966), and hydrodynamic phenomena when studying thermal flashes

(Henyey and Ulrich 1972). These time scales arise by recasting the

various relationships into non-dimensional forms and defining the

factors multiplying the time derivatives as position dependent time

scales. This produces a good approximation for the discussion of some

attributes of pulsations.

However, we have been using asymptotic formulations to study

stellar pulsations (Buchler 1978), and when we attempted Lo include Lhe

outer layers (Buchler and Regev 1982a), it was found necessary to

understand the thermal structure in a more quantitative way. Ocher

authors have found modes that are dominated by thermal eftects (Wood


1976; King 1980; and Saio and Wheeler 1982) which have been called

"strange modes." In this work, a mathematical formalism for studying

modes with these characteristics is given.

We shall clarify what the thermal response of a star is and show

effects in pulsations that result. A Sturm-Liouville operator, embedded

in the linear, non-adiabatic radial pulsation equations, is used to give

what will be called the thermal eigenfunctions of a star. It will be

shown that under physically realistic assumptions, the eigenvalues are

real and represent decays back to the initial state, here assumed to be

one of thermal equilibrium.

Asymptotic forms for the eigenvalues and eigenfunctions are shown,

and several realistic models have their spectra calculated. A discussion

of the physical requirements for the quasi-adiabatic approximation for

the stability coefficient is given, along with a way of evaluating it

without the introduction of cutoffs in the integration variable. The

lowest order pulsation operator is shown to have a more complicated

structure than the simple adiabatic operator, when the outer layers are

included. This operator, and the accompanying stability coefficient, are

derived for a general star and compared to earlier work. A new class of

modes, dominated by thermal effects, is shown to exist.



The theoretical study of stellar pulsations can be divided into

several general areas. First is the adiabatic wave equation and the

discovery that the periods predicted by this theory agree with

observations (Ledoux and Walraven 1958; Cox, J.P. 1980). Here the entire

star is assumed to react faster via an acoustic process than thermally,

implying that stellar pulsations are basically sound waves propagating

in the envelope. In much of a star's mass, adiabaticity is a very good

assumption and the periods calculated will be correct as long as this is

true. The stability analysis of an adiabatic mode uses the quasi-

adiabatic approximation (Ledoux and Walraven 1958; Ledoux 1965),

although the integrals must have a cutoff introduced to prevent the

outer layers from dominating the result. These calculations demonstrate

that the quasi-adiabatic interior would damp any oscillations unless an

active source of driving was present. Several authors (Cox 1960 and

1958; Zhevakin 1963 and his earlier papers) showed that the

thermodynamic work derivable from the ionization of once ionized helium

could provide the necessary driving for the classical Cepheid variables.

This driving mechanism requires that some of the star be non-adiabatic

and a finely tuned relationship between the location in the star .:. che


ionization zone and the transition from adiabatic to non-adiabatic

behavior be satisfied. Predictions of this theory are narrow instability

strips in the H-R diagram (Cox, J.P. 1980), with sharply defined blue

edges where the evolving star first becomes unstable.

Other classes of stars are driven by similar mechanisms. Both

radial and non-radial oscillations in the ZZ Ceti variables are driven

bY hydrogen ionization i(Sarrfield et al. 1982; Winger, Saio and

Robinion ]9:). Hot wnite dwarfs, in particular PG1159-035, have been

.hown to be jntable to radial pulsations in Starrfield, Cox and Hodson

(19,0), arn this analysis has been extended in Starrfield et al. (1983),

to include non-radial motion and to show the driving is due to the

ionization of carbon and *:-.ygen. Variables such as the RR Lyrae, BL

Herculae, :. Scuti and W VirgLnis classes, are driven primarily by helium

ionization although hydrogen does participate in the destabilization

(,XC J.P. 19~,!.

imnce crie des5tablilzinrg echanism was known, computer codes were

developed that followed tre nonlinear behavior of pulsating stars

.Cnhristy 19o-; King et al. 1964), and reproduced many observed features

of the Cepheid and KR Lyrae variables. Details such as the asymmetry of

the light and velocity cur'.es and the phase shift of the luminosity

ma3.imurn ith respect to the radius minimum were explained using large

land expensive) codes. A number of new problems arose from this work,

two major ones being tne mass discrepancy in Cepheids (Cox, A. N. 1980a)

and the lack or drivlini in beta Cephei variables (Osaki 1982; Lesh and

Aizcrnman 19v'). in the auttept to understand these difficulties, fully

non-adiabatlc, but linear, stability analyses were developed (Baker and


Kippenhahn 1962; Castor 1971; and Saio, Winget, and Robinson 1983).

These analyses have the advantages of being much easier to use and less

expensive than the nonlinear codes. They are easily modified to include

more physics, such as convection (see, for example, Baker and Gough


The linear analysis shows only whether a given model is self-

exciting or damping. An unstable model in the linear regime is a

candidate for further study in a hydrodynamic code. However, the use of

the linear, non-adiabatic eigenvectors in a perturbation expansion is

limited due to their non-periodic nature. To understand how non-

linearities affect an oscillation, without using a hydrodynamic computer

code, a purely real spectrum of eigenvalues is required in first order.

The adiabatic eigenfrequencies satisfy this requirement, but treat the

outer layers in an extremely poor fashion.

As the bulk of the star is quasi-adiabatic, Buchler (1978)

developed a formalism using two-time asymptotic techniques modeled after

Cole (1968) and Kuzmak (1959). This method links the linear analysis

with the nonlinear calculations by providing a consistent expansion in

powers of the amplitude (Buchler and Perdang 1979). The requirement of

periodicity is inherent in an expansion of this type, but the payoff is

the capability to study a pulsation amplitude as it starts at an IniLtal

value and evolves to a final state. Additional problems that can be

studied include modal selection, mode switching, and the feedback of the

pulsation on the thermal evolution of the star, all of which are open

questions in this field (Regev and Buchler 1981; Simon, Cox and ho~d:n


In the simplest version of the theory, the quasi-adiabatic results

are recovered in the zero amplitude limit, but at the cost of

introducing a cutoff in the integrals. When this was done, it became

obvious that these layers were being incorrectly treated and that a more

reliable method needed to be developed. A straightforward way of

including these layers was presented in Buchler and Regev (1982a), where

the structure of the eigenvectors was determined by the condition of

constant luminosity above the cutoff point in the model. While giving

good results for stars of the Beta Cephei class (Pesnell, Regev and

Buchler 1982), cooler stars, like Cepheids, were not treated much better

than with the simpler adiabatic approximation. The present work was

undertaken to quantify the thermal behavior of pulsations. The results

presented are encouraging and explain or clarify several phenomena found

in stellar pulsations.



The structure and evolution of stars is governed by a set of four

equations: the conservation of momentum, energy, and mass, and a flux

equation. Other relationships are constituitive, such as the equation of

state, opacity (see Appendix A), and nuclear energy generation rates.

Written in a Lagrangian form which uses the interior mass as the

independent variable, these equations are (see Clayton 1968, ch. 6;

Appendix B below)

b r3
= 4ip (3-1)

d2 r Gm 4r2 g(r,S) (3-2)
d -t2 r2 4rr2 g(rS) (3-2)
d t2 r2 a

TdS q(p,T) T k(r,S) (3-3)
d t )8 m

L(m) = (4r22 ac T(34)
3< s m3-

The symbols are defined here to be

m = mass beneath the point of interest


r = radius from the center of the star to this point

p = density

T = temperature

P = pressure, P = P(p,T)

S = entropy, S = S(p,T)

L = luminous flux

K = opacity, K = K(p,T)

q = nuclear energy generation rate, q = q(p,T).

We have assumed radiative heat transfer in the diffusion

approximation and that all constitutive quantities are given by analytic


To study this system, the equations are solved with all time

derivatives set to zero to find an initial model in hydrostatic and

thermal equilibrium. The equations are then linearized in r and S and

combined to l w.'e trw operator equations with two unknowns, 6r and 6S.

Written in operator form these are

i- *.:r + "6-fS = A.6r B*6S (3-5)
d k \

S (--i + +r6S C.6r + D6S. (3-6)

it i~ well inoJn that, with the appropriate boundary conditions,

the A operator i1 of the Sturm-Liouville class and yields a symmetric

matrix operator (Ca(tor 1971; Buchler and Regev 1982b).

The operator A is the adiabatic pulsation operator (6S=0) which can

be written as (see Cox, J. P. 1980)

d Ir Ppr ]} + {r3d [(3r1-4)P] pr } 0
dr r r dr dt2 r

and in the Schrodinger form of this equation (Morse and Feshbach 1953,

ch. 6), the new independent variable has the dimensions of time and

represents the acoustic travel time from a point in the star to the

surface and is given by

R, Fr (r) P(r)/2 R
tacous (r)= r dr p(r) = r dr/c(r) (3-7)

where c(r) is the local adiabatic sound speed and R* is the equilibrium

radius of the star. The fundamental radial pulsation period is roughly

given by t acous(Rcore), the sound travel time from the core to the

surface. As the acoustic time corresponds to the time for a dynamic

response to occur, t acous(r) can be identified as the dynamic time scale
as a function of position in a star.

The nice properties of A come from the fact that it is of second

order in the spatial derivatives (the time dependence is assumed to be

exp( iWt )). Looking at the D operator, the same second order

derivatives are present. We now show that there is a representation

where D has the same properties as A, and a natural definition of the

thermal time scale will result.



For an optically thick medium, the thermal operator D can be

written in a continuum form as

d6S 6S dS 1 d 6S/
d 6d 6S d
T = q(pe ) qT C d m m)(4 T) + Lj) d-I d


/I d L
Here q (ln) d-eq = q(p ,T ), i.e. an initial model in thermal
a la\e nT/ 'd m nq eq
balance, T = T1MT C =(-) and all derivatives are evaluated at
\ P v \ Tp
the equilibrium densities and temperature. This form for D has been

linearized at constant density (radius) and is strictly valid at points

deep in the star so that T > 1.

We asume a time dependence of the form exp( -ot ), and note that

the operator Is of second order in the spatial derivative. Therefore, an

irLtegracirin fctor, p exists that transforms D into a self-adjoint

operator. liil factor can be written as

m0 d
,im' = e (4-2)

witerc r, Is .a reference mass, here chosen to be more, the mass of the



core where nuclear burning occurs. Defining a new eigenvector as
d In L
= 6S/C and X = ___ eq, equation 4-1 becomes
v dm

T d cT d In-r d I
T Co = [X(4 KT qT) ] + [ mdl m ]

Appropriate boundary conditions (which make 4-3 self-adjoint) are

d 6L/L b d
d m a m(m,) + b (m,) = 0 (4-4a)


F(m ) = 0. (4-4b)

The surface boundary condition (4-4a) is equivalent to the assumption

that at the surface there is no radiation incident from the outside

(Castor 1971). In this form, 4-4a shows that the thermal readjustment

time of the outer surface is infinitely short. As the weight function on

the left-hand side is positive definite, and the function multiplying

the derivative on the right-hand side is negative definite, 4-4 is of

the Sturm-Liouville class of operators (Boyce and DiPrima 1969, ch. 11).

Comparing 4-3 to the standard Sturm-Liouville equation, since
In 1- > 0 and pTC /L > 0 all eigenvalues are real and
d m v eq
positive. With each eigenvalue on is associated an eigenvector Fn whose

number of nodes corresponds to the numeric position of the {On} in an

increasing ordering of the eigenvalues starting with a0. In addition,


the eigenvectors are non-degenerate, orthogonal with respect to the

weight function pCvT/ Le and can be normalized so that ( m, is the

total mass of the star)

J p T C Si dm = 6...
0 v 1 J L 1J

Due to their Sturm-Liouville character, these eigenvectors form a

complete set for describing the 6S readjustment in a star. In

combination with the adiabatic radial wave functions (Cox, J. P. 1980),

a complete set for describing spherically symmetric disturbances of a

sAir is found, provided that equations 3-1 through 3-4 are always valid.

To understand the structure of the eigenvectors, the approximation

of neglecting nuclear burning can be made. This reduces the equation to

a version suitable for discussing thermal effects in a stellar envelope.

Without nuclear burning the initial luminosity is constant and the

analysis is halted at a mass ( more ) outside the region where this

assumpction i not valid. is the structure of the n 's will show, this is

equi.alent co as5dmanf wee are looking at "short" thermal time scales

chat do not probe the core. Radial stellar pulsations are an envelope

phenomenon, io this aisumption is justified.

I:-electcin nuclear burning and defining L, to be the static

lurninoit c entcri&, the envelope from below, equation 4-3 becomes

ST- -L T n d + (4-5)
Fllo and Pekers d m d 4

Fciloiare Ledou.: and Pekeris (1941), we multiply eq. 4-5 by n and

integrate over the interval m = more m m.. A self-adjoint system

satisfies a minimum principle and any well behaved eigenvector can be

used to estimate the lowest eigenvalue, a0. Choosing S0(m) = 1 for all

m, this estimate is
m, d IT m* d T
J LI(m) dm J dm
m0 dmm m0 d m
S= L L
0 m
a p TC dm
T dm (4-6)

If the thermal time is defined as the reciprocal of o0 this gives the

overall thermal time of the envelope as

t (4-7)
th L m d KT
pI dm
m0 dm
There is an apparent mathematical inconsistency with this

expression. When the opacity derivative, KT, is constant in the

envelope, the thermal time is undefined. However, comparing 4-6 to the

expression for the adiabatic radial pulsation fundamental frequency

(Ledoux 1965),

S 4xr2 dm [ ( 3 rl- 4 ) P ] dm
WO m=

0 r2 dm
they have the same form only if = 0 so that the formula
dP d 1
for tth is incomplete. If --= 0 and -d = 0, then w0= 0, but the
th d m dm 0
variations of P are large in a star, varying many orders of magnitud=,

while rF can vary only from 1 to 5/3 and the neglect of the pressure


derivative would not be correct. The inclusion of nuclear burning in 4-7

remedies this situation and the net thermal time for the star is then

C =< (4-7')
th,nuc m d KT
L* 0f [(4 T- q) d ] dm

which allows for the luminosity variations and shows that the apparent

infinite value of tth from 4-7 when KT is constant in the star is not

real. The term allowing for nuclear burning will always give a finite

value for tth,nuc When KT is constant in the envelope, the assumed

form of & is a poor approximation and gives an unrealistic result. The

posicI.'lt of t is difficult to guarantee, but a simple case would
th.nuc d K
d InL T
require tht - qT > 0, --- m> 0, and -- < 0 all of which can

t-e satisfed in most of a star's mass, although in some nuclear burning

regions, qT can be much greater than zero, and the opacity variations

are .i,:j3-ssed below. The inclusion of nuclear burning as it relates to

scability of Ftellar models is discussed in Hansen (1978) and Gabriel

(19;-'). in this work, this area will not be discussed as the

simplii.:acion that result from the neglect of the core make the

interpretation of the results easier.

the opacity is of the Kramers' form ( K T-3.5 < 0 ) near the

surface but inside the ionization zones. Near the core, electron

sca terling -ill dominate ( KT -->0), so that KT decreases going

outward in the s r. Therefore, except in a small region near the

ionizatlion one=, tth will have only positive contributions in the

inregral. The mass of the ionization zones is small compared to the

total mass of a scar, so that the negative contributions to the integral


will not cause a negative overall time scale. In addition, each point in

the star is weighted by the ratio

m 4 Tm) T()
P(m) T(mo

which gives the interior regions far more importance than the outer

layers, further reducing the influence of the negative portions of the

integrand. Physically tth represents the time necessary to remove the

p weighted internal energy by the radiative luminosity present,

allowing for the position (and therefore temperature) dependence of the

opacity derivative.

Another thermal time can be found from equation 3-3 as the time

necessary to remove the material thermal energy at the static luminosity

present in the envelope. Writing this in a integral form (Saio, Winget

and Robinson 1983) yields

t ,p m J T C dm/L,. (4-8)

Comparing the two definitions, there are several differences. The value

tth,p is counting only the thermal energy in the matter that makes up

the star. In equation 4-7, the formula emphasizes the radiation

( p ~ T4 ) and the method of energy transport, as well as giving the

interior stellar regions more weight than the outer layers. These

differences are minor and are overshadowed by a major shortcoming

of tth. The definition in 4-8 can be generalized to calculate a thermal

time as a function of mass by replacing the lower limit in the Incegral

by the desired mass. Values for tth,p found are always positive and can

be interpreted as time scales. However, when the same technique is done

with tth, the value returned is not positive definite due to the

ionization zones near the surface. Mathematically, this generalization

means that ((m') = 1, m'>m and F(m') = 0, m'
function at m. The derivatives of S in 4-5 will yield delta functions at

m, adding contributions to tth that are not easily understood. To

alleviate these difficulties, tth shall be redefined in terms of the

independent variable that results from analyzing eq. 4-3 by standard

asymptotic methods.

Asymptotic Expansions

Following the technique outlined in Morse and Feshbach (1953, ch.

6) we make the substitutions:

y(x) = [ 2 (dl n

m TC /2
1 _v d InT 2 d
J m6 L. d m

mj TTCv d InT f/2dm
mO m [- d m

k- = J2a ,
2 a

1 d2 d v /4 2 L d KT
'-d---m TCv2 ]
(d 1 1/TC] dx2 C2 TCv dm
\d m/ v

to transform eq. 4-3 to

d -- + ( k2 w(x) ) y = 0. (4-9)

The units of J are (time)/2 so a natural definition of the overall

thermal time of a stellar envelope is

tth= 2. (4-10)

We note that the integrating factor p does not appear in this

definition. As p is an artificial quantity introduced for mathematical

convenience, this is physically appealing.

From the equilibrium model (B-4 in Appendix B below) there is a

relationship between L, and dlnT/dm. Inserting this into eq. 4-10 gives

m, 3 C r
th [ J ( -) 2dm ]2 (4-11a)
t 0 (4tr2)2 ac T3

or, using the continuity equation (3-1), RO = r(m0) and R* = r(m,),

R, 3 C K 1,
tth R ( v )1/2 pdr ]2. (4-11b)
0 a c T3

This definition shows the effect of both the internal energy that must

be dissipated and the method of dissipation. Rewriting 4-11(b) in a

slightly different way

R 3 K T C 1
tth =[ R CT )2 2p dr]2 (-1lc)
R0 c aT4

or, defining Ug = T C ,

R U /r
t = [ I (Kp--gas )dr ]2 (4-12)
th c R0 Urad

we have a interesting definition for tth. The photon diffusion length is

present ( 1/Kp ), so increasing the optical thickness increases the

relaxation time. The 3/c term can be interpreted as the diffusion

velocity of the photon; however, this is discussed below. From the ratio

of internal energies, it is seen that in the central regions where this

ratio is large, the contributions to tth increase. Near the surface,

where the density and opacity approach zero and the temperature becomes

constant, the radiation dominates and the integrand becomes small.

If the approximation of a monatomic ideal gas with black body

radiation is made, a simple relationship exists for the ratio

U / U Defining n to be the reciprocal of the mean molecular weight
gas rad
(or the molar density in moles/gr), the pressure of the gas can be

written as

a 4
P = P + P = nRpT + a T,
gas rad 3

or, introducing p (Chandrasekhar 1967,ch. 6),

P = pP and P = (l-P)P.
rP rad


The internal energy of the gas can be written as

3 4
U = U + U = nRT + aT /p,
gas rad 2

so that

3 3
Ug = P /p B P/
gas 2 gas 2


Urad = 3Prad/p = (1-P)P/p.
"rad rad

With these assumptions, the ratio becomes

U gas/Urad =1/2 P/(0 ).
gas rad 2P1 -

This mixture has a thermal time given by

tth = R [ p P/(1-p) 2dr }2. (4-13)

If P is constant in the envelope

t =-3 R (--/2dr ]2.
th 2c 1 p Rf


Here we see that the thermal time is a diffusion length

1= [ R (Kp) 2dr ]2 (4-15)

divided by a diffusion velocity

v = c (1 p)/, (4-16)

where c is the speed of light in vacuo.

The diffusion velocity must be less than c; therefore, the value

of 3 must be greater than 0.4. This formula shows that in stars where

0 i' of order I, the thermal times are large as it takes a long time for

Spnoton to mo.e through the material. In a star that is radiation

dominated, the thermal time is short as the photons can move freely

[lirroujh the en-elope. These interpretations are qualitative, depending

on the i5sumed nature of U They are suitable for estimating whether
a tcar ujrnergoe; adiabatic pulsations or moves in some way that shows

large thermal effects. A star with a very tenuous envelope, and

cliereforc mall p, would be expected to show large thermal effects.

Star. of this cpe include the R Corona Borealis class (King 1980; Saio

and wheeler 1982', as well as a Cygnus (Lucy 1976; Pesnell 1981).

RecturnLir. to the original form of tth, a position dependent time

scii e car be defined by replacing the lower limit by the desired radius.

The res;llr nt formula represents the thermal travel time from that point

tIo cti Surface and can be written as

th(r R, 3 [ j C /2 p dr ]2. (4-17)
r a c T3

In figure 4-1, we compare this thermal time with tth,p, also

generalized to a position dependent form as discussed above (Cox, J.P.

1980; Saio, Winget, and Robinson 1983). The time scales are plotted

versus the logarithm of temperature (with the surface of the star at the

left of the diagram) for a model of BW Vulpeculae and in figure 4-2, for

a generic Cepheid, discussed in chapter 5. Except for the constant

separation, it is readily seen that the two functions agree quite well

in the interior but diverge from one another at the surface. This effect

is a consequence of tth,p containing only the material thermal energy,

which for a hot star is quite considerable and the drop to zero at the

surface is due only to the integral formulation for tth,p. In the

equation for tth, the method of energy transport reduces to the grey

atmosphere approximation at the surface, and the thermal time must

therefore approach zero more smoothly for any class of star.

This disagreement does not change the normal use of the thermal

time. Several classes of variable stars are driven by ionization zones,

including Cepheids (Cox, J.P. 1980) and ZZ Ceti (Starrfield et al. 1982;

Winget, Saio, and Robinson 1982). For this driving mechanism to

function, the ionization zone must coincide with the point in the star

where the transition from (quasi-)adiabatic (constant entropy) motion to

highly non-adiabatic (constant luminosity) motion occurs (Cox and

Whitney 1958). This transition zone is found at the radius, Rtr'

where tth,p(Rtr) ~ IO, the thermal time is of the same order as tIe
th,p tr 0


period of the pulsation. The expression for tth is derived from a

frequency, and the location of this region would be where
o(Rtr ) = tth(Rtr)- 0= 2n/H0, which is analagous to critical damping.

The constant separation in the figures is approximately 2n; so both

formulations predict the same location for the transition zone, provided

this region is below the atmosphere.

Asymptotic Thermal Eigenfunctions

Once the normalized coordinate has been defined, the asymptotic

form- for the eigenvalues and eigenfunctions of eq. 4-9 can be

3calculated. Here, asymptotic implies that k2 = J2 o > w(x), for all

values cof .:; conSidEred. Due to the surface effects, this is not a good

assu-ptcion, but ai- Important feature in the eigenfunctions can be seen

in trhi appro-:imation; so, it is discussed.

RTe asyivpcotic form of eq. 4-9 is the harmonic oscillator and has

boundary, conditions

Yk'll) + b' d-(l) = 0 (4-18a)

,0) = 0 (4-18b)

where a' and b' are a and b from eq. 4-4 transformed to the x coordinate

jatem. a tmuming ,I.;i = A sin( kx + 6 ) replaces the outer boundary


condition by a transcendental equation of the form

tan( k + 6 ) = b'/a' k (4-19)

In the limit k > 1 this equation has roots that correspond to half-

integer multiples of n (Butkov 1968, ch. 9), k = ( 2 + 1 ) ,

S ( ( + ) n )2/J2, and we shall let 6 = 0 .

The physical eigenvector can be found from the asymptotic

expression as

(x) (x) [12 TC v/d inT II
n(x) = Yn(X) / [ 2 dln (4-20)
n Yn L:)( Ld m ) ( )

In Figure 4-3, we plot the x variations for the Beta Cephei model.

For log(T)<5.7, the shape of the eigenvector is determined solely by the

denominator in 4-20. In this star the integrating factor decreases

outward and the net effect is an eigenvector that peaks at the surface.

(Figures 4-4 and 4-5). The minimum in i at log(T)~4.6 is due to a

maximum in KT at that temperature caused by the second ionization stage

of helium. The shoulder in y at this temperature is also from this


Moving to a Cepheid model, the situation changes considerably. In

Figure 4-6, the thermal coordinate is plotted, and we see that x~l

until log(T) >4.5, well interior to the hydrogen ionization zone. When

the asymptotic (n=10) eigenvector is plotted (Figure 4-7), a large peak

occurs at log(T)=3.9. This variation cannot be due to the cosine term in

4-20 as the argument is essentially 1.0; so, this feature must be due to


the denominator, whose major variations come from the integrating

factor. A large minimum in p can be seen in Figure 4-8, where the

logarithm of p is plotted. This minimum corresponds to the maximum

in and is caused by the large positive value for KT in the cool side

of the hydrogen ionization zone. The integrating factor has a
4 T
temperature dependence of roughly T and in regions where

KT > 4, p increases outward instead of decreasing as it does elsewhere.

This structure corresponds to the peak in the temperature eigenvector in

a standard linear, non-adiabatic pulsational stability analysis. As the

numeric results will demonstrate, any ionization zone will cause a local

peaking of the eigenvector, but a region where KT > 4 causes an effect

that is best removed when calculating thermal effects in these stars,

and plotting the resultant eigenvectors. There are secondary minima

in p at log(T)~4.2 and 4.7 due to the ionization stages of helium. Any

effects these have on the eigenvector y (Figure 4-7) are completely

hidden by the peak from the hydrogen.

However, it is not true that k2 = J2o > w(x) everywhere. At the

surface, T-->TO and < -->0, so that w(x) increases outward. While the

interpretation of the "bumps" in the eigenfunctions as maxima in KT

will be .alid in the actual numeric results, the frequencies and

eigenfiJnct ions are not accurately represented in these expressions. With

crts knowliege in hand we now move on to solving the equations

nIuaericall, and show the thermal modes for several stars.






bo >
^ -I

p r

4J' CO
*ll -

60'6 ILi E S SCE




* H



o >O


4- -0






-h ----------- -- "L-- t-4- -4
4 *t' *.' '44 i,,C

~~_~ ~

r-_ -H-







06_ __ 0 9 : a ho r
/ "'-'


-+ -~--. --


0 0

*3 '-


00:1 06:L 00: L)LO 0:0 0!0 0:0 0:*
31Ui~lOH00 ~UHYFI a




b 0


-O "9 02 "L- on -'


3,^ 4 C friCL C CC


4 r



SJ -,ll ] AHr 3.]. i LS

LL 09 0 l1 rh 0Il ? 1 9


ji. I 1ii








' 12 H






l o

tT s"0a Go'61 00" 1 00 6 0 n" 0"


I ---










-A -'----------
0 o7 0L- 02 09 oh' s0
--^---^-i--c--I- 4




In this section the eigenvalues and eigenvectors for two stellar

models will be discussed. In each case, the initial model was found by

using the model builder described in Appendix B, with mass, luminosity,

effective temperature, and composition appropriate to the class of star

the model represents. One model is of BW Vulpeculae, a hot supergiant of

the Beta Cephei class of variables, and the other is a generic Cepheid.

To find the matrix operator, equations 3-3 and 3-4 are

differentiated at constant radius (density) and written in a form

resembling that of Castor (1971):

d 6S
T d 6S = AK2(T6S). (5-1)
d t

Although the matrix AK2 is not symmetric, it is tridiagonal and a

numeric analog of the integrating factor does e
Applying this transformation to AK2 yields the D operator as a real,

symmetric, tridiagonal matrix. The diagonal transformation m*riy is

defined as

1/2 (1) = 1

2 (1) P 2 (,-1) [ AK2(I-1,3)/AK2(I,I) 2, I=2,N,

and the matrix D is found by using /2 as a similarity transformation

D = 4/2 AK2 L '2

so that the spectrum of D is identical to that of AK2. In a convective

region of a star, one procedure is to neglect the Lagrangian

perturbation of the convective flux (Saio, Winget and Robinson 1983).
This will have no direct effect on the definition of p/2 as the two

matrix elements, AK2(I-1,3) and AK2(I,1), are multiplied by the same

fraction of the luminosity being carried by the radiation. Any gross

features of 1/2 such as peaks in ionization regions, will still be

present in a model where convection is permitted, although they may be

extended over a larger region due to the smaller temperature gradient.

The transformation has one side benefit, the matrices analogous to AG2

and AK1 in Castor (1971), are much closer in magnitude to the elements

of the symmetric matrices A and D and have smaller variations as a

function of position in the model, resulting in better numeric accuracy.

Real eigenvalues are guaranteed by the symmetric, tridiagonal form

of the matrix operator (Ralston 1967). For negative definite eigenvalues

(representing decays) the diagonal elements must be negative. This can

be seen in a simple example.

A negative diagonal element means that increasing the entropy

temperaturee) in a zone increases the luminous flux out of the zone


more than it decreases thi flux entering it. This stabilizes the zone as

the additional energy is lost and the temperature decays to its initial

value. A simple example can be used if the luminosity gradient had a

local relationship. Equation 5-1 would then be written as

-= -a(S S0

so that

S(t) SO = Sat

If a>0, the entropy decays, but if a<0 the entropy increases. The latter

situation is possible in a nuclear burning core but not in the envelope.

To preserve the Sturm-Liouville character of the continuum operator

in 4-5 in the matrix operator, the off-diagonal elements must be

positive. This condition is composed of two separate requirements. The

luminosity exiting a zone must increase if the temperature of the zone

above is decreased, and the luminosity entering a zone must also

increase if the temperature of the zone below is increased. Both of

these are easily satisfied if the opacity is independent of temperature

or decreases with temperature, ( KT < 0 ). From the form of the

continuum operator (4-5), only if KT > 4 can violations of thEes

conditions occur and these regions are precisely where the interpoiatlon

formulae for the luminosity are essential. In the initial model the

temperature ratio is not allowed to exceed a given nuriber in che

ionization zones. This process limits the variation of [he opacity 3a


well and raises the upper bound for positive matrix elements. A first

order expansion of the super-diagonal matrix element shows that for

KT < 1/ln(T(I-1)/T(I))~17 ,

these matrix elements are always positive. We have never found any

negative off-diagonal matrix elements in any calculated envelope model.

This process also limits the variations of the matrix elements so that

the sum of the off-diagonal elements is less than the diagonal elements,

permitting the use of Geschgorin's Theorem (Wilkinson 1965).

The boundary conditions are incorporated into the matrix operator

by giving certain elements of the matrix special functional forms.

Assuming the outer surface to be radiative and given by the grey

atmosphere formula is equivalent to assuming that this surface has an

infinitely short thermal response time. Therefore, the analytic form of

this boundary condition is that of eq. 4-4, even though the diffusion

approximation is not very accurate here. For numeric computations this

expression is unsuitable, as it is valid at a point exterior to the

actual surface of a model, and we have assumed that the exiting

luminosity is given by

L(N+I) = 4tI R2(N-+) T4 + 2
eff 4 3

Whin this formula is linearized and inserted into the operator equation,

the outer equation reduces to one of first order and closes the system.

The inner boundary condition is 6S/C --> 0 and is included by zeroing

any elements in the matrix that refer to 6S/C (0). Again, this reduces

the appropriate equation (here the innermost zone) to one of first


As the thermal time of the outer surface is zero, and the overall
thermal time ( ~ 00 ) is normally longer than a pulsation period, the

spectrum of D must contain eigenvalues which bracket the radial

fundamental adiabatic pulsation frequency. Here is the information that

we seek regarding the transition zone. All parts of a star that can

relax faster through a thermal process than acoustically will move in a

highly non-adiabatic fashion. These regions of a star will be called the

"sudden" regions and correspond, in some stars (see below), to the outer

layers where d L/ = 0 in a fully non-adiabatic stability analysis.
d m

The inner regions of a star have thermal response times longer than the

fundamental period ( o < 0 ) and are referred to as quasi-adiabatic.

Eigenvalues for this matrix can be found by a variety of means. For

the low order modes iterating on the outer boundary condition is

possible (Castor 1971), but the structure of the high order modes is

such that this method fails. A general eigenvalue routine (Smith et al.

1976) is used to find all eigenvalues of D and the necessary

eigenvectors are found by standard Gaussian back substitution. The

eigenvectors of D and AK2 are related by the same similarity

transformation as the matrices and the more desirable set can be

calculated, depending on the effective temperature i T ). of the star

under study. If Teff > 15,000 K, then the hydrogen lonlzal1,,i, zone is

absent and the En's give a good representation of cne rhermal modes.

When Teff < 15,000 K, this ionization zone causes a spike to be pre.-ent


in Fn as shown in the section on asymptotic eigenvectors in Chapter 4.

To remove this feature, the eigenvectors of D are calculated and

displayed. This results in a much clearer picture of what is occurring

in the surface regions of the star.

We cannot look at the long decay rate modes as we have removed the

core and the effect it has on the eigenvectors. The lower modes can have

significant amplitude in the core, but the higher modes are effectively

isolated from the core even more so than the radial pulsation modes.

mhit .zauieo the dlftfcultv roced earlier, the number of nodes in a

caliculated eienv.'ector is no longer an indicator of it's position in the

pectrum: the i:ze of tr-e nuoDrs causes the computer to underflow and

the interior p.:rtion of the elaen.'ector is set to zero! Physically, this

impilee thltc cre e'.ter.r canr react to perturbations on much shorter

time icjles than the interior. Ho.weer, the outer layers are linked

quite strjonell, t- tcer interior through the low order modes. When a heat

pulse is lnIerted at the baje of the envelope, the entire model is

affected. The outer la.ers react to maintain thermal balance with the

neu lui;n.: .itL from the heat pholie. The eigenfunctions of the diffusion

equJatli.n e>xhibst an infincte propagation velocity for perturbations

Iliore and fel.hbach 1953, ch. 7), causing this behavior. The asymptotic

cl.genfInctiuonl are trigir)onetric wjth an amplitude that increases

outward and we can therefore speak of the thermal time of a point in the

.tar a- the reciprocal of the largest eigenvalue whose eigenvector has

Aigntilcant ampittude at that p,-ition. This may seem ambiguous, but as

thermal mjadei canr be used to e pand the pulsation equations in the

therral tiCesZ, tniS identtficatrion will enable us to study the thermal


effects in an oscillation. In particular, the quasi-adiabatic stability

coefficient and pulsations that are dominated by thermal effects can be


Concentrating on those modes with have eigenvalues comparable to

that of the radial fundamental pulsation frequency, we now examine two

stellar models.

BW Vulpeculae

The physical parameters for this star are plotted in Figure 5-1 and

the eigenfunctions 6r/r, 6p/p and 6L/L for the radial, adiabatic,

fundamental pulsation mode are shown in Figure 5-2. The results of a

linear, non-adiabatic radial pulsation analysis (LNA) are shown in Table

5-1. To find a realistic model for this star, it was necessary to

integrate into the nuclear burning core slightly, but the thermal modes

discussed below are not altered by this.
-4 -1
The fundamental frequency (w0) is 3.9 10 sec and the four

thermal eigenvectors with eigenfrequencies closest to this are plotted

in Figure 5-3.

Important features are the cutoff near log(T)=5.2 and the two

peaks. The peak near log(T)=4.6 is due to the second ionization of

helium causing a local maximum in KT. At log(T)~5.1, a -ecoriar,- peak

associated with the first due to the use of the Rosseland .eain opacity.

The weighting function in this mean has a maximum at a temperature

roughly four times that where the ionization occurs and cause a3 ma limum

in KT at this temperature. This peak has been studied by tellingwerf


(1978) as well as Saio and Cox (1980), as a possible driving mechanism

for these stars.

The luminosity variations ( 6L/L ) for each mode are shown in
d 6L/L
Figure 5-4. Near the surface dm = 0 showing the domination of the

short thermal times in this region. Moving inward (increasing the

temperature), the luminosity variation starts to decrease at log(T)=5.0,

oscillates due to the nodal structure of , and drops to zero by

log(T)=6.0. These modes are completely isolated from the nuclear burning

core (log(T)>6.5) as stated earlier.

Figures 5-5 and 5-6 show the real and imaginary parts of the

entr pc.p variations ai calculated from the LNA pulsational stability

inali t. The pcakj at log(T)=4.6 and 5.1 as well as a sharp cutoff

abofv 1lo,(Ti=5.3 car. be seen in both components. Each feature has an

iana Ioa in the Lhcrml L eigenvectors in Figure 5-3. The modulus of the

luLmnosit .arijatcion from this analysis is shown in Figure 5-7. The
d ,:-L L
change from d = 0 is present, albeit at a slightly higher

cemperatire chan it, Figure 5-4. Some of this is due to contamination

wicn low order cnearmil modes, but the significant amplitude in 6L/L (up

to iot>T,1=1.)l ii due to the adiabatic density fluctuations in this

rec ion.

Ti-s i.model siho-s that the shape of 6S/C in an LNA pulsational

stability analvyis 1i determined by the thermal modes of that star. As

always, [he model is stable and, at present, we do not find any new ways

of dcscablitzinr IL. The structure of the eigenvectors serves to

demonrrace and clArlaify E[e physical reasons behind a transition zone

ana the spilttins of the itar in two regions, adiabatic (log(T)>5.4) and

instantaneous thermal balance (or sudden) at log(T)<5.4. More

interesting behavior can be seen when a model of a Cepheid variable is


Cepheid model

This star is a generic Cepheid whose characteristics were obtained

from A. N. Cox (1980a). The results of the LNA pulsational analysis and

summary of the input data are found in Table 5-2. A graph of the

physical parameters is found in Figure 5-8 and a graph of the adiabatic

eigenvectors 6r/r, 6p/p and 6L/L is shown in Figure 5-9. The model is

unstable in the first three pulsation modes and corresponds to a model

of the star DL Cassiopelae. The results presented here are illustrated

very nicely by this star and are similar to those of other Cepheids


In Figure 5-10 and 5-11 the four thermal modes that bracket the

pulsation frequency 0 are plotted. These are the eigenvectors of AK2,

and we see almost no detail (or the individual modes) due to the spike

in the hydrogen ionization zone. The origin of this spike was discussed

in the section on asymptotic eigenvectors, and as it eliminates any

other features present, we choose to remove it from the eigenvectors.

This is accomplished by using the eigenvectors of D, as shown in Figure

5-12. These vectors have the same eigenfrequencies as those in Figure 5-

10, but much more detail can be seen. The bumps at log(T)=3.82 and 4.12

are due to the hydrogen ionization zone. They are the lower slopes of

the spike, and the minimum at log(T)=4.0 is a result of modulation


effect of the integrating factor being a maximum at this point. The

helium ionization zone is visible at log(T)=4.6 as an inverted peak, but

the second peak at log(T)=5.1 is absent. From Figure 4-6, where the

thermal coordinate for this star is plotted, we see that the value for x

starts to drop away from one at this point. In the interpretation of

features for log(T')<-.7. the neglect of the trigonometric function

in r 1t valid. Dut interior to this, the actual nodal structure

of r will start to d>rmin3ie.
Examining the luminosity variations which correspond to Figure 5-12

(Figure 5-1). the iriv.irig mechanism responsible for destabilizing this

tlar is visible .se below). To accentuate this, the variations

for o < h.'- are sh,;n In Figure 5-14 and 5-15 and a > w0 in Figure 5-16

and 5-!;. In the low frequency results, the luminosity looks much like

that of BW v*l ( Figure 5--); increasing the frequency to those in Figure

5-lr, ani 5-17, the dccrea3e in 6L/L near log(T)=4.6 has become larger

than any other part of the vector. Here is the physical reason for the

drn ing zone of this jtar.

For a scar to be unstable in a pulsation mode, a region where
d -UL
d < mJiit exitr (i :a, J. P. 1980). This is seen if the work

intteral is wricr In itr e form (Aizenman and Cox 1975; Buchler and Regev

1982 b,).

[ d m. 6L/
) f 1- U) () (d L/L dm .
3 P d m

The quan.,tity r3- I is always positive and for the fundamental mode

< (jas -- iricreae. outward) (see Figure 5-2); thus, when
p r

d 8L/L > 0 the star is stabilized and for d L/ < 0 the star is
dm dm
d 8L/L
unstable. In the BW Vul model (Figure 5-7) dm is zero in the

atmosphere and 6L/L generally decreases with increasing temperature. The

thermal modes have the same basic luminosity behavior (Figure 5-4),

showing that this star is stable. In the Cepheid model, the pulsational

6L/L (Figure 5-26) shows a definite region near log(T)~4.6 where
d 6L/L < 0 This outward decrease of 6L/L cannot be due to the
adiabatic density fluctuations as they will always have 6L/L increasing

outward. The thermal modes show a change in structure as the decay rate

is increased through the value a~". The change in the shape of 6L/L in

this region is the cause of the destabilization of the CGpheid.

Calculations of the BW Vul model with o/u < 1 and w/o < 1 do not show

the same change. Examining Figure 5-18 and 5-19 where the thermal modes

with /l~1lO, and Figures 5-20 and 5-21 where /o~l10, we do not see the

same change in the structure of the thermal modes. This leads to the

conclusion that the star's thermal structure prevents this type of

driving mechanism from functioning.

The location of the transition zone in each model should be

discussed as well. The BW Vul model shows a definite change from

adiabaticity to sudden behavior at log(T)=5.2. The Cepheid shows a less

dramatic change, as the thermal eigenvectors do not show an extremely

rapid cutoff at the transition. However, unlike BW Vul, the Cepheid has
d 6L/L
a region where d L < 0 that begins near the transition. Also, it
would seem more difficult to establish a transition zone for this model
d 6L/L
as d is not zero until well above the helium ionization. In fact,
the luminosity variations from the pulsation analysis (Figure 1-2bi do


not have a vanishing derivative until above the hydrogen ionization

zone. It is difficult to use the timescale arguments to locate the

transition zone for such a model. The use of a w0 to partition the

thermal modes is, in this case, a fairly reliable way to locate this

region. For o < wo the thermal mode has a 6L/L vector that is, in

general, increasing outward and has only a small dip in the helium

ionization zone. When a > w0, 6L/L is dominated by this region,

until on > w when the amplitude of the eigenvector (6S/C ) becomes

small enough that the mode is isolated from that region of the star.

To see the benefits of the integrating factor, the real and

imaginary parts of the LNA pulsational analysis are shown in Figures 5-

22, 5-23, 5-24 and 5-25. In Figures 5-22 and 5-23, the eigenvectors of

the AK2 matrix are shown, with the "spike" due to the hydrogen

ionization zone dominating the graphs. When the eigenvectors of D are

examined in Figures 5-24 and 5-25, the variations have been greatly

reduced. Lb 5-22, (6S/C )real varies by two orders of magnitude and the

imajinar', part changes by almost three. The eigenvectors of D vary by

less than one order of magnitude and the details of the helium

iont at3ion an other structures can be seen.

The eigenvectors shown in this section have some numeric noise that

*jc3cr'c: coouient. The system that is solved here is identical to that

used in tre L1IA stability analysis, and any noise in the thermal modes

is also present in that analysis. The spikinesss" of the eigenvectors

pre-enc in Figure 5-14 can be eliminated by putting more zones in the

model here tr250) but this puts the eigenvalue routine at a loss for

aolutlons. The envelope mass can be decreased, but at the expense of the

accuracy in the pulsation modes. The basic features of the thermal modes

are well represented and accurately calculated. Interior of the

transition zone, the adiabatic density fluctatiuons will overwhelm the

small luminosity changes due to the entropy variations.

The optimal numeric solution to an eigenvalue problem is found when

the matrix elements are approximately equal in magnitude. The

transformation of the mechanical variable introduced by Castor (1971),

6r- XO = /DM2 6r does this for the adiabatic wave equation. In a model

where the shell mass is a geometric progression increasing inward, this

transformation is equivalent to each zone having the same acoustic

traversal time. The sound speed is roughly proportional to the square

root of the temperature and the thickness of the shell is roughly

proportional to the mass; so, the integrand in 3-7 has equal

contributions from each zone.

This transformation has the unfortunate side effect of decreasing

the numeric niceness of the thermal operator, D. The thermal diffusion

velocity (4-16) decreases as the sound speed increases. To optimize the

matrix D, the mass of each zone should decrease inward, the opposite

direction required of A. This tradeoff between the adiabatic and thermal

operators will be important only if the mode is dominated by thermal

effects. For quasi-adiabatic modes, the interior is dominated by the

adiabatic density fluctuations and the inaccuracy in the thermal modes

is unimportant. We are working on ways of circumventing this difficulty

so that both operators are represented as accurately as possible.

/ /








., ) H



n 3 r


I :I"IL I WIl.lnl I II '. 'I" l IL II



( t

oo inu


/ 2 2VT LSD '0- C6,0- 9911- AD 2- 61 --
i C **"

f m ,0
I / M)^

I I =r *r
1 / ..CT U -

1 / T-l ^

1~, /u nl
1 / U t
\ / ,

I / a

I~ / 1

S~OL3 nN31" 1 t-JI0


Oh'o LEO L2iU

., C *

r- o

*S C

m HI


o m

m I

Cd -'



uv I ~u ;



S l-I
ul <



) l

"-------------------- ------------I------------ tt0 El

". 01

l 0)
1 1 1- -+-------- 2 0H
0'"h [r4 0' hy2 00'? ts!i 6 09!0 20!) L 0 0- r 0
1' ig

*ri H




-4 U


0 ,
u a

00 0





Ul 4

0 3H
i- *
a a)




0 7

a *




,- -,

'I. l I II


:3 4-J





0 U

t C
0 1




11 901 ONU






U 1


, M
u) H
n -


CL (

1111 1





w 0
C) *H
z0 0
o cC


C )r


EI"2eh C5-ELt OCCCE Ch'LEE LZ'PE8 COCIh Ee'C)I 29'Chi OCC'96 BI'Ch *' '1


L n





., 3"


w) C


0 ,0

0 1



1]v 0:0


h;E LI E OE : hhf2 1012 OL:T EE I

09*a ES2O E['


T ,


Elu i

o 1


i 0

U ca

) U

El t>
II- l
n -
a) o
O; 4-1


-- I

i 1. I I I i~

, '

.,',1 -*
1 I

r ~
.L .


~J 1
? ii i



L" --


,--;- '-

CuoE S92'E ;3 C6UI h'l

------ H ~- t----- -
OL O 90O0 o0'- 22'-
A l.'q



, u

C*'-1 '


0- >r i
.. I


o C
Ca) .



K r-1

.T _

--3; ~--
~ ..:;. ~~. ~. .I: I
__ I. r
-_ --~~r ----
;-----. .

---~--; ----I--






m I
O w

,| "




+ -4

CP 4 1cI

s ^ .; or
I~oe I-'

i'ii- n
i *r
(" a g 5

/ '' 5
i a n
/ -~---~--- t---------- ^

[ r~in n'in o




0 u
7 .0
1-4 C*

. .,

0- -- __
-: "

C ) I

6h' V 0" 29'2 912 1 Ou I 99"0 F.'o 1'0- 6 414'a- F.E
--CT*H tU
is s-i
Z 3 t
I M p
(1; 01
U u *

(U -0 /-
in e !

6h'E~~~~~~. (0' I9Eh' EI 19L h[ Od hd LI
h3 >l^


ai -o

P c


S 0



S .

, wI O-.9 aT iq-sg- g r i- a)Qn}- )s*Lr?- ra "a rven I
^ a 73
( -


r 1

*=G__^1;'.. 1

"~"'" ^s'3s^"w ^^ ^.^ ^, __ a j >

T3> ")^' 0; 4-' n

C.o U
U a

C,- ,- Ir

W-U *


: Ul,
La 0

a) 0 l

11 1 Ira

El" U -

Y-^s~f'^^" -' '

""l/'n "
; *'

L .-

G'O z['- E9'0-


T [f;n



d) M
w I


0 0

0 !


0 -4

uT! r

r- *
* ir-l


- .-. -

i-;--- -I------- -H----1--

Lh '* -o~L

o n


S o0
14 0
> In

S 0 -C

4 -

S0 '
H )
- '-^ H *

Lh'O g9cD 92 a'

so' 9 I '=



-I ,

c .

S '
4- "oa




0 0






0 n

9C 1 bL*6- W12- 9 4(h's- woi0- -L'21 f
C )

0) .



w 0




a o

I u 0

a um

A 1/SeJ -t s u 1 2 W cc 'EL 2 a



a, o



0 0

C :

o pc



S 0 Ifl On "0 [j 90" 0E0 '00 I0
f (
m a~
\- S m

(Z~ ~ ~ ~ QS OIY [^-?1 l


?0'0 20'0-
I -n ll I

0 0


0 u




O "

0 r,


m w
E o

H .
* H C

C1 O


a-' 1-


f- 4---- t--'--c--- _______
0/C' )C'C $4- Ci/' f-C)



0 0
0 o


0 -

I 0

h L 2 ..,[ :
// "' L

4J 44V
N.' c6

<^ -- ft


Table 5-1

BW Vulpeculae

M=11.0 M

Tff= 25120 K

oI= 0 193 ~n=-1.9 10-4

nIl= 0d147 '1=-2.3 10-3

II2= 0 118 ,2=-1.5 10-2

Composition X=0.70, Y=0.28, Z=0.02

Envelope mass = 5.02 M

Table 5-2

DL Cassiopelae

M=5.69 M

log(L/Lsun )=3.57

T = 5848 K

110= 7 87 0 = 6.8 102

nI= 5 64 n1= 1.9 0-1

nl2= 4d28 2= 6.2 10-2

Composition X=0.70, Y=0.28, Z0.02

Envelope mass = 2.76 M



To understand the utility of the thermal eigenvectors, we shall

look at various forms that the linearized pulsation equations (3-5 and

3-6) acquire when the n's are used to expand them in the ratio of

thermal time to pulsation period. Eliminating 6S one can rewrite 3-5 and

3-6 as an integro-differential equation:

w2 6r = ( A + B i1--- C ) 6r. (6-1)
i) D

we have an eifenvalue problem for w2 and 6r. Eliminating 6r 3-5 and 3-

, can b. written in a form suitable for studying the secular behavior of

i.^: = (D + C ---B) 6S. (6-la)
W2 A

The i,. 0 cerm is a resolvent that can be expanded in the complete set

or eieenveccors of D. Adopting the Dirac (or bra-ket, see Davydov 1965,

ch. 5) notation of quantum mechanics, this set is { -0 ,I n > }

.n.d equaijon r -l becomes

[n >< Snl
02 6r = ( A + B C ) 6r .

Next, the mechanical eigenvector is

and eigenvectors of A, { 21 yi >

adiabatic effects on a mode that is

expansion will be written as

expanded in the set of eigenvalues

}. We would like to examine the non-

primarily adiabatic, and this

6r E I y > = I i > + ji j I yj >

w = ). + Ai


Inserting this expansion and multiplying eq. 6-2 by < y give

(1 + Yjiajl 2 ) 02 = < yilAl Yi > +

< Y, IBI 5n >< n ICI yi>
n i) + o

< Yj BI Sn >
j n(j i + o


"i2 I yi > = A I yi > .

< yi IBI n >< n ICI Yj >
i + 0 j

< Yk IBI >< n CI Yj >
kj' iw + a


This expression is rather cumbersome, but following the normal quantum

mechanical perturbation scheme, we evaluate the change in the

eigenfrequency using the initial eigenvector. Choosing the adiabatic

mode as the starting point, the quasi-adiabatic correction can be

written as

< IBI n >< nICI Yi >
t2-_a002 = n -(6-4)
1 n id) + o

which is Identical cit equation (26) in Castor (1971). The normal quasi-

adiabtiL: approximation is a /w. / 1 for all possible n, reducing 6-4
0 1

-'- = < Yi |BC y. > (6-5)

I I = I.

-II iI'n

Tri; e:rpresio:.n fr lth quasi-adiabatic correction can be further

as3lzed it .-1V C I. Then w2 = w + 2w.Aw and eq. 6-5 can be rewritten
1 1

AW = I1 iBCl yi >
2 W2 1

a result identical to the Ti (equation 8, in Buchler 1983). Here, Aw is

purely imaginary, as all quantities on the right are real, and generally

shows that the star is stable (iAw < 0) as the diagonal elements of B

are positive and the diagonal elements of C are negative.

This recapitulation of the quasi-adiabatic analysis serves to point

out that by knowing the structure of the thermal modes, a better

approximation can be made. Partitioning the thermal modes into two sets,

S = { o n n = 1, NT : O i < I }


S2 = { ,n n = NT+ : 0 ia < 1 },

eq. 6-2 can be rewritten as an eigenvalue problem assuming an adiabatic

starting point only to determine NT

2 I y > = A y > +

BT B n >< n n ICI y >

N BI u >< I CI Y >
l iu + 0
n=NT+1 n

The second term on the right is the quasi-adiabatic correction. rno .iLh

an appropriate cutoff. However, the cutoff is not in the spatial


variable, as the integrals implied by the notation extend over the

entire mass of the star. The correct cutoff is to stop the sum at the

thermal eigenvalue whose magnitude is just below the pulsation

frequency, as the sum over the set Sl is doing. For this reason, the

second term should be more correct in the calculation of this part of

the stability coefficient.

The third term produces a correction in the lowest order, and

therefore cannot be treated as a perturbation. To see this effect, the

denominators of the sums are expanded in a binomial series and eq. 6-6


w2 I y > = A y > +

1 T BIn Sn IC y >
it I + a /io +
n=l n

N BI &n >< X n I CI y >
a (1 + ) (6-7)
n=N +1 n n+

)2 1 y > = Al y > +

1 nT a
i T n >(1- )< n IC y > +

I BI n >(1 )< X n IC! y >. (6-8)
n=NT+1 n

The first term is the quasi-adiabatic correction already discussed

anld cr.e second is the effect of the outer layers on the pulsation. The

expansion is in /a n, and in the limit l/on -->0, or very short thermal

times; this sum becomes

I B Cn >-o n IC| y >.
n=NT+1 n

This term represents a zero-order quantity and should not be neglected

when the initial eigenvalue is calculated. If the thermal modes had a

theta-function behavior at the transition zone, the set S2 would form a

complete set in and of itself for describing the thermal response of the
outer layers of the star on time scales shorter than w The actual

eigenvectors have a fairly sharp fall at this point; so, a first try at
including these terms is to write the sum over states as D' where D'

is the outer [(N NT)*(N NT)] submatrix of D that includes only those

zones with I > NT Substituting D' for the second sum over states gives

2 I y > = A y > +

B n >< n ICI y > -

BI(D')-IIC y > + iwBI(D')-2CI y >. (6-9)

This equation can be expanded in the ratio of the proper time

scales; the second term on the right al= o/o0 and the last two

in a2= 0/0 Writing a=al=a2 to give the simplest expression and

arranging the equations in a hierachy of terms by order in a g i:

0 Y > = ( A- B(D') C ) Y0 > (6-10a)

Sy1 > + 2w0w Y0 > = ( A B(D')-C)I y > +

-l 0Bl ><^ |C| yg>+
%=T B n n CI YO > +
0 n=l

iwoB(D')-2 C y >; (6-10b)

I y > = YO > + a I Yl > + a2 y2 > +

w = 0W + aw + a2w2+ O(a3).

To eliminate secular terms in the equation for yl >, the scalar

product < yO I y, > is set to zero. For convenience we shall set

=1. This gives the initial eigenvalue and the first order
ztabililt c.,efficient as

..= < yO IA B(D')-1CI Yo > (6-11a)

1 NT n ICI y > +
2w6 n=l

1/2< Y IB(D')-2CI YO>. (6-11b)

At any point the matrix D' can be replaced by the sum over states that

it appro ii,.actes

0 I(A + X BI n > < n C) yO >, (6-12a)
n=NT+1 n

S xNT < Y0 IBI n >n X n 1CI y0 > +
2w20 n=l1

1/2 < Y IB > --< n |C| yo>. (6-12b)
n=NT+1 02

Equations 6-11a and 6-11b agree with equations 25 and 39 respectively in

Buchler and Regev (1982b), who derived them from physical arguments

regarding the constant luminosity of the outer layers. It should be
d 6L/L
noted that the assumption of = 0 is not valid over the entire
section of the Cepheid model above the transition zone. The sum over

states is valid, and comparing the thermal mode 6L/L to that of the LNA

analysis (Figure 5-16), we see that the functional dependence of the

LNA 6L/L is reproduced by the thermal modes.

Using the infinitely sharp cutoff limit, the eigenvectors of 6-10a

can be calculated in a straight forward way. Defining the error in yO >


E = y(I) r(I) 122.

where 6r is the real part of the mechanical eigenvector given b7 the LNA

stability analysis and both eigenvectors are normalized :o units in

their respective scalar products. The error in the eigenfrewuencv is

defined as


These quantities are plotted in Figure 6-1 against the location of the

transition zone, log( Ttrans ). A sharp minimum is seen at

log(Ttrans ) = 5.23, corresponding to the location of the actual
transition zone as predicted by the luminosity eigenvector plotted in

Figure 5-6. In this model, the change is rather abrupt, falling a factor

of 4 when log( Ttrans ) changes by -0.13 or +0.04 The change w0 is not

as dramatic, reflecting the "minimum" principle that a quasi-adiabatic

oscillation will follow.

When the same quantities are plotted for the Cepheid (Figure 6-2),

it is obvious that the agreement is not much better than using the

adiabatic eigenvectors. The error is reduced at a point where the

transition zone is located, but not by a factor of 10 as in BW Vul. This

results from the more complicated behavior of 6L/L in the Cepheid (see

Figure 5-16). The thermal modes have a large peak in the ionization

zones due to the integrating factor's temperature dependence but do not

have the rapid decrease inward of those in BW Vul. Therefore, the

approcimation of replacing the sum over states by the submatrix D' is

not a rooa one for this class of star. As the sum over states contains

the necesr ary information about the luminosity variations in the

eitericr, including the integrals as in 6-12a should give much closer

agree.~ienc. This problem is discussed in Pesnell and Buchler (1983).

Comparing the initial eigensystem 6-12a to the adiabatic wave

equation produces three differences:


1) There is, at present, no guarantee of either the reality or

positivity of the eigenvalues of either 6-11a or 6-12a.

2) The eigenvectors < yo I and I yo > are not necessarily

related by a transposition, much less being equal. The matrix

A- B(D') C is not symmetric but is of an upper Hessenberg

form whose only advantage is being diagonally dominant. The

eigenvectors come as right and left hand vectors and must be

calculated separately.

3) There is no proven relationship between the number of nodes

in I Yo > and the numeric position of w02 in an increasing
ordering of the spectrum of A B(D') C .

However, for stars which are predominantly adiabatic, as most

pulsators are, the eigenvalues are close to the adiabatic frequencies,

the second part of the operator having little or no effect on A. The

transposed eigenvector has a serious numeric problem. In the region of

the transition zone, the transposed matrix has a structure such that a

discontinuity appears in < y0 The use of eq. 6-12a instead of 6-11a

will spread this feature over several zones, removing the discontinuity.

When the transposed eigenvector is calculated using the LNA analysis,

the feature is present but is much smoother. A change in the method of

calculation is indicated, not a revision of 6-12a.

A final observation regarding 6-12a concerns the limit Ni -->i.. A

a normal star is almost adiabatic ( NT -->N), the opposice liiEc must

describe an abnormal star. For the condition NT --.>0 o be satisfied,

the non-adiabatic regions of the star must be ectensl'c, coverlri, all of

the star that is oscillating. These conditions may be found in the R


Corona Borealis stars or other highly evolved stars where the envelopes

are tenuous and have low gas to total pressure ratios. In these stars,

trying to recover the adiabatic limit in one obvious way will lead to

some rather strange behavior of w2.

Introducing an arbitrary scale factor, X into equations 3-5 and

3-6, which represents the ratio of the pulsation period to some overall

thermal time, they become

w26r = A'6r + B.6S (6-13a)

iw6S = X(C-6r + D-6S). (6-13b)

For most stars, when -+O the adiabatic limit is recovered as 6S is

forced to 0. However, in the integro-differential form, splitting the

thermal modes as before,

M2 I y > = A y > +

-w 1 + Xo /i +
n=1 n

N B >< I Y > (6-14)
n=N+1 ni + Xn

Tne .alue of f:_ may depend on X for most stars, but for the radiation

domninsce. .acmospheres, the thermal response time is very short for all

mode= who=e cutott in amplitude is above the normal transition zone.

I.Ising the ideal eas as before, the atmosphere should have a low density


gradient and opacity gradient so that the diffusion length (4-15) is

large. If the outer modes have o n/ > 1 for all n > N then N is

almost constant during this discussion.

In the limit of X -->0 ( NT constant), eq. 6-14 is expanded using

the binomial theorem as before,

I2 y > = A y > +

BI n > < n I|C y > + O(X) (6-15)
n=NT+1 n

the B(D')- C term does not disappear. For highly non-adiabatic

envelopes, this term can dominate and change the spectrum so that it no

longer reduces to the adiabatic case. This would, in part, explain the

modes that have been named "strange" modes (King 1980; Saio and Wheeler


This explanation of the strange modes does not require that a mode

be near (in some sense) to an adiabatic mode. There are no theorems that

prove that all of the eigenvalues of 6-12a are real, and starting from

an adiabatic mode, the non-adiabatic corrections could move the

frequency far away from the initial value. It is IncerestLrni chat

if NT can be considered constant, the adiabatic limlt is never recovered

in some models. This property is just what a strange mode show.

Modes that possess the property of NT being independent of shall

be called "sudden" modes. They are characterized by 3a ondiciori of

thermal balance existing over a large fraction of the mechanical

eigenvector's amplitude. The quasi-adiabatic scabilic, coefficient would


be a poor approximation for these modes, as the second term in 6-12a

will dominate w1. We must stress that the condition of thermal balance

does not imply that the luminosity has a zero derivative in these modes.









0 c
S' 1------ -

on o h o L- 't- o 2- O 2 0 -
^^. a) a

i., .i : .. *., ... . . ,.: *4 os co
i1, .i II I 1 ,111




--) P.

w C
1 I I

) 0


0 >

X C.
-4 44


1 0


s e

09-6- iPgo- L.d- h ii- Z o- 1 '- 9;
I' Jr




a ,.

) *


c 0





0 11
\ 4J (U


0 O
'., 4 1
*i-l o

H 4-
n) n




5 A


Q) 0


0, a)-
4 FJl -

^ 4-1

^ (1)


We have shown that the thermal structure of stellar pulsations can

be represented by the eigenfunctions of the diffusion equation. As this

equation is of second order in the spatial derivative, it can be mapped

onto a Sturm-Liouville operator, with the concomitant theorems regarding

the reality of the eigenvalues and the form of the eigenfunctions. This

transformed operator defines a thermal time scale through a new

independent variable with the units of (time) 2 This time scale agrees

with that of earlier authors in the interior and predicts the same

location for the transition zone in a star.

The thermal modes that arise from the eigenvalue problem can be

used to analyze the thermal effects in a first order perturbation

expansion. The traditional adiabatic wave operator is replaced by an

operator that includes the knowledge that the outer layers are in
d 6L/L
thermal balance, and that this does not imply that d 6 = 0 in these


Partial spectra of the thermal modes are presented for two

realistic stellar models, along with diagrams of the eigenvectors.

Mathematically consistent expansions of the pulsation equations are

presented that give the quasi-adiabatic stability coefficient and the


analogous expression when the outer layers are correctly accounted for.

It is shown that a new type of mode, the "sudden", mode can be defined

using these expressions.

The driving mechanism in Ocpheids has features in the thermal modes

that do not show up in a star of the Beta Cephei class. This casts doubt

on the viability of ionization driving being responsible for these

stars' variability.

In an aesthetic vein, a way of displaying the entropy or

temperature variations for stars that normally have a "spike" in the

hydrogen ionization zone is given, showing much of the detail in this

region without the introduction of an arbitrary cutoff. This

transformation also increases the numeric accuracy of solutions in the

linear non-adiabatic pulsation stability analysis by limiting the

variation of the elements of the coupling matrices.



For the purposes of these calculations, an equation of state has

been developed, based in part on a routine supplied by Arthur Cox (Cox,

A. N. 1980b). In this appendix the various formulae and constants used

are discussed.

As our analysis requires the entropy of an ionizing gas, we first

calculate the free energy (F) as a function of the specific volume (v),

temperature (T), and composition (X,Y,Z). All thermodynamic functions

can be evaluated by derivatives of F (Zel'dovich and Razier 1966, ch.


entropy: S =- ( ,N

pressure: P = --

internal energy: E = -T2 ("F)

In order to calculate F, the approximation of an ideal, lrnizing

gas is adopted with the characteristics of the five species listed in

Table A-I. Values for the mass weighting of the elements has been



adapted from Cox and Tabor (1976), and the atomic masses and ionization

potentials are from Novotny (1973). For completeness, a term is added to

include the effects of black body radiation. This is assumed to be that

appropriate to isotropic radiation so that

F v
rad A

The 1inized fraction of ezac element is determined by minimizing

the free energy- at the gi.en temperature and specific volume. This

procedure leads to a Saha equation that is solved using a Newton Method

to a relatl.e accuracy of i part in 102 for the electron density. When

the electron density iis no,-n, the reciprocal of the mean molecular

weight fin moles'gram) 1i

n = 1,'. = :,I. 1.0'797 + Y.'4.0026 + Z/21.02447 + ne.

The relative abunrdance of each ionization stage, ni, is calculated from

the ratios of the appropriate partition functions, and the free energy

15 wuritte-i as a .su o.er these ,pecies.

F = I r n(n uI ,n + n| -2.5048 + ln(vT3/2)


where the constant term is:

1 + In(2xm k) ln(NA) -3 ln(h),
2 p A

mi is the species gram molecular weight, and ui is the temperature

independent partition function that is assumed to the ground state


The expressions for the entropy (ergs/gram/K), P (dynes/cm2), and E

(ergs/gram) are

S = R{ nl n(m32ui/n ) + n[ -1.0048 + ln(vT3/2) ] + vT3

P = nRT/v + T

E = nRT + avT4 + ni

Derivatives of the pressure and the internal energy are found by

analytic differentiation of the above formulae.

The subroutine NEOS performs the calculations required to find the

thermodynamic quantities as outlined above. A listing of this routine


e n>

w0 P
-'&-- E-

1F-> U

CP U -

.n W _
C E- P
Lj.JI E 0 I 0
> r I 0* 0
F. lI H O1

1 x W; W Q wz pq 11 O

zI JU -I I C) | Il a |
-4 ) HI m :l I g* 0= H W )W
4- 4ail -4 z r-ZCNJ
I/fllllu l/] (1 IC PTWL
3 E- -i0 I H nr-
-- S.->iiZ-. H I N ulcr-
L 3-j H EI-.-4-.--.- 3 oiN
N E -- I rrc
u >.>-C-- I . -o
2 11 I a rn w inu M r w m < 0
r C u w 1 ratsc-i-.xr :j . > <* I-n :F u
I H i 3 nO (COP5 (-
* > aS z::z...! 1- Hi HE-i II* 0 US ua3 e-i nc
>-I ~lLrJJJ- E g E(IEDa4l |IY H H. H E-4I=

* V It UHUC )iU -O4E
1 ~ Cc toC* C "* e-umC c x T rU -H
1-I ,2..ccCa.-. ij C) Q NrCariU 0. .3Z
i C '5-l Ow fiyllui g N nr M EN (0O = =C =i- viH -E a ......n.*. .... i 1I Ez EU;.

,-" -,,n H I- I QNNNc* 0 01. a*0 oa- z 0
- i*-A ^' ^^.- -J r-l j C NEUC .11 (/2 * '-C a'^CLS .2zlC -^ t I HI-li- NN N I E(oC
r* C 4 i--' 1 -,-E----** ( flo Ul C(J

* v Je- *T .-'. *. H *C .-) = I VN a 4 OW
* 0. "t czi.LuLj

3: 0 INN.(CHN* H II1w HU liii Q 11

XLLa %_- 11 i i 11 1i MI H H i iH C U NHNaC ZN 50CC

Y Cli us~ 4 iC 4 4i a P I C.H E24-W4. 00
I *
>u I *
I *

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs