2 Modelling method

The method of modelling stellar H₂O maser emission at a single stellar phase of the pulsation cycle of a LPV star is analogous to that used in the stellar SiO maser model of Humphreys et al.(1996) (H96). The hydrodynamic stellar model employed here is the pulsating CE of a M-Mira variable. It is identical to that detailed in Sect. 3 of H96, except that in the present work we have extrapolated the model beyond its outer radial boundary of 53 AU to 134 AU using the radial dependency derived from the Bowen data. We therefore do not describe the stellar model in detail again here. The model star has the characteristics given in Table 2. We note that a stellar radius of 244 $~R_{\odot}$ (1.1 AU) is rather small when compared with values typically measured for M-Mira variables, at around 2 AU. However, our model star is similar in size and in ${\dot{M}}$ to measurements for R Cas (see Sect. 5.1). The physical conditions in the spherically symmetric CE at the epoch of our calculations are shown in Fig. 1. We outline the overall procedure for coupling the H₂O maser model to the CE model, and give details of the H₂O maser model, in Sects. 2.1 and 2.2 respectively.

 

 

Table 2: Characteristics of the model star used to compute the physical conditions in the CE of a typical M-Mira long-period variable.

Parameters of the model M-Mira variable
Mass	1 $M_{\odot}$
Fundamental period	332 days
Stellar radius	244 $R_{\odot }$ ( $1.7 \times 10^{11}$ m)
Effective temperature	3002.2 K
Maximum inner boundary speed	3.93 kms-1
Mass loss rate	$1.8 \times 10^{-7}~M_{\odot }$ yr-1

Figure 1: The physical conditions in the model CE at the epoch of our calculations. a) Gas velocity as a function of radial distance; b) kinetic temperature, $T_{\rm k}$; c) hydrogen number density, n(H2). 1 R* corresponds to 1.1 AU at this epoch.

   
2.1 Overview of the single epoch model

It is evident from the observational data at 22GHz that only certain regions at the same radius in the CE are suitable for the formation of bright masers. In order to mimic the clumpy nature of the CE, the cause of which is unknown, we must break the spherical symmetry of our CE model. At the epoch of the stellar cycle chosen for our calculations, we follow the procedure described below.

(i) Sites of potential maser action are distributed in the model CE by a standard Monte Carlo method. We require a distribution of sites which has a lower filling factor with increasing distance from the star, in order to imitate a scenario in which water maser clumps are produced near to the star and then travel out to occupy shells of larger volume. This is based on observational data in Richards et al. (1999), and is realised by a uniform distribution of sites over radial distance, r, a distribution of polar angle $\theta'$ which is related to a randomly chosen value of $\theta$ between 0 and $\pi$ of $\theta'=\cos^{-1} ( 1 - 2 \theta / \pi)$, and a uniform distribution of azimuthal angle, $\phi$. H₂O maser emission is calculated at these Monte Carlo sites only, at which the abundance of H₂O is assumed to be n(H₂O)/n(H ₂)=10^-4. On the basis of trial calculations, 3000 sites are used in the current model between the radial limits of 1-50 R_* (1.13-56.5 AU) from the centre of the star. 3000 sites were chosen in order to yield the number of bright components observed towards evolved stars. The outer radial boundary was chosen such that no inversions for any H₂O maser transitions were found to occur beyond this distance.

(ii) The physical conditions at each maser site are derived from the CE model. r defines the physical conditions at the site, whereas $\theta'$ and $\phi$ are necessary (i) to calculate the line-of-sight velocity gradient, $\alpha _{\rm los}$, at r(Eq. (4) in H96); and (ii) to define the position, in the plane of the sky, of any maser component which may develop.

 

 

Table 3: Ranges of physical conditions which yield the five brightest components at each frequency. r is the coordinate of the maser site; $\alpha _{\rm los}$ is the velocity gradient along the direction of maser propagation in the line-of-sight; $\partial V$/ $\partial R$ is the bulk radial velocity gradient at r; V / r is the tangential velocity gradient at r.

$\nu$	r	$T_{\rm k}$	n(H2)	$\mid{\alpha_{\rm los}}\mid$	$\mid\partial V$/ $\partial R\mid$	$\mid V / r \mid$
(GHz)	(cm)	(K)	(cm-3)	(kms-1 pc-1)	(kms-1 pc-1)	(kms-1 pc-1)
22	$\rm 3.5E13\rightarrow4.6E13$	$\rm 1819\rightarrow4213$	$\rm 9.0E8\rightarrow3.4E9$	$\rm 1.1E5\rightarrow7.1E5$	$\rm 9.0E5\rightarrow2.6E6$	$\rm 1.5E5\rightarrow2.5E5$
321	$\rm 3.5E13\rightarrow5.8E13$	$\rm 1819\rightarrow4007$	$\rm 1.9E8\rightarrow3.4E9$	$\rm 1.9E4\rightarrow1.2E5$	$\rm 5.6E5\rightarrow1.1E7$	$\rm 1.8E4\rightarrow2.8E5$
325	$\rm 3.7E13\rightarrow6.1E13$	$\rm 2070\rightarrow2552$	$\rm 1.7E8\rightarrow2.5E9$	$\rm 1.9E5\rightarrow8.8E5$	$\rm 4.4E5\rightarrow1.8E6$	$\rm 1.3E4\rightarrow2.0E5$
183	$\rm 1.7E14\rightarrow2.2E14$	$\rm 1343\rightarrow1583$	$\rm 1.1E7\rightarrow1.7E7$	$\rm 5.7E4\rightarrow7.4E4$	$\rm 6.5E3\rightarrow3.5E5$	$\rm 7.8E4\rightarrow9.3E4$

(iii) At each site, the population distribution of H₂O molecules over energy states is calculated using the Large Velocity Gradient (LVG) approximation (see Sect. 2.2.2). Where population inversions result for maser transitions, maser radiation is propagated in the line-of-sight.

(iv) The emission from sites which yield masers (we refer to emitting sites as "components'' hereafter) is combined to produce single-dish lineshapes and VLBI-type images. Each component produces a spectral emission line corresponding to each maser transition from which emission is calculated to occur. For any given maser transition, the spectral lines from all components are combined at their appropriate line-of-sight velocity to form single-dish lineshapes. In order to produce VLBI-type images, for any given maser transition, the spectral line from each component is averaged over all velocity channels and plotted at its projected coordinates. We emphasise that the success of any maser site for becoming a component depends on the physical conditions found at that site in the CE model and is not constrained to occur in any way. For example, less than half of the 3000 Monte Carlo sites distributed give rise to 22GHz emission in the line-of-sight (see Sect. 5.1).

   
2.2 The $\mathsfsl{H_{2}}$O maser model

The H₂O maser model consists of 100 rotational energy levels of the ground vibrational states of both ortho and para-H₂O. It is a 1-D radiation transport code in which the radiation field is treated classically and the molecular response is quantum mechanical. The effects of saturation and competitive gain are included in the propagation of maser radiation through a gas which contains population inversions, based on Field & Gray(1988).

2.2.1 Inputs to the model

The r coordinate of a maser site places it in a radial zone of the CE model. Linear interpolation between the boundary values of the zone provides values of the gas kinetic temperature ($T_{\rm k}$), density and bulk radial velocity at the site. The radial and line-of-sight velocity gradients at r can be calculated directly, using the $\theta'$and $\phi$ coordinates of the site, see Sect. 2.2.3. We assume the gas to be composed of molecular hydrogen, with number density n(H₂). Our value of the H₂O abundance (10^-4) is based on the models of stellar H₂O abundance by González-Alfonso & Cernicharo (1999). It is also consistent with the range of values of n(H₂O) found to result in population inversions in Yates et al. (1997; Y97 hereafter). We assume a thermal equilibrium ortho-to-para water abundance ratio of 3. The radiation field in the model is described in Sect. 2.2.2. The maser propagation distance is derived from the CE model, see Sect. 2.2.3.

   
2.2.2 Calculating the populations

Calculation of the H₂O level populations involves the solution of the kinetic master equations for 100 rotational energy levels in the ground vibrational state of both ortho and para-H₂O. Einstein A-values for the 422 dipole-allowed transitions involved are calculated using an algorithm described in Bayley (1985). Rate coefficients for inelastic collisions between H₂O and He involving the first 45 levels, up to level 7₇₀ for ortho and 7₇₁ for para-H₂O, are taken from Green et al. (1993). As we assume the collision partner to be molecular hydrogen rather than helium, the rate coefficients in Green et al. (1993) are scaled to account for the difference in reduced mass. For collisional rate coefficients for which either level lies above 7₇₀ and 7₇₁, the rate coefficients are estimated using the parametrisation adopted in Neufeld & Melnick (1991). The radiation field in the model consists of the stellar continuum, suitably spatially diluted according to the radial distance of each maser site from the star, and a dust radiation field. The stellar continuum radiation field is that of a black-body at the effective temperature of the model star. The dust local to the H₂O maser zone is assumed to be at its frequency-dependent radiative equilibrium temperature. The dust temperature is computed using a spherical radiation transport model, using dust opacity data from Laor & Draine (1993). The maser amplifies this background radiation field, which is composed of the dust and the stellar continuum.

To calculate self-consistent populations and line and continuum radiation fields, solutions of the kinetic master equations for the populations of the rotational energy levels of H₂O are obtained using the Sobolev or LVG approximation. We have employed this method, rather than use an accurate Accelerated Lambda Iteration (ALI) method, due to the time constraints imposed by the large number of calculations required to produce this model. In the present work, the LVG method has a further advantage that it provides a size for the maser region which is independent of other parts of the CE, the analogue of the observational clump size. We note that the use of the LVG approximation is rarely appropriate. It is a useful qualitative method, however, when the maser pump scheme is dominated by collisional excitations. The H₂O maser at 22GHz is well-known to be pumped by such a scheme (see e.g. Y97; Cooke & Elitzur 1985). However, if the pumping requires predominantly radiative excitations, the LVG approximation is no longer an effective method. The drawback here is that the maser lines at 437, 439 and 471GHz, which are found to have a predominantly radiative pump in Y97 using an ALI method, will be poorly represented by an LVG model and cannot be modelled here. Of the 100 levels included in the calculations for each of ortho- and para-H₂O, we exclude the results from the energy levels higher than 60. The truncation of energy levels may lead to inaccurate line optical depths for the higher levels.

 

 

Table 4: Data ranges which result for the five brightest components at each frequency. Fractional saturation distance = saturation distance/propagation distance.

$\nu$	Propagation	Fractional	FWHM
	Distance	Saturation
(GHz)	(m)	Distance	(kms-1)
22	$\rm 2.98E11\rightarrow1.72E12$	$0.33\rightarrow0.62$	$\rm0.98\rightarrow1.47$
321	$\rm 1.75E12\rightarrow3.00E12$	$\rm0.47\rightarrow0.90$	$\rm0.92\rightarrow2.67$
325	$\rm 2.42E11\rightarrow1.27E12$	$\rm0.27\rightarrow0.52$	$\rm0.94\rightarrow2.71$
183	$\rm 2.50E12\rightarrow3.00E12$	$\rm0.40\rightarrow0.47$	$\rm0.84\rightarrow1.00$

Figure 2: Velocity-intensity-distance surfaces at 22, 325, 321 and 183 GHz for emission arising from a single masing component. For clarity, intensity has been normalised to the peak value in the 22 GHz surface. This component was chosen as it is particularly bright at 325 GHz. The physical conditions at this site in the CE are: $T_{\rm k} = 2070$ K; hydrogen number density, n(H $_{2})= 2.5\times 10^{9}$ cm-3; line-of-sight velocity gradient, $\alpha _{\rm los}= 4.5\times 10^{5}$ kms-1 pc-1.

   
2.2.3 Maser propagation

Maser radiation is propagated in the line-of-sight for maser transitions in which population inversions are calculated. We identify the line-of-sight by calculating the velocity gradient in the line-of-sight at r, according to Eq. (4) of H96. The length of gas traversed during the amplification process is constrained to have a maximum value of $3\times 10^{12}$ m, a distance which has been empirically derived from resolved 22GHz maser features in RT Vir (Richards, private communication), or a value such that the velocity gradient in our line-of-sight causes a shift of 3 Doppler widths, which ever is the smaller. Field et al. (1994) showed that, in the presence of a velocity gradient along the path of maser propagation, significant amplification may take place for velocity shifts exceeding 1 Doppler width, i.e. for propagation distances greater than $\Delta v_{\rm th}$/$\mid$ $\alpha _{\rm los}$$\mid$, where $\Delta v_{\rm th}$ is the thermal linewidth and $\mid$ $\alpha _{\rm los}$$\mid$ is the velocity gradient. We note that simulations show that, when the velocity shift exceeds 3 Doppler widths, amplification typically becomes negligible.

 

 

Table 5: Break down of the results by temperature. Total summed maser output lies in the ratio 5:3:4:1 for 22, 321, 325 and 183GHz masers respectively.

Temperature Range (K)	<999	1000-1999	2000-2999	3000-3999	4000-4999	5000-5999	6000-6999
Frequency (GHz)	$\%$ Summed maser output at each frequency (number of contributing sites)
22	<0.001 (245)	24 (752)	48 (250)	18 (40)	8 (14)	<2 (20)	<0.02 (1)
321	0	16 (183)	46 (231)	<24 (40)	9 (14)	5 (20)	0.2 (1)
325	<0.2 (430)	36 (716)	54 (249)	9 (37)	<0.5 (12)	<0.6 (19)	<0.05 (1)
183	19 (1608)	65 (716)	13 (244)	2 (37)	<0.5 (11)	0.5 (19)	<0.03 (1)

The following method is used to treat the development of the component spectra over the propagation distance in the line-of-sight (see Eq. (1) of H96). At each numerical integration point in the propagation of masers through a masing zone, the molecular velocity distribution is divided into 101 bins covering 15 Doppler widths. The distribution is shifted appropriately in frequency at each integration step to take account of the local velocity field. Complete velocity redistribution (CVR) is assumed throughout the present calculations. This is achieved by summing the populations of all bins at each propagation step, taking account of saturation. The summed, possibly saturated, populations are then redistributed among the velocity bins according to a Gaussian profile. We note that the line centre of a component's spectral line does not correspond to the rest frequency of the maser transition in the presence of a velocity field. Microturbulence of velocity 1 kms^-1 is added in quadrature to the Doppler linewidth.