THERMAL EFFECTS IN STELLAR PULSATIONS
B '
WILLIAl [DEAJ PESF.ELL
.A, ISSlF.R'aTION t R'LSLITEL' TO rhE iCrADLTATE SCHOOL
OF THFE lNIVE'ii' ,F FitLORIDL'A
id P[A.TIAiL ILFILLfLtJT OF THE PFKELilR.IHE NTS
FOF THE [EFGR.EE Of( L)CTO'n Of Pti lLOSu'piu
LObfjLLKSIT'i OF FLORIDA
ACKNOWLEDGEMENTS
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
discussions.
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 were d.,e at the University of Florida, primarily under the
IUSIC sIs ,'teu, for which I ad grateful.
TABLE OF CONTENTS
ACKNOWLEDGMENTS ... ............ ........ ................... i
ABSTRACT ... ........ ...... ....... .... ..... .............. iv
CHAPTER
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
APPENDICES
A THE ANALYTIC EQUATION OF STATE ..................93
B THE INITIAL MODEL INTEGRATOR ...................103
BIBLIOGRAPHY ...... .......................................... 125
BIOGRAPHICAL SKETCH ...................................... 129
iii
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
THERMAL EFFECTS IN STELLAR PULSATIONS
By
WILLIAM DEAN PESNELL
AUGUST 1983
Chairman: Dr. Jean Robert Buchler
Major Department: Physics
A SriurmLi,':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 et5 Crcphel variable are given showing the differences in the two
ilasse of varlt3rle'. TrIe tie scales introduced agree with earlier ones
in L eir ci m1i:n ue 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
introduced.
CHAPTER I
INTRODUCTION
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 nondimensional 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
2
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 SturmLiouville operator, embedded
in the linear, nonadiabatic 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 quasiadiabatic 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.
CHAPTER II
BACKGROUND
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 quasiadiabatic 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 nonadiabatic
and a finely tuned relationship between the location in the star .:. che
4
ionization zone and the transition from adiabatic to nonadiabatic
behavior be satisfied. Predictions of this theory are narrow instability
strips in the HR 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 nonradial 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 PG1159035, 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 nonradial 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
nonadiabatlc, but linear, stability analyses were developed (Baker and
5
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
1979).
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, nonadiabatic eigenvectors in a perturbation expansion is
limited due to their nonperiodic 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 quasiadiabatic, Buchler (1978)
developed a formalism using twotime 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
1980).
In the simplest version of the theory, the quasiadiabatic 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.
CHAPTER III
THEORY
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 (31)
Sm
d2 r Gm 4r2 g(r,S) (32)
d t2 r2 4rr2 g(rS) (32)
d t2 r2 a
TdS q(p,T) T k(r,S) (33)
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
8
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
expressions.
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 + "6fS = A.6r B*6S (35)
d k \
S (i + +r6S C.6r + D6S. (36)
it i~ well inoJn that, with the appropriate boundary conditions,
the A operator i1 of the SturmLiouville 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 [(3r14)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) (37)
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
acous
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.
CHAPTER IV
THE THERMAL OPERATOR
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) dI d
(41)
/I d L
Here q (ln) deq = 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 selfadjoint
operator. liil factor can be written as
m0 d
,im' = e (42)
witerc r, Is .a reference mass, here chosen to be more, the mass of the
10
11
core where nuclear burning occurs. Defining a new eigenvector as
d In L
= 6S/C and X = ___ eq, equation 41 becomes
v dm
T d cT d Inr d I
T Co = [X(4 KT qT) ] + [ mdl m ]
eq
(43)
Appropriate boundary conditions (which make 43 selfadjoint) are
d 6L/L b d
d m a m(m,) + b (m,) = 0 (44a)
and
F(m ) = 0. (44b)
The surface boundary condition (44a) is equivalent to the assumption
that at the surface there is no radiation incident from the outside
(Castor 1971). In this form, 44a shows that the thermal readjustment
time of the outer surface is infinitely short. As the weight function on
the lefthand side is positive definite, and the function multiplying
the derivative on the righthand side is negative definite, 44 is of
the SturmLiouville class of operators (Boyce and DiPrima 1969, ch. 11).
Comparing 43 to the standard SturmLiouville 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,
12
the eigenvectors are nondegenerate, 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
eq
Due to their SturmLiouville 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 31 through 34 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 43 becomes
ST L T n d + (45)
Fllo and Pekers d m d 4
Fciloiare Ledou.: and Pekeris (1941), we multiply eq. 45 by n and
integrate over the interval m = more m m.. A selfadjoint 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 (46)
If the thermal time is defined as the reciprocal of o0 this gives the
overall thermal time of the envelope as
t (47)
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 46 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
dP
they have the same form only if = 0 so that the formula
dm
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
14
derivative would not be correct. The inclusion of nuclear burning in 47
remedies this situation and the net thermal time for the star is then
C =< (47')
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 47 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
te 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,:j3ssed 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 T3.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
15
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 33 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,. (48)
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 47, 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 48 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 45 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. 43 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
'dm TCv2 ]
(d 1 1/TC] dx2 C2 TCv dm
\d m/ v
to transform eq. 43 to
2
d  + ( k2 w(x) ) y = 0. (49)
dx
The units of J are (time)/2 so a natural definition of the overall
thermal time of a stellar envelope is
tth= 2. (410)
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 (B4 in Appendix B below) there is a
relationship between L, and dlnT/dm. Inserting this into eq. 410 gives
m, 3 C r
th [ J ( ) 2dm ]2 (411a)
t 0 (4tr2)2 ac T3
or, using the continuity equation (31), RO = r(m0) and R* = r(m,),
R, 3 C K 1,
tth R ( v )1/2 pdr ]2. (411b)
0 a c T3
This definition shows the effect of both the internal energy that must
be dissipated and the method of dissipation. Rewriting 411(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 (Kpgas )dr ]2 (412)
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 = (lP)P.
rP rad
19
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
and
Urad = 3Prad/p = (1P)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/(1p) 2dr }2. (413)
If P is constant in the envelope
t =3 R (/2dr ]2.
th 2c 1 p Rf
20
Here we see that the thermal time is a diffusion length
1= [ R (Kp) 2dr ]2 (415)
0
divided by a diffusion velocity
v = c (1 p)/, (416)
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 enelope. These interpretations are qualitative, depending
on the i5sumed nature of U They are suitable for estimating whether
gas
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. (417)
r a c T3
In figure 41, 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 42, 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 nonadiabatic (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
22
period of the pulsation. The expression for tth is derived from a
frequency, and the location of this region would be where
1,
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. 49 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
assuptcion, but ai Important feature in the eigenfunctions can be seen
in trhi appro:imation; so, it is discussed.
RTe asyivpcotic form of eq. 49 is the harmonic oscillator and has
boundary, conditions
Yk'll) + b' d(l) = 0 (418a)
d
and
,0) = 0 (418b)
where a' and b' are a and b from eq. 44 transformed to the x coordinate
jatem. a tmuming ,I.;i = A sin( kx + 6 ) replaces the outer boundary
23
condition by a transcendental equation of the form
tan( k + 6 ) = b'/a' k (419)
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 (420)
n Yn L:)( Ld m ) ( )
In Figure 43, 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 420. In this star the integrating factor decreases
outward and the net effect is an eigenvector that peaks at the surface.
(Figures 44 and 45). 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
minimum.
Moving to a Cepheid model, the situation changes considerably. In
Figure 46, 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 47), a large peak
occurs at log(T)=3.9. This variation cannot be due to the cosine term in
420 as the argument is essentially 1.0; so, this feature must be due to
24
the denominator, whose major variations come from the integrating
factor. A large minimum in p can be seen in Figure 48, 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, nonadiabatic 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 47) 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.
m
oa
4O
0
N
t)
P,
bo >
^ I
HO
p r
,c
4J' CO
*ll 
Hr
60'6 ILi E S SCE
*H
C,
0
a)
* H
0
0
H
04
o >O
l
0
*1
a)
4 0
.e4
.C
0
aJc
n)
Fr
h   "L t4 4
4 *t' *.' '44 i,,C
~~_~ ~
(1U
w
r_ H
tc
C
g;
a
3
0
S0
r1
U
C
A
06_ __ 0 9 : a ho r
C
JJ
C
/ "''
N
+ ~. 
G I
U
r)
C
0 0
*3 '
z
a
00:1 06:L 00: L)LO 0:0 0!0 0:0 0:*
31Ui~lOH00 ~UHYFI a
28
a
0
bo
b 0
>
4
O "9 02 "L on '
o3
t0
4
3,^ 4 C friCL C CC
o
4 r
m
aC
Ul
C
SJ ,ll ] AHr 3.]. i LS
a
LL 09 0 l1 rh 0Il ? 1 9
I
ji. I 1ii
10
ci
0
c0
H
'0
a)
E9
'0
JU
;i
' 12 H
0
oJ
0
o
e
*r
.
1
.i
l o
tT s"0a Go'61 00" 1 00 6 0 n" 0"
l)
I;
I 
*n
Cr
0
',
a
0
r1
I U
ca
o
0
a)
E
A '
0 o7 0L 02 09 oh' s0
^^icI 4
I
CHAPTER V
NUMERIC RESULTS
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 33 and 34 are
differentiated at constant radius (density) and written in a form
resembling that of Castor (1971):
d 6S
T d 6S = AK2(T6S). (51)
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
2 (1) P 2 (,1) [ AK2(I1,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).
1/2
This will have no direct effect on the definition of p/2 as the two
matrix elements, AK2(I1,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
35
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 51 would then be written as
dS
= 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 SturmLiouville character of the continuum operator
in 45 in the matrix operator, the offdiagonal 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 (45), 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
36
well and raises the upper bound for positive matrix elements. A first
order expansion of the superdiagonal matrix element shows that for
KT < 1/ln(T(I1)/T(I))~17 ,
these matrix elements are always positive. We have never found any
negative offdiagonal matrix elements in any calculated envelope model.
This process also limits the variations of the matrix elements so that
the sum of the offdiagonal 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. 44, 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
V
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
order.
As the thermal time of the outer surface is zero, and the overall
1
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 nonadiabatic 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 nonadiabatic 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 quasiadiabatic.
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
eff
38
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 tre 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
39
effects in an oscillation. In particular, the quasiadiabatic stability
coefficient and pulsations that are dominated by thermal effects can be
defined.
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 51 and
the eigenfunctions 6r/r, 6p/p and 6L/L for the radial, adiabatic,
fundamental pulsation mode are shown in Figure 52. The results of a
linear, nonadiabatic radial pulsation analysis (LNA) are shown in Table
51. 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 53.
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
40
(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 54. 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 55 and 56 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 53. The modulus of the
luLmnosit .arijatcion from this analysis is shown in Figure 57. The
d ,:L L
change from d = 0 is present, albeit at a slightly higher
cemperatire chan it, Figure 54. 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.
Tis i.model sihos 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
examined.
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 52. A graph of the
physical parameters is found in Figure 58 and a graph of the adiabatic
eigenvectors 6r/r, 6p/p and 6L/L is shown in Figure 59. 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
calculated.
In Figure 510 and 511 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
512. 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
42
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 46, 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.
n
Examining the luminosity variations which correspond to Figure 512
(Figure 51). 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 514 and 515 and a > w0 in Figure 516
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
5lr, ani 517, 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
d
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 52); 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 57) dm is zero in the
atmosphere and 6L/L generally decreases with increasing temperature. The
thermal modes have the same basic luminosity behavior (Figure 54),
showing that this star is stable. In the Cepheid model, the pulsational
6L/L (Figure 526) 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
dm
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 518 and 519 where the thermal modes
with /l~1lO, and Figures 520 and 521 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
dm
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,
dm
the luminosity variations from the pulsation analysis (Figure 12bi do
44
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, 523, 524 and 525. In Figures 522 and 523, 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 524 and 525, the variations have been greatly
reduced. Lb 522, (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
preenc in Figure 514 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 37 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 (416) 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 quasiadiabatic 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.
I
/ /
/
/
S/
I.
I
ILI
OO
o
J)
0
ua
41
., ) H
IO
60
ca.
u0
n 3 r
.H
r
I :I"IL I WIl.lnl I II '. 'I" l IL II
47
0U
( t
HU
cu
oo inu
cc
/ 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 / Tl ^
1~, /u nl
1 / U t
\ / ,
I / a
I~ / 1
S~OL3 nN31" 1 tJI0
,n
Oh'o LEO L2iU
A3/SU
nC
., C *
r o
0C
*S C
m HI
1
ucc
Crd
o m
~3
m I
nl
Cd '
a)
BOl
uv I ~u ;
49
0
EO
M)
S lI
ul <
00
c'a
) l
SI
0
O EI

" I tt0 El
ri
". 01
H U
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
>g
0
41
4
00
4 U
o
0 ,
u a
00 0
02
tHl
0
S0
041
'r
aP
Ul 4
0 3H
i *
a a)
0
CO
so
(u
0 7
0
0
a *
,8
m
o
0g,
C
,i
, ,
5.
14
'I. l I II
0a
0
*u
00
:3 4J
OC
00
0r
oa
o
41
O H
0 U
t C
0 1
a,
40
0''
.HB
r1
11 901 ONU
~
SC,
1
o
0
Ul
vO
.H
,
U 1
4
0
, M
u) H
IJ0
U)
n 
0U
'4
CL (
1111 1
01
.,
ar
m
Cl
40
w 0
C) *H
z0 0
o cC
0''0
A/ CO
C )r
C.c
EI"2eh C5ELt OCCCE Ch'LEE LZ'PE8 COCIh Ee'C)I 29'Chi OCC'96 BI'Ch *' '1
A/SO La
C
L n
0
to
S0i
C
., 3"
)E
w) C
4i
0 ,0
0 1
o CGO
.I
CC
1]v 0:0
WEC
h;E LI E OE : hhf2 1012 OL:T EE I
A3!/s
09*a ES2O E['
ij..
T ,
.';,
Elu i
C
la)
o 1
10'
i 0
co
Sca
U ca
) U
El t>
II l
n 
a) o
O; 41
<"".
C;
 I
i 1. I I I i~
, '
"'*
.,',1 *
1 I
r ~
.L .
1
~J 1
? ii i
rr
I'
L" 
L'
;2I
,; '
_~nrc
CuoE S92'E ;3 C6UI h'l
 H ~ t 
OL O 90O0 o0' 22'
A l.'q
1r
U)
, u
41
C*'1 '
o
0 >r i
El
.. I
44
C)O
o C
Ca) .
.ci
oaC
K r1
.T _
3; ~
IT~7~ii
~ ..:;. ~~. ~. .I: I
I
__ I. r
_ ~~r 
;. .
;e
T
~~?;,
~;~
~; I
=I~
~i~.r;r.~~_
c
~~~~
___~_~_~
~
60
10
!D
*I
m I
O w
*H1
, "
C C
I.I.
S
ct
+ 4
CP 4 1cI
cc
4.Icc
'C
m
s ^ .; or
I~oe I'
i'ii n
i *r
(" a g 5
/ '' 5
i a n
/ ~~ t ^
[ r~in n'in o
*d,
61
en
u
0 u
7 .0
14 C*
I_
. .,
a
0  __
: "
C ) I
6h' V 0" 29'2 912 1 Ou I 99"0 F.'o 1'0 6 414'a F.E
A31SO
LOa)
CT*H tU
is si
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^
62
ai o
P c
00
SO
'V
S 0
*H
a)
o
S .
, wI O.9 aT iqsg g r i a)Qn} )s*Lr? ra "a rven I
a
^ a 73
( 
U~
r 1
*=G__^1;'.. 1
"~"'" ^s'3s^"w ^^ ^.^ ^, __ a j >
T3> ")^' 0; 4' n
C.o U
U a
C, , Ir
WU *
4J
3
't
: Ul,
La 0
SI
a) 0 l
0,II
11 1 Ira
El" U 
Y^s~f'^^" ' '
""l/'n "
; *'
L .
G'O z[' E9'0
K
T [f;n
1/C1
CO
1:
i
d) M
w I
0_
0 0
0
*H'
0 !
cI
Sv
0 4
U
4C
Cj
uT! r
r *
* irl
~tt~
.0
o~~i
 .. 
i; I H1
Lh '* o~L
o n
>
0
S o0
14 0
> In
S 0 C
4 
I
S0 '
H )
 '^ H *
Lh'O g9cD 92 a'
Al/SO
so' 9 I '=
66
0
o
Sr
co
I ,
c .
S '
JI
C"C
rl
4 "oa
0
0
a)
0H
0)0
co
0 0
4.j
00

4
,,
CC
04(0
0 n
Cv
9C 1 bL*6 W12 9 4(h's woi0 L'21 f
C )
0) .
S4,
68
0a
w 0
00
a
.Ja
a
o
a o
I u 0
00
a um
o"3
A 1/SeJ t s u 1 2 W cc 'EL 2 a
*rl
no
69
a, o
0
mi
0 0
a
C :
Sa0)
O
o pc
C,'
ca
al,
a,>
S 0 Ifl On "0 [j 90" 0E0 '00 I0
f (
m a~
\ S m
(Z~ ~ ~ ~ QS OIY [^?1 l
02*0
?0'0 20'0
I n ll I
0 0
a)
a)
0 u
CO
CO
C
50
J C
am
O "
0 r,
C
La
m w
E o
Q
'Cq
I
H .
0
* H C
C1 O
HI
0"
ajc
a' 1
CZ
LO'O
f 4 t'c _______
0/C' )C'C $4 Ci/' fC)
71
(o
0 0
e
0 o
0
0 
I 0
h L 2 ..,[ :
_00
// "' L
4J 44V
N.' c6
<^  ft
:
Table 51
BW Vulpeculae
M=11.0 M
sun
log(L/Lsun)=4.23
Tff= 25120 K
oI= 0 193 ~n=1.9 104
nIl= 0d147 '1=2.3 103
II2= 0 118 ,2=1.5 102
Composition X=0.70, Y=0.28, Z=0.02
Envelope mass = 5.02 M
sun
Table 52
DL Cassiopelae
M=5.69 M
sun
log(L/Lsun )=3.57
T = 5848 K
110= 7 87 0 = 6.8 102
nI= 5 64 n1= 1.9 01
nl2= 4d28 2= 6.2 102
Composition X=0.70, Y=0.28, Z0.02
Envelope mass = 2.76 M
sun
CHAPTER VI
DISCUSSION
To understand the utility of the thermal eigenvectors, we shall
look at various forms that the linearized pulsation equations (35 and
36) acquire when the n's are used to expand them in the ratio of
thermal time to pulsation period. Eliminating 6S one can rewrite 35 and
36 as an integrodifferential equation:
w2 6r = ( A + B i1 C ) 6r. (61)
i) D
we have an eifenvalue problem for w2 and 6r. Eliminating 6r 35 and 3
, can b. written in a form suitable for studying the secular behavior of
i.^: = (D + C B) 6S. (6la)
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 braket, 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 .
n
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
I
where
Inserting this expansion and multiplying eq. 62 by < y give
(1 + Yjiajl 2 ) 02 = < yilAl Yi > +
< Y, IBI 5n >< n ICI yi>
n i) + o
n
< Yj BI Sn >
j n(j i + o
n
(62)
"i2 I yi > = A I yi > .
< yi IBI n >< n ICI Yj >
i + 0 j
n
< Yk IBI >< n CI Yj >
kj' iw + a
(63)
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 quasiadiabatic correction can be
written as
< IBI n >< nICI Yi >
t2_a002 = n (64)
1 n id) + o
n
which is Identical cit equation (26) in Castor (1971). The normal quasi
adiabtiL: approximation is a /w. / 1 for all possible n, reducing 64
0 1
' = < Yi BC y. > (65)
I I = I.
L
II iI'n
Tri; e:rpresio:.n fr lth quasiadiabatic correction can be further
as3lzed it .1V C I. Then w2 = w + 2w.Aw and eq. 65 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 quasiadiabatic 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 }
and
S2 = { ,n n = NT+ : 0 ia < 1 },
eq. 62 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 quasiadiabatic correction. rno .iLh
an appropriate cutoff. However, the cutoff is not in the spatial
78
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. 66
becomes
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 + ) (67)
n=N +1 n n+
or
)2 1 y > = Al y > +
1 nT a
i T n >(1 )< n IC y > +
n=l
N
I BI n >(1 )< X n IC! y >. (68)
n=NT+1 n
The first term is the quasiadiabatic 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
NI
I B Cn >o n IC y >.
n=NT+1 n
This term represents a zeroorder quantity and should not be neglected
when the initial eigenvalue is calculated. If the thermal modes had a
thetafunction behavior at the transition zone, the set S2 would form a
complete set in and of itself for describing the thermal response of the
1
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
1
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 > 
n=1
BI(D')IIC y > + iwBI(D')2CI y >. (69)
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 > (610a)
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 >; (610b)
where
I y > = YO > + a I Yl > + a2 y2 > +
and
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 > (611a)
1 NT n ICI y > +
2w6 n=l
1/2< Y IB(D')2CI YO>. (611b)
At any point the matrix D' can be replaced by the sum over states that
it appro ii,.actes
N
0 I(A + X BI n > < n C) yO >, (612a)
n=NT+1 n
S xNT < Y0 IBI n >n X n 1CI y0 > +
2w20 n=l1
N
1/2 < Y IB > < n C yo>. (612b)
n=NT+1 02
Equations 611a and 611b 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
dm
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 516), 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 610a
can be calculated in a straight forward way. Defining the error in yO >
as
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
E = 1W0 LNAI'LNA
These quantities are plotted in Figure 61 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
trans
transition zone as predicted by the luminosity eigenvector plotted in
Figure 56. 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 quasiadiabatic
oscillation will follow.
When the same quantities are plotted for the Cepheid (Figure 62),
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 516). 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 612a should give much closer
agree.~ienc. This problem is discussed in Pesnell and Buchler (1983).
Comparing the initial eigensystem 612a to the adiabatic wave
equation produces three differences:
83
1) There is, at present, no guarantee of either the reality or
positivity of the eigenvalues of either 611a or 612a.
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
1
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. 612a instead of 611a
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 612a.
A final observation regarding 612a 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 nonadiabatic regions of the star must be ectensl'c, coverlri, all of
the star that is oscillating. These conditions may be found in the R
84
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 35 and
36, which represents the ratio of the pulsation period to some overall
thermal time, they become
w26r = A'6r + B.6S (613a)
iw6S = X(C6r + D6S). (613b)
For most stars, when +O the adiabatic limit is recovered as 6S is
forced to 0. However, in the integrodifferential form, splitting the
thermal modes as before,
M2 I y > = A y > +
NT BI >
w 1 + Xo /i +
n=1 n
N B >< I Y > (614)
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
85
gradient and opacity gradient so that the diffusion length (415) 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. 614 is expanded using
the binomial theorem as before,
I2 y > = A y > +
N
BI n > < n IC y > + O(X) (615)
n=NT+1 n
the B(D') C term does not disappear. For highly nonadiabatic
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
1982).
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 612a are real, and starting from
an adiabatic mode, the nonadiabatic 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 quasiadiabatic scabilic, coefficient would
86
be a poor approximation for these modes, as the second term in 612a
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
o
60
0
0
0
04
Tl
>4
a)
0 c
o
S' 1 
on o h o L 't o 2 O 2 0 
^^. a) a
,l~^
i., .i : .. *., ... . . ,.: *4 os co
i1, .i II I 1 ,111
0
r4
m
) P.
0c
w C
1 I I
) 0
CON
Cra
0 >
CN
X C.
4 44
o)
1 0
JC
4O
,za)
v4
s e
096 iPgo L.d h ii Z o 1 ' 9;
I' Jr
14
II
to
0
a ,.
0
) *
H
a)
41
BO
c 0
001
0
C Z
H
0
01
tw
0 11
\ 4J (U
,7u
*44
0 O
0.
'., 4 1
*il o
H 4
n) n
90
cl
to
0
5 A
ro
Q) 0
4j
mC
OC
0
m0o
0, a)
4 FJl 
4c
^ 41
0il
^ (1)
CONCLUSIONS
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 SturmLiouville 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
dm
layers.
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 quasiadiabatic stability coefficient and the
92
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 nonadiabatic pulsation stability analysis by limiting the
variation of the elements of the coupling matrices.
APPENDIX A
THE ANALYTIC EQUATION OF STATE
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.
III):
entropy: S = ( ,N
pressure: P = 
internal energy: E = T2 ("F)
N
In order to calculate F, the approximation of an ideal, lrnizing
gas is adopted with the characteristics of the five species listed in
Table AI. Values for the mass weighting of the elements has been
93
94
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 wurittei 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
degeneracy.
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
follows.
r)
e n>
w0 P
'& E
UHU< 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 rZCNJ
I/fllllu l/] (1 IC PTWL
3 E i0 I H nr
 S.>iiZ. H I N ulcr
L 3j 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 ratsci.xr :j . > <* In :F u
I H i 3 nO (COP5 (
* > aS z::z...! 1 Hi HEi II* 0 US ua3 ei nc
>I ~lLrJJJ E g E(IEDa4l IY H H. H E4I=
* V It UHUC )iU O4E
1 ~ Cc toC* C "* eumC c x T rU H
1I ,2..ccCa.. ij C) Q NrCariU 0. .3Z
i C '5l 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 rl j C NEUC .11 (/2
* 'C a'^CLS .2zlC ^ t I HIli 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 E24W4. 00
I *
>u I *
I *
I*
I*
.JUUUtJUU'jUUUtt)JUU'JUUUUUUUUUUUUC0 ) UUUUUU JUUUUU
