text:
# Numerical simulations of stellar SiO maser variability

## 1 Introduction

It is well-known from VLBI experiments that bright $`v=1`$ $`J=10`$ SiO maser components in M-Mira and Supergiant stars lie in approximate ring-type structures of several AU in diameter (in the range 1.5 – 4.0 R), which are assumed to be centered on the stellar position (Colomer et al. colomer92 (1992); Diamond et al. diamond94 (1994); Greenhill et al.greenhill95 (1995); Miyoshi et al. miyoshi (1994); Boboltz et al. boboltz (1997); Doeleman et al. doele (1998) at 86-GHz; Desmurs et al. desmurs (2000); Yi et al. yi (2001)). The width of the projected ring of clumpy emission is relatively narrow. The ring diameter changes with stellar pulsation phase, as emission brightens and dims, during the cycles of M-Mira stars such as TX Cam and R Aqr (Diamond & Kemball diamond99 (1999); Boboltz et al. boboltz (1997)). M-Miras and other classes of star on the thermally pulsing AGB (TP-AGB) are of interest due to their high mass loss rates (in the range 10-7 – 10-4 Myr-1) which replenish the interstellar medium with heavy elements and dust grains. SiO masers, which may now be studied at 43-GHz with an angular resolution of 200$`\mu `$as and a velocity resolution of 0.1 kms-1 using the VLBA, provide excellent probes of the dynamics and time-varying physical conditions of the extended atmosphere/inner CE of AGB stars. A detailed VLBA study of $`v=1`$ $`J=10`$ (43-GHz) masers in the Mira TX Cam provides data on the inner CE at 50 epochs over 1.25 stellar cycles (Diamond & Kemball diamond99 (1999)).

An important role in the mass loss mechanism of TP-AGB stars is thought to be played by a combination of stellar pulsation-driven shock waves and radiation pressure on dust. In this scenario, a succession of shocks greatly extends the Mira atmosphere, and material is accelerated, overcoming the gravity of the star with accompanying mass loss. Observational evidence for shock waves is provided by large amplitudes ($`\mathrm{\Delta }v`$ = 20 – 30 kms-1) in the velocity curves derived from photospheric CO lines. These indicate that an outwardly propagating shock wave is associated with the stellar pulsation (Hinkle et al. hinkle84 (1984); hinkle97 (1997)). However models are unable to reproduce the relatively high mass loss rates measured for TP-AGB stars without the inclusion of radiation pressure on dust, which must be formed close to the star (Bowen bowen88 (1988)). This is possible in the shocked model in which relatively dense gas is levitated to regions of lower temperature by the shocks, allowing dust condensation to take place within a few stellar radii from the photosphere. Radiation pressure from stellar photons then accelerates the grains, driving a slow, cool wind through frictional coupling with the circumstellar gas. Infrared interferometry measurements support the formation of dust shells with inner radii of 3 – 5 R towards a number of Miras and Supergiants (Danchi et al. danchi (1994)). Towards the Supergiant VX Sgr, coordinated SiO maser and infrared interferometry observations indicate that the $`v=1`$ $`J=10`$ (43-GHz) ring of bright SiO masers and the inner radius of the dust shell are well-separated, the masers lying at $``$3 R inside the dust region of the supergiant star (Greenhill et al. greenhill95 (1995)).

Much remains to be understood about the complex processes in these stars. For example, recent multi-wavelength observations of the S-type Mira, $`\chi `$ Cygni, show that the star is largest at minimum light, whereas current non-linear pulsation models predict the maximum physical size to occur close to maximum light (Young et al. young (2000)). With respect to the inner CE region, measurements of the 8-GHz “radio photospheres” of several stars indicate that shocks must have been significantly damped by the inner radius of the SiO maser zone and propagate outwards with velocities of $`<`$5 kms-1 (Reid & Menten reid (1997)). Yet proper motion measurements for the SiO masers in TX Cam show that some components are outflowing at velocities of $``$10 kms-1 (Diamond, private communication). In addition, magnetic fields in the SiO maser zone may be of the order of several gauss (Kemball & Diamond kemball (1997)), and will therefore play a crucial role in the gas dynamics of the inner CE, or may be only a few tens of milligauss (Wiebe & Watson wiebe (1998)). Finally, recent observations appear to have detected rotation in the SiO maser zone of some stars (e.g. Boboltz & Marvel boboltz00 (2000)).

As a first essay towards understanding some of these observations, we combine a model of a pulsating TP-AGB star, in which mass loss is driven by shocks and radiation pressure on dust, with an SiO maser model. Humphreys et al. (humphreys (1996); hereafter H96) showed how, at a single stellar phase, the key features of SiO maser emission are reproduced by such a model. In the present study, the crude representation of the stellar IR radiation field results in a decoupling of the effects on SiO masers of a varying IR radiation field and shocks in the inner CE. The aim of the present study is to perform simulations of the time evolution of SiO masers, both in intensity and in spatial distibution, throughout a stellar pulsation cycle. Although constant radiative pumping is present, we effectively determine whether the effect of shocks alone is sufficient to reproduce the observed features of SiO maser variability.

## 2 SiO maser temporal variability

Three types of variability may be identified. Firstly, there is variability from cycle to cycle, random in behaviour in the sense that SiO maser spectra are quite different in appearance from one cycle to the next. Secondly, there may be relatively orderly long-term variability within a cycle, in which maser flux passes through a maximum at an optical phase not far removed from maximum light, as described in Sect. 2.1. Thirdly there is rapid variability over a period of a day to a few tens of days, as described in Sect. 2.2. We are concerned mainly in the present paper with the second and third types of variability, which we refer to as long term and short term (or rapid) variability. We do not address variability from cycle to cycle in this work, except for an investigation of the stellar phase at which maser emission is severely disrupted (Sect. 4.3).

The long term variability of SiO masers has been investigated in monitoring programs towards many Mira variables, for example $`o`$ Ceti (Mira), R Aqr, TX Cam, U Her, R LMi, IK Tau (=NML Tau), R Leo, W Hya, $`\chi `$ Cygni, R Cas, U Ori, R Aql, R Cnc, X Hya and T Cep. The most complete sets of observations performed to date are those of Martínez et al. (martinez88 (1988); hereafter MBA88) and Alcolea et al. (alcolea99 (1999)) for $`v=1`$ 43-GHz masers and Nyman & Olofsson (nyman (1986); NO86 hereafter) for $`v=1`$ 86-GHz masers. These studies demonstrate long and short term variability. A statistical approach to the study of long term variability has been taken by papers involving Cho. The same maser transition was observed towards a large sample of stars, yielding statistical information on that maser as a function of stellar phase. Very short term variability on a timescale of 12 to 24 hours measured over a small part of a stellar cycle has been investigated in one detailed study of R Cas and R Leo (Pijpers et al. pijpers94 (1994)) at 43-GHz $`v=1`$ $`J=10`$. We consider first the observational data for long term variability.

### 2.1 Long term variability: correlations with optical phase of the host star

The following general characteristics can be derived from observations of Miras.

(i) The period is stable from cycle to cycle, but spectra vary greatly (NO86). Indeed, so do photospheric sizes, as measured by speckle interferometry (Bonneau et al. bonneau (1982)).

(ii) Optical, infrared and low frequency (43 and 86-GHz) SiO maser fluxes are correlated. There is firm evidence for an average phase lag of $``$0.2 of a period between optical and maser maxima, but this may vary between objects, and for the same object, between cycles (MBA88, Cho et al. cho96b (1996); Alcolea et al. alcolea99 (1999)). From the 43-GHz data in MBA88 the phase lag had a mean of 0.18, standard deviation 0.1.

(iii) SiO maser emission appears to vary in phase with near and mid-IR lightcurves. Alcolea et al. (alcolea99 (1999)) show that the SiO maser maxima are closely related to the light curve at 3.79 and 4.64 $`\mu `$m for R Aqr and IK Tau. However, the amplitudes of the IR and maser lightcurves are not correlated. A good correlation of the maxima of 86 GHz masers and the IR continuum at 1.04 $`\mu `$m was found in 6 out of 8 miras studied in NO86, but failed for two objects, R Leo and R Aql.

(iv) SiO masers often have strong linear polarization (Barvainis & Predmore barvainis (1985)) which may introduce bias into observations made with a single linearly polarized feed. NO86 and MBA88 argue the bias is small however. Papers involving Clark are the best in this respect, since they measure the Stokes parameters I, Q and U.

(v) The contrast between the highest and lowest peak lineshape intensity (and lineshape area) during a cycle is very variable from cycle to cycle and star to star. For the Mira-type stars studied at 43 GHz in Alcolea et al. (alcolea99 (1999)), this value varied between 1 up to as high as 20, in the case of $`o`$ Ceti. However, at 86 GHz the dynamic range during one cycle of $`o`$ Ceti was at least 100, and could have been considerably larger, since masers in this source fell below the 3$`\sigma `$ detection limit of 5 Jy near minimum.

(vi) The velocity extent of emission is typically $``$15 kms-1, with the cycle-averaged peaks of the spectra close to the stellar velocity (NO86). Cho et al.(cho96b (1996)) find the mean velocity of their sample was red-shifted with respect to V for optical phase 0.3 to 0.8 with blue-shifted emission appearing from $``$0.85 and dominating between 0.0 and 0.2. The cycle-averaged spectral peak at $`v=1`$, $`J=10`$ was 0.3 kms-1 to the red of V, compared to the blueshift of $``$1 kms-1 found by NO86. In addition, SiO profiles sometimes have unusually broad, low-intensity linewings which can exceed the expansion velocity of the circumstellar envelope as traced by thermal CO measurements, at least at 86-GHz (Herpin et al.herpin98 (1998)).

### 2.2 Short term variability

Balister et al.(balister (1977)) noted that a number of SiO maser sources (in $`v=1`$, $`J=10`$ at 43-GHz) showed variability over a period of a few days. Spectra are shown of the semiregular variable R Dor and the supergiant AH Sco. Both stars showed a marked increase in maser flux, especially AH Sco, over a period of only 4 days. This remarkable phenomenon was not however studied in any detail until the work of Pijpers et al. (pijpers94 (1994)) who observed R Cas and R Leo over phases 0.96 to 0.09 and 0.14 to 0.23 respectively. With respect to total flux, R Cas showed a slow decline and R Leo went through a marked maximum over this period of 40 days. In both cases these variations probably represented correlation with the stellar phase, as discussed in Sect. 2.1, rather than short term variability. Line profiles however changed within a period of days. Channel by channel variations ranged between 10% and 30% of the peak flux, with typical timescales of variation of 10 to 20 days. The velocity centroid of the dominant maser peaks apparently shifted by $``$1 kms-1 on the same timescale. No clear periodicity was observed in these short term variations.

### 2.3 Survival of masers during a stellar cycle

Clark et al. (clark84a (1984)) and Clark et al. (clark85 (1985)) decomposed their 86-GHz maser spectra into a number of overlapping gaussian components. They find for R Leo and W Hya that individual gaussian components in the Stokes I spectra often persist in successive spectra over the greater part of an optical cycle. Miller et al. (miller84 (1984)) reach the same conclusion for $`o`$-Ceti. The consistency of polarization position angle in these individual components is an important factor in reaching this conclusion. Given that the centre velocity, velocity width, polarization position angle and fractional polarization of individual features persist throughout much of the optical cycle, Clark et al. (clark85 (1985)) argue that individual maser zones are stable and discrete physical entities. The continuity of maser features was reported to be broken at, or just prior to, maximum light for W Hya (Clark et al. clark85 (1985)), for R Cas (Clark et al. 1982a ), for R Leo (Clark et al. 1982b ; Clark et al. clark84a (1984)) and for $`o`$-Ceti (Miller et al. miller84 (1984)). NO86 however found no such result for R Cas, R Leo and $`o`$-Ceti, but note that a lack of polarization information (and line blending) could have caused the effect to be missed. According to MBA88 and Alcolea et al. (alcolea99 (1999)) however, continuity appears to be broken where SiO maser brightness is at minimum, at around an optical phase of 0.7, rather than near optical maximum.

### 2.4 Proper motions during a stellar cycle

VLBI studies have been made towards several M-Miras: R Aqr (Boboltz et al. boboltz (1997); Hollis et al. hollis (2001) for $`v=1`$ 43-GHz masers); TX Cam (Diamond & Kemball diamond99 (1999) for $`v=1`$ 43-GHz masers; Desmurs et al. desmurs (2000) for $`v=1`$ and $`v=2`$ 43-GHz masers; Yi et al. yi (2001) for $`v=1`$ and $`v=2`$ 43-GHz masers); $`o`$ Ceti (Gardner et al. gardner (2000) for $`v=1`$ and $`v=2`$ 43-GHz masers); R Cas (Yi et al. yi (2001) for $`v=1`$ and $`v=2`$ 43-GHz masers). Boboltz et al.(boboltz (1997)) made the first proper motion measurements of circumstellar SiO masers towards R Aqr. They observed an infalling mean component proper motion of 1 mas in 98 days, corresponding to 0.22 AU over 0.25 stellar periods or 4 kms-1, for their assumed distance to R Aqr of 220 pc. In the long-term monitoring of TX Cam, more complex motions are evident. Maser components are generally expanding outwards from the star at a typical velocity of 3.65 kms-1 (Diamond & Kemball (diamond99 (1999)), although some have constant outflow velocities as high as $``$10 kms-1. In one quadrant of the images, components decelerate smoothly as they appear to run into denser material. Some of the ring also shows features sliding in non-radial directions. Masers in one region of the ring fall inwards towards the star and then appear to “bounce back”, an effect which could be due to a new shockwave moving through the maser zone (Diamond, private communication). Disruption to the maser ring appears to occur at maser minimum light at around optical phase 0.67 (Diamond, private communication). A new maser emission ring, of smaller angular extent than the disrupted ring, is then observed to appear.

## 3 The model of SiO masers in the circumstellar environment

The technical means employed in this paper have been described in some detail in H96. A hydrodynamic pulsation model for M-type Miras, based upon that in Bowen (bowen88 (1988), bowen89 (1989)), is coupled to a model for SiO maser action based on that described in Doel et al. (doel95 (1995); D95 hereafter) and Gray et al. (1995). H96 considered only a single phase of the stellar pulsation and showed how, at this phase, the physical conditions in the CE were such that a ring of SiO masers should form about the star at a distance from the photosphere of $``$1 R (1 R = 1.1 AU). The synthetic map produced showed an encouraging resemblance to data acquired using the VLBA (Diamond et al. diamond94 (1994) and new data presented in H96). The lineshapes of masers in various transitions in v=1 and v=2 ranging from $`J=10`$ to $`J=76`$ were also presented and were found in form and width to resemble lines observed in the very large number of sources in the literature. H96 allowed certain conclusions to be drawn with respect to the location of SiO masers, that is, that they do not in general reside in the stellar wind, defining the stellar wind as that part of the flow in which the radial velocity of the CE is always directed away from the photosphere with no further episode of infall. The agreement between the calculated and observed position of SiO masers lends support to the shock driven pulsation model, since it is this model which determines spatial variation of number density, temperature and velocity field which in turn dictates where the SiO masers should be found.

The present work is an extension of H96 involving calculations at 20 phases of stellar pulsation, carrying the evolution of the CE through a complete stellar cycle. As the cycle develops, the number densities, kinetic temperatures and bulk velocity fields change as a shock propagates through the photosphere and into the CE. This in turn causes populations of rovibrational levels of SiO to be modified. Inversions between rotational levels therefore increase, decrease or disappear and maser emission is modified accordingly. The remainder of this section is devoted to a description of the hydrodynamic pulsation and the kinetic and radiative SiO maser models.

### 3.1 The hydrodynamic pulsation model

The model of the time-varying CE is based upon that described in Bowen (bowen88 (1988),bowen89 (1989)) and is identical to that described in H96. The models of Bowen have been chosen over more recent models, which include time-dependent dust formation, such as the work of Bessell, Höfner and Fleischer, as they have tended to concentrate on carbon-rich variables (Fleischer et al. fleischer91 (1991), fleischer92 (1992); Bessell et al.bessell (1996); Höfner et al.hofner96 (1996); Höfner hofner99 (1999)), which clearly do not support SiO masers since the greater part of the oxygen budget is bound up in CO. The models by Bowen represent oxygen-rich Miras, which comprise the vast majority of the AGB star population, although dust in these models is introduced as a parametrised opacity.

Briefly, the stellar oscillation is driven by a sinusoidally varying force situated below the photosphere of the star. Radiation pressure is assumed to act on dust, and on H2O molecules. The quantity of dust present depends on the grain temperature, with no dust for Trad = 2000 K and a maximum achieved for Trad $`<`$ 1000 K. The outflow of material from the star becomes a steady wind at 40 to 50 R, the outer boundary of the model. The same stellar parameters are used as in H96 and for ease of reference these are shown again in Table 1. The stellar phase in the model (the ‘model phase’) is defined as zero when the sinusoidally varying inner boundary is moving outwards with maximum velocity. We refer to the problem of relating the model phase to stellar optical phase in Sect. 4.

In the CE, a parcel of gas for example close to the photosphere, initially experiences a brief but very strong outward acceleration, and thereafter decelerates, reversing direction. After a period of infall it is then driven outwards again and this cycle is repeated with successively weaker impulses as the material finds itself increasingly further from the photosphere. Thus if one were arbitrarily to choose some initial radial position within the CE, and ride on that parcel of gas in the co-moving Lagrangian frame, then one would experience a continuously changing set of number densities, temperatures and velocity fields. The nature of these changes would depend on the initial choice of radius. This is illustrated in Fig. 1 a,b and c and in Fig. 2. Fig. 1a shows the variation in velocity at four different model phases, Fig. 1b the variation in temperature and Fig. 1c the variation in number density. The model phases chosen correspond to zero, 83 days, 166 days and 249 days, the stellar pulsational period being 332.0 days. Values are plotted out to 10 AU only, in order to show the zone of interest for SiO masers. Any vertical line drawn through the four graphs in each figure illustrates how the chosen parameter would vary with stellar phase. Fig. 2 illustrates the density and temperature history, in a Lagrangian frame, for four masing parcels of gas during a stellar cycle i.e. Fig. 1 gives an Eulerian view, and Fig. 2 the Lagrangian view following a particle. Points in Fig. 2 are separated by 16.6 days. Similar phase plots may be made for velocity-temperature and velocity-density.

### 3.2 Overview of the computational procedure

The calculation of maser spectra and of maps proceeds as follows. Sites of maser action in the CE are first chosen on a random basis, for one epoch of the calculations only. At present, it is assumed for simplicity that maser zones maintain their integrity throughout the cycle involving model phase zero to 1. The trajectories of these initially chosen sites are followed in terms of their radial distances from the star as the cycle progresses, using data from the stellar pulsation model. At any chosen phase, for each radial distance a number density, kinetic temperature and velocity field may be assigned, as is apparent from Fig. 1a,b,c and Fig. 2. For each such distance, SiO rovibrational level populations are calculated. If these populations lead to inversions, the intensities of the resulting masers are calculated for paths in our line-of-sight. In this way, maps of SiO maser emission may be calculated for as many phases of the stellar cycle as desired and for any masing transition. The series of steps involved in synthesising maser images is described in more detail in the succeeding sections. Maser spectra are generated at any phase by summing over the contributions of individual maser spots, with line-of-sight velocity shifts appropriate for their position.

### 3.3 Choice of masing zones: proper motions in the stellar cycle

VLBI images of SiO emission at 43-GHz and 86-GHz clearly show that at some chosen radial distance from the photosphere only certain portions of gas in the CE are capable of sustaining maser action. This suggests that there exists small-scale, non-radial velocity structure superimposed on the large-scale velocity structure calculated by the stellar pulsation model and/or inhomogeneity in the density and temperature structure of the CE, such as clump formation via thermal instabilities (Cuntz & Muchmore cuntz (1994)). Since the nature of such small-scale structure is unknown, sites of maser action were chosen randomly in H96. This turned out to be a successful means to reproduce the appearance of VLBI maps (of TX Cam) and we adopt the same procedure here. Thus, as in H96, we choose sites of maser action in the CE at random radial distances from the star within a specified range of radius, $`r`$, and in fact we start with the same set as in H96. Only directions of maser propagation in our line-of-sight are considered (see Sect. 3.5). The range of radius over which values are randomly chosen extends from an inner boundary rin = R = 1.7x1011 m to an outer boundary rout = 5 R. Tests with larger values of rout showed that this range of radius was sufficient to include all masers in all transitions. The number of random values of $`r`$ used is a free parameter of the model. We found in H96 that 1500 values, of which $``$25$`\%`$ yielded inversions, was appropriate to reproduce the appearance of VLBI maps of TX Cam in terms of the number of bright maser spots around the host star. The same value was chosen in the present study. Maser propagation is discussed in more detail in Sect. 3.5.

The CE is seen in projection in the plane of the sky and the location of masing zones within the spherical CE must be fully specified in order to reproduce VLBI maps. We use spherical polar coordinates r, $`\theta `$ and $`\varphi `$, choosing random points in such a way as to provide uniform sampling by volume. Thus a set of random r, $`\theta `$ and $`\varphi `$ is chosen and weighted as described in detail in H96, giving r’, $`\theta `$’ and $`\varphi `$’, where primed quantities refer to weighted values. In summary, r’ defines the physical conditions and $`\theta `$’ and $`\varphi `$’ are used to locate the position, in the plane of the sky, of any maser spot that may develop. The next stage is to advance the value of radius r’ to that corresponding to the next phase for which the physical conditions and maser location are required. The distance $`\delta `$r’ that any mass of gas of velocity V travels in time $`t`$ is given by

$$\delta \mathrm{r}^{}=_{\mathrm{t}_0}^{\mathrm{t}_0+\delta \mathrm{t}}\mathrm{dt}\mathrm{V}(\mathrm{r}^{},\mathrm{t}),$$
(1)
which is computed using $`\delta r^{}=\delta tV(t=t_0)+\delta t^2/2dV/dt+\delta t^3/3!d^2V/dt^2`$ etc where the differentials are evaluated at $`t=t_0`$. The total differentials $`dV/dt`$ and $`d^2V/dt^2`$ are obtained from the output of the stellar model as follows. The pulsation model describes how the velocity at a given r’ changes as phase advances and how for any one phase the velocity changes with radial distance in the neighbourhood of the radial position r’. These data may then be used to obtain $`dV/dr^{}`$ and $`d^2V/dr^2`$ through an expansion in partial differentials, where differentials up to third order have been included. Variations of r’ with phase trace out the proper motions of potential maser sites in the course of the stellar cycle. The CE model generates 520 values of number density, temperature and velocity over a range of radial distance 1.5 x 1011 m to 7.3 x 1012 m from the centre of the star. The physical conditions associated with any position are obtained through linear interpolation. Conditions are generated for each of the 1500 values of r’, $`\theta `$’ and $`\varphi `$’ for each of the 20 phases studied, making in all 30,000 sets of conditions. A certain proportion of the 1500 sets of conditions at any phase yield population inversions, and excluding the very small proportion which lie behind the disk of the host star, one may trace the proper motions of brighter persistent masers in the CE, as described in Sect. 4.3. In this connection, SiO rovibrational populations assume values corresponding to the prevalent physical conditions on a timescale which is very rapid compared with the timescale of variation of the physical conditions. Collisional and radiative events determine the SiO populations, the rate-determining steps being collisional. A typical rate coefficient for the transfer of rovibrational energy through collisions of H2 with SiO may be of the order of 10-12 cm3s-1 or greater. The number density lies between 108 and 1010 cm-3 in a typical maser zone, and thus the timescale for pumping maser inversions does not exceed a few x 104 seconds and is typically rather less than 104 seconds. Thus the SiO populations effectively respond instantaneously to changes in the environment.

### 3.4 Calculation of rovibrational level populations

The methods used for the calculation of the populations of SiO rovibrational levels have been described in detail in D95 with modifications described in H96. Populations of levels are initially calculated ignoring the presence of inversion in the medium or the influence of maser radiation. In this connection, population inversions involve rotational transitions whose line emissivities are negligible by comparison with those of rovibrational transitions (H96; Yates et al. yates97 (1997)). The effects of the presence of masers are discussed in Sect. 3.5. A discussion of the work of other groups involved in modelling SiO maser emission (e.g. Langer & Watson langer (1984); Lockett & Elitzur lockett (1992); Bujarrabal 1994a , 1994b ) is given in D95. Calculation of the SiO populations involves solution of the master equations for 200 rovibrational energy levels of SiO, involving $`v`$ = 0 – 4 and $`J`$ = 0 – 39 in each vibrational state. Stellar continuum radiation is included with a geometrical dilution factor of 0.1 (D95). As the model has a fixed photospheric temperature, with no phase variation (see Table 1.), and is represented as a black-body, the radiative part of the pumping scheme is represented only crudely. We note that the present hydrodynamic model provides no means of tracing either the true radius or temperature of the photosphere in the near to mid-IR wavelength range. Consequently, the present combined model is not capable of testing the observed phase correlation between the IR continuum and low-frequency SiO maser peaks. The current work is therefore a useful test, in that it provides information about those aspects of maser time variability which can and cannot be reproduced by a pump in which collisions provide that part of the pumping scheme which is dependent on stellar phase. An infrared dust radiation field is also present. The dust continuum is represented in a crude manner by a black-body function at a single dust temperature of 200 K modified by a wavelength dependent factor, $`(\lambda /\lambda _0)^p`$, where $`\lambda _0`$ is 80$`\mu `$m and $`p=1.1`$ (Rowan-Robinson et al. 1986), and the shortest value of $`\lambda `$ is $``$ 8 $`\mu `$m. D95 shows that the only significant property of the dust radiation field (or any other external field) is that it should be weak. Thus our very simple and inaccurate representation of the dust field is not a significant source of error in computation of the SiO energy level populations except for high $`v`$-states and minor isotopomers. The dust radiation field is spatially diluted through a factor of 0.01, since it is assumed that significant quantities of dust are not found in regions occupied by SiO masers (Greenhill et al. greenhill95 (1995)). The dust may also be patchy, with a covering factor significantly $`<`$1 (D95). Collisional rate coefficients for the transfer of rotational and rovibrational energy between SiO and its collision partners (see below) have been taken from Bieniek & Green (1983a ,1983b ; hereafter BG83) where values have been extrapolated, as described in D95 and Doel (doel90 (1990)), to include all $`v`$ and $`J`$ states in the maser model. We note that Lockett & Elitzur (lockett (1992)) also employed the rate data of BG83, and extended the range of vibrational transitions treated up to $`\mathrm{\Delta }`$v = 4. Langer & Watson (langer (1984)) used the BG83 data for collisions of molecular hydrogen with SiO, and additionally estimated rates for collisions of SiO with atomic hydrogen. The Sobolev or large velocity gradient (LVG) approximation is used to calculate self-consistent populations and line and continuum radiation fields, as in Lockett & Elitzur (lockett (1992)). In support of the use of this approximation, CE models clearly indicate the presence of supersonic velocity shifts over distances smaller than or comparable to the dimensions of the maser zone. However the physical conditions in the CE close to the photosphere change markedly over a typical Sobolev length, whilst in the simple form of LVG used here conditions are assumed constant over a Sobolev length (= Doppler width divided by the local velocity gradient). The LVG approximation evidently introduces a coarse-grained interpretation of observational data. As described in H96, the number density may change by as much as an order of magnitude over the Sobolev length, representing a rather severe spatial averaging. As in H96 we use a general expression for the angle-averaged photon escape probability in a spherically symmetric system, using values of the radial velocity gradient $`dv/dr^{}`$ and the tangential velocity gradient $`v/r^{}`$ taken from the stellar pulsation model for each value of radius. As we discuss in H96 we ignore the influence of non-monotonic velocity gradients, which may cause a ray of light to encounter additional surfaces in the medium at equal velocity to an emitting point (Rybicki & Hummer rybicki (1978)). D95 concluded that masers are pumped by collisions at typically 1500 K and a number density of $``$5 x 109 cm-3, in an environment containing large velocity gradients, perhaps in excess of 105 kms-1pc-1, and experiencing only a weak dust continuum radiation field. D95 pointed out that these conditions were consistent with those predicted by CE models (e.g. Willson willson (1987); Bowen bowen88 (1988); Bowen bowen89 (1989); Bowen & Willson bowen91 (1991); Fleischer et al. fleischer92 (1992)) for a zone lying within about a stellar radius from the photosphere.

The chief sources of error in the calculation of the populations of SiO energy levels are:

(i) the use of inaccurate rate coefficients for energy transferring collisions between SiO and collision partners. Calculated values of BG83 are for He. Collision partners in the CE are H2 and H, in unknown proportions, with a relatively small contribution from He.

(ii) use of the LVG approximation. Exact methods, omitting velocity gradients, have been used in Bujarrabal (1994a ,1994b ) for SiO masers. Exact methods involving Accelerated Lambda Iteration have been developed for treating radiative transfer in molecular systems (Jones et al. jones (1994); Randell et al. randell (1995); Yates et al. yates97 (1997)) and these could be used to advantage here (see Sect. 5). However, they would be inconsistent with the symmetry-breaking assumption used to simulate clumping in the present work.

(iii) a lack of any independent chemical model to calculate the SiO abundance. The ratio of SiO is fixed at 10-4 of the total particle number density, where the number density is represented in terms of H2 molecules. This value was chosen on the basis of the estimation given by Doel et al. (doel95 (1995)), in which compositions for circumstellar gas and dust were assumed and then a gas to dust mass ratio argument was used to yield an upper limit to the SiO abundance. We note that this value is somewhat higher than the generally accepted value of $``$5 10-5 for the inner CE, see for example Bujarrabal et al. (burra89 (1989)). The effect of different SiO abundances upon SiO maser emission was investigated by Doel et al. (doel95 (1995)). In the range 5 10-5 – 2 10-4, the unsaturated maser gain coefficients were found to be essentially proportional to n(SiO). In the absence of any chemical model, this abundance remains unchanged throughout these calculations for all conditions encountered in the SiO maser zones for all phases of the star, excepting when kinetic temperatures exceed 5700 K. On the basis of thermodynamic calculations, which indicate that SiO will be largely dissociated above this temperature in the inner CE (Field, private communication), we include a crude cut-off assuming negligible abundance of SiO at those sites for which Tk $`>`$ 5700 K.

### 3.5 Maser propagation in the circumstellar envelope

There is good evidence that bright SiO maser spots are strongly saturated in M-Miras, as discussed in detail in D95. Populations are therefore substantially modified in the volume of gas occupied by bright masers. Observations indicate that masers fill only a very small proportion of the total volume of the maser zone (e.g. Greenhill et al. greenhill95 (1995)) and the assumption is made here, as in H96, that saturating masers are sufficiently spatially confined that they do not affect the general pumping cycle elsewhere in the maser zone. Hence the omission of masers from the calculation of level populations described in Sect. 3.4. Maser propagation occurs through exponential growth followed by a region involving saturation, if rays achieve sufficient intensity. Maser polarization (Sect. 2.1 and 2.3; McIntosh et al. mcintosh94 (1994); Nedoluha & Watson ned90 (1990), ned94 (1994)) is not included in our calculations. Maser saturation and coupling of the maser radiation to the kinetic master equations are treated with a semi-classical formalism, in which the radiation is treated according to Maxwell’s equations but the response of the molecular ensemble is calculated using quantum mechanical density matrix theory (Field & Richardson field84 (1984); Field field85 (1985); Field & Gray field88 (1988)). All effects of saturation and competitive gain are included in the propagation of a maser beam through the gas containing population inversions. The masing zone is assumed to amplify a black-body background at the appropriate wavelength at the local kinetic temperature of the maser zone. Our calculations amplify this background as a function of frequency within the maser line.

Propagation of maser rays is performed by numerical integration of the set of coupled equations

$$\frac{\mathrm{dI}_{\mathrm{ji}}^\nu }{\mathrm{dz}}=\mathrm{A}_{\mathrm{ji}}\frac{\lambda _{\mathrm{ji}}^2}{8\pi }(\rho _\mathrm{j}\rho _\mathrm{i})\mathrm{I}_{\mathrm{ji}}^\nu $$
(2)
where $`I_{ji}^\nu `$ is a specific intensity function. Maser rays are assumed to be contained in a vanishingly small solid angle. Thus, in the standard expression relating the angle averaged intensity to specific intensity, the function is a $`\delta `$-function in the angles. The upper and lower level populations per sublevel, $`j`$ and $`i`$ respectively in Eqn. 2, are those which respond at a certain frequency, $`\nu `$, within the inhomogeneously broadened line. Populations as a function of frequency in the presence of saturating maser radiation are calculated using the expression

$$\rho _\mathrm{p}=\rho _{0\mathrm{p}}\mathrm{exp}\left\{\mathrm{\Sigma }_\mathrm{r}\left\{\mathrm{T}_\mathrm{r}^\mathrm{p}(\mathrm{L},\mathrm{S})\mathrm{I}_\nu \right\}\right\}$$
(3)
originally derived in Field & Gray (field88 (1988)) and used extensively elsewhere, for example in D95. In Eqn. 3, the $`r`$ subscript refers to a transition and the $`p`$ subscript or superscript to an energy level. $`T`$ is a maser intensity independent function of $`L`$ and $`S`$, which in turn are functions of all the non-maser radiative and kinetic events built into the inversion pumping scheme. The possibility that saturating maser beams may intersect one another is ignored. Counter-propagating beams (or ‘streams’) are not included in our calculations (Elitzur elitzur (1993)). The length of material traversed by the maser, the ‘gain length’, is constrained to have a maximum value equal to the stellar diameter of 3.4 x 1011 m or a value such that the velocity gradient in our line-of-sight causes a Doppler shift of three Doppler widths, whichever is the smaller. When the velocity shift exceeds three Doppler widths, amplification almost invariably becomes negligible. In integrating the coupled Eqns. 2 and 3, the line-of-sight is identified by calculating the velocity gradient in the line-of-sight at $`r^{}`$ and using this calculated value for maser propagation. The velocity gradient in the line-of-sight, $`\alpha _{los}`$, may be shown to be related to the radial and tangential velocity gradients by

$$\alpha _{\mathrm{los}}=\left(\mathrm{dV}/\mathrm{dr}^{}\mathrm{V}/\mathrm{r}^{}\right)\mathrm{sin}^2\theta ^{}\mathrm{cos}^2\varphi ^{}+\mathrm{V}/\mathrm{r}^{}$$
(4)
Propagation of SiO masers in a velocity field may readily be performed since, in Eqn. 2, maser propagation is treated as a function of frequency within the gain profile. At each numerical integration point in the propagation of masers through a masing zone, the molecular velocity distribution is divided into 61 bins covering 12 Doppler widths. The distribution is shifted appropriately in frequency at each integration step to take account of the local velocity field. Complete velocity redistribution 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 populations are then redistributed among the velocity bins according to a Gaussian profile, as described in Field et al. (field94 (1994)).

## 4 Results of calculations

Data are presented in a manner as far as possible similar to that of observational work described in Sect. 2. The appearance of synthetic SiO maser lineshapes and images at 43 and 86-GHz in the $`v=1`$ state is discussed, and the predicted proper motions of bright maser components are described. Unfortunately, there is uncertainty in establishing the relation between the model phase of our simulated data with stellar optical phase. This comes about as the exact optical phase at which a pulsation-driven shock front emerges from the stellar photosphere (model phase 0.0) is unknown. In order to compare the results of our simulations with variability observations, we calibrate our results against the data for TX Cam in Sect. 4.1, point (v).

### 4.1 Spectra of SiO masers at 43 and 86-GHz: long term variability

We now consider the points in Sect. 2.1 relating to long term variability and discuss how observed and calculated features correspond. Fig. 3 shows, using the $`v=1`$ $`J=10`$ (43-GHz) maser as an example, how time series of maser spectra can be output from the variability simulations for transitions in the range $`v=1`$$`3`$, $`J=10`$$`J=109`$. These data may be combined to produce cycle-averaged spectra, as used by NO86, and these are shown in Fig. 4 for $`v=1`$, 43 and 86-GHz lineshape calculations. The variation of peak maser lineshape intensity as a function of stellar phase for 43 and 86 GHz masers is shown in Fig. 5.

(i) The general forms of the data in Figs. 3 and 4 are similar to those found in the most extensive data available (NO86, MBA88).

(ii) For much of the cycle, Fig. 5 shows that the peak intensity of the 86-GHz masers is greater than that of the 43-GHz masers at the corresponding epoch, as typically observed for Miras (Pardo et al. pardo (1998)). During epochs of bright maser emission (Epochs 13 - 20), the ratio of the peak lineshape intensities, $`v=1`$ $`J=21`$/$`v=1`$ $`J=10`$ $``$ 1 - 3. Cho et al. (cho98 (1998)) find a comparable typical value of $``$ 1 - 2.

(iii) In Alcolea et al. (alcolea99 (1999)), maser minimum typically passes to maximum over $``$ 0.5 periods (i.e. the asymmetry factor $`f`$ is 0.5) for the $`v=1`$ 43-GHz data. Our actual values are 0.35 and 0.3 stellar periods for $`v=1`$ 43-GHz and $`v=1`$ 86-GHz masers respectively. However, within the typical statistical uncertainties, given the number of emitting components contributing to the lineshape at each phase, our data is consistent with observations. In Fig. 5, maser minimum may occur between Epochs 7 – 11 and the maximum between Epochs 15 – 17.

(iv) The FWHM of periods of bright simulated maser intensity are 0.2 periods for the 43-GHz masers and 0.25 for 86-GHz masers. These values underestimate observed values which lie in the approximate range 0.3 - 0.7 for the Miras in Alcolea et al. (alcolea99 (1999)).

(v) By equating the model phase of maser minimum light with the observed optical phase of 0.67 (Sect. 2.4) for the observations of $`v=1`$ 43-GHz masers in TX Cam, we derive the optical maximum of our model star to correspond to a model phase of 0.78, i.e. Epoch 16 - 17. Maser maximum for $`v=1`$ $`J=10`$ emission, which occurs at Epoch 17 (phase 0.85) in our simulations, therefore peaks somewhere between the optical and IR peaks if the IR peak is delayed by 0.1 periods with respect to its optical counterpart. The current model is therefore not inconsistent with a small ($`<`$0.1 period) phase difference between low-frequency maser and IR peaks. Although the mean observed phase lag between optical and maser maximum light is $``$0.1 - 0.2, we noted in Sect. 2.1 (ii) that this value is very variable between stars and between different cycles of the same star. Indeed, for some cycles of some of the stars studied in Alcolea et al. (alcolea99 (1999)), for example R Cas and U Her, there is also negligible phase lag between stellar and maser maxima.

(vi) The contrast (ratio of maximum to minimum peak lineshape intensity during the cycle) in our simulations is 175 at 43-GHz and 850 at 86-GHz. As the observed contrasts are typically less than $``$10 $`\%`$ of these values, this is a clear failing of our simulations, though we note that very large contrasts have been observed in $`o`$ Cet (NO86): see Section 2.1.(v). In Fig. 3 it is evident that maser emission is almost eradicated between Epochs 6 to 11 (0.25 stellar periods). We discuss the reasons why the model fails in this respect in Sect. 5.

(vii) Fig. 4 shows that the bulk of emission originates from velocity extents of $``$12 kms-1, with broader low intensity wings extending from around -15 to +15 kms-1, as observed by Herpin et al. (1998). Herpin et al. (1998) consider a number of processes which may give rise to this emission, such as turbulence and asymmetric mass loss. Since this CE model does not include turbulence and is spherically symmetric, it appears that the wings can arise simply through weak emission originating from regions of the CE which are infalling or outflowing at relatively high velocity.

(viii) The cycle-averaged peaks of both spectra in Fig. 4 occur close to the stellar velocity, peaking 0.3 kms-1 to the blue of V. The cycle-averaged spectral peak at $`v=1`$, $`J=10`$ was 0.3 kms-1 to the red of V in Cho et al.(cho96b (1996)), compared to the blueshift of $``$1 kms-1 found by NO86. Adopting the data calibration in (v), we find that in our simulated data $`v=1`$ $`J=10`$ emission is redshifted with respect to V for an optical phase range of 0.73 – 0.98 and blue-shifted from 0.23 – 0.63. This does not agree with the statistical result of Cho et al. (cho96b (1996)) in Sect. 2 (vi).

### 4.2 Spectra of SiO masers at 43 and 86-GHz: short term variability

The present model cannot deal with the very short term variability of a few days or less reported in Pijpers et al. (pijpers94 (1994)), which the authors attribute to short wavelength sound waves generated through convection, a phenomenon not included in the present stellar pulsation model (see Pijpers & Hearn pijpers89 (1989)). However conditions in the CE do change significantly between model epochs, on timescales of $``$17 days. This is clearly illustrated by Fig. 3 and 5.

### 4.3 Synthetic images at 43-GHz and 86-GHz: survival and proper motion of maser components

H96 showed that masers at model phase zero (Epoch 1) were disposed in an approximate ring around the host star at a radius of $``$1 R from the photosphere. The present work shows that this ring expands, contracts, undergoes severe disruption and reappears according to the stellar phase. Fig. 6 shows the time series of synthetic images resulting from our simulations at 43-GHz, which can be displayed for any transition. The interval between epochs in Fig. 6 is 16.6 days, corresponding to 0.05 stellar periods. The choice of 1500 maser sites in order to generate these images was based on numerical experiments at Epoch 1. It would now appear that a smaller number might have been more appropriate, since the ring structures in Epochs 14 – 18 are rather more complete than observations (of TX Cam) currently suggest. Fig. 7 shows a VLBA image of the 43-GHz masers in TX Cam, for comparison with our simulated data. Fig. 8 shows maser ring radii for 43 and 86-GHz masers as a function of the stellar pulsation phase. At all events a number of useful conclusions can be drawn from these data.

(i) In Fig. 6, it is evident that tangential amplification is the norm, but that an occasional feature may form over the disk of the star, as at Epoch 17 at 43-GHz. Data for TX Cam indeed show that an occasional maser spot may form over the disk of the star, as shown in Fig. 7.

(ii) The radius of the 43-GHz ”ring” of masers varies between 1.93 – 2.3 AU in Fig. 8. The average radius of the maser ring declines between Epochs 1 – 6, at a mean rate of 0.2 AU in 83 days, corresponding to 4.2 kms-1 (cf. infall of 4 kms-1 in R Aqr over a similar duration of the stellar cycle). Contraction of the ring is followed by a larger shock-driven expansion at a rate of 3.2 kms-1 (cf. dominant expansion in the TX Cam observations, at a typical velocity of 3.65 kms-1). In our simulations, expansion of the ring is accompanied by a period of bright maser emission, induced by the temperature and density enhancements of the post-shock gas.

(iii) The radius of the $`v=1`$ 86-GHz ring of masers is consistently similar, to within 0.1 AU, to that of the $`v=1`$ 43-GHz ring throughout the stellar cycle. Gray & Humphreys (gray2000 (2000)) show that this is not the case for masers in different vibrational states. The radial motions of the rings are highly coupled as bright $`v=1`$, 43 and 86-GHz maser emission often originates from shared components in the CE. This is supported by observational data of Colomer et al.(colomer96 (1996)) and Doeleman et al.(doele (1998)), although their multi-transition results were not made simultaneously.

(iv) More sites in the CE produce maser emission at 43-GHz (16 $`\%`$ of the 1500 sites at Epoch 15) than at 86-GHz (7.5 $`\%`$ at the same epoch). However, of the components which do emit, $`v=1`$ $`J=21`$ emission is typically significantly stronger than that of at $`v=1`$ $`J=10`$, leading to the more intense images and lineshapes in our data. (There is a general trend in our data for fewer emitting components as a function of increasing $`J`$-value for the transition.) It is not clear why the $`v=1`$ $`J=21`$ should be the stronger. Results in D95 also suggested that 43 and 86-GHz emission in $`v=1`$ form under similar conditions. As discussed in D95, competition betwen masers and cycling of populations between the $`v=1`$ and $`2`$ and $`J=1`$ and $`2`$ states result in a complicated interplay whose outcome determines whether the 43 or 86-GHz line will emerge as a strong maser.

(v) In our simulated data, the passage of a shock front through the maser zone causes minimum maser light. The combination of physical conditions in the zone are unsuitable for yielding bright maser emission, in particular, the high temperatures which dissociate SiO. Although temperatures in real stars are unlikely to be as high as in this model, we believe that the general picture provided here is the correct one. That is to say a shock front disrupts the existing maser ring and new features then form in the gas in the wake of the shock i.e. a new maser ring will appear to form which has a smaller angular extent than the previous ring.

## 5 Concluding remarks

These simulations have some rather severe limitations, which we discuss below. However, a significant conclusion of the present work is that coupling a SiO maser model to a pulsating AGB stellar model, keeping the stellar IR radiation field constant throughout the cycle, reproduces qualitatively much of the available observational data. This is particularly evident in the location, tangential amplification and proper motion of SiO masers (although clearly our spherically symmetric model will not reproduce asymmetrical effects) in the CE. The essential point is that shocks in the inner CE have an effect on SiO maser emission which, given the lack of time-dependent chemistry and efficient molecular cooling in our present simulations, could explain much of the SiO maser variability phenomenon. Of course the inclusion of a varying stellar IR radiation field, in addition to the improvements outlined above, is essential for modelling fully the environments provided by such stars. We note that it is not yet established whether SiO masers are predominantly radiatively or collisionally pumped.

A number of important quantitative shortcomings are evident in our simulations, which may be summarised as follows. The duration of the period of bright maser emission is underestimated, the contrast in the maser lightcurve is overestimated, maser emission is very weak for a significant portion of the stellar cycle and in observed data, red-shifted emission dominates for the greater portion of the stellar cycle, whereas in the simulated data blue-shifted emission dominates. In this discussion we should recall that the model parameters set out in Table 1 are appropriate to $`o`$-Ceti, and these represent an arbitrary choice which may not necessarily yield the behaviour which may be characterised as “typical” of M-Miras.

An important problem remains the difficulty in relating accurately the phase of our model star to stellar optical phase. As we have stated, the current hydrodynamic model is not really capable of deciding this issue because it has a fixed-temperature photosphere. VLBI observations, however, are at present consistent with the view that the shock wave in the envelope arrives at the maser zone in a reasonably well-constrained phase range for a group of objects with very different stellar parameters. Assuming that the observed maser ring in a VLBI experiment corresponds well to the radius of the shock at the time the shock impacts the maser zone - as is the case in our model - we deduce from the five stars in Table 2, that the mean phase, $`\overline{\varphi }`$ of shock arrival is 0.71 in the modellers’ definition, with a standard deviation on the mean, $`(_{j=1}^n(\varphi _j\overline{\varphi })^2/(n(n1)))^{1/2}`$, equal to only 0.029 periods. In the case of R Cas, Phillips et al. (phillips2001 (2001)) give sufficient data to make a reasonable estimate of the root-mean-square uncertainty in the time required for the shock to reach the maser zone. The error in the angular measure is dominated by the uncertainty in the stellar diameter, quoted as 30$`\pm `$6 mas. The fractional error in the angular distance to be crossed by a shock travelling to the maser ring, of diameter 57 mas, is therefore 6/27 = 0.22. To this we add the fractional error in the distance to R Cas, about 0.2 from the authors’ figures, to get 0.44. This level of uncertainty is certainly large, but we are basing our results on a group of five, not on a single source. Assuming that 0.44 is a reasonable fractional error for all five objects, we obtain a root-mean-square uncertainty on the mean for the five objects in Table 2 of 0.15 periods. This is certainly larger than the standard deviation on the mean, but is still usefully restrictive. This information is quite independent of any assumed pumping scheme, as it is based only on observation. Given that most objects have a roughly fixed delay from shock impact to maser maximum, and a further constant phase shift links the modellers’ and optical phase definitions (see Section 4.1.5) then we expect a fairly constant phase shift between the peaks of the optical (or IR) light curves and the maser peak. Observationally, this shift is close to 0.0 in the IR case for most, but not all objects. This link would be quite natural if the shock and IR photosphere are dynamically linked in many objects. Such a link is plausible because much of the opacity in the 3-10 $`\mu `$m range is derived from water molecules. The shock would provide a dense layer which could be optically thick for part of the stellar cycle. Observational support for this view comes from observations by Yamamura, de Jong & Cami (yama (1999)) who found that the warmer of two water layers, located at around two nominal stellar radii in $`o`$ Cet was optically thick. The high temperature suggests the layer is shock-heated, and the mid-IR photosphere would therefore lie somewhere in the post-shock gas.

A significant source of error in these calculations is likely to be the lack of molecular coolants such as CO in this model. The shocks in the present work (see Fig. 1a) are quite close to being adiabatic. The effect of coolants, such as CO is not included and cooling of the shocks in the present model takes place broadly on the hydrodynamic timescale. This timescale depends on the phase of the star and lies in the approximate range of 60-100 days. If the effects of cooling by CO are included (Cuntz & Muchmore cuntz (1994); Woitke et al. woitke (1996)), the shocks cool on a radiative timescale of only 1 – 2 days and the shocks are much closer to isothermality than those used in the present work. This means that SiO maser emission would not be severely weakened for a large portion of the stellar cycle, as it is in the current work, partly due to the dissociation of SiO. We would therefore expect more accurate shock conditions, and possibly time-varying IR radiation, to reduce the maser dynamic range from the high levels found in the simulations. Time-dependent chemistry is also a requirement for providing the abundance of SiO and of other coolant species.

A number of the discrepancies mentioned above may also arise from the use of a spherical pulsation model. This yields too clean-cut behaviour of the CE. Variations in the physical conditions in the CE almost certainly do not take place in the same manner equally in all directions. Moreover clumping of material (Sect. 4.1) due to thermal instabilities is inherently spatially inhomogeneous. Given non-spherical behaviour (Sect. 2.4), the CE at any one phase of stellar pulsation will exhibit a spread of physical conditions for any given radius, rather than the unique set used here. This will increase the duration of bright emission and reduce the contrast between maximum and minimum maser brightness. In this connection the VLBA data on TX Cam at a variety of phases show that the entire maser ring never fades away in its entirety as our present models suggest (Diamond: private communication). The rings wax and wane in intensity but $``$50% of the ring generally remains visible. An example of a prediction of the model that is verified by observations of maps at different stellar phases is that maser emission is not exclusively tangential and that an occasional maser spot may indeed be found over the disk of the host star, as we predict.

To improve our understanding of the extended atmospheres of AGB stars, the following are necessary. The relationship between the optical and model phases should be reliably fixed by a combination of interferometry observations and theory. The hydrodynamic pulsation model should in future incorporate a varying IR radiation field, which in turn requires the hydrodynamic model to incorporate an accurate prediction of the radius of the IR photosphere as any phase. Additionally, we require some chemical modelling of SiO formation, of the H, H2 balance and of coolant molecules such as CO and H2O. As noted above, this is significant in determining the radiative timescale for cooling of the shocks in the CE. Chemical modeling also has implications both for the relative amplification factors at different phases and for the nature of the collisional pumping mechanism, since H and H2 have significantly different rate coefficients for rotationally and rovibrationally inelastic collisions with SiO in the new calculations by Jimeno et al. (jim (1999)). In addition, a time-dependent description of silicate, that is, olivine dust formation should be introduced, noting that dust formation also affects the local thermal balance. At present, the importance of including magnetic field effects in the CE is unclear. With respect to the maser model, accelerated lambda iteration methods should be used to replace LVG (Jones et al. jones (1994); Randell et al. randell (1995); Yates et al. yates97 (1997)). Accelerated lambda iteration methods have the important characteristic that they can incorporate the velocity, temperature and number density structure of the CE explicitly into the calculation of the populations of SiO rovibrational levels. These methods require development for a spherical geometry for molecular systems and are computationally expensive, though not prohibitively. New rate coefficients for collisional energy transfer between H2 and SiO must also be computed. The outstanding observational requirement is for more sets of VLBI images disposed at equal intervals of $``$0.1 in phase (near-)simultaneously for several SiO maser transitions towards suitable Mira long-period variables.

###### Acknowledgements.

EMLH thanks the Swedish Foundation for International Co-operation in Science (STINT) for the financial support of this project.


format: text

text:
# Search for New Physics in Photon-Lepton Events in 𝑝𝑝̄ Collisions at √𝑠=1.8TeV

## Abstract

We present the results of a search in $`p\overline{p}`$ collisions at $`\sqrt{s}=1.8\mathrm{TeV}`$ for anomalous production of events containing a photon and a lepton ($`e`$ or $`\mu `$), both with large transverse energy, using 86 pb-1 of data collected with the Collider Detector at Fermilab during the 1994-95 collider run at the Fermilab Tevatron. The presence of large missing transverse energy ($`\overline{)}E_T`$), additional photons, or additional leptons in these events is also analyzed. The results are consistent with standard model expectations, with the possible exception of photon-lepton events with large $`\overline{)}E_T`$, for which the observed total is 16 events and the expected mean total is $`7.6\pm 0.7`$ events.

preprint: FERMILAB-PUB-02/031-Epreprint: EFI-02-66

An important test of the standard model (SM) of particle physics SM  is to measure and understand the properties of the highest-energy particle collisions. The observation of an anomalous production rate of any combination of the fundamental particles of the SM would be a clear indication of a new physical process. This Letter summarizes an analysis of the inclusive production of a photon and a lepton ($`e`$ or $`\mu +\gamma +X`$), including searches for additional photons, leptons, and large missing transverse energy, using 86 pb-1 of data from proton-antiproton collisions collected with the Collider Detector at Fermilab (CDF) during the 1994-95 run of the Fermilab Tevatron jeffPRD .

Production of these particular combinations of particles is of interest for several reasons. Events with photons and leptons are potentially related to the puzzling “$`ee\gamma \gamma \overline{)}E_T`$” event recorded by CDF eegg . A supersymmetric model susy  designed to explain the $`ee\gamma \gamma \overline{)}E_T`$ event predicts the production of photons from the radiative decay of the $`\stackrel{~}{\chi }_2^0`$ neutralino, and leptons through the decay of charginos, indicating $`\mathrm{}\gamma \overline{)}E_T`$ events as a signal for the production of a chargino-neutralino pair. Other hypothetical, massive particles could subsequently decay to SM electroweak gauge bosons, one of which could be a photon and the other a $`W`$ or $`Z^0`$ boson that decays leptonically. In addition, photon-lepton studies complement similarly-motivated inclusive searches for new physics in diphoton gg , photon-jet gj , and photon-$`b`$-quark events gb .

The CDF detector cdfblurb  is a cylindrically symmetric spectrometer designed to study $`\overline{p}p`$ collisions at the Fermilab Tevatron. A superconducting solenoid of length 4.8 m and radius 1.5 m generates a magnetic field of 1.4 T and contains tracking chambers used to measure the momenta of charged particles. A set of vertex time projection chambers is used to find the $`z`$ position CDFcoo  of the $`\overline{p}p`$ interaction. The 3.5-m-long central tracking chamber (CTC) is a wire drift chamber which provides up to 84 measurements between the radii of 31.0 cm and 132.5 cm in the region $`|\eta |<1.0`$. Sampling calorimeters, used to measure the electromagnetic and hadronic energy deposited by electrons, photons, and jets of hadrons, surround the solenoid. Each tower of the central ($`|\eta |<1.1`$) electromagnetic calorimeter (CEM) has an embedded strip chamber for the measurement of the 2-D transverse profile of electromagnetic showers. Muons are detected with three systems of muon chambers, each consisting of four layers of drift chambers. The central muon (CMU) system is located directly outside the central hadronic calorimeter, and covers $`|\eta |<0.6`$. Outside of the CMU is 0.6 m of steel shielding, followed by the central muon upgrade system. The central muon extension system provides muon detection for $`0.6<|\eta |<1.0`$.

Events with a high-transverse momentum ($`p_T`$EtPt  photon or lepton are selected by a three-level trigger jeffPRD , which requires an event to have either a high-$`E_T`$ photon or a high-$`p_T`$ lepton ($`e`$ or $`\mu `$) within the central region, $`|\eta |<1.0`$. Photon and electron candidates are chosen from clusters of energy in adjacent CEM towers; electrons are then further separated from photons by requiring the presence of a CTC track pointing at the cluster. Muons are identified by requiring CTC tracks to extrapolate to a reconstructed track segment in the muon drift chambers.

To reduce background from the decays of hadrons produced in jets, both the photon and the lepton in each event are required to be ‘isolated’. The $`E_T`$ deposited in the calorimeters in a cone in $`\eta \varphi `$ space of radius $`R=0.4`$ around the photon or lepton position is summed, and the $`E_T`$ due to the photon or lepton is subtracted. The remaining $`E_T`$ in the cone, $`E_{cone}^{iso}`$, is required to be less than 2$`\mathrm{GeV}`$ for a photon, or less than 10% of the lepton transverse momentum. In addition, for photons the sum of the $`p_T`$ of all tracks in the cone must be less than 5$`\mathrm{GeV}`$.

Inclusive photon-lepton events are selected by requiring an isolated central photon with $`\mathrm{E}_\mathrm{T}^\gamma >25`$ GeV and an isolated central lepton ($`e`$ or $`\mu `$) with $`\mathrm{E}_\mathrm{T}^{\mathrm{}}>25`$ GeV. The technical criteria used to identify leptons and photons are very similar to those of Refs. gg ; gj ; gb , and are described in detail in Ref. jeffPRD . A total of 77 events pass this selection: 29 photon-muon and 48 photon-electron candidates.

The production of pairs of new heavy states that decay via cascade decays can lead to final states with multiple photons or leptons; in contrast, the dominant SM background processes lead to signatures with only one photon and one lepton observed in the detector, as discussed in detail below. The inclusive sample is consequently analyzed as two subsamples: a “two-body inclusive photon-lepton sample” typical of a two-particle final state, and a “multi-body inclusive photon-lepton sample” typical of three or more particles in the final state. The two-body sample selection requires exactly one photon and exactly one lepton, with an azimuthal separation $`\mathrm{\Delta }\phi _\mathrm{}\gamma >150^{}`$, but excludes those events for which the invariant mass of the photon and electron, $`M_{e\gamma }`$, is within 5$`\mathrm{GeV}`$ of the mass of the Z boson, $`M_Z`$ (these events are used as a control sample, as described below). The multi-body sample is composed of the remaining inclusive photon-lepton events. The multi-body sample is then further analyzed for the presence of large $`\overline{)}E_T`$, and additional isolated leptons and photons. The $`\overline{)}E_T`$ threshold of 25$`\mathrm{GeV}`$ was determined a priori in previous analyses r\_papers  as a significant indicator of a neutrino arising from leptonic decays of the $`W`$ boson. Figure 1 shows the breakdown of the inclusive sample into the final categories.

The dominant source of photon-lepton events at the Tevatron is electroweak diboson production, in which a $`W`$ or $`Z^0`$ boson decays leptonically ($`\mathrm{}\nu `$ or $`\mathrm{}\mathrm{}`$) and a photon is radiated from either an initial-state quark, a $`W`$, or a charged final-state lepton. The number of such events is estimated using leading-order (LO) matrix element calculations diboson matrix  for which the computational code diboson mc  was then embedded into the general-purpose event generator program pythia pythia , followed by a full simulation of the detector. The uncertainty in this number has roughly equal contributions from higher-order processes, simulation systematics, luminosity, proton structure, and generator statistics.

A jet can contain mesons such as the $`\pi ^0`$ or $`\eta `$ that decay to photons, which then may satisfy the photon selection criteria. The number of lepton-plus-misidentified-jet events is determined by counting the number of jets in a sample of events with a lepton and then multiplying by the probability of a jet being misidentified as a photon, $`P_\gamma ^{jet}`$. The factor $`P_\gamma ^{jet}`$ is determined from samples of jets and photons in events with a lepton trigger, using the distribution in $`E_{cone}^{iso}`$. By fitting to a sum of the expected distribution for prompt photons and that measured for jets, the misidentification rate is found to be $`P_\gamma ^{jet}=(3.8\pm 0.7)\times 10^4`$ jeffPRD .

The dominant source of misidentified photon-electron events is $`Z^0e^+e^{}`$ production, where one of the electrons undergoes hard photon bremsstrahlung or a track fails to be reconstructed. We assume that photon-electron events consistent with $`Z^0`$ production are not a significant source of new physics, and use them to estimate the probability $`P_\gamma ^e`$ that an electron is reconstructed as a photon. The number of misidentified photon-electron events in the control sample divided by the number of electron-electron events with the same kinematics gives $`P_\gamma ^e=(1.28\pm 0.35)`$%. For any other subset of central electron pairs, the contribution to the corresponding photon-electron sample is the product of $`P_\gamma ^e`$ and the number of central electron pairs.

Other, smaller, backgrounds are due to hadrons faking muons, and to leptons from the decay of bottom and charm quarks. Charged hadrons may penetrate the calorimeters into the muon chambers, or may decay to a muon before reaching the calorimeters. These contributions are determined by identifying isolated, high-momentum tracks in the inclusive photon sample, applying the probability of each track of being misidentified as a muon, and summing this probability over all tracks in the sample jeffPRD . The contribution to photon-lepton candidates from heavy-flavor produced in association with a prompt photon is estimated using Monte Carlo event generation pythia  and detector simulation, and found to be negligible.

New physics in small samples of events would most likely manifest itself as an excess of observed events over expected events. The significance of an observed excess is computed from the likelihood of obtaining at least the observed number of events, $`N_0`$, assuming that the null hypothesis (the SM) is correct. The “observation likelihood”, $`P(NN_0|\mu _{SM})`$, is defined as the fraction of the Poisson distribution for the number of expected events from SM sources, with a mean $`\mu _{SM}`$, that yields outcomes $`NN_0`$ delta\_mu . The likelihood $`P(NN_0|\mu _{SM})`$ is computed from a large ensemble of calculations in which each quantity used to compute photon-lepton event sources varies randomly as a Gaussian distribution, and the resulting mean event total is used to randomly generate a Poisson-distributed outcome $`N`$. The fraction of calculations in the ensemble with outcomes $`NN_0`$ gives $`P(NN_0|\mu _{SM})`$.

The predicted and observed totals for two-body photon-lepton events are compared in Table 1. Half of the predicted total originates from $`Z^0\gamma `$ production where one of the charged leptons has escaped identification; the other half originates from roughly equal contributions of $`W\gamma `$ production, misidentified jets, misidentified electrons, and misidentified charged hadrons. The likelihood of the observed total is 9.3%.

The predicted and observed totals for inclusive multi-body photon-lepton events are also compared in Table 1. About half of the predicted total originates from $`Z^0\gamma `$ production, a quarter from $`W\gamma `$ production, and the remaining quarter from particles misidentified as photons or leptons. The likelihood of the observed inclusive multi-body total is 10%. The predicted and observed kinematic distributions for these events are compared in Figure 2. The difference between the observed and predicted totals can be entirely attributed to events with $`\overline{)}E_T>25\mathrm{GeV}`$. Figure 2 also shows the distribution in $`H_T`$, the scalar sum of the $`E_T`$ of all objects in the event plus the magnitude of $`\overline{)}E_T`$, a variable correlated with the production of massive particles jeffPRD .

The predicted and observed totals for multi-body $`\mathrm{}\gamma \overline{)}E_T`$ events are also compared in Table 1. For photon-electron events, requiring $`\overline{)}E_T>25\mathrm{GeV}`$ suppresses the contributions from $`Z^0\gamma `$ production and from electrons misidentified as photons, which have no intrinsic $`\overline{)}E_T`$, while preserving the contribution from $`W\gamma `$ production. As a result, 57% of the predicted $`e\gamma \overline{)}E_T`$ total arises from $`W\gamma `$ production, 31% from jets misidentified as photons, only 3% from $`Z^0\gamma `$ production, and the remaining 9% from other particles misidentified as photons. The observed $`e\gamma \overline{)}E_T`$ total agrees with the predicted total, with a 25% probability that the predicted mean of $`3.4`$ events yields 5 observed events. One of these 5 is the $`ee\gamma \gamma \overline{)}E_T`$ event eegg .

For photon-muon events, requiring $`\overline{)}E_T>25\mathrm{GeV}`$ does not completely eliminate the contribution from $`Z^0\gamma `$, for a second muon with $`|\eta |>1.2`$ and $`p_T>25\mathrm{GeV}`$ can escape detection and induce the necessary amount of $`\overline{)}E_T`$. The rate at which this occurs is modeled well by the $`Z^0\gamma `$ event simulation, however, since it is largely a function of the CDF detector geometric acceptance. Of the 4.6 multi-body events predicted to originate from $`Z^0\gamma `$ production, 2.2 events are predicted to contain a second visible muon, 1.0 events are predicted to have less than 25$`\mathrm{GeV}`$ of $`\overline{)}E_T`$, and only 1.0 events are predicted to pass the 25$`\mathrm{GeV}`$ $`\overline{)}E_T`$ selection. One event is observed with a second muon, in agreement with the $`Z^0\gamma `$ prediction. The predicted total for multi-body $`\mu \gamma \overline{)}E_T`$ events consists of 47% $`W\gamma `$ production, 24% events with jets misidentified as photons, 23% $`Z^0\gamma `$ production, and the remaining 7% from particles misidentified as muons.

The $`\mu \gamma \overline{)}E_T`$ event total is higher than predicted (11 observed vs. 4 expected), with an observation likelihood of 0.54%; the observation likelihood of the $`\mathrm{}\gamma \overline{)}E_T`$ total is only slightly higher at 0.72% e\_mu\_diff . The predicted and observed distributions of the kinematic properties of multi-body $`\mathrm{}\gamma \overline{)}E_T`$ events are compared in Figure 2. The observed photon $`E_T`$, lepton $`E_T`$, $`\overline{)}E_T`$, and $`H_T`$ distributions are within the range expected from the SM referee\_motivation .

The predicted and observed totals of multi-lepton events are compared in Table 1. Nearly all of the predicted total is expected from $`Z^0\gamma `$ production. Approximately 6 events are expected; 5 events are observed, including the $`ee\gamma \gamma \overline{)}E_T`$ event. No $`e\mu \gamma `$ events were expected, and none were observed.

The predicted number of multi-photon events is dominated by $`Z\gamma `$ production, for which only 0.01 events are expected. The single event observed is the $`ee\gamma \gamma \overline{)}E_T`$ event, whose (un)likelihood is described in Ref. eegg .

In conclusion, we have made an a priori search for inclusive photon+lepton production. We find that subsamples of this data set agree well with their SM prediction, with the possible exception of $`\gamma \mathrm{}\overline{)}E_T`$. However, an excess of events with 0.7% likelihood (equivalent to 2.7 standard deviations for a Gaussian distribution) in one subsample among the five studied is an interesting result, but is not a compelling observation of new physics. We look forward to more data in the upcoming run of the Fermilab Tevatron.

###### Acknowledgements.

We thank the Fermilab staff and the technical staffs of the participating institutions for their contributions. We thank U. Baur and S. Mrenna for their critical calculations of the SM $`W\gamma `$ and $`Z\gamma `$ backgrounds used in this analysis. This work was supported by the U.S. Department of Energy and National Science Foundation; the Italian Istituto Nazionale di Fisica Nucleare; the Ministry of Education, Culture, Sports, Science, and Technology of Japan; the Natural Sciences and Engineering Research Council of Canada; the National Science Council of the Republic of China; the Swiss National Science Foundation; the A. P. Sloan Foundation; the Bundesministerium fuer Bildung und Forschung, Germany; the Korea Science and Engineering Foundation (KoSEF); the Korea Research Foundation; and the Comision Interministerial de Ciencia y Tecnologia, Spain.


text:
# Comment on “Exact Spectral Functions of a Non Fermi Liquid in 1 Dimension”

## Abstract

We point out the incorrect statement of the recent manuscript by K. Penc and B. S. Shastry.

In a recent manuscript  K. Penc and B. S. Shastry wrote in the first paragraph:

“Schulz and Shastry  have introduced a new class of gauge-coupled one-dimensional (1D) Fermi systems that are non Fermi liquid in the sense that the momentum distribution function has a cusp at the Fermi momentum $`k_F`$ rather than a jump as in a Fermi liquid . This behavior is of the sort first found by Luttinger in the context of his study of a one-dimensional model that is popularly known as the Luttinger model . The model introduced by Schulz and Shastry (SS) is in fact intimately connected to the Luttinger model, and is best viewed as a reinterpretation of Luttinger’s original model as a gauge theory. Particles of different species exert a mutual gauge potential on each other, and this is sufficient to destroy the Fermi liquid. This model has the added property that the charge and spin correlations are unaffected by the interaction, owing to the ‘gauge’ nature of interaction.”

I would like to point out that the main statement of this paragraph is incorrect. In fact, this can be seen from the second reference from  (which is the Reply to my previous Comment ) and that Comment itself. It turns out that all the properties, mentioned in the cited paragraph, namely —

* The gauge coupling between 1D Fermi systems, which produces the non Fermi liquid behavior (the Luttinger liquid behavior);
* Particles of different species in that class exert a mutual gauge potential on each other (which is sufficient to destroy the Fermi liquid);
* The charge and spin correlations are unaffected by the interaction, owing to the “gauge” nature of interaction

— had been already known for the class of models, introduced by us in Refs.  in 1992. This is why the priority of introducing the models with these properties belongs to us. One can check that particles in  of different species are connected with each other via ‘gauge’ potentials (reminiscent of the Peierls phase factors), similar to the ‘gauge’ potentials in , and that we emphasized on the Luttinger liquid behavior of our 1D models with this ‘gauge’ couplings already in 1992, much earlier than . Even the title of  is “Exactly solvable models of an effectively two-dimensional Luttinger liquid”. Notice that in  we interpreted the additional index, which distinguishes species of particles, as a number, which enumerates 1D chains, coupled with each other via gauge potentials. By the way, SS in their first paper of  never used this interpretation, but now, in , the authors already imply that the class of models, introduced by SS described coupled 1D chains.

According to the Reply , SS introduced some class of models, which had all mentioned above properties of the class, introduced in , but with the additional constraint: Some recursion relations for ‘gauge’ potentials have to be satisfied. SS especially pointed out that difference, see, please, the footnote  of the Reply , where they wrote: “We use the term ‘class of models’ in the specific sense that the members share a common method of solution, rather than a vague sense in which many models share certain physical properties.” Obviously, only the properties, mentioned in the first paragraph of , cannot properly define the class of models, introduced by SS, because they belong to both classes:  and .

Actually, there are two alternatives:

* (1) Either SS introduced a subclass (with some specific, additional properties) of the class of models  (which had been earlier introduced in our papers) with the common properties, mentioned in the disputed paragraph. Certainly, one cannot introduce any new class of models in 1998 with the same properties as the models, introduced in 1992.
* (2) Or SS (cf. their Reply ) introduced some new class of models, different from ours. However in this case the authors of the manuscript  could clearly define the properties, which determine only the models, introduced by SS, but which are not present in the class, introduced earlier in . At least they could carefully distinguish between those two classes, and not to emphasize on the properties, which belong to the other class of models, introduced in  in the definition of the models, introduced by SS.

In both cases the statements of the disputed paragraph  are wrong. I do not imply that the model, studied in  belongs to the class of models, introduced in  in the sense (2). However, clearly, when writing about the properties of the models, introduced in , the authors of  could properly write about the features, which pertain only to the class introduced by SS, but not about the ones, which had been known for the other class of models, introduced earlier by us in . It turns out that according to the definition of the Reply , all mentioned in the disputed paragraph of  properties have namely “a vague sense in which many models share certain physical properties”.

This is why, immediately after the manuscript  appeared in the ArXive, I asked the authors of  to correct the statement of the disputed paragraph. However the authors did not agree, i.e., they insist that the class of models with the properties, mentioned in the cited paragraph of their work, was introduced in . But then a contradiction exists: If the statements of the Reply  are correct, then the statements of the first paragraph of the manuscript  are obviously wrong.

It is the goal of this my Comment to point out this contradiction.