\part{Modelling of radio line emission from circumstellar molecules in
a sample of M-type AGB-variables}

\chapter{Observations and modelling of CO \& SiO circumstellar line emission}    

Radio observations of CO and SiO thermal line emission towards a
number of M-type AGB-variables have been used for studying the
characteristics of their stellar wind and the properties of their
circumstellar medium.

\section{The molecules under study}
\markboth{CHAPTER 4. OBSERVATIONS AND MODELLING}{4.1. THE OBSERVATIONS}


The first detected circumstellar molecule was maser emission from OH.
While maser lines owe their existence to rather specific causes,
thermally excited lines are more probable and more common.  The first
detection of such a line was the CO (v=0; $J$=1$\rightarrow$0) line at
115 GHz in IRC+10216 (Solomon et al. 1971).  Thermally
excited lines allow a more detailed exploration of the envelopes and
reveal information that cannot be derived from the maser lines.  But
these lines are not amplified and therefore they are more difficult to
detect.

\subsection{CO}

The two most important excitation processes for CO molecules in the
winds of AGB stars are the excitation of the rotational levels by
collisions with H$_{2}$ molecules and the photo-excitation of the
vibrational levels by infrared pumping.  The strongest thermal lines
are in general the rotational transitions in the ground vibrational
state of \co.  Because the molecule in its ground electronic and
vibrational state has a small permanent dipole moment, the
$J$=1$\rightarrow$0 transition of CO has a relatively modest Einstein
value and it takes only a relatively small critical density of H$_{2}$
to maintain collisionally a thermal population of this molecule in the
first rotational excited state.  There exist measurements of thermally
excited CO in more than 300 circumstellar envelopes of which about a
half are oxygen-rich.  The highest transition detected so far in
ground based telescopes is the $J$=7$\rightarrow$6 line in the carbon
star IRC+10216 but much higher transitions have been detected with KAO/ISO.

In the excitation analysis of the CO molecule we have used 40
rotational levels in each of the ground and first excited states. 
Excitation to the $v$$=$2 state has been ignored, since the $v$$=$1
state is not well populated.  We neglected collisional transitions
between the vibrational states because of the low densities and the
relatively low temperatures.

\subsection{SiO}

Circumstellar lines of SiO have been detected in vibrationally excited
lines but they are of maser nature.  SiO lines in the ground
vibrational level are, however, thermally excited.  Most observations
of thermal SiO emission have been made in the $J$=2$\rightarrow$1 and
$J$=1$\rightarrow$0 transitions, but transitions as high as
$J$=10$\rightarrow$9 have been detected in some O-rich stars.

It is commonly thought that the silicon in the extended atmospheres of
AGB-stars is incorporated into dust grains. SiO is expected to be a
major component of grain nucleation and growth in M-type AGB stars.
As a result, the SiO gas phase abundance should fall off with increasing
distance from the star as SiO molecules in the outflowing stellar wind
are incorporated into the grains, a quite uncertain process since it
does not proceed at thermal equilibrium. However, the nucleation rate
only peaks around 600K, which is significantly lower than the typical
temperatures at the inner edge of the dust shells.

In the radiative excitation analysis of SiO 50 rotational levels in
both the ground and first excited vibrational state are considered. 
The energy levels of this linear rotor are calculated using the
molecular constants from Mollaaghababa et al. (1991). 

The radiative rates are calculated using the dipole moment from
Raymonda et al. (1970).  Collisional deexcitation rates has been
calculated by Turner et al. (1992) in the temperature range 20--300\,K and
up to $J$$=$20.  The original data set has been extrapolated both in
temperature and to include levels up to $J$$=$50.

                   
\section{The sample}

The sample used in this work contains mainly semiregular (SRVs) and
irregular (IRVs) M-type AGB-variables. A few O-rich AGB Miras where
included for extending the study of SiO depletion onto dust grains to
objects with higher mass loss rates.

The original sample of semiregulars was used in a previous work for a
comparison of pulsational properties, effective temperatures, mass
loss rates, luminosities, scale heights and volumen densities in the
galactic plane with the corresponding characteristics of optical Mira
variables.  A total of 1442 SRa and SRb variables with reliable
periods were selected from the General Catalogue of Variable Stars.
The source positions were cross-correlated with the entries in the
IRAS-PSC in order to get positional associations.

A subsample of these SRVs containing 91 objects were selected for
observations of CO circumstellar radio lines on the basis of their
IRAS-60\micron~flux densities (typically S$_{60}$ $\geq$ 5 Jy).  Only
some of them ($\approx$ 40\%) had already been detected in these
lines.  The observations [presented in a previous work] yielded a total of
44 detections.

A total of 31 well studied O-rich IRVs were selected also on the basis
of their IRAS-60\micron~flux densities with the purpose of detect them
in CO circumstellar emission for the first time.  The
IRAS-60\micron~flux threshold for selection was, however, lowered to
about 2 Jy in this type of stars in order to get a comparable number of
IRVs and SRVs. A number of 22 IRVs were detected in CO radio lines. The
detections were also presented in a previous work.

Only stars redder than 1.2 mag in the IRAS [12]--[25] colour were
observed, thus biasing the sample towards stars with detectable
circumstellar dust envelopes.  There is a possibility that stars with
detectable gas mass loss rates, but with very little circumstellar
dust, were missed due to this.  Nevertheless, the detection rate was
rather high, about 60\%, and the conclusions drawn in this work should
be representative for the stars in the sample.

The CO sample used in the modelling work presented in {\bf Paper~III}
consisted of 69 M-type irregular (the 22 objects with CO detections)
and semiregular (the 44 objects detected in CO plus new detections in
3 other stars) AGB-variables. This sample is believed to be
representative for the mass losing stars of this type.

The modelling work presented in {\bf Paper~IV} was performed in a
subsample of the CO sample for which circumstellar SiO thermal radio
line emission was detected. A total of 60 stars were observed. The
detection rate was 60\%: 8 irregulars (detection rate 44\%) and 28
semiregulars (detection rate 68\%). In this paper, 12 M-type AGB Miras
were added to the sample in order to extend the study to objects with
higher mass loss rates. First, we modelled their CO circumstellar
radio line emission from which we derived the mass loss rates. Then,
we studied their SiO thermal emission together with the irregulars and
semiregulars.

\section{The observations}

\subsection{CO observations}

The CO ($J$=1$\rightarrow$0, 2$\rightarrow$1, 3$\rightarrow$2, and
4$\rightarrow$3) data which are used as the observational constraints
for the mass loss rate determinations were collected over the years
1993 to 1998.  They were obtained using the 20\,m telescope at Onsala
Space Observatory (OSO), Sweden, the 15\,m Swedish-ESO Submillimetre
Telescope (SEST) on La Silla, Chile, the IRAM 30\,m telescope on Pico
Veleta, Spain, and the James Clerk Maxwell Telescope on Mauna Kea,
Hawaii.  Three additional SRVs detected in CO were observed at OSO in
May 2000.

\subsection{SiO observations}

The SiO (2$\rightarrow$1) observations were obtained using the 20\,m
telescope at Onsala Space Observatory (OSO), Sweden, and the 15\,m
Swedish-ESO Submillimetre Telescope (SEST) on La Silla, Chile. At
SEST, a sizable fraction of the stars were observed also in the SiO
$J$=3$\rightarrow$2 line.  The higher-frequency lines,
$J$=5$\rightarrow$4 and $J$=6$\rightarrow$5, were observed towards
only 3 stars. The data were collected over the years 1992 to 2001.

\section{Model description}

\subsection{Model purpose}

The goal of modelling the CO circumstellar radio line emission was to
derive the mass loss rates for the semiregular and irregular M-type
AGB-variables in the sample. The SiO thermal radio line emission was
modelled to determine circumstellar SiO abundances, and thus, to study
the gas depletion onto dust grains.

The CO emission stems from the outer parts of the CSEs, where the gas
has already reached the terminal velocity. By fitting the CO line
profiles, it is possible to derive the gas terminal velocities for the
stars in the sample. However, the SiO emission comes from inner
regions of the CSEs, where very likely the gas is still accelerating
and dust condensation is taking place. Therefore, a comparison between
the gas expansion velocities derived from the CO and SiO line emission
was expected to be a direct probe of the CSEs dynamics.
 	
	
\subsection{The Monte Carlo code}

A non-LTE Monte Carlo code was used to model the circumstellar
emission and derived the characteristics of the circumstellar
envelopes.  An accurate treatment of the molecular excitation is
needed for a correct description of the radiative transfer in an
expanding circumstellar envelope.  A Sobolev approximation has been
shown, for example, to introduce significant errors.  The Monte Carlo
method has the advantage of being very flexible.  Thus, it is well
suited for solving multi-level radiative transfer problems in complex
geometries and velocity fields.  The main drawback is its slow
convergence ($\approx \sqrt{N_{iter}}$).

The modelling of the radiative transfer aims at obtaining the
steady-state level populations of the molecule under study, using the
statistical equilibrium equations.  The line photons are simulated by
a number of model photons, each of which represents a large quantity
of real photons assigned by a weight, {\it W}.  These model photons
are followed through the envelope which is divided in a number of
concentric shells with their own set of physical parameters.  The
number of absorption events caused by the radiation that is
represented by a model photon is calculated and stored.  As the model
photon travels in the envelope, its weight is continuously modified in
order to account for these absorptions and for the stimulated
emissions that also take place.  The total number of induced
de-excitations (u$\rightarrow$l) that the model photon gives rise to
when it travels a distance {\it s$_{1}$} in a shell {\it m} is

\begin{equation}
    N_{ul,m} = \frac{h\nu_{ul}}{4\pi} \phi_{\nu} n_{u,m} B_{ul} 
    \int_{0}^{s_{1}} 
    W(x)\,dx \, .
\end{equation}


The weight of the model photon will be modified according to 

\begin{equation}
    W(s_{1}) = W_{0}\,e^{-\tau_1}\, ,
\end{equation}
%
where $W_{0}$ is the weight at {\it x}=0, and the optical depth is given by

\begin{equation}
   \tau_1 = \frac{h\nu_{ul}}{4\pi} \phi_{\nu}
   \left(n_{l,m} B_{lu} - n_{u,m} B_{ul} \right) s_{1} \, .
\end{equation}


The number of radiative de-excitations per molecule in the upper 
state, {\it u}, in shell {\it m} is

\begin{equation}
    S_{ul,m} = \frac{N_{ul,m}}{n_{u,m}V_{m}} =
    \frac{h\nu_{ul}}{4\pi} \,\phi_{\nu}\,B_{ul}\,\frac{W_0\,
    s_{1}}{V_{m}\,\tau_{1}} \left(1 - e^{-\tau_1} \right) \, .
\end{equation}

When all photons have been emitted and followed through the cloud, 
new level populations are calculated from the equation of statistical 
equilibrium according to the Monte Carlo formalism. 

\begin{equation}
    n_{l,m} \left(\frac{g_u}{g_l} \, \sum S_{ul,m} + C_{lu} \right) =
    n_{u,m} \left(A_{ul} + \sum S_{ul,m} + C_{ul} \right) \, .
\end{equation}

The whole process is then repeated until some criterion for
convergence is fulfilled.


\subsection{The circumstellar model}

For modelling the observed molecular line emission, we have made used
of a standard circumstellar model. A density law consistent with an
isotropic constant mass loss and expansion velocity was assumed. It
can be derived from the conservation of mass,

\begin{equation}                                             
   \rho_{H_2} = n_{H_2}m_{H_2} = \frac{\dot M}{4\pi r^2v_e} \, .
\end{equation}

This is certainly an appropriate scenario for modelling the CO line
emission. However, there is a possibility that neither the mass loss
rate nor the expansion velocity are constant in the regions were the
SiO thermal emission is originated.

The envelope inner radius was taken as 3R$_\star$. A turbulent
velocity of 0.5 Km/s, which contributes to the Doppler broadening of
the local line widths, was considered. In principle, it is expected
that for low mass loss rate objects the choice of these two parameters
has an effect on the results. The model results, however, indicate
that this is not the case. A 50\% change of the turbulent velocity
affects the line intensities by less than 20\%. The dependence on the
envelope inner radius is of the same magnitude.
    
A radiation field consistent of a central blackbody that characterises
the star, and another simulating the background microwave radiation at
2.3 K which injects photons from outside, were considered. For the
sample of Mira variables we also included a second central blackbody
to represent the dust-processed radiation, since for these high mass
loss rate stars the excitation of the SiO molecules may be affected by
dust emission.

In order to determine the temperature structure in the CSEs, the energy
balance equation

\begin{equation}                                             
    \frac{dT}{dr} =
    (2-2\gamma)\frac{T}{r}+\frac{\gamma-1}{n_{H_2}kv_e}(H-C)\, ,
\end{equation}

%
was solved including a number of cooling and heating processes. This
equation was resolved in the modelling of the CO emission for each
CSEs, and the derived temperature structures were used in the
modelling of their SiO emission as well. This is reasonable since the
SiO line emission contributes very little to the cooling of the gas,
and in most cases the SiO molecules are mainly radiatively excited.

$\gamma$ is the adiabatic index and $k$ is the Boltzmann constant. The
first term on the right hand side is the cooling due to the adiabatic
expansion of the gas. Additional cooling is provided by molecular line
emission from CO and H$_2$. The molecular cooling due to CO was
calculated from the derived level populations using the expression of
Crosas \& Menten (1997). For the H2 cooling we used the approach by
Groenewegen (1994). 

Radiation pressure on the dust grains 
gives them a drift velocity, v$_{dr}$, relative to the gas

\begin{equation}                                             
    v_{dr} =
    \left(\frac{Lv_eQ}{\dot Mc}\right)^{1/2}\, ,
\end{equation}
%
where L is the luminosity, Q is the averaged momentum transfer
efficiency, and c is the speed of light. As a result of the dust-gas
drift, kinetic energy of the order of $\frac{1}{2}m_{H_2}v_{dr}^2$
is transfered to the gas each time a particle collides with a
dust grain. This was assumed to provide the dominating heating of the
gas (Goldreich \& Scoville 1976),
\begin{equation}                                             
   H_{dg} = (n_d\sigma_dv_{dr})\frac{1}{2}\rho_{H_2}v_{dr}^2\, ,
\end{equation}
%
where $n_d$ is the number density of dust grains, and $\sigma_d$ the
geometrical cross section of a dust grain. Introducing the dust-to-gas
mass ratio this equation can be written

\begin{equation}                                             
   H_{dg} = \frac{3}{8}m_{H_2}^2n_{H_2}^2\frac{\Psi}{(a_d\rho_d)} \frac{v_{dr}^3}{\left (1+\frac{v_{dr}}{v_e}\right )}\, ,
\end{equation}
%
where $a_d$ and $\rho_d$ are the average size and density of a dust
grain, respectively.

When solving the energy balance equation free parameters describing
the dust, i.e., via the dust-gas collisional heating, were
introduced. These are highly uncertain, but affect the derived line
intensities. In this work we assumed that the Q-parameter, i.e., the
efficiency of momentum transfer, to be equal to 0.03 (Habing et
al. 1994), and defined a new parameter that contains the other dust
parameters,

\begin{equation}                                             
   h = \left(\frac{\Psi}{0.01}\right)\left(\frac{2.0~g~cm^{-3}}{\rho_d}\right)\left(\frac{0.05~\mu m}{a_d}\right)\, ,
\end{equation}
%
where the normalized values are the ones we used to fit the CO line
emission of the well-studied C-rich AGB star CW Leo (IRC+10216) using
our model, i.e., h = 1 for this object.

Additional heating is provided by the photoelectric effect,

\begin{equation}                                             
   H_{pe} = k_{pe}n_{H_2}\, ,
\end{equation}
%
where $k_{pe}$ is a constant that depends on the properties of the dust
grains. Following Huggins et al. (1988) we adopted $k_{pe}$ = 1 10-26 erg
s-1. This heating mechanism is important in the cool, tenuous, outer
parts of CSEs around high mass loss rate stars.

The CO fractional abundance in the CSEs was assumed to be
$f_0$=2$\times$10$^{-4}$ (Kahane \& Jura 1994). This is essentially a
free parameter in the modelling, although its upper limit is given by
the abundance of C, i.e., 7$\times$10$^{-4}$ for solar abundances. The
SiO fractional abundance was determined through the modelling of the
circumstellar SiO thermal emission.

\subsection{Envelope sizes}

\subsubsection{The CO envelope size}
The spatial extent of the molecular envelope is generally an important
input parameter, and the derived mass loss rates depend on this. The
radial distribution of CO in the CSE was estimated using the modelling
presented in Mamon et al. (1988). It includes photodissociation,
taking into account the effects of dust-, self- and H$_2$-shielding, and
chemical exchange reactions.  In their model the radial abundance
distribution of CO with respect to H$_2$, f(r), can be described by

\begin{equation}                                             
   f(r) = f_0\,exp\left[-ln2\left(\frac{r}{r_p}\right)^\alpha\right] \, ,
\end{equation}
%
where $f_0$ is the initial (photospheric) abundance, $r_p$ is the
photodissociation radius (where the abundance has dropped to $f_0$/2),
and is a parameter describing the rate at which the abundance
declines. Both $r_p$ and $\alpha$ depend on the mass loss rate ($\dot
M$), the expansion velocity ($v_e$), and $f_0$. In our modelling we
used analytical fits to these results from Stanek et al. (1995)
\begin{equation}                                             
  r_p = 5.4\times10^{16}\left(\frac{\dot M}{10^{-6}}\right)^{0.65}\left(\frac{v_e}{15}\right)^{-0.55}\left(\frac{f_0}{8\times10^{-4}}\right)^{0.55}+7.5\times10^{15}\left(\frac{v_e}{15}\right)\,cm  \, ,
\end{equation}
% 
and Kwan \& Webster (1993) 
\begin{equation}                                             
   \alpha = 2.79\left(\frac{\dot M}{10^{-6}}\right)^{0.09}\left(\frac{15}{v_e}\right)^{0.09} \, ,
\end{equation}
%
where the units of $\dot M$ and $v_e$ are M$_\odot$yr$^{-1}$ and
kms$^{-1}$, respectively. We note here that a somewhat more
sophisticated photodissociation model was developed by Doty \& Leung
(1998) in which they include scattering by dust. For the CSE of CW
Leo they derive a CO envelope size which is 30\% smaller than what is
obtained using the Mamon et al. model. It is however not clear by how
much the envelope size changes for an object with a significantly
lower mass loss rate, and hence we have used the results of Mamon et
al.

\subsubsection{The SiO envelope size}

The size of the SiO envelope is not easily observationally determined
nor theoretically estimated.  Early work assumed that the whole CSE
contributes to the observed SiO thermal line emission (e.g, Morris et
al. 1979).  The mostly Gaussian single-dish SiO profiles found by
Bujarrabal et al. (1986; 1989) towards O-rich CSEs suggested that this
is not case. The generally small size of the SiO thermal line emitting
region requires interferometric observations in order to resolve it.

The model code used in this work allow us to derive SiO envelope sizes
provided that multi-line SiO observations are available.  The emission
from higher-$J$ lines (e.g., $J$=5$\rightarrow$4 and
$J$=6$\rightarrow$5) comes very likely from the warmer inner regions
of the SiO envelope.  Therefore, the intensities of these lines were
fitted by varying only the SiO abundance, i.e., they are insensitive
to the outer radius of the SiO envelope. Once the SiO abundance was
found, the lower-$J$ lines were used as constraints to derive the
size of the SiO envelopes, since their emission is photodissociation
limited (i.e., not excitation limited).

Only five sources in our sample were observed in higher-$J$ SiO
lines. The resulting $r_{\rm e}$:s obtained from the modelling were
combined with the interferometric observations of SiO thermal line
emission by Lucas et al. (1992) on the six stars in common with our
samples. It provided us with a relation of the mass loss rate
dependence of the SiO envelope size.  In this way, we derived envelope
sizes for the bulk of stars in the sample.

We determined the characteristic radii for the two processes that
compete for the decrease in circumstellar abundance of SiO molecules:
condensation onto dust grains and photodissociation by interstellar UV
photons. The former was estimated to play only an important role for
high mass loss rate stars.  It was reassuring that the calculated
photodissociation radii and the estimated SiO envelope sizes were
rather equal.

                                                                  
\subsection{The fitting strategy}

The radiative transfer analysis produces model brightness
distributions.  These are convolved with the appropriate beams to
allow a direct comparison with the observed line intensities and to
search for the best fit model.  In the modelling of the circumstellar
CO radio line emission there are two free parameters in the fitting
procedure, the mass loss rate $\dot{M}$ and the $h$-parameter (which
contains the dust properties). For stars observed in several CO lines
these two parameters were varied until the best-fit model was found.
In the modelling of the SiO circumstellar radio line emission the SiO
abundance and the SiO envelope size are allowed to vary in the five
sources with multi-line observations until the best fit model is
found.

For stars observed in only one or two CO transitions the h-parameter
was fixed in the CO modelling to a median value derived from those
objects with multi-transition observations. Likewise, for stars
observed in only one or two SiO transitions the SiO envelope size was
derived using an indirect method (see Sect. 5.1.2). In the modelling
of both molecules, the quality of a fit was quantified using the
chi-square statistic,
%
\begin{equation}
\label{chi2_eq}
\chi^2_{\rm red} = \frac{1}{\rm N-p}\,\sum^{\rm N}_{i=1} \frac{[I_{\mathrm{mod},i} -
I_{\mathrm{obs},i})]^2}{\sigma^2_{i}},
\end{equation}
%
where $I$ is the total integrated line intensity, $\sigma_i$ the
uncertainty in observation $i$, p the number of free parameters, and
the summation is done over all independent observations N. The errors
in the observed intensities are always larger than the calibration
uncertainty of $\sim$20\%.  We have chosen to adopt
$\sigma_i$\,=\,0.2$I_{\mathrm{obs},i}$ to put equal weight on all
lines, irrespective of the S/N-ratio.
