RADIATIVE TRANSFER IN CIRCUMSTELLAR DUST
By
CHRISTOPHER ALVIN HARVEL
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1974
to my parents
____ ___ j
ACKNOWLEDGMENTS
I would like to thank Professor Sabatino Sofia for
suggesting this problem to me and for the inspiration and
guidance he has given me toward a solution. I would also
like to thank Professors Frank Bradshaw Wood, KwanYu Chen,
Howard L. Cohen, and James W. Dufty for serving on my com
mittee.
I have enjoyed, greatly appreciated, and profited
by many valuable discussions with Andy Endal, Wayne McClain,
Joe Mullen, Liz Mullen, Al Rust, and Wilbur Schneider.
The curve plotting done by Andy Endal, the drafting
done by Charley Bottoms and Joe Mullen, and. the typing of
the entire manuscript done by Nadia Scheffer are also
greatly appreciated.
Computing time for the project was donated by the
Northeast Regional Data Center of the University of Florida
and the Central Florida Regional Data Center of the
University of South Florida.
___
TABLE OF CONTENTS
Section Page
ACKNOWLEDGMENTS . . . . . . . . . iii
ABSTRACT . . . . . . . . . . . v
I. INTRODUCTION . . . . . . . . 1
II. BASIC RELATIONS . . . . . . .. 30
III. NUMERICAL PROCEDURE . . . . . ... 47
IV. RESULTS OF MODEL CALCULATIONS . . . . 56
LIST OF REFERENCES . . . . . . . .. 87
BIOGRAPHICAL SKETCH . . . . . . ... 90
Abstract of Dissertation Presented to the Graduate
Council of the University of Florida in Partial
Fulfillment of the Requirements for the
Degree of Doctor of Philosophy
RADIATIVE TRANSFER IN CIRCUMSTELLAR DUST
By
Christopher Alvin Harvel
August, 1974
Chairman: Frank Bradshaw Wood
Major Department: Astronomy
A procedure is developed to calculate the radiation
field as a function of position, direction, and wavelength
within a spherically symmetric circumstellar dust shell.
The dust shell is assumed to consist of grey isotropically
scattering dust particles in thermal equilibrium with the
radiation field and to be characterized by seven parameters
(1) the radius of the central star, (2) the inner radius of
the shell, (3) the outer radius of the shell, (4) the total
optical depth of the shell, (5) an index which specifies
the density distribution, (6) the albedo of the dust parti
cles, and (7) the temperature of the central star.
The temperature distribution and the radiation field
within several model dust shells are determined, and used to
calculate for each shell the spectralenergy distribution of
the radiation emitted. These models are compared to those
of other authors. It is found that fitting the visual
spectralenergy distribution of an observed object does not
necessarily mean that the model parameters so obtained are
consistent with the infrared part of the objects'spectral
energy distribution, and that the shape of the spectral
energy distribution can be strongly dependent on the index
in the assumed power law density distribution.
The procedure is applied to the two infrared objects
HD 45677 and R Monocerotis.
SECTION I
INTRODUCTION
This paper will consider the problem of radiative
transfer in a spherical shell of dust particles illuminated
by a central star. In such a shell a dust particle will
absorb, reemit and scatter the radiation from both the cen
tral star and the other particles of the shell. Since the
star will, in general, be much hotter than the particles,
the radiation field will consist of a visual part, due pri
marily to the star and an infrared part due primarily to
the dust. Therefore, the solution of the transfer problem
must take into account absorption, reemission and scattering
of both these parts of the radiation field.
The solution of this transfer problem will yield the
quantities which characterize the radiation field at any
point in the shell and at any frequency the intensity and
the flux for example. Using these quantities we can de
scribe the flow of energy from the star to the outer edge
of the shell and predict the observational characteristics
of the shell. Also, from a theoretical point of view the
radiation field within the dust shell is of considerable
interest, since it determines a number of physical quanti
ties related to the matter in the shell. The temperature
of the dust particles, the radiation pressure on the dust
and gas and the excitation and ionization states of the gas
are all dependent on the radiation field.
From a knowledge of the dust temperature within the
shell we can limit the number of possible materials of which
the dust particles can be composed. According to Gilman
(1969), at temperatures above 15000K, refractory silicates
such as Al2SiO5 will be the most important condensates
around oxygenrich stars while carbon will dominate around
carbonrich stars below this temperature a number of other
compounds can begin to condense. A knowledge of the radia
tion pressure as a function of distance from the star will
allow us to compare its effect on the material there with
that of other forces such as the gravitational force and the
force due to high velocity particles which might be ejected
from the central star. The ionization and excitation state
of the gas associated with a dust shell can be of consider
able importance in calculating models for collapsing, rotat
ing photostars interacting with an external magnetic field
(Prentice and ter Haar, 1971).
The observational characteristics of the shell which
the:solution of the transfer problem would yield are of
considerable interest because of the recent availability of
a great deal of observational information at infrared wave
lengths. In the past few years a large number of objects
have been discovered which emit more infrared radiation
than expected. Many of these infrared objects emit most of
their radiation in the region to the red of 1i and have
gone undetected in the past because of the low sensitivity
of detectors in that region. Since 1965 the region between
1 and 20p has been surveyed extensively and several thousand
infrared sources have been discovered. A partial list of
these surveys follows: the California Institute of Technology
survey at 2.2p which discovered approximately 5600 sources
(Neugebauer and Leighton, 1969), surveys for red stars using
IN plates by Hetzler (1937), Chavira (1967), Ackermann et al.
(1968), and Ackermann (1970), and surveys at wavelengths
longer than 10P (e.g., Low 1970).
Many of the objects discovered in these surveys have
an excess of radiation within a given infrared wavelength
interval, as compared with a blackbody with an effective
temperature appropriate to the objects observed spectral
type. This excess of infrared radiation is usually referred
to as an "infrared excess." Most of the 5600 objects dis
covered in the 2.2p survey have been identified with late
type giant stars (Grasdalen and Gaustad, 1971) but about 50
of the reddest of these objects could not be identified
with cataloged visual sources. These very red sources are
manifestations of a wide variety of physical objects, in
cluding according to Neugebauer et al. (1971), HII regions,
planetary nebulae, M supergiant stars, latetype giant
stars, RV Tauri stars, novae, Be stars, and T Tauri stars
and related objects. Most of the infrared excesses asso
ciated with the T Tauri stars and related objects are as
sumed to result from thermal radiation from a circumstellar
dust shell. This assumption is also made for Be stars but
in their case it has been argued by some authors, e.g.
Woolf et al. (1970), that the observed infrared emission is
consistent with freefree emission from ionized hydrogen.
For those objects other than the Be stars and the T Tauri
stars and related objects thermal radiation by dust is prob
ably present but may not be the principal cause of the
observed infrared excess.
Since a general solution of the transfer problem for
circumstellar dust shells would be readily applicable to the
T Tauri stars and related objects a more detailed descrip
tion of these objects will be given at this point.
The T Tauri stars and related objects form an especial
ly interesting group from both an observational and a
theoretical point of view. Poveda (1965a) pointed out in
1965 that T Tauri stars, which for the last 20 years have
been assumed to be young gravitationally contracting objects,
might have thick circumstellar dust shells. He argued that
such a shell would be a remnant of the contracting cloud
from which the star formed and that it would be bright in
the infrared. Shortly after this Mendoza (1966) found that
most T Tauri stars do have infrared excesses and later
Mendoza (1968) found that the bolometric luminosity of these
objects was from 1.3 to 6.6 times that expected from the
visual observations. For R Monocerotis, an object closely
related to the T Tauri stars, Mendoza (1968) found the
bolometric luminosity to be 58 times that expected from the
visual data. Figure 9 shows the spectralenergy distribu
tion of R Monocerotis and it can be seen that most of the
observed luminosity is emitted beyond li.
The T Tauri stars are Ha emission stars, and in the
early 1940's all such emission stars were studied spectro
scopically because of the intrinsic value of their spectra.
It was not until 1945 that Joy (1945) recognized T Tauri
stars as a distinct class of emissionline variables asso
ciated with nebulosity.
From a study of eleven stars, Joy established the cri
terion to be used in assigning stars to this new class. He
describes in detail each of the eleven stars and says that
the spectra in general are F5G5, but that there is a small
variation in spectral type with the brightness. Also,
absorption lines are usually lacking but many bright lines
of low excitation are present. He points out the similarity
between the spectra of these stars and the solar chromo
sphere. For five of the stars he finds color excesses, which
he attributes to strong selective absorption of the nebulosi
ty. He finds that usually the emission lines seen are dis
placed to the violet of the absorption lines (if any absorp
tion lines are present), and that it is very difficult to
say anymore about the radial velocities of these stars.
In a review article on T Tauri stars and related
objects, Herbig (1960) points out that the only conclusive
test for membership in the T Tauri class is provided by the
spectroscopic data.
The classification system in use today is described in
and employed by the Third Edition of The General Catalogue
of Variable Stars (Kukarkin et al., 1969). The criteria
given in the catalogue are those adopted by the IAU Commis
sion 27. These criteria are spectroscopic with one exception,
which is, that if the variable is connected with diffuse
nebulae it is designated InT whereas if it is not it is
designated IT. The prototype is stated to be T Tauri.
The catalogue breaks the irregular nebular variables
into several classes, one of which, InT, represents the T
Tauri stars. Its class Is is the same as Hoffmeister's RW
Aur class, and differs little from the InT class. Three
more of the classes defined by the catalogue, In, Ins, and
Inb, contain young objects of spectral class FM with ir
regular light variations (these have been called the Orion
variables). Of its nebular variable class the catalogue
class Ina (early spectral type Orion variables) differs the
most from class InT.
In some of the current literature stars falling into
all these classes are referred to as T Tauri stars and
related objects, although usually only in the sense that
the star is in an early evolutionary stage can it be grouped
with the T Tauri variables.
R Mon, the unusual object observed by Mendoza, is not
included in class InT because it is assigned by a Russian
author to spectral class AFep, which is too early a
spectral class for an InT variable. Joy (1945) originally
classified it as approximately G, and Mendoza (1968) classi
fies it as K. Thus this star has been assigned spectral
classifications from A through K. This points out an impor
tant fact; namely, the spectral classification of T Tauri
stars is very difficult because of the lack of a strong
absorption spectrum.
Most of these stars are rather faint; the brightest of
the class is about my = + 9 at maximum light. The magnitude
usually quoted for T Tauri variables is the mean of the
observed magnitude, even though this may not be the most
meaningful value.
For the last twenty years T Tauri stars have been
assumed to be premainsequence (PMS) stars and recently PMS
models have become available which might help to verify this
assumption. The best protostar models at present are those
of Hayashi (1966) and Larson (1969b). Both of these models
are directly applicable to the problem of the T Tauri
phenomena.
For his model Larson assumes that the initial proto
stellar cloud is spherical, not rotating, and free of
magnetic fields. He does not assume, as Hayashi does, that
the collapse is homologous and that the density follows a
polytropic law. Initially the cloud has the Jeans radius
(R = vs[Gp]1/2, where vs = the velocity of sound, G = the
gravitational constant, p = the density), and a density of
1.10 X 1019 gm/cm3. The time scale for collapse in this
model goes as 1//p, so that the density gradient increases
with time, resulting in a dense core, which finally gains the
size and mass of a star. The core then evolves toward the
Main Sequence (MS) while accreting matter from the lower den
sity cloud about it.
If the mass of the cloud is 2.5M.or less, the proto
star will become visible somewhere along the lower part of
the fully convective Hayashi track; if it is somewhat more
massive it becomes visible at some point on the radiative
track running horizontally toward the MS; if it is very mas
sive (as massive as an 0 or B MS star) it does not become
visible until reaching the MS.
The time required in Larson's model for the accretion
of the envelope depends on the mass of the protostar, and
ranges from 105 to 106 years.
Hayashi's computations start with a protostar of much
higher density, and in his model the accretion of the cloud,
which brings the protostar to the beginning of the fully
convective phase, takes place in a few decades.
The two models differ in their treatment of the shock
phenomena at the interface between the shell and the dense
core, but the latter parts of their evolutionary tracks are
in agreement. In both models we have a dense core and a
circumstellar shell, and Strom (1972) points out that we
should expect some young stellar objects to show evidence of
the remnants of their shells. These protostar models
neither explain nor predict the variability, line emission,
and mass ejection observed for the T Tauri stars.
There are a number of reasons for believing that the T
Tauri stars are PMS objects: (1) their large red and infrared
excesses and deficiencies; (2) their position on the HR
diagram (above and below the MS); (3) their association with
dense interstellar matter; (4) their association with O and
B stars; (5) their apparent rotational velocities; (6) their
anomalous Li abundance; (7) their polarization.
As noted earlier, Joy (1945) found that almost half of
his T Tauri stars had color excesses. Since then other
authors have found large infrared excesses for some stars,
and infrared deficiencies for others.
Mendoza (1966) did photometry in the UBVRIJKLM bands
from .361 to 51 for 26 stars, including T Tauri, RW Aur, R
Mon and Lk Ha 120. His observations show red excesses for
most of the stars, a deficiency (EVR 0.2) for DF Tau most
likely caused by the Ha line in emission in the V band, and
infrared excesses (from 0.9 to 5V) for all the stars
(EVM = 8.5 for R Mon).
Two years later Mendoza (1968) published another set
of photometric observations on Johnson's nine color system;
this time for 33 stars, some of which he had observed
previously (Mendoza, 1966). His results were similar to
those of 1966. He states that the two stars with the largest
infrared excesses are R Mon and R CrA. Most of the objects
observed by Mendoza (1966, 1968) were T Tauri stars.
In making the observations described above Mendoza
used the 60" infrared telescope at the Catalina observatory
in Arizona. He notes that making observations of these
objects is made difficult by their faintness and their
close association with nebulosity. He tried various dia
phram sizes from 36 to 13 seconds of arc in diameter, and
found that in many cases a diaphram smaller than 13 seconds
was actually needed.
From his observed fluxes, Mendoza (1968) computes
bolometric corrections, and then absolute magnitudes for
several of these objects. He finds that they have high lu
minosities, and speculates that Joy (1945) classified them
as dwarfs because of their rather wide absorption lines.
To explain his 1966 observations Mendoza (1966)
advanced two theories: (1) a multiple star system in which
one of the components is a very cool infrared star, and (2)
a coreenvelope model in which the UBV photometry pertains
to the core and the infrared photometry to the envelope. In
this paper (Mendoza, 1966) he does not favor one theory over
the other.
In the second paper Mendoza (1968) states that his
previous coreenvelope model is probably the explanation for
the large infrared excesses observed for objects such as R
Mon and R CrA. For the objects with small infrared excesses
(infrared just a few times the visual) he suggests that the
cause is heating of the stellar wind by cosmic rays.
He points out that his observations indicate that the
envelope around the core behaves like a neutral filter in
the visual (UBV) region. Assuming that all the infrared
radiation comes from this envelope, he finds that it would
have a radius of about 10 AU and a mass of about 103 M
Low and Smith (1966) published the results of their
photometry of R Mon at 20p. They found that at that wave
length the star was about one fourth as bright as Mondoza
(1966) found it to be at 3.4p.
To explain their observations they computed the
spectralenergydistribution to be expected from an optical
ly thin circumstellar dust shell in which scattering and
reabsorption of infrared radiation emitted by the dust can
be neglected. The theoretical distribution they obtained
matched the available observations fairly well.
Strom et al. (1972) observed 42 stars in NGC 2264 at
1.6, 2.2, and 3.4p,and found infrared excesses and some
infrared deficiencies. The authors state that these obser
vations are an indication that the stars have circumstellar
shells. They used profiles of the observed hydrogen lines
to estimate the surface gravities of some of the PMS A stars
in their sample of 42 NGC 2264 stars. From these surface
gravities they computed luminosities, and found that those
stars that seemed the faintest, as compared with their
computed luminosities, all had infrared excesses. They sug
gest that these infrared excesses are due to a dust shell,
radiating in the infrared, energy it has absorbed from the
star in the visual part of the spectrum. The stars that
they found to have the largest infrared excesses, were those
that lie below the MS.
Strom et al. (1972) also found that these shells have
a very high ratio of total to selective absorption, which
would suggest large dust particles. The dust shells of
lowest temperature were around the latest stars.
They speculate that the infrared deficiencies observed
could be due to a very cool dust shell at a great distance
from the central star. Of course, observations in the far
infrared should find the missing radiation.
Martin Cohen (1973b)has observed stars in NGC 2264, IC
5146, VI Cyg and NGC 7000 (The North American Nebulae). NGC
7000 contains the T Tauri star Lk Ha 190, which was discover
ed by Herbig in 1958, and has recently brightened visually
by 6 magnitudes.
Cohen used the 60" Catalina infrared telescope, and
made observations in all wavelengths in a single night (this
is very important since these objects vary rapidly). He
concluded that for most of these stars (many of them not T
Tauri stars) thermal emission from circumstellar dust is the
cause of the infrared excesses he observed. In a few cases,
he says that freefree emission from a circumstellar shell
may be important out to about 11 .
Breger (1972), like Cohen, finds evidence for circum
stellar shells about PMS stars in NGC 2264. As well as
infrared excesses, he finds light variability, which he at
tributes to circumstellar shells.
Strom et al. (1971) found a small amount of continuous
and Ha emission, similar to that found for T Tauri stars, in
the spectra of some AF stars in young clusters. They attri
bute this to very thin circumstellar shells. They also say
that the strength of these features increases as the
circumstellar absorption does.
According to Strom (1972) we should expect to see free
free (electronH ion) and freebound emission from envelopes
in the infrared, if they have sufficient size and density.
Dust envelopes, as stated previously, would show thermal
infrared reradiation; thus, Strom (1972) ascribes the cir
cumstellar emission in the optical region of T Tauri spectra
to gas (this is the emission that washes out the normal
absorption spectrum). He says that, in a joint paper, Kuhi
and the Stroms indicate that gas shell emission may also do
minate the infrared region.
In support of his idea that gas and not dust shell
emission is dominant in the infrared, Strom (1972) offers
the following argument and observational evidence. He
finds that for some young stars (especially AF stars) with
circumstellar shells, Ha emission and infrared excesses
there exists a correlation between log (Ha intensity) and
the L magnitude (3.4p). In an optically thin case of the
dust model the dust emission is proportional to the optical
depth (T = p x length), whereas the gas emission (Ha)
depends on p2 x length. Therefore, unless we believe that
all circumstellar shells have the same densities, the dust
theory is difficult to defend, because we would not get the
observed correlation.
According to Herbig, in the Taurus cloud, Orion
nebulae, NGC 2264, NGC 6530, and IC 5146, we find that if we
fit the brightest T Tauri stars to a MS with a distance modu
lus of 6 the Ml to M3 stars observed are too bright. Thus,
the observed MS in these regions appears to be turned up on
the red end; however, the observed MS does not curve upward
more and more as we observe later and later stars, but runs
parallel to the Zero Age Main Sequence (ZAMS) after sharply
breaking off from it at about AO. This observed MS forms a
band with its upper edge several magnitudes above the ZAMS.
Herbig (1960) believes that the uncertainties in the
observations and theory are large enough so that we need not
be concerned by the fact that the observed MS does not curve
upward. He also believes that part of the spread in the
observed lower MS is caused by a spread in the times of
stellar formation in the cluster.
Hayashi's evolutionary tracks for protostars would
lead us to expect young stars to lie above the ZAMS.
Poveda (1965a) first brought attention to objects
falling below the ZAMS in NGC 2264. He said that their lo
cation below the ZAMS was caused by thick circumstellar
dust shells.
Strom et al. (1972), as we have seen, were able to
show that those stars, in NGC 2264, lying below the MS would
be above it, if their circumstellar shells were removed.
Presumably, the stars below the MS are younger than
those above it, and have not had time to flatten their
circumstellar shells. According to Mendoza (1968), the T
Tauri stars with the least infrared excess are the oldest,
and have the most flattened disklike shells. If we are
looking at a disk edge on, we should see a large discrepan
cy between the luminosity as determined from the spectral
classification, and the bolometric luminosity observed.
Strom (1972) finds, for a sample of 10 Herbig Ae and Be
stars, that for three, Lgp > > Lbolometric
According to Strom et al. (1972), the presence of
latetype stars near the MS in young clusters is caused by
circumstellar obscuration, and not by a large spread
(several x 107 years) in times of stellar formation, as
argued by Herbig (1960).
T Tauri stars are found directly associated with the
clouds in which they presumably formed. They can not be
interlopers from the general field around the region in
which they are found, since their space density averages 1
or 2 orders of magnitude above that of the field stars of
the same luminosity (Herbig, 1960).
Cohen (1973c) has observed the T Tauri type objects
at the tips of several cometary nebulae. Mendoza (1966)
states that T Tauri and R Mon illuminate bright nebulosities.
Both of these nebulosities are variable cometary nebulae.
R Mon is in NGC 2261, Hubbles variable nebulae, and T Tauri
is in NGC 1555, Hind's nebulae.
Poveda (1965c) describes the cometary nebulae
(elephant trunks, bright rims, etc.) as HII regions expand
II
ing into a surrounding HI region, or as a HII region expand
ing past an area of high density, compressing and deforming
it. If either of these explanations is correct, the tip of
such a nebulae will be a region of high density, suitable
for star formation. Thus, cometary nebulae are very young
phenomenon illuminated by even younger stars.
If we assume the stars have the space velocities of
the eddy of gases from which they formed (about 2 km/sec),
and if we estimate the size in km of the tip of the nebulae,
we can estimate how long it would take the star to drift
out of the tip, and thereby arrive at an upper estimate for
the age of the star. Poveda does this for several nebulae,
and finds ages running from 5 x 103 to 50 x 103 years.
Mendoza (1968) also makes a calculation to obtain the
ages of some T Tauri stars. He uses the PMS evolutionary
tracks of Iben to determine the age of a star from its
position on the HR diagram. The ages Mendoza obtains in
this way agree with those of Poveda for some stars and
disagree for others.
For R Mon the two methods agree rather well, but
for the other stars Poveda's results are, in general, lower
than Mendoza's.
Herbig (1960)states that in regions such as Orion, we
find T Tauri stars associated with 0 and B stars, which are
considered to be young objects. It seems logical then, to
consider the T Tauri stars to be of similar age (hence
young), but of lower mass, and thus not as greatly evolved.
Joy (1945) classified the T Tauri stars as MS dwarfs,
which we now know to be incorrect. Those T Tauri stars
that do show absorption lines, have very wide lines, which
might indicate that they were luminosity class V objects.
However, it is now believed that these abnormally wide
absorption lines are caused by rotational broadening (Herbig,
1960).
The rotational velocities found for these stars
(20 km.sec < v sin i < 65 km/sec) are usually high for
stars of their spectral class. Normal GK stars have
v sin i < 15 km/sec.
If we compute the radii of the T Tauri stars from
their observed luminosity and colortemperature, and then as
sume they contract without loss of angular momentum along a
Hayashi track, we can find their rotational velocities on the
MS. When this is done, the computed values are in close
agreement with observed values of v sin i for MS stars of
comparable mass. It may of course be true that the line
broadening is not caused by rotation.
The T Tauri stars as a group appear to be overabundant
in Li, as compared with the Sun, by 1 or 2 orders of magni
tude, which is taken as an indication of their youth (Herbig,
1960). Herbig compared the Li abundance of the T Tauri
stars to that of chondritic meteorites, which presumably
have the Li abundance of the presolar nebulae.
It is possible that the large abundance of Li is being
caused by Li production in the photospheric layers of the T
Tauri stars, in which case it may not be an indication of
their youth. It seems unlikely though that such Li produc
tion is actually occurring, since the most plausible
mechanisms suggested for Li production require large, strong
magnetic fields, which have not been observed.
Zellner (1970), Capps and Dyck (1972) and Breger and
Dyck (1972) observe polarization for some stars in the near
and far infrared. Zellner's measurements for R Mon indicate
polarization typical of large grains (radius u 0.1l) while
the other authors cited indicate that their measurements
favor electron scattering similar to that observed for Be
stars. In a more recent paper Zellner and Serkowski (1972)
state that R CrA exhibits polarization similar to R Mon.
They suggest that in objects such as R Mon and R CrA the
dust cloud producing the infrared excess is in the form of
an equatorial disk and the light seen at shorter wavelengths
is scattered light which has filtered out the poles of this
disk.
Breger and Dyck find that between 3000 and 7500 A,
four out of a sample of 35 stars in NGC 2264 show intrinsic
polarization. Most of the stars they observed were not T
Tauri stars, but of the four showing polarization one is a
T Tauri star. Breger and Dyck comment that its polarization
has a different X dependence from that they observe for the
other three polarized stars in their sample, and that its
polarization is also different from that of T Tauri, as
given by Serkowski (1971).
These investigations indicated the presence of circum
stellar shells about T Tauri stars and stars in clusters as
sociated with T Tauri stars, such as NGC 2264.
Some other important characteristics exhibited by
young stars and T Tauri stars are the following: (1) mass
loss and the P Cyg effect, (2) mass accretion and the inverse
P Cyg effect, (3) ultraviolet excess, and (4) correlations
of the visual magnitude and the infrared emission, of the Ha
UV and IR emissions, and between mass infall and the V
magnitude.
Most T Tauri stars show a P Cyg type spectrum, that is,
their emission lines havebleshifted absorption features
superimposed. According to Strom (1972) and Herbig (1960),
the theory that this effect is caused by an expanding envel
op with self absorption is widely believed.
Herbig (1960) states that the star T Tauri has a set
of two absorption lines superimposed on its H and Ca emis
sion lines, corresponding to radial velocities of 70 and
170 km/sec, and that RY Tau has one line corresponding to a
radial velocity of 90 km/sec.
From an analysis of the forbiddenline radiation, com
ing from the small forbiddenline region surrounding T Tauri,
Herbig says that Varsavsky, in 1960, finds a mass loss of
about .5 x 105 Me/yr for T Tauri, assuming that the
forbiddenline region is 750 AU in diameter, and that the
region is sustained by material rising at 170 km/sec.
The accretion of matter by some young stars is indicated
by the presence of an inverse P Cyg spectrum (sometimes
called a YY Orionis spectrum). Out of a sample of 23 UV
excess stars observed by Walker (1969) in NGC 2264, he finds
10 show an inverse P Cyg effect. None of these are T Tauri
stars, but they are all closely related to T Tauri stars. In
1972, Walker (1972) reports finding two more young stars
with inverse P Cyg effect, RR Tauri and 01C Orionis.
T Tauri spectra show a strong continuum in the visual,
which is sometimes called the 'blue continuum', and a
continuum starting in the XX37003800 region and rapidly
rising toward the shorter wavelengths. Unlike the blue
continuum, this UV continuum adds an appreciable amount of
energy to the luminosity of the star. According to Herbig
(1960), this UV continuum has been attributed to synchrotron
emission like that from some solar grains, and to the running
together of the higher Balmer lines. The latter of these
hypotheses he refers to as B6hm's hypothesis. The higher
Balmer lines blend together because of turbulent broadening,
according to this hypothesis.
Mendoza (1966) finds UV excesses for the majority of
the 26 stars he observes. Anderson and Kuhi (1969) give
detailed observations of the T Tauri star AS 209 (1900: 16h
43m.6, 140 13'), and report that it shows a large UV excess.
They say they favor B6hm's hypothesis, since Kuhi's
observed line width for the HB line, indicates a turbulent
velocity as large as 100 km/sec, which would be more than
enough to blend the higher Balmer lines.
In general, Ha emission and UV excess in T Tauri stars
is correlated, but for individual stars there does not appear
to be a correlation between the two as the star varies in
brightness. In the case of AS 209, Anderson and Kuhi looked
for such a correlation, and found none.
Various correlations have been noted between observed
T Tauri phenomena. Herbig finds that as RW Aur gets fainter
in the V band, it gets redder. Anderson and Kuhi (1969) find
that AS 209 increases its UV excess as it gets fainter in the
V band. Cohen (1973a) finds that as T Tauri gets fainter in
the visual it gets brighter in the near infrared.
Walker (1969) does a detailed study of SU Ori, an Inb
star with inverse P Cyg spectrum and UV excess, and finds
that when the redward displaced absorption feature is most
prominent, the star is brightest. He shows five spectro
grams of SU Ori taken when the star was at different magni
tude, and the redward displaced absorption line disappears
as the star gets dimmer.
Poveda (1965b, 1965c) suggests that the circumstellar
shells about T Tauri stars and related objects are pre
planetary systems. He says that some of the variability of
these stars must be due to eclipses by protoplanets (large
clouds along the orbits of the planetstobe), and that
phenomena such as FU Orionis occur because the circumstellar
material can form into large particles very rapidly (a few
years); thus, rapidly dropping the optical depth of the
material, and brightening the star.
In light of the observed characteristics of the T
Tauri stars and related objects described above it is not
surprising that there have been a number of papers published
in the past few years dealing with the problem of radiative
transfer in circumstellar envelopes. In each of these
papers the author makes a number of assumptions about the
object or objects he is attempting to model in order to sim
plify the radiative transfer problem. These simplifying as
sumptions tend to reduce the accuracy and the versatility of
the resulting model and should therefore be avoided if
possible.
In all of these papers it has been assumed that the
envelope is spherically symmetric about the central star and
that the effects of the spherical symmetry are significant.
If the entire shell were metrically thin and distant from
the central star it would be possible to simplify the prob
lem by assuming the shell to be stratified into parallel
planes. This simplification could also be made if the outer
layers of the envelope were optically thick. However, it
appears from the observations that the circumstellar clouds
about the T Tauri stars and related objects can not be
adequately represented by a model based on a plane parallel
approximation. The assumption of spherical symmetry while
much better than the simple plane parallel assumption is
still an idealization since it has been pointed out that
visually these objects appear to have a very complicated
outer structure (e.g., Herbig, 1968).
In all these papers it is assumed that a power law of
the form K(r) = [C/rn] governs either the volume absorption
coefficient or the number density of particles. Here r is
the distance from the central star, C is a constant, and n
is a variable index (Chandrasekhar, 1934). Some assumption
as to the functional relationship between K and r is neces
sary (see Section II), and this one seems to allow reason
able flexibility with the adjustment of only one parameter
(the constant is determined by the geometry and total
optical depth of the shell). In a few of these papers it is
further assumed that n = 0; that is, that K(r) is constant
with r
Several authors (Stein, 1966; Low and Smith, 1966;
Krishna Swamy, 1970) have considered the less general prob
lem of radiative transfer in an optically thin circumstellar
shell in which scattering and reabsorption of infrared radia
tion emitted by the dust can be neglected. Stein's paper
pertains to the observed infrared excesses of certain B
stars; Low and Smith's, to the much larger infrared excess
of R Monocerotis and Krishna Swamy's, to the infrared excesses
of several infrared objects, among them, VY Canis Majoris.
Stein in his paper states that the total optical depth
of the shells he is considering is about 103. For a shell
this optically thin it is quite reasonable to assume that
the star is the only significant energy source in the system.
In the paper by Swamy and in the the one by Low and Smith,
the objects modeled (e.g., VY CMa and R Mon) show an amount
of infrared radiation as compared with the optical which
indicates that these objects have optically thick shells
(Mendoza, 1966; Low et al., 1970).
Huang in one paper (Huang, 1969a) considers the two
extreme cases of very small and very large optical thickness
and then in a second paper (Huang, 1969b) considers the more
general problem of intermediate optical thickness. In both
these papers Huang makes the assumption that the opacity of
the dust is proportional to the geometrical cross section
of the grains.
In the first paper he solves the transfer problem for
an optically thin shell in much the same manner as Stein
(1966). His approach to the problem of the optically thick
shell is analogous to that of Chandrasekhar (1934) for
extended stellar atmospheres.
In the second paper he solves the transfer equation
for a spherical shell having its inner boundary in contact
with the surface of the central star, assuming grey
particles, the Eddington approximation, isotropic scattering
and Local Thermodynamic Equilibrium. He does not actually
state that the shell is in contact with the surface of the
star but that the radiation at the inner boundary of the
shell is completely diffuse. For all practical purposes as
long as the inner edge of the shell is within a few stellar
radii of the surface of the star the radiation field can be
considered diffuse.
Huang (1971) published a third paper dealing specifi
cally with shells distant enough from their central stars
such that the stars could be considered point sources. He
points out in this paper that the assumption of a complete
ly diffuse radiation field throughout the shell represents
an idealized case of one extreme while the assumption that
the star is a point source represents the other extreme.
In his second paper Huang (1969b) argues that the as
sumption that the particles scatter radiation as greybodies
is valid. He states that in the case of interstellar red
dening blue light is preferentially scattered out of the
line of sight, whereas in the case of a circumstellar dust
cloud we observe scattered radiation coming in all direc
tions; such that, the amount of radiation of a certain wave
length which is lost in direct transmission along the line
of sight is regained due to the overall scattering process.
For his intermediate case Huang (1969b) assumes that
the radiation field can be broken into two parts a visual
part and an infrared part. He then proceeds to solve the
overall transfer problem as though it were two less diffi
cult problems: (1) a nonconservative scattering problem in
volving only the stellar radiation in the visual (the
infrared emission of the star is assumed to be negligible)
and (2) the problem of an atmosphere in local thermo
dynamical equilibrium for the infrared radiation emitted by
the dust, with an energy source, namely, the visual radia
tion absorbed in the nonconservative optical region. There
fore, the results for his intermediate case described in his
second and third papers (Huang, 1969b, 1971) refer to these
two broad wavelength regions.
In his results Huang (1971) states that the radiation
field and the temperature distribution do not critically
depend upon the index n in the expression assumed for the
density distribution and indirectly that the ratio of the
inner to the outer radii of the shell is an unimportant fac
tor. In Section IV it is shown that these results are not
in general correct.
Huang does not apply this method to any observed
object, but in 1970 Herbig (1970) uses it successfully in
the development of a model for VY Canis Majoris, and recent
ly Apruzese (1974) has used it to model the infrared object
HD 45677 (a Be star with a large infrared excess). The
approach followed by Apruzese parallels Huang's for the in
termediate case and diffuse radiation except that he
modifies Huang's technique slightly to allow the determina
tion of the radiation field in ten spectral regions to the
blue side of 1w For each of these spectral regions he
has in effect a problem of nonconservative scattering.
Apruzese's model will be considered further in Section IV.
Larson (1969b) in studying the emission of collapsing
protostars finds an approximate solution for the radiative
transfer in a spherical dust shell in which the particles
scatter none of the radiation but absorb an amount propor
tional to AP, where A is the wavelength and p is a varia
ble index. He derives an approximate expression for the
temperature distribution in the cloud and using this solves
for the emitted spectrum. The weighting function integra
tion method he develops to compute the emitted spectrum
from the known temperature distribution is also used by
Herbig in the paper mentioned above. For the index n which
is used to establish the density distribution in the cloud
Larson uses 1.5 and states that this is the value determined
by his dynamical calculations and published in a previous
paper (Larson, 1969a).
Hummer and Rybicki (1971) have a procedure involving
iteration on the Eddington factor, K/J, which they use to
solve the transfer problem, under the assumption of the con
servative grey case, for a spherical cloud extending from a
point of radiation at its center out to some finite radius,
R.
None of the papers mentioned above adequately deal
with the problem of radiative transfer in circumstellar en
velopes of intermediate optical depth. The only two proce
dures that come close to solving this case are those due to
Huang (1969b, 1971), Larson (1969b), and Apruzese (1974).
They all assume that the density in the shell is proportion
al to 1/rn and Huang and Apruzese assume that the particles
are isotropic scatters Larson neglects scattering. The
first of these assumptions seems as good as any other for
the density distribution and is not critical to the solution
of the problem. The second assumption is necessary in
order to take advantage of the spherical symmetry of the
shell. It allows us to reduce the problem to one involving
only two spacial coordinates. If the particles are assumed
to scatter according to an arbitrary phase function the
problem becomes much more complex. Huang assumes that the
opacity of the cloud is due to the geometrical cross section
of the particles and therefore that it is grey, whereas,
Apruzese and Larson assume that the opacity is wavelength
dependent. The procedure described by Huang is restricted
to cases where the central star can be considered a point
source or the radiation is completely diffuse. It is
further restricted to cases where the Eddington approxima
tion is valid. Most importantly his procedure does not de
scribe the radiation field at a particular wavelength but
merely breaks it into two regions visual and infrared. The
procedure described by Apruzese is based on Huang's proce
dure and involves one major improvement his procedure can
be used to determine the radiation field at a number of
wavelengths in the visual. Larson's procedure, as well as
neglecting scattering, uses a rather crude approximation
for the temperature distribution in the shell. His proce
dure does however allow the determination of the spectral
energy distribution emergent from the shell in the visual
and infrared.
In Sections II and III a technique is developed
capable of solving the transfer problem for the radiation
field at any point in a circumstellar shell, for any wave
length, and any optical depth.
This technique is applicable to cases where the
central star can be considered neither a point source nor a
diffuse source of radiation. To simplify the problem the
following assumptions are made: (1) the cloud is spherically
symmetric about the central star, (2) the cloud is made of
small particles which have no preferential orientation, (3)
the particles reemit absorbed radiation isotropically and
scatter radiation isotropically, (4) the particles have an
albedo which is independent of wavelength, (5) the particles
radiate as greybodies at a temperature T where T is a func
tion of the radiation field about the particle, (6) the
volume absorption coefficient as a function of the distance
from the central star, K(r), is constant with wavelength and
proportional to the density p(r), and (7) the optical depth
T and the distance from the central star are related by the
usual expression (Chandrasekhar, 1934),
dt r
dT = b n t =
where b is a constant.
This technique characterizes each dust shell by a set
of parameters related to the size of the shell and the
density distribution within the shell, the size and tempera
ture of the central star, and the albedo of the dust
particles within the shell. The equation of transfer is
solved in integral form by a computer code using an itera
tive procedure.
SECTION II
BASIC RELATIONS
In this section the physical and mathematical relations
necessary to solve the transfer equation for the radiation
field in a circumstellar dust cloud will be developed. These
relations will be in terms of a number of physical parameters
that characterize the dust cloud.
The parameters to be used can be placed into three
groups, one that specifies the geometry of the cloud, one that
specifies the albedo of the particles, and one that specifies
the temperature of the central star. The first group consists
of five parameters, the radius of the central star, the inner
radius of the cloud, the outer radius of the cloud, the total
optical depth of the cloud, and an index which specifies the
density distribution. These five quantities will be desig
nated respectively as R,, RI, R, TI, and n. The second and
third groups have only one member each: the albedo of the
particles, designated o in the case of the second group and
the temperature of the central star, designated T*, in the
case of the third.
For some set of these seven parameters the dust cloud
model should be capable of reproducing the observed spectral
energy distribution of an object such that that set of param
eters can be considered to characterize the observed object.
Consistent with the simplifying assumptions set forth
in Section I, relations can be derived for a number of useful
quantities. Using the expression for dT given in Section I
we can derive the relation for the absorption per unit volume
at a distance r from the central star, K(r). We have
dr = Rdt
d bRnldr
rn
but
dT = K(r)dr,
and therefore,
bRn1
K(r) =
rn
or
B
K(r) = (1)
n
r
where B is a constant. We can also derive expressions for
the optical depth at r, T(r), and for B in terms of R R,
TI, and n. In the following n $ 1 unless specified otherwise.
We have
r
T(r) = f K(r)dr
o
or
r
T (r) = Bf rndr
0
and, therefore,
ln
T (r) = B( r ) + C (2)
1n
where C is a constant of integration. Since T(RI) = TI and
T(R) = 0 we have
1n
and
Rln
0 = B (  ) + C
1n
from which we obtain,
C = B ( Rn).
1n
Substituting this into (3) we have for Ti,
R 1n
TI = B( )
l1n
R1n
 B( )
1n
B (R 1n R1n
I 1n
So
I (ln)
B =
(R 1n RI )
and we obtain from equation (1)
K(r) = (1n)
(R1n R 1n)rn
I
if R < r < R
= 0, if r > R, or r < R
and from equation (2)
(r) = T(R1 rln)
Tir) = R
(R1n R 1n
I
In the case with constant density we have n = 1 and
the relations for K(r) and T(r) become
1
K(r) = TI[Zn(R/RI)r] RI 5 r < R
= 0, r > R, r > R (7a)
T(r) = in(R/r) [ n(R/Ri)] (7b)
Equations (6) and (7b) apply only for RI r < R. For
r < RI, T(r) = TI, and for r > R, T(r) = 0.
The spectralenergy distribution is observed in terms
of flux received at the surface of the earth per unit wave
length, so this is the quantity we want to calculate with
our model. If we define F(A,r) as the flux at wavelenth A
through a unit area of the shell of radius r about the
central star, where F(A,r) is in units of ergs/cm2 
str sec, we have
F(X,D) = F(X,R)R2/D
where D is the distance between the object and the earth.
Note that if D is not a known quantity we can still find a
scaled flux and thus determine at least the shape of the
spectralenergy distribution.
We also have
F(A,r) = f I(A,r,u)cos0 dw (8)
where I(A, r, w) is the specific intensity of radiation at
wavelength X, at any point a distance r from the central
star, and coming from some direction w radians from the cen
tral star. The integration is over 47 steradians. Because
of our assumption of spherical symmetry we can define I
uniquely with reference to only two of the three spherical
coordinates, r and 6; there will be no 4 dependence. It is
also true that I(A, r, 6) = I(X, r, 6), so the integration
over 6 need only be carried out from 0 to i. Thus,
Tr
F(X,r) = 2 /f I(A,r,6)sin6 cosO d6
0
where 6 is an angle in a single plane.
In order to determine F(A,r) and the F(A,D) we must be
able to calculate I(A,r,6) numerically. This can be
expressed as
I(A,r,6) = m et S(X,re,6e)dt, (9)
0
where S(X,r ,e ) is the source function for radiation of
wavelength X emitted in a direction 6e radians from the
central star, at any point a distance re from the central
star. The integration is carried out over optical depth
along the line of sight, from a point at a distance r from
the central star and in a direction 6 radians from the cen
tral star. The upper limit for the integration, Tm,
represents the maximum optical depth attainable along the
line of sight and will always be less than 2TI. The quan
tities re and 8e are functions of r, 8 and t. Here t is
optical depth along the line of sight and the variable of
integration. Thus we have re = re(O,r,t) and e = e(6,r,t).
The expression for I(A,r,6), giveniin equation (9),
does not help us unless we can evaluate S(A,re,e ). The
source function can be expressed as the ratio of the emission
to the absorption coefficient, that is,
S(X,r ,ee) = j(X,r ee)/K(r ),
where j(X,re,8e) gives the emission per unit volume at a
distance re from the central star, at wavelength X, and in a
direction ee radians from the central star. Since the dust
model is based on the idea that the dust scatters part of the
radiation it intercepts and absorbs and reradiates the rest,
we can write S(A,re, e) in the following form:
S(X,re,Oe) = [jB(A,r, ee) + js(X,re,Oe)] / K(re) (10)
where k(A,re,9e) has been broken into two parts, scattered
emission and thermal or greybody emission, designated
respectively js and jg .
It is evident that both emission terms in equation (10)
are related to the radiation field, I(X,re,6), at re. There
fore, we need to determine the exact functional relationships
expressing jB(X,re,6e) and js(A,r e) in terms of I(X,re,6).
In doing this we will first need to derive expressions for
re = re(O,r,t) and 0e = 8e(8,r,t).
From consideration of the geometry in Figure 1, we
have
re(8,r,t) = (r2 + s2 2rs cosS)1/2, (11)
where s is the metric distance along the line of sight to the
point of emission and
C~
e = sin1 [(r/re)sin8].
The variable of integration along the line of sight, t,
can be evaluated from the expression
s
t = f K(r) ds (12)
o
once we have assigned a value to n in the formula for K(r).
The amount of radiation scattered by a particle will
be that amount which is intercepted by the cross sectional
area of the particle and not absorbed; that is; the amount
intercepted times o The amount of this radiation scattered
in a direction De will, in general, be given by
js(X,re'e)=(kpo /4T) I P(a)I(X,re,O)du, (13a)
where kp is the cross sectional area of the particle and
P(a) is a phase function defined such that
(1/4T)/ P(a)dw = 1
The quantity a is the angle through which the radiation is
scattered and is thus a function of 0 at re and of e (see
e
Figure 1).
Since we are assuming isotropic scattering, where P(a)
= 1 for all a, we can drop ee from our notation and rewrite
equation (13a) as
js(X,re) = (kpwo/4T) I I(A,re,6) dw (13b)
Equation (13b) gives the radiation scattered per particle
and we want the radiation scattered per cm3. To derive the
needed form we must replace k in equation (13b) with the
cross sectional area of all the particles in a cm3, which we
will denote k. Then,
js(A,re) = (kwo/4r) f I(A,re,O)dw
which because of spherical symmetry we can rewrite as
7T
js(A,re) = (kwo/2) I I(A,re,6)sine dO. (14)
0
Two additional equations for js are formed by separat
ing I(X,re,6) into a term due to the cloud and a term due to
the star. Let y(O) be defined such that y(6) = 1 if
0 < < *(re) and y(6) = 0 if 6 > 8*(re), where e*(re) =
1 R.
sin (1 ) is the angular radius of the central star as seen
e
from a distance re; then
I(A,re ) = Ic(A,re,e) + y(6) I*(A)
In the above equation I the term due to the cloud, is given
c
by
TC t
Ic(A,re,9) = f et S(A,r) dt
where r = r(O,re,t) is obtained analogously to re (see equa
tion 11),
c = Tm y() (Tb Tm)
and Tb represents the optical depth along the line of sight
to the inner boundary of the cloud. The term due to the star
is I.(A) = eTbB(X,T*), where B(X, T*) gives the Planckian
emission of the star.
From equation (14) and the above expression for I we
have
7T
js(X,re) = (kow/2)[ / Ic(A,re,O)sine dO +
o
8*
f I,(A)sinO de ], (15a)
0
and if we are only interested in the bolometric emission, we
have
7T
js(re) =(kow/2)[ / I (re,6)sinO dO +
0
6* TeT I oT
I e ( )sin6 de ], (15b)
0
where
0
Ic(re,6) = / I(A,re,e)dX,and a (15c)
o
is the StefanBoltzmann constant.
When 6, is very small we can use the following
approximate versions of equations (15a) and (15b):
Tr
js(X,re) = (kmo/2) I Ic(X,re,O)sin6 dO
0
+ (kmo/4n)exp(TerI)R B(A,T,)/r2 (16a)
and
7T
js(re) = (kwo/2) I Ic(re,e)sin6 dO
+ (kw/4r)exp(eT )R2 oaT/r2 (16b)
where T is the optical depth from the surface of the cloud
to re and the other quantities used have been previously
defined.
We now turn to the problem of determining the grey
body emission of the particles, jB(A,re,6). As long as any
elemental volume is small enough for every particle in it
to experience nearly the same radiation field, the amount of
energy absorbed will be independent of the shape of the
element and its orientation. This absorbed energy will be a
function of I and the effective cross sectional area of the
particles per unit volume, k. Thus, we can derive an expres
sion for k by considering the energy absorbed by an idealized
unit volume; namely a small cylinder of unit height and unit
cross sectional area oriented such that the incident radia
tion is perpendicular to one end.
For such a cylinder the amount of the intensity lost
due to absorption, IL, will be IL = I(1 e ), where T here
is the total optical depth along the axis of the cylinder,
or if we assume the medium to be thin enough so that occulta
tion of one particle by another are rare, we have
IL = (Apart/Acy)I, where Apart and Acy respectively denote
the cross sectional areas of the particles and the cylinder.
Now we can write
part = (1 e)Acy
Apart cy
and since T < < 1, we have
A = [1 (1 T)]A = AC T
part cy cy
or
Apart = AcyK(r) h
where h is a length over which we have an optical depth T and
K(r) has been assumed constant over this length. For a cyl
inder of unit volume we can take Ay = 1 and h = 1; then
Apart/unit vol = K(r), and k(r) = K(r).
We now need to derive a relation between the effective
cross sectional area of the particles in a volume and the
effective surface area enclosing that volume. Let A be the
sum of the effective surface areas of all the individual
particles in a volume element V and let k be the sum of their
effective cross sectional areas. Define the effective cross
sectional area of a single particle, to be its geometrical
cross sectional area averaged over all possible orientations.
Now if the number of particles in V is large, and they have
no preferential orientation; then k will not vary with the
direction of observation w, and we will be able to derive a
relation for A in terms of k; thus defining A.
For the volume, V, the energy absorbed, Ea, is given
by the relation
Ea = I I(C)k d ,
if the density is low enough such that occultations of one
particle by another can be neglected. The energy emitted,
Ee, is given by Ee = OT A. Therefore,
f I(w)k dw = oT4A ,
if the particles are in thermal equilibrium with the radia
tion field. This can be rewritten as
T4 = (k/Ao) f I() dw (17)
and if we let W = k/A, we have
T4 = (W/a) f I(w) dw (18)
This equation is true for any volume of particles and any
radiation field.
If we apply equation (18) to a volume situated in a
blackbody cavity at temperature TB, we have
TB = (W/o) I B(TB) dw,
and since B(TB) is constant with w, we have
4
TB = [W47B(TB)]/a (19)
Since B(TB) is not a function of X, it represents the
Planckian integrated over A. Thus,
B(TB) = i B(X,TB)dX = (aT)/T ,
0
and equation (19) becomes T4 = 4WT4 from which we obtain
B B
W= 1/4, a value indentical to that obtained for a sphere.
Substituting this value for W into equation (18) gives
us
1/4
T =[(1/4a)f I(w)dw]/4
This equation can be used to find the temperature of any
greybody due to an arbitrary radiation field, I(w),
provided that the body is tumbling randomly.
Since W is the ratio of the effective cross sectional
area to the effective surface area of a particle, A = 4kp
for any particle.
From our assumptions that the particles have no pref
erential orientation and that they are in thermal equilib
rium with the radiation field, we have that A = 4kp and
that
F(A,T) = (1 w )iB(B,T)
where A is the effective surface area of a particle, T
is the equilibrium temperature of the particle, B(A,T) is
the Planckian and F(A,T) is the flux per unit effective
area of the particle. The thermal isotropic emission of a
single particle can now be written as
jB(X) = (A/4T)F(X,T)
or
jB(x) = k (lwo)B(A,T)
The corresponding emission for a unit volume of such
particles will be
jB(A,re) = K(re) (W )B(A,T)
Note that the temperature T, like K, is a function of the
distance from the central star, r.
If this emission is considered only bolometrically,
F(T) = (1o) aT4
jB(r ) = [K(re)/F] (lWo)oT4 .(20)
The condition of thermal equilibrium gives us the
relation
Ea(r) = Ee(r) (21)
where E (r) and E (r) denote,respectively,the energy absorbed
a e
and emitted by a unit volume of particles at a distance r
from the central star. Now from equations (20) and (21) we
have
Ea(r) = I (1oo) (0/)K(r)T4 (r)dw ,
and since the emission is isotropic
Ea(r) = 4K(r) (lw )oT4(r)
From this formula for E we can obtain the following
expression for the temperature,
T(r) = {[1/4a(lmo)] [E (r)/K(r)]}1~/4
For the nonbolometric case we have,
E (r) = K(r) (10 ) ff I(A,r,0)dXdw ,
a mA
giving us
1/4
T(r) = [(1/40) ff I(X,r,0)dXdw]
Because of spherical symmetry we can write this as
iT 1/4
T(r) = [(T/20) /f I(A,r,6)sin6 dX dO]
oA
If we divide the radiation field into two parts, as
we did before, we have
7T
T(r) = {(1/40)[2r f I (r,0)sin0 dO +
o
OT,4 TTI 1/4
/ e sin6 do]} (22)
o
where Ic(r,O) is the bolometric intensity defined by equation
(15c). When the central star can be considered a point
source we have the additional equation,
T(r) = {(1/4a) [27 / I (r,O)sin8 dO +
o
(TTI) 2 4 2 1/4
e (R*oT/r2)]}
Note that in equation (9) I is a functional of S and
that in equation (10a) S is a function of j. We have seen
that j is a functional of I; therefore, I is a functional of
S and S is a functional of I. Thus, equations (9) and (10a)
can be rewritten as a recurrence formula for I(A,re,w), which
can be used in an iterative solution for the radiation field
throughout the cloud and hence for the desired spectral
energy distribution.
The recurrence formula is
In(A,r,w) = fm et S[ I In_1(A,re,w)dw]dt, (23)
o o
where the subscripts n and n+l refer to iterations n and
n+l. The true radiation field consistent with a given set
of parameters will be given by
I(A,r,w) = lim T (,r,w)
nt
This completes the definition and derivation of the
basic quantities needed to implement an iterative solution
of equation (11).
Once the proper set of the parameters, TI, R R and
n, is determined, we can obtain values for the density with
in the cloud, p(r), and the total mass of the cloud M.
If we let D equal the number of particles per cm and
d equal the average diameter of a particle, where d < < R,
then
D = CSA/cm3 = particles/cm3
CSA/particle
The cross sectional area of the particles per cubic centi
meter is a function of r and is given by K(r); thus,
D(r) = K(r) = 4 K(r)d2
(2 2
For the volume and mass of a particle, Vp and m respective
ly, we have
V= 4 d(d)3 cm3/particle,
p 3 2
and
m p ( ) grams/particle ,
where p is the specific gravity of a grain. The density of
the cloud at a distance r from the central star will then be
given by
p(r) = D(r) mp
4 d 3 K(r) or
p(r) n( ) p d or
3 2 P (d) 2
p(r) = p K(r) (d)
Substituting for K(r), we find
2 T (ln)
p(r) p d [ ] 1 (24)
3 P RlnRln
I
From the above expression for p(r) we can derive an
expression for the total mass of the cloud. We have
R
M = I (47r2)p(r)dr ,
RI
from which we obtain,
M = 8ppTI(1n) (R3n R 3n)d
3(3n)(Rl l(25)
3(3n)(Rln Ri1n)
SECTION III
NUMERICAL PROCEDURE
In this section a procedure based on the equations of
Section II is developed, to effect a numericsl solution for
the radiation field in the dust shell. This procedure will
involve iterative numerical integration of equation (23)
to obtain a set of convergent series of n(A, ri,w) Here i
specifies the ith such set and ri, therefore specifies the
distance from the central star at which the I(A,ri,w)n for
that set are to be evaluated. The n signifies that this is
a value obtained in the nth iteration. All the integration
will be performed using GaussLegendre formulas.
In order to accomplish the necessary integration it
will first be necessary to establish an integration grid in
keeping with the boundary conditions imposed by the geometry
of the shell and the central star. Consider the spherical
shell of dust extending from RI to R and divided into N
thinner concentric shells each with the same optical thick
ness, TI/N (see Figure 2 ). The inner and outer radii of
these shells are a function of T and can thus be found from
the equation
r(T) = [R1n (T/T) (Rln r n) 1(26)
ln (26)
For each of the N points lying on the semidiameter of
the cloud, SR, (see Figure 2) at distance ri from the central
star, we assign M values of the angle 6, which will define M
lines radiating from each point. We then break these
radial lines into segments of approximately equal optical
depth, where the end points of the segments are all the
points at which the radial lines intersect the N circles of
radii ri.
Thus, about each of the N points on SR we have M
radial lines, and each of these lines is broken into from 1
to 2N segments, where the actual number of segments for any
particular line depends on how many of the N circles of
radii ri that line intersects.
The M values of 8 used at each of the N points along
SR would in general be the GaussLegendre points for the
interval 0 to i; however, it is usually advantageous to
break this interval into two intervals with M/2 Gaussian
points each one from 0 to 9m(ri) and the other from 8m to ir,
where 6m(ri) is an angle that must be determined for each of
the N points ri. If ri is small enough so that the central
star shows an appreciable disk, then 6m(ri)=9*(ri), but if ri
is large enough so that the star can be considered a point
source Om(ri) will be an angle between 0 and T/2. The value
of 6m(ri) in this case should be determined, such that, it
lies between u/2 and the value of 9 at which I(ri,6) has its
greatest value. Defining 9m(ri) in this manner allows us to
concentrate our Gaussian points on the part of the interval
0 to a where I(0) is varying most rapidly. In one special
case the M Gaussian points must be determined in a different
manner. When ri = rn = R the M O's are restricted to the
interval 0 to f/2 since for i/2 < 6 < i, I(R,6) obviously
vanishes.
We have broken each of the M radial lines into segments
as a first step in the determination of the locations of L
Gaussian points to be used in the numerical evaluation of
equation (23). Figure (2) shows a few representative radial
lines for the ith of the N points on SR. Also, the radial
lines are shown broken into segments with end points a
distance sk from R, where sk is measured along the radial
line and k is an index ranging from 1 to 2N.
Let x and y be the rectangular coordinates of a point
in a righthanded coordinate system, where the central star
is the origin, and the semidiameter SR is the positive
yaxis. Also, let the slope and yintercept of a radial line
be denoted by m and b respectively; then, each pair of
values for m and b defines a particular radial line. In
terms of previously defined quantities, m = Cot (6) and
b = ri When b and m are specified, the rectangular coordi
nates of all the points where a particular radial line inter
sects the N circles of radii r. can be found by solving, for
x and y, the system of equations
x2 + y2 ri2 = 0
y mx b = 0, for i = 1 to N.
We are only interested in points to one side of SR, so we
exclude all solutions with x > 0. Points which are occulted
by the star are also excluded.
The starting point of a radial line will lie on SR,
and have coordinates x = 0 and y = ri; the kth point of
intersection on that line will have coordinates xk and yk'
Thus, the distance along the line from the starting point to
the kth point, which is denoted sk in Figure 2 is
sk = [xk2 + (rik) 2]/2
To accomplish the integration of equation (23) we must
first perform the integration indicated by equation (12)
along each of the radial lines using the end points of the
segments on these lines as the limits of integration.
Equation (11) is used here to find the distance from the
central star to a point on the line and hence the values of
K(r) needed by (12) from (5) or (7a).
From (12) we will obtain for each of the N x M lines of
sight (LOS) a set of accurate optical depths, tk, and cor
responding distances, sk, measured along the LOS. Because
of the discontinuity caused by the central region of the
dust cloud (r < RNR) it is helpful to break the LOS into two
intervals to consist of L/2 Gaussian points.
Unless the LOS intersects the central region of the
cloud (r = rl = RNR), the midpoint of each LOS, at which we
can break the LOS into two segments, is roughly determined
by dividing the total number of points by two. However, if
the shell rl is crossed two midpointss" are established
(say, MID1 and MID2) where the first such point is the point
of intersection of the LOS and rl closest to the starting
point and the second is the next point further along the LOS.
The segments from the beginning of the LOS to MID1 and from
MID2 to the intersection of the LOS with rl are to be
considered the two halves of the LOS (see Figure 2 ). One
additional case remains; if the LOS crosses the shell rl and
then intersects the surface of the star (r = R*) the LOS is
to be considered terminated at rI This will obviously
occur whenever 6(ri) < e*(ri) for a LOS.
Next each half of the LOS is broken into 19 segments
(this is the specific number used in the present version of
the computer code it can be varied if necessary) as nearly
the same length in optical depth as computing time will
allow (the computer code is given a parameter which sets the
precision desired). This is done by an iterative procedure
in which the computer code uses an AitkenLagrange inter
polation routine (Hildebrand, 1956; Mullen, 1974) to inter
polate in the table of sk and tk values for the set of sk
values corresponding to a new set of evenly spaced tk values.
A Gaussian integration routine is then used to evaluate
equation (12),thus redetermining the values of tk If the
values of tk found by this integration are not equally
spaced (to the precision requested) another interpolation
and integration are performed, and this process is continued
until the desired accuracy in spacing is achieved.
If the original set of sk and tk values contains only
two points (the case where the LOS starts at the outer
boundary rn and then intersects rn without intersecting
rnl) the procedure described above breaks down since a set
of two points can have no midpoint. In this special case
the code sets up a more extensive table containing 20 points
(midpoint at 10) by breaking the interval s2 to sl into 19
equal segments thus obtaining 20 sk values and then finds
the corresponding tk values by integration. The resulting
tables of 20 sk values (equally spaced) and tk values (un
equally spaced unless n = 0 in density relation) are then
used in the iterative process to find the set of tk values
at equal spacing in optical depth.
Since we want to evaluate equation (23) by Guassian
quadrature we must be able to evaluate the integral at the
Gaussian points, tkg, on the interval 0 to 1,; therefore,
we next interpolate in our table of sk and tk values for
the sk values, designated skg corresponding to the tkg
values.
We have that S(t) = S(T(t)), where is the optical
depth measured from the outer boundary of the cloud to the
point specified by t; therefore, we want to find those
values of T which correspond to the tkg. From the skg
determined by interpolation we can calculate, using the law
of cosines, a set of corresponding rkg values and from them
a set of Tkg values. We can then find the values of S(Tkg)
by interpolating in a table of S(Ti), where Ti is the optical
depth at the ith of the N shells. All the interpolations
mentioned above are performed by the same AitkinLagrange
interpolation routine.
The table of N values of S(Ti) is of course based on
the previous iteration of the code.
The computation to determine sets of T(ri), K(ri),
.S(ri,6,X), I(ri,6,X) and ultimately F(rN,X) starts from a
set of initial values of T(ri), which should be as close to
the expected final values as possible, and a set of the
presumed physical parameters. The set of K(ri) can be found
immediately from the given parameters using either equation
(5) or (7a). Using the K(ri), the initial set of T(ri) and
their associated optical depths, we can evaluate the integral
in equation (23) along each of the M lines radiating from
each of our N points. Then from the M values of I(ri,6,X)
obtained about each of the N points, improved sets of T and
S can be calculated.
This process is repeated until successive sets of
T(ri) converge to within a few tenths of a degree. The sets
of T, K, S, and I obtained in the last iteration are taken
as the true values for the cloud. A set of F(rN,A) is
obtained from the I(rN,6,X) using equation (8).
To save computation time a bolometric case with =
0.0 is run to determine the set of T(ri). If in the numer
ical procedure, we are performing the calculations for L
values of A, then the computation time will be cut by a
factor of approximately 1/L when we drop the A dependence
and run a bolometric computation. When we use the bolometric
expressions instead of their wavelength dependent analogs,
equation (23) becomes the bolometric expression
0
I(r,9)n+1 = I et S[I(r,6)n] dt.
Like equation (23) this can be solved iteratively for the
specific intensity. The set of T(ri) is then found from
equation (22) and used as the set of initial values in the
desired wavelenth dependent solution.
Note that we can assume that the temperatures deter
mined by a bolometric computation are the same as those
that would result from a detailed many wavelength computa
tion only if the particles in the cloud are greybodies. If
we want to consider cases where wo = wo(X) then we can not
use a bolometric computation to determine the temperature
distribution.
The basic equations in Section II which involve inte
gration over 8 can be rewritten such that the integration
are over p = cosO by making the change of variable
a = cosl (); dB = dp/sin6 The code has been tested
with these equations programed with 6 as the variable of
integration and with cos6 as the variable of integration 
the results do not differ significantly.
When the central star can not be considered a point
source it is necessary to integrate over the angular diameter
of the star. Since the specific intensity of the stars'sur
face can be considered constant (we can ignore the small
effects of limb darkening) an integration from 6 = 0 to
8 = ,*(ri) need only be done once for each ri and used in
place of the relevant point source approximation.
For a shell very close to the star the correct contri
bution by the star to the mean intensity and the approximate
value can differ by as much as a factor of 2.0 the point 
source approximation giving the lower value. The mean
intensity, J, the Eddington flux, H, and the socalled K
integral, K, are all moments of the mean intensity I (Mihalas,
1970) that is, M(n) = 1 I()PndVJ,
21
where I = M(0), H = M(1), and K = M(2). The flux F should
not be confused with H; they are related such that F = 4wH.
When RNR >> R* we need not worry about inaccuracies due to
approximating the star as a point source. For example, when
R, % 0.2 x 1013 and RNR 4 0.15 x 1014 we find an error due
to the approximation of only 0.46%.
Since the code is capable of evaluating the specific
intensity as a function of 6 and X at any point in the shell
it can be used to find values of Jx(T), H (T), Kx(T) and the
Eddington factor, KX(T)/JX(T) The Eddington factor can
vary from 0.0 to 1.0 but is usually assumed to be 1/3 (the
Eddington approximation, K = 1/3 J) since this is the value
it would have in an isotropic radiation field.
SECTION IV
RESULTS OF MODEL CALCULATIONS
The computer code used to determine the radiation
field within a circumstellar dust shell is broken into
three main programs (program 1, program 2, and program 3).
The first is a program used to set up the integration grid
and store it in the computer (usually on disk storage), the
second is a program used to determine the bolometric radia
tion field within a dust shell and to find the temperature
distribution, and the third program is used to calculate
the monochromatic radiation fields for a number of wave
lengths using the temperature distribution determined by
the second program. The three programs must initially be
run in the order given but program 1 need only be run once
for each set of the geometric parameters R., RNR, R, TI,
and N. Therefore it is possible to run a number of differ
ent models with the same integration grid by varying the
two nongeometrical parameters T, and o .
The models that we want to consider here are charac
terized by the parameters in Table 1.
Table 1. Values of the seven parameters used for each of
the models computed (cgs units; T* in oK)
Model R, RNR R TI n Wo T,
I .202x1013 .15x1014 .75x1015 5.0 1.5 0.1 6000
II 0.0 "
III 1.0 "
IV 0.6 "
V 7.0 0.0 "
VI 0.1
VII .202x1013 .202202x1013 0.0
VIII .625x1015 .75x1015 2.0 2.0 0.0
IX .15x1014 5.0 1.3 0.1
X .8352x1012 .8352x1015 3.97 0.0 0.8 12500
One additional model, model XI, will be described later in
this section.
From Table 1 we can see that the most common values of
Wo used are 0.0 and 0.1. This is because the amount of
computing time necessary to determine the radiation field
is a function of o For small values of o()o<.2) program
3 requires approximately six iterations to determine values
of the flux, for example, converged to three significant
figures, whereas, if Mo = 1 the program will need at least
15 iterations to achieve the same accuracy. Because of the
high cost of models with large values of wo only two are
computed, models IV and X. Model III was computed for only
one wavelength and therefore is not comparable in cost to
models IV and X which were computed for 11 wavelength. When
to = 1.0 we can save time by running the program for only
one wavelength, Al say, and then multiplying the quantities
obtained for that wavelength by the factor B(12,T,)/B(AI,T*)
to obtain the quantities appropriate to X2. This is a
valid procedure when wo = 1.0 since under that condition all
the radiation leaving the star must reach the surface of the
shell unchanged in its spectralenergy distribution.
When program 2 is run the total outgoing bolometric
flux through each shell is computed at each iteration. Since
the only source of radiation is the central star, the total
outgoing bolometric flux at each shell should be equal to
the bolometric luminosity of the star. This is true to within
0.5% for some of the models in Table 1 and to within 4% for
all the models in Table 1 except model VII. The flux
(bolometric) emitted by the outer boundary of the shell in
model VII is too low by 23% principally because of the very
severe geometry used in model VII. This model was run as a
test of the code to determine its accuracy under extreme
conditions.
To determine the temperature distribution in the shell
program 2 is run. The temperature distributions for the
models in Table 1 are shown in Figure 3. Note that the
curve for model IX is very nearly a straight line. This is
due to the geometry of the shell and the fact that, as we
shall see, the radiation field in the shell is highly
anisotropic. For a highly anisotropic radiation field all
moments of the intensity are equal, that is J = H = K. Thus,
1/4 1/4
J F and we know that T a J/ so T F/4. Because of
conservation of flux in a spherical shell we have that
F l/r2; therefore, T (1/r2)1/4 or T l/r05. There
fore, for n = 1.5 we have from equation (6) that T(r) =
C1 + C2/r05, where C1 and C2 are constants. This means
that
T(l/r0.5) = C1 + C2 (l/r05)
and, therefore we can see that T will be a linear relation
of T.
The temperatures of the dust grains as mentioned in
Section I are consistent with the assumption that they are
carbon grains where the temperature is below about 21000K
and consistent with the assumption that they are refractory
silicates for temperatures less than about 1500K. Note
that the spread in temperatures at T = 0 is small even
though some of the models have very different physical
parameters.
Program 3 is run to determine a variety of quantities
for the cloud. These quantities are: I(X,T,6), J(A,T),
K(X,T)/J(X,T), and F(A,T). Figures 4 and 5 show some
representative values of these quantities (computed for
model IX). The conversion of visual to infrared radiation
can be seen on Figure 4 by following the relative intensi
ties of a few wavelengths from T = 7 to T = 0. Figure 5
also shows this by way of the graphs for J and F. The graphs
for K/J show that the radiation field is generally anisotropic
although at depth the infrared radiation field does approach
the Eddengton approximation, K/J = 1/3. Thus, the use by
Huang and Apruzese of the Eddington approximation can be ex
pected to decrease the accuracy of their results significant
ly.
In Figure 6 the run of temperature with optical depth
is shown for two models calculated with this code and for two
models calculated by Huang (1971). The two models calculated
with this code are models VIII and XI. Model XI has the fol
lowing parameters: I = 2.32, R* = .202 x 1013 cm,
RNR = .609 x 1015 cm, R = .75 x 1015 cm, n = 2.0 and
T* = 60000K. For optical depths larger than about 2/3 the
curves for models VIII and XI are very similar the curves for
Huang's corresponding models. From optical depth 2/3 to 0
the curves for models VIII and XI drop off much more rapidly
than Huang's curves. This is just what we would expect to
happen due to the erroneous value of the Eddington factor
used by Huang for the outer part of the cloud.
Figure 7 shows the spectralenergy distribution for
HD 45677 as given by Swings and Allen (1971) and Low et al.
(1970). On the same graph the spectralenergy distribution
for model X is also plotted. Model X is based on a set of
parameters derived by Apruzese (1974) for HD 45677. He
determined these parameters by fitting the observations in
the visual since his model is not capable of directly
fitting the infrared observations (actually his model does
use the infrared data in that he attempts to make the total
infrared emission calculated by his model equal to the
total infrared emission given by the observations). The
graphs show that his parameters do give a very good fit in
the visual region and that the total emission in the
infrared predicted by his parameters is approximately cor
rect; however, his parameters do not give a good fit to the
actual infrared spectralenergy distribution. A different
set of parameters could be found that would fit both the
visual and the infrared regions.
Figures 8 and 9 show the spectral energy distributions
for a number of the models listed in Table 1. On Figure 8
the graphs for models II, V, VIII and X are shown. The
curve for model VI would fall between those for models II
and V. We can see the effect of changing the various model
parameters by studying these curves.
The curves for models II and V show that increasing
the total optical depth of the shell causes the peak of the
distribution to move toward the red (Figure 3 shows that at
the same time the temperature of the inner part of the shell
goes up) and causes the amount of visual radiation filtering
through to decrease by a large amount.
Models VIII and X show, as we would expect, that when
the optical depth of the shell is decreased by a large
amount the central star causes the distribution to peak in
the visual.
Figure 9 shows the spectralenergy distributions for
models I, IV, IX and the means of the observations of R Mon
due to Mendoza (1966, 1968), Low and Smith (1966) and Low
et al. (1970), where vertical bars show the amount of
variability from those means. The vertical error bars on
the observations at 20 and 22V indicate the rootmean
square error of a single observation, approximately 20%.
The U through M values plotted are intensities
calculated from Mendoza's observations (corrected for inter
stellar extinction), using Johnson's (1966) calibration of
the UBVRIJKLM system. The interstellar extinction correc
tion was made, assuming that R Monocerotis is located in
NGC 2264, and using the following adopted data for that
cluster and hence R Monocerotis: EBV = 0.06 (Strom et al.
1971) and a reddening curve of the NGC 2244 type with R = 5.4
(Johnson, 1968).
Using these values we have AV = 0.324, and from the
reddening curve we obtain the extinction in magnitudes, AM,
for each of the observed colors (see Table 2).
The rootmeansquare errors of a single observation
are shown for the 201 and the 22p observations in Figure 9,
but for all the other bands the rootmeansquare errors of
a single observation are negligible, considering the
intrinsic variability of the object. The largest such er
ror is for the M band, and it amounts to less than 0.1
magnitudes.
Table 2. Interstellar extinction in magnitudes, Am, for R
Monocerotis in each of the colors observed
Band Am
0.424
0.383
0.324
0.277
0.224
0.189
0.141
0.094
0.083
The curve for model IV shows that a large increase
in the value of wo causes the peak of the distribution to
narrow and move toward the blue, it also shows that the
direct contribution of the star in the visual increases
greatly. The curves for models I and IX show that a small
increase in the index n causes the peak of the distribu
tion to move toward the blue, while the ends of the curve
in the visual and far infrared drop.
Model I is the best fit to the observations of R Mon
but it looks as though Model IX with a slightly larger wo
would fit somewhat better.
The following discussion pertains to the development
of the model for R Mon.
The spectral class of R Monocerotis, roughly G (Joy,
1945), and luminosity of 870 L (Low and Smith, 1966) lead
to the values of T* and R, chosen for these values since
4
4IRR*oT* 870L The value for TI was determined from a
cursory examination of the observations in the following
manner. We have the ratio of the visual flux, FV, to the
bolometric flux, F for R Monocerotis; namely, Fv/FB .006
(Low and Smith, 1966). From this and the relation
TI = Zn(FV/FB) we obtain for T(RI); TI ~ 5.
It should be noted that this method gives only a
lower limit to TI, since we cannot say what portion of FV
is actually coherently scattered radiation. Thus, wo and
T cannot be determined separately with any ease. The
values of w0, RI, R and n used were just rough initial
approximations.
We must consider the following observational data
for R Monocerotis: (a) the spectrum does not show selective
absorption (Joy, 1945), (b) strong variable optical
polarization exists indicating scattering by large grains
(Zellner, 1970), and (c) irregular light variations occur
of several magnitudes (Joy, 1945; Mendoza, 1966, 1968; Low
et al. 1970). In addition we have the UBVRIJKLMQ observa
tions of Mendoza, Low and Smith, and Low et al. mentioned
previously.
Points (a) and (b) above are consistent with the as
sumption that the cloud is made up of solid particles which
act as grey absorbers and emitters. For a period of
variation of the polarization which is fairly short, the as
sumption that the particles are randomly oriented should be
reasonable. Because of the variability mentioned in point
(c), those observations in several colors, which were made
on the same day are of the greatest interest, and where such
observations are lacking, we are forced to compare the
computed fluxes to mean observed fluxes.
Using the parameters from Case I given by
TI = 5 ,
RI = .15 x 1014 cm
R = .75 x 1015 cm,
n =1.5,
and assuming the particles to have the density of carbon or
silicon dioxide; that is, pp ~ 2.0, we have from equation
(25) for the mass of the cloud.
M = 13dM
This gives us for particles about twenty microns in
diameter a cloud mass of about 102 M .
The set of parameters for Case I also gives us for
equation (24) the following form:
d
p(r) = (1.5 x 109) r1.5
If d = 20p and r ranges from 1013 to 1015 cm, then p(r)
will range from approximately 1013 to 1016 grams/cm.
In summary, the best fit to the observations was
obtained by a model with the following parameters:
RI = .15 x 1014 cm (lAu) ,
__
66
R = .75 x 1015 cm (50 Au)
T =5,
= 0.1
n = 1.5
M~ 102M
1016 < p(r) < 1013 gm/cm3
2250K < T < 2,1000K .
~
Radiation is shown falling on point 1 at
distance r from the central star. It is
incident there at an angle 01 after being
emitted from point 2 at an angle 0 .Part
of this emitted radiation is radiation
from point 3 which was incident on point
2 at an angle e2, and scattered through
an angle a. Point 2 is a distance re
from the central star and a distance s(t)
from point 1. Analogous triangles, with
vertices at the star, point 1 and any
point on the line directed from point 1,
at the angle 81, can be constructed.
Note that the angle 61 is designated
simply 6 in the text.
Figure 1.
/
/
/
/
/
/
/
/
/

f 
STAR
A schematic diagram of the integration grid
for the case when N = 4. Four lines of sight
are shown to indicate the manner in which
the boundary conditions at rl and rn are
handled. The filled circles on the lines
represent the initial grid of intersection
points. The tick marks on the lines
represent the points found at equal inter
vals of T by an iterative technique. The
open circles represent the GaussLegendre
points for a three point quadrature formula.
Figure 2.
U)
0
H
0)
i,
a)
0
'iU
> El
tI
.10
>
F 0
Qo
.rl
CO
H a
0o
 m
OCi
ta
I
re
00
Hi
Os
.5La
n
0
0
H
Log [I(x,6,T)] vs. 6 for model IX. The angle
0 is in radians. The three panels from top
to bottom are for T = 7.0; 3.68; 0.0. The
curve are labeled with letters designating
their respective wavelengths on the
UBVRIJKLMNQ system. Note that in the case
for T = 7.0 the star is not considered a
point.
Figure 4.
v
M
r 7.0
13
12 K M
10
11 : 3.6 8
1M M
10 .
t 2
Figure 5. Log [J(A,T)], K(X,T)/J(X,r), and Log [F(X,T)]
vs. T for model IX. The curves are labeled
with letters designating their respective wave
lengths on the UBVRIJKLMNQ system.
km
 a)HOa)H r
rd x ? a) o a c a)
0 O4O a 0a4 n
0 H r I 4 I 0 ,I .r4 0
. a) a) 0 R a) P E 4
$ 44rO O 9 0 ., E z ; >
P 4J 0 0 H4 00 0f a 
a) OO ) a H o GtH 0 a
4 ar d 0>axP c X
E 'Q I 4 HO N *HC
4J 0 o ' 0 0 H O
4' o at) *000W
14 k0 I O, ( 0 ,
a 0 00 OWOOd
P a) n tM. 4 2 S 5 4
OXO (d H rid rd rd 0
0 O rl P r ib U 0
*Ha r l O 4J H H
.,I C 4.)a IQ I 4 n a)H , 0
Srl rd X (d  > H Eo H
NO lOH 0 >H .
() UHl H H :J EC
SH a)H a)r O mH
H ) >M H > \
o>.d 9> 0 p r  C)N
4o. J W ) O O 0 r1
'd CW 0 OMH Oa 'a
SOU()ri 4J R > 4)
* rI a) O4 P 0 C Nl
(a Q) 0 a) 3 a) ri 0 En rd 0
4 Q'0 m(0 4 0 rO Q ,'. 4r.4 CN
w 0 4 E2 ) rd
> H Hr 4) )
O H rd 4 rd> c
4HN kO U UNE Hr
O4 0 > W L*H W H (d 3
E1 4 ) 0 ,C 4 U044
N\ H rmn 04 0
H 0 0n 0 *H 0 H 0a N
H C0 4 a r 0 H OO
H 0 OHOE 2Ca'H >H O'
O4 O P 0 H 00 4 OH
O) O H H w r: H
4J *a oa x l OnN (' OHr
10 a rl Icl 0H ( d H Win a, ^
SI orig oEra i l
aN W )k a (d W E 140
2 *HO i(d.Ct N 2 H 0
) O r0 41 a) > (3 d ad
 0 4 > NQA rd a
4a I ) a) (1) p >'44 ric
rWWEn Z0 > 04 00 0
Oa) if U >0VCN 0 ) U
N P d ) 4N 0H 4 U0 > r0
rH H (Z 12O44 l0) 44 a lJ
N 400 (0 O :I Q4
(d U aWr) PE 4 H 3rH a OO
ERl 124 a) d > P 0U Q 43
N 4J r4 Vra N V rH HO V 4 0 N
0 04 0 rd 9 Z nd ni 0 0 X 0
Z 04P 04' d u 04F 2caJt4
I
*dl
h
oJ c
Im
ra 0 a)
$1 >
U *.
04 > 4
4n) a)4
r4 r 
u 141 >
> Pi
0 4 34
cd m 0
41 (L > ,
J 0 0 ,
S0 'a .c
O H1
rd c i
0 0 H
*H )
) $4 ) W 4V
tP 0 40 03
w 4
aV Q 4 4
e0 Ol ,
X24)
H 0i '
0 Ord
> 0 00
W ) d .0
0 04 *
H0 4) 0I
$02 04)V
)HH c a) a)
SWH 0, a
r0
H
6>1
4 8 a 3 3
r.
H
a
o
0
0
0.I
0)
'a
H
co
H
in
CO ~
11
nF
82
> O0
,l d
o c (o (/ (\ o
. o 6 L< o 6 o
___
4J 4 )
'O
a) ) 0 n
I I a
H M n *
0 I 0
s 1 0)H (
a >7 on
0 U 0 0*
acao Sa
a) 0U 4)
4o 0 
5c 044 1 41
0)0 l 0)
o o
0< .. C 4)
IJ 0 0 rI
0 0) o
4O H $ I c
H 0) >
rO
(a r4 V .0
to 0 > a
44 > z)
0 c 1
0 0 (n 0)
<0 H *()
S 0 E41 ,A
140 r 00
tUL2 I
0)
4) CIL 4
*: 0 1 aS)
Un H (a rc
4, 414
.1
l
0 CDO
6 6 6
O
0
40
0 0 4"
0 H 0
On 
0 *0
) a)o
0o i 1
m
00
m m
01 H1
0o
OM
0 0
ONa
0 N I
u0 0
0
O m
0 > 0
0i
H 0
Oo
0 0
4> W
*H o H
0' 0 C
o I u
4 )i cI
.q 0 rU
aC
0'
cH
** u
u)i *
* () +
n 4 1 n
N )
I I I I IT
L I_ I I
0
CO
6
0
6
10
or
00
6
LIST OF REFERENCES
Ackermann, G., Fugmann, F., Hermann, W., and Voelcker, K.
1968, Z. Ap., 69, 130.
1970, Astron, Ap., 8, 315.
Anderson, L., and Kuhi, L.V. 1969, NonPeriodic Phenomena
in Variable Stars, ed. L. Detre (DordrechtHolland:
D. Reidel Publishing Co.), p. 93.
Apruzese, J.P. 1974, Ap. J., 188, 539.
Breger, M. 1972, Ap. J., 171, 539.
Breger, M., and Dyck, H.M. 1972, Ap. J., 175, 127.
Capps, R.W., and Dyck, H.M. 1972, Ap. J., 175, 693.
Chandrasekhar, S. 1934, M.N.R.A.S., 94, 444.
Chavira, E. 1967, Bol. Obs. Tonantzintla y Tacubaya, 4, 29.
Cohen, M. 1973a, M.N.R.A.S., 161, 85.
1973b, ibid., p. 97.
1973c, ibid., p. 105.
Gilman, R.C. 1969, Ap. J., 155, L185.
Grasdalen, G.L., and Gaustad, J.E. 1971, A.J., 76, 231.
Hayashi, C. 1966, Ann. Rev. Astr. and Ap., 4, 171.
Herbig, G.H. 1960.
1968, Ap.J., 152, 439.
1970, ibid, 162, 557.
Hetzler, C. 1937, ibid., 86, 509.
Hildebrand, F.B. 1956, Introduction to Numerical Analysis
(New York/Toronto/London: McGrawHill), p. 49.
Huang, SuShu 1969a, Ap. J., 157, 835.
1969b, ibid., 157, 843.
Huang, SuShu, 1971, ibid., 164, 91.
Hummer, D.G., and Rybicki, G.B. 1971, M.N.R.A.S., 152, 1.
Hyland, A.R., Becklin, E.E., Neugebauer, G., and Wallersten,
G. 1969, Ap. J., 158, 619.
Johnson, H.L. 1966, Ann. Rev. Astr. and Ap., 4, 193.
1968, Nebulae and Interstellar Matter, ed. B.M.
Middlehurst and L.H. Aller (Chicago: University of
Chicago Press), p. 202.
Joy, A.H. 1945, Ap. J., 102, 168.
Krishna Swamy, K.S. 1970, Ap. and Space Science, 9, 123.
Kukarkin, B.V., Kholopov, P.N., Efremov, Yu, N., Kukarkina,
N.P., Kurochkin, N.E. Medvedea, G.I., Perova, N.B.,
Federovich, M.S., and Frolov, M.S. 1969, General
Catalogue of Variable Stars, Third Edition (Moscow).
Larson, R.B. 1969a, M.N.R.A.S., 145, 271.
1969b, ibid., 145, 297.
Low, F.J., and Smith, B.J. 1966, Nature, 212, 675.
1970, SemiAnn. Tech. Rep. Contract No. F 19628
70C0046 Project No. 5130 (ARPA).
Low, F.J., Johnson, H.L., Kleinmann, D.E., Latham, A.S.,
and Geisel, S.L. 1970, Ap. J., 160, 531.
Mendoza, V.E.E. 1966, Ap. J., 143, 1010.
1968, Ap. J., 151, 977.
Mihalas, D. 1970, Stellar Atmospheres (San Francisco: W.H.
Freeman and Company), p. 5.
Mullen, J. 1974 (Private Communication).
Neugebauer, G., Becklin, E., and Hyland, A.R. 1971, Ann.
Rev. Astr. and Ap., 9, 67.
Neugebauer, G., and Leighton, R.B. 1969, TwoMicron Sky
Survey a Preliminary Catalog (NASA SP3047).
Poveda, A. 1965a, Bol. Obs. Tonantzintla y Tacubaya, 4, 15.
1965b, ibid., 26, 15.
1965c, ibid., 26, 22.
Prentice, A.F.R., and terHaar, D. 1971, M.N.R.A.S., 151, 177.
Serkowski, K. 1971, Proceedings of the Conference on Late
Type Stars, ed. G.W. Lockwood and H.M. Dyck (Kitt Peak
Contr. No. 554), p. 107.
Stein, W. 1966, Ap. J., 145, 101.
Strom, K.M., Strom, S.E., and Yost, J. 1971, Ap. J., 165,
479.
Strom, S.E. 1972, PASP, 84, 745.
Strom, S.E., Strom, K.M., Brooke, A.L., Breger, J., and Yost,
J. 1972, Ap. J., 171, 267.
Swings, J.P., and Allen, D.A. 1971, Ap. J., 167, L41.
Walker, M.F. 1969, NonPeriodic Phenomena in Variable Stars,
ed. L. Detre (DordrechtHolland: Reidel Publishing Co.),
p. 103.
Walker, M.F. 1972, Ap. J., 175, 89.
Woolf, N.J., Stein, W.A., and Strittmatter, P.A. 1970,
Astron. Ap., 9, 252.
Zellner, B. 1970, A.J., 75, 182.
Zellner, B.H., and Serkowski, K. 1972, PASP, 84, 619.
BIOGRAPHICAL SKETCH
Christopher Alvin Harvel was born on December 7, 1944
in Columbia, South Carolina. He attended public schools in
North Carolina, California, New Jersey, and Virginia, gradu
ating from McLean High School in McLean Virginia in 1963.
Four years later, in 1967, he received his Bachelor of
Science degree with a major in Astronomy from Georgetown
University in Washington, D.C. and started graduate school
at the University of South Florida.
In December of 1968 he was drafted into the United
States Army and in December of 1969 he began a oneyear tour
of duty in the Republic of South Vietnam.
At the end of his military service in December of
1970, he returned to the University of South Florida and
received a Master of Arts degree in Astronomy in 1972. In
that same year he started graduate school at the University
of Florida where receipt of the Ph.D. degree should occur in
August, 1974.
I certify that I have read this study and that in my
opinion it conforms to acceptable standards of scholarly
presentation and is fully adequate, in scope and quality,
as a dissertation for the degree of Doctor of Philosophy.
Frank Bradshaw Wood, Chairman
Professor of Physics and
Astronomy
I certify that I have read this study and that in my
opinion it conforms to acceptable standards of scholarly
presentation and is fully adequate, in scope and quality,
as a dissertation for the degree of Doctor of Philosophy.
Sabatino S. Sofia
Professor of Astronomy
University of South Florida
I certify that I have read this study and that in my
opinion it conforms to acceptable standards of scholarly
presentation and is fully adequate, in scope and quality,
as a dissertation for the degree of Doctor of Philosophy.
K hen
KwahYu en
Associate Professor of Physics
and Astronomy
I certify that I have read this study and that in my
opinion it conforms to acceptable standards of scholarly
presentation and is fully adequate, in scope and quality,
as a dissertation for the degree of Doctor of Philosophy.
Howard L. Cohen
Associate Professor of Physical
Sciences and Astronomy
I certify that I have read this study and that in my
opinion it conforms to acceptable standards of scholarly
presentation and is fully adequate, in scope and quality,
as a dissertation for the degree of Doctor of Philosophy.
James W. Dufty 1P
Assistant Professor of Pnysics
This dissertation was submitted to the Graduate Faculty of
the Department of Physics and Astronomy in the College of
Arts and Sciences and to the Graduate Council, and was
accepted as partial fulfillment of the requirements for the
degree of Doctor of Philosophy.
August, 1974
Dean, Graduate School
