© ESO, 2012

1. Introduction

Stars of less than eight solar masses evolve onto the asymptotic giant branch (AGB) when they have exhausted much of the available H. They become very large (R_⋆ > 1 AU) and cool (T_⋆ ≈ 2500 K) and lose mass at the rate of up to an Earth mass per year. Red supergiants (RSG) had main sequence masses  > 8 M_⊙ and swell to ten times the size of AGB stars, with ten times the mass loss rate. Although the internal stellar processes are mass-dependent, the winds are similar, being cool, dusty and molecular. O-rich stars produce SiO, H₂O and OH masers with spatially very compact, spectrally narrow lines which can be imaged at milli-arcsec resolution using radio interferometry. A combination of pulsations, radiation pressure and possibly other e.g. magnetic forces drives material away from the stellar surface. The exact mechanisms are still under debate (e.g. Woitke 2006) and at a few stellar radii (R_⋆), SiO masers show both infall and outflow (e.g. Assaf et al. 2011). Dust nucleation is likely to start close to the star (Wittkowski et al. 2007) and by 4–6 R_⋆, radiation pressure on dust is the main force driving the wind. The H₂O masers around these stars are found in approximately spherical shells extending from about 5–50 R_⋆. They provide a very high-resolution probe of the structure, kinematics and conditions in the circumstellar envelope (CSE) in the region where the outflow is most strongly accelerated and passes through the escape velocity, see e.g. Bowers & Johnston (1994), Yates & Cohen (1994) and also Habing (1996) for a general review.

MERLIN¹ has been used to monitor H₂O and OH masers around a number of evolved stars with no known companions. We analyse 2–7 epochs of images of 22-GHz H₂O masers around seven objects, listed in Table 1. The observations of the RSG VX Sgr and S Per (first epoch) were reported in Murakawa et al. (2003) (M03) and Richards et al. (1999) (R99), respectively. The first epochs of the AGB stars IK Tau, U Her, U Ori and RT Vir were reported in Bains et al. (2003) (B03). W Hya MERLIN 22-GHz observations have not been published previously.

These stars are close enough and unobscured enough to have AAVSO optical monitoring data, Hipparcos or other parallax measurements and in many cases optical or IR interferometric meaurements of the stellar diameter. They are among the objects which have been monitored by the Pushchino radio telescope at 22 GHz for many decades. Component fitting allows the size of AU-scale clouds to be measured from MERLIN data. The H₂O maser properties of the sources above declination 0° were investigated in detail by Richards et al. (2011) (R11), who found that they are predominantly unsaturated. The maser appearance is mostly characteristic of spherical clouds but, in some objects, amplification appears to be occuring along the long axis of flattened clouds.

The H₂O masers have also been imaged using VLBI and the VLA, see Table 1. However, a large fraction of the flux (half is not unusual) is resolved-out by VLBI, whilst the VLA cannot resolve individual maser clumps. Allowing for this, our observational results are consistent with other images. Comparisons with previous interpretations are included in the relevant sections of this paper. For completeness, we note that we observed several other objects with MERLIN at 22 GHz. The first epochs of NML Cyg and VY CMa were published by Richards et al. (1996, 1998); further monitoring of these two exceptional RSG will be presented by Yates et al., as will obervations of the heavily-obscured SRb, R Crt. We also observed R Cas (20000404, 20010511, 20020405) and R Leo (20000429) but these were non-detections at a conservative upper limit of 200 mJy. This is consistent with non-detections for several years around this period (upper limit 10 Jy) using Pushchino (Pashchenko & Rudnitskii 2004, R Cas; Esipov et al. 1999, R Leo).

This paper presents a large-scale analysis of multiple observations (between two and seven epochs per star, spanning up to 7 years). Data acquisition and analysis is summarised in Sect. 2, including estimating the centre of expansion, aligning epochs and assessing the astrometric accuracy. We examine proper motions and cloud survival in Sect. 3 and compare images with Pushchino single dish monitoring in Sect. 4, to investigate the survival of masers with respect to their parent clouds. We analyse the kinematics and estimate the mass concentration in maser clouds in Sect. 5. Local and global maser variability is discussed in Sect. 6. The size of the maser clumps and their relationships with stellar properties is discussed in Sect. 7 and Sect. 8 presents the summary and future work.

2. Data acquisition and analysis

Radio interferometry observations were made using 5 antennas of MERLIN, longest baseline 217 km, centred on the 22.235079 GHz H₂O maser line, corrected in the correlator to the V_LSR of the source. 2 GHz of bandwidth, providing 25 km s^-1 velocity coverage in 0.105 km s^-1 channels, was used for most epochs. This was found to be insufficient for IK Tau and S Per and so 4 GHz bandwidth, providing 50 km s^-1 in 0.21 km s^-1 channels, was used for these sources in 1999 and later epochs. VX Sgr was observed at both epochs in 4 GHz bandwidth and 0.42 km s^-1 channels were used for analysis. Data were recorded in two hands of circular polarization, calibrated independently (apart from amplitude normalisation). High spatial and spectral resolution is required to detect the weak polarization of H₂O (see e.g. Vlemmings et al. 2005 for results using VLBA data), and circular polarization is neglible for these sources at MERLIN resolution, so we made and analysed images in total intensity only.

The reduction of the 1994 data (and 1999 data for VX Sgr) is described in M03, R99 and B03. Similar techniques were used at later epochs, with the addition of phase referencing in some cases, see Sect. 2.2. All data reduction and initial analysis was performed using aips (Greisen 1994). Observational parameters and measurements are given in Table 2. The epochs of observation are superimposed (red lines) on AAVSO (American Association of Variable Star Observers) light curves in Figs. 1 to 3. The noise levels achieved depended on the source elevation, the weather and the total time on target. The weighting used during imaging was chosen to provide the best resolution without increasing the noise by more than 20%. In all cases, the dynamic range was limited to a few hundred, or about one hundred for low-declination sources with poor visibility-plane coverage.

The flux scale was derived from either 3C 84 or 3C 273, using 22-GHz monitoring data kindly provided by H. Teräsranta from the Metsähovi antenna. These sources may be slightly resolved by MERLIN at some epochs and vary during the few weeks or months between the Metsahovi and MERLIN observations. This dominates the uncertainty in the flux scales, which is about 10%.

In Appendix A we present contour maps of the maser emission averaged every few channels.

The position and peak intensity of each patch of maser emission brighter than 4.5–6σ_rms (depending on visibility plane coverage) in each channel was measured by fitting a 2-D Gaussian component. We rejected components fainter than the the local sidelobe level close to bright peaks. Their size (with the restoring beam deconvolved) and integrated intensity was also measured for all sources at positive declinations i.e. excluding VX Sgr and W Hya, for which the resolution N-S is too coarse. The accuracy of these measurements is proportional to (beam size)/(signal-to-noise ratio) (Condon 1997, adapted for sparse arrays by Condon et al. 1998; Richards et al. 1999) so the different sensitivies are taken account of via the uncertainties used in analysis.

The maser component positions are shown in Figs. 4 to 14. Maser components separated by a few milli-arcsec form series in adjacent channels which we term features. Our selection criteria imply that valid features must have a broader minimum total velocity span, if they were observed with broader channels. In the case of S Per and VX Sgr, the number of features remains similar despite changes in channel width, probably because the broader channels provided greater sensitivity. However, the sensitivity was reduced (due to worse weather and shorter observing times) for the 2000–2001 IK Tau epochs, despite the doubling of channel width compared with 1994. This may explain why significantly fewer features were detected at the later epochs.

The measured parameters for all sources are given in Tables B.1 to B.23. All but the faintest have Gaussian velocity profiles and R11 summarise the evidence that these are discrete, water vapour clouds. We rejected isolated components which did not form a series of at least 3 components in successive channels spanning the average feature size or less. The width of an individual component represents the beamed maser size, but R99 showed that the largest angular separation of components comprising a feature is at least 80% of the physical size of the cloud (as long as the emission is not resolved-out). Feature parameters are given in Tables B.1 to B.23.

Note that the individual components represent the line of sight through the cloud where masing is strongest in the velocity interval sampled by the channel. The arangement of series of components represents the velocity gradient with position but it does not necessarily show the full shape of the cloud. A linear or curved series of components can arise from a spherical cloud, especially if the masers are tightly beamed. This was discussed further in R11, who selected well-resolved components and features for beaming analysis whereas here we include all data. In addition, in aligning the epochs for U Her, we found a few misallocated components and revised a few feature properties, and we also improved the fitting of shell limits slightly for some sources/epochs. Hence, there are small discrepancies between their Table 2 and our Table 2. This does not affect any of the results in R11 except to reduce the scatter slightly in some of the analysis involving feature properties.

2.1. Finding the centre of expansion and aligning epochs

We initially assume that the winds are dominated by radial expansion from the central star. Yates & Cohen (1994) describe finding the centre of expansion by fitting a thick shell in 3 dimensions to angular radius and velocity. This gives well-constrained results for shells which are symmetric along any axis passing through the centre, although not necessarily spherical. This was used for objects with well-filled CSEs, i.e. all epochs of VX Sgr, S Per and RT Vir and early epochs of IK Tau and U Her. Many of the individual epochs of AGB stars are poorly filled and asymmetric, and for these we assumed that emission at close to V_⋆ should appear within a thick ring centred on the position which maximises the average angular separation from the centre of expansion. Bains et al. (2003) found that the two methods gave similar results for the 1994 observations of RT Vir, IK Tau and U Her. We used the second method for all epochs of U Ori and W Hya and for the later epochs of U Her and IK Tau.

We then aligned successive epochs of the same star using the initial estimates of the centres of expansion and looked for matching features, allowing for proper motions consistent with the spectral extent of the H₂O maser shell. Matches were found for those epochs listed in Table 2 as having N_prev features in common with the previous epoch. Components belonging to matched features are outlined in black in Figs. 4 to 14. The criteria used are described in more detail and the results are analysed further in Sect. 3.

Having identified potential matches, we refined the alignment by minimising the sum of the offsets of matched features (Richards et al. 1998) and if necessary applied additional shifts to the epoch with the greater uncertainties. These were in the range 1–5 milli-arcsec for all sources except W Hya, where the uncertainties are greater in declination and adjustments up to 10 milli-arcsec were needed. In the cases of U Ori and IK Tau, where phase-referenced positions were available at two epochs each, the adjustments were within the uncertainties of the astrometry (see below). The component positions in Figs. 4 to 14 are given relative to the matched centres of expansion for each star.

2.2. Phase referencing results

All positions are given in J2000 ICRF. Phase referencing was used at the epochs marked φ in Table 2. For U Ori and W Hya, we also used a more distant source to provide initial phase-rate and amplitude solutions, allowing fainter sources within 2° to be used for phase solutions. In all cases, the maser peaks were bright enough to give better signal-to-noise in a single spectral channel, than were the phase-reference sources observed in the maximum available 13-MHz bandwiths. There are 4 factors which contribute to the position uncertainty:

• 1.

  The noise-based error in measuring the position by fitting aGaussian component, negligible ( ≲ 1 mas) for the brightest masers and for the phase reference sources unless otherwise stated.

• 2.

  The uncertainty in transfering the phase solutions between reference source and maser. This is estimated from the sky separation in minutes of time, e.g. 2° = 8 min. We examine the phase drift in the equivalent time on the raw data for a bright source observed continuously for many minutes, usually the bandpass calibrator. A typical phase drift of half a turn corresponds to half a beam position uncertainty. This assumes that the atmospheric effects are the same throughout the observations; if the weather or the elevation changed significantly we added an allowance for this.

• 3.

  The phase reference position uncertainty. The target positions have been adjusted to the most accurate positions available.

• 4.

  Uncertainty in the telescope positions. This is about 1 cm,  ~1 wavelength at 22 GHz, on the longest MERLIN baseline, giving an uncertainty of about one beam (see Table 2).

The naturally-weighted beam size must be used for these estimates. Only 1. affects data for different components in a single source at a single epoch. 3. and 4. do not affect data taken using the same array and the same phase-reference at different epochs. The uncertainties for each target and phase reference source are given in the paragraphs following. Where possible, we initially transfered the phase reference solutions to the maser and then repeated the process in the opposite direction. This provides two independent estimates of the position of maser emission in the reference channel. The discrepancy is used to refine the estimate of error 2. In other cases, the phase reference source was too weak to detect until after applying the maser solutions.

We used Hipparcos positions and proper motions (van Leeuwen 2007), to predict the stellar positions of U Her, U Ori and W Hya at our epochs of observation. The formal errors, derived from the uncertainties in position at epoch 1991.25 and in the proper motions, are 10–20 milli-arcsec. The optical position of IK Tau is less well known due to its highly obscured location, having an uncertainty of about 100 milli-arcsec in either coordinate (e.g. Zacharias et al. 2005). We compare these positions with our estimates of the centres of expansion derived from the masers, assumed to be the stellar positions.

U Ori

(Fig. 7). The phase reference source J0552+1913 at 05:52:25.8858 +19:13:40.275, uncertainty (10, 10) milli-arcsec (Browne et al. 1998), just over 1° from U Ori, was used in 2000 and 2001. The brightest maser channel was at  − 38 km s^-1, at both epochs. This is close to V_⋆ and thus the maser feature will have a significant proper motion in the plane of the sky due to the expansion of the CSE, as well as to stellar proper motion. We therefore compared the positions of the centre of expansion. These were 05:55:49.1699 +20:10:30.624 and 05:55:49.1697 +20:10:30.619 in 2000 and 2001, respectively, fitting uncertainty (5, 5) milli-arcsec, total astrometric uncertainty (23, 32) milli-arcsec per epoch. This corresponds to a proper motion of (–3, –5) milli-arcsec, which, although close to the fitting uncertainties, is in a direction similar to the Hipparcos proper motion of ( − 12.12,  − 6.13), uncertainty (1.52, 0.73) milli-arcsec yr^-1. The positions of U Ori derived from Hipparcos were 05:55:49.1697 +20:10:30.685, uncertainty (14, 7) milli-arcsec, and 05:55:49.1688 +10:30.679 (15, 7) milli-arcsec, in 2000 and 2001, respectively. The positions derived from the masers are offset by (3, –61) and (13, –60) milli-arcsec from the estimated Hipparcos positions. The combined uncertainties are (27, 35) milli-arcsec. The good agreement between the expansion centres estimated independently from the maser emission at the two epochs suggests that our alignment is accurate. If the star really was at the location calculated using the Hipparcos data, it would be at the outer edge of the maser shell (opposite to the maser peak), which suggests that there is an unrecognised source of astrometric error of order tens of milli-arcsec in one or both sets of data.

U Her

(Fig. 9). The phase reference source J1635+1831 at 16:35:39.1588 +18:31:03.723, uncertainty (7, 12) milli-arcsec (Browne et al. 1998) at 2$\stackrel{\mathrm{°}}{.}$4 separation from U Her, was used in 2001. We obtained a position for the centre of expansion of U Her of 16:25:47.4710 +18:53:32.887, with a fitting uncertainty of (4, 4) milli-arcsec and a total astrometric uncertainty of (20, 32) milli-arcsec. The Hipparcos position at this epoch is 16:25:47.4702 18:53:32.841, uncertainty (6, 7) milli-arcsec. The position derived from masers is offset by (11, 46) milli-arcsec, close to the combined uncertainties. The brightest U Her maser feature observed in 2001 is at  − 15.88 km s^-1, position 16:25:47.4668 +18:53:32.879, total extent 29 milli-arcsec, 2.1 km s^-1. Vlemmings et al. (2002) carried out MERLIN observations less than 4 weeks later, optimised for astrometry by using 3 phase reference sources and a velocity resolution of 0.42 km s^-1. They found the brightest feature at  − 15.7 km s^-1, position 16:25:47.468 +18:53:32.849 (10, 10). Our position is offset from this by (–17, 30) milli-arcsec, within the combined uncertainties for observations with the same array but different reference sources.

IK Tau

(Fig. 11). The phase reference source J0409+1217 at 04:09:22.0087 +12:17:39.845, uncertainty (0.2, 0.2) milli-arcsec (Titov 2004), 4° separation, was used in 2000 and 2001. The bright, blue-shifted feature used for self-calibration is within 1 milli-arcsec of the centre of expansion at 03:53:28.8919 +11:24:21.940, and 03:53:28.8924 +11:24:21.932 in 2000 and 2001, respectively, fitting uncertainty (3, 3) milli-arcsec, total astrometric uncertainties (30, 42) milli-arcsec and (33, 43) milli-arcsec. The 2001 position is offset from the 2000 position by (5, –8) milli-arcsec. This is less than expected from the optically measured proper motions, of (23.5, –26.1) or (8, –80) milli-arcsec yr^-1 from NOMAD (Zacharias et al. 2005) or the PM2000 Bordeaux Proper Motion catalogue (Ducourant et al. 2006), respectively. The 2001 maser-derived position is (137, –86) milli-arcsec from the Bordeaux position at the same epoch, 03:53:28.8825 +11:24:22.018, uncertainty (27, 27) milli-arcsec. The NOMAD position differs from both by a total of 0$\stackrel{\mathrm{″}}{.}$75 (0$\stackrel{\mathrm{″}}{.}$15) (with greater discrepancies at the earlier epoch).

W Hya

(Fig. 14). The phase reference source J1342-2900 at 13:42:15.345598 –29:00:41.83163, uncertainty (0.3, 0.8) milli-arcsec (Petrov et al. 2011), under 2° separation, was used at all epochs. W Hya has a high proper motion, (149, –60) milli-arcsec and the phase-referenced maser positions were offset in the expected direction from the pointing position (by 2–3 arcsec) but the positions at each of the 4 epochs differed by several tenths of an arcsec or more from the positions predicted from Hipparcos data, in a non-systematic fashion. Phase transfer at such low elevation was not accurate enough to give astrometrically useful results.

In conclusion, the phase referencing results for U Ori and IK Tau vindicate the methods used to align the epochs. However, the uncertainties in our measurements and the Hipparcos data are too great for sufficiently accurate comparison to confirm that the star is located at the centre of expansion. Vlemmings et al. (2002) suggest that the brightest feature in U Her is in the line of sight to the star itself. Our model places it at a separation of about 40–60 milli-arcsec from the centre of expansion. Unfortunately, the uncertainties are just too great to distinguish between these models astrometrically. Such tests will have to wait for imaging of the radio star and masers simultaneously with e-MERLIN, as was done for W Hya by Reid & Menten (1990) using the VLA.

2.3. Maser shells

Figures 15 to 21 show the maser separations from the centre of expansion projected against the plane of the sky, $\mathit{a}\mathrm{=}\sqrt{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}}$, as a function of V_LSR, centred at V_⋆. We fitted ellipses by eye to the inner and outer limits of the 22-GHz maser emission, giving the inner and outer radii and expansion velocities of spherical shells, r_i, r_o, v_i and v_o. The velocity gradient is given by ${\mathit{K}}_{\mathrm{grad}}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{o}}\mathrm{-}{\mathit{v}}_{\mathrm{i}}}{{\mathit{r}}_{\mathrm{o}}\mathrm{-}{\mathit{r}}_{\mathrm{i}}}$(1)and the logarithmic velocity gradient by $\mathit{ϵ}\mathrm{=}\frac{\log \mathrm{(}{\mathit{v}}_{\mathrm{o}}\mathit{/}{\mathit{v}}_{\mathrm{i}}\mathrm{)}}{\log \mathrm{(}{\mathit{r}}_{\mathrm{o}}\mathit{/}{\mathit{r}}_{\mathrm{i}}\mathrm{)}}\mathrm{·}$(2)These parameters are given in Table 2 and used to investigate the mass loss as a whole in Sect. 5.1.

3. Feature persistance and proper motions

We matched features at successive epochs as described in Sect. 2.1. The numbers matched with the previous epoch are given in Table 2. We required the mean feature V_LSR velocities to coincide to within a thermal line width ΔV_th ≈ 1.4 km s^-1 for H₂O in the maser shell. We allowed potential feature matches between epochs if the corresponding proper motions did not exceed the maximum expansion velocity  | V_LSR − V_⋆ | , plus an allowance for position inaccuracy and feature size. For reference, an expansion velocity of 5 km s^-1 corresponds to a proper motion of  ≈1 AU yr^-1 (4–10 milli-arcsec yr^-1 for the AGB stars at 266–99 pc, and 2–3 milli-arcsec in 5 yr for the RSG at 2.3–1.7 pc, using the stellar distances given in Table 1). We tested the the angular and velocity search regions by decreasing or increasing them by 50%. We also tested whether matches were random coincidences by offsetting one epoch by the match radius in x and y before looking for matches. Since one of the goals of this project was to investigate the CSE kinematics, we did not want to impose conditions on the directions of proper motions.

We found 5 or more matched pairs of features between the RSG epochs 5 years apart and between most AGB epochs separated by 1 year or less (up to 3/4 of all features in successive epochs). The number of closest matches are listed in Table 2. Changing the match radius led to fewer matches or more multiple matches but few if any new unique matches. Applying shifts reduced the number of matches to 25% or fewer, suggesting that the majority of the listed matches are correct.

We only analyse reliable matches, using the above criteria and requiring the majority of matches to be unique. In some cases we found possible matches over longer AGB epoch separations but, as the potential proper motion and hence search radius increased, there were more multiple matches and as many matches were found after applying a shift as before, so these cannot be used.

More details are given for each source:

VX Sgr.

We used a maximum allowed separation of 23 milli-arcsec in 5 yr, and about half of all features were matched.

S Per.

We used a maximum allowed separation of 13 milli-arcsec in 5 yr, and about half of all features were matched.

U Ori.

We used a maximum allowed separation of 12 milli-arcsec in 1 yr, or 30 milli-arcsec in 5 yr. We only found reliable matches between 2000–2001, for about 1/3 of the features.

U Her.

We used a maximum allowed separation of 12 milli-arcsec in 1 yr, or 45 milli-arcsec in 6 yr. About 1/2 of the features were matched from 2000–2001. Some possible matches were seen between 1994–2000 but none were unique and a similar number of matches were seen after applying the test shifts. After comparing the evolution of the H₂O maser shells as a whole in Sect. 6.1, we found that the ring of emission seen in U Her in 1994 could have expanded to the location of the outer ring seen at later epochs, suggesting the persistence of a large-scale structure.

IK Tau.

We used a maximum allowed separation of 14 milli-arcsec in 1 yr, or 50 milli-arcsec in 6 yr. About 1/2 of the features were matched from 2000–2001. Some possible matches were seen between 1994–2000 but none were unique and a similar number of matches were seen after applying the test shifts. No obvious patterns appear to survive.

RT Vir.

We used a maximum allowed separation of 9 milli-arcsec over 10 weeks in 1996, and 25 milli-arcsec for 1994–1996. Between 1/2 and 3/4 of features were traced from epoch to epoch in 1996, with eleven (about 1/4) being identified at all six epochs. No matches were found between 1994–1996.

W Hya.

We used a maximum allowed separation of 20 milli-arcsec in 1 yr; 1/3 to 1/2 of the features could be matched from one year to the next, including six at three successive epochs. No features were matched for all four epochs.

3.1. Identification of features

We are confident that in most of the matched features, we are seeing masers from the same cloud, for two main reasons. Firstly, the 22-GHz H₂O masers clouds are much denser than the average wind density (R99, B03, M03, Sect. 5.2). Secondly. the more extended features (better-resolved in the higher-declination sources) often show a comparable internal structure. This can be seen for some, although not all of the matched features in Figs. 5, 9, 11, and 13. We zoom in on two examples from S Per and one from RT Vir.

Figure 5 shows a series of features labelled A, around (30,35) milli-arcsec, peaking at  − 46 –  − 47 km s^-1, which appears to have been rotated counterclockwise by about 60° between epochs (as well as bulk motion to the S). Another pair of features labelled B, at around ( − 37,  − 25) milli-arcsec,  − 30 km s^-1 are made up of a radial spoke of components. seen more clearly in the enlargement in Fig. 22.

What effects produce these variations within features? Bulk rotation could be responsible for the changes shown in Fig. 22 (left). This consists of a series of 4 matched features (separate spectral peaks) whose components form a spatially and spectrally continuous series, taken to be a single, intrinsically asymmetric (amœboid) cloud, with internal density variations. After removing the mean proper motion of the cloud (3 mas), it is rotating anticlockwise at  ≲6 ± 4 km s^-1. The feature in Fig. 22 (right) shows an almost radial gradient of velocity with position, with two slight wiggles seen at both epochs (with greater position uncertainty in 1999 as it has become fainter). The proper motion corresponds to an expansion velocity  ≈10 km s^-1.

One of the matched RT Vir features is shown in detail in Fig. 23. The bright head of the “tadpole” has a similar appearance at all epochs, but it has two faint tails; the more redshifted one is seen in 960 405, 960 429 and (a short piece) 960 612, and the less redshifted (orange-yellow) tail in 960 421, 960 429, 960 515, 969 524 and 960 612. Turbulence inside the feature could affect the beaming intensity and/or direction, for example causing the higher-velocity (redder) tail to fade and to be replaced by a lower-velocity tail. Note that the most sensitive observations were made in 960 429 and 960 524 (Table 2), hence the remnant of the red tail and the longest extension of the new tail seen at these epochs, respectively. The tails have extents  ≳ 1 AU and a separation of about 0.5 AU. If the change is due to a single disturbance, this must be broad enough to affect both regions almost simultaneously, since a localised, non-dissociative shock could not traverse the distance in the short time intervals of a few weeks between observations. Alternatively, the entire cloud has rotated about an axis close to the plane of the sky such that beaming from different angles was favoured.

If the clouds are approximately spherical and the masers are amplification-bounded, the brightest emission at the line centre is tightly beamed along the line of greatest amplification. Small changes in the cloud will have little effect on the apparent location of this emission although they may affect the appearance of weaker, less strongly beamed emission from the line wings. On the other hand, flattened or shocked clouds may possess matter-bounded masers which are less tightly beamed. R11 found that U Ori, U Her and to a lesser extent IK Tau possess some matter-bounded masers, whilst the other sources studied are mostly amplification-bounded.

Overlapping clouds along the line of sight can produce strong amplification of a single feature, see Sect. 6.2. The survival of maser features is discussed in Sect. 4.2.

3.2. Expansion proper motions

We aligned the results from successive epochs for each source and matched features as outlined in Sect. 2.1. Table 3 gives the error-weighted proper motions; for multi-epoch matches, we performed least squares fits to position as a function of time. The mean total proper motion (and formal uncertainty) is  ⟨ xẏ ⟩ . We independently estimated mean proper motions in the radial and tangential directions,  ⟨ ȧ ⟩  and .

In a few cases the direction of  ⟨ xẏ ⟩  has a large scatter and and  ⟨ ȧ ⟩  are consequently small and of low significance. This is most noticeable at the earlier, noisiest epochs of W Hya, probably due to measurement errors. Since its values of change sign from epoch to epoch, they are probably not significant. The U Ori proper motions also show more random behaviour. Coupled with the failure to match features in the one-year interval between 1999–2000, this suggests that the U Ori CSE is particularly subject to turbulence or other disruption. Many of its H₂O masers appear matter-bounded, hence less tightly beamed, so that the apparent centroid may be shifted easily by small disturbances (see Sect. 3.1 and R11).

In all other cases the radial motions  ⟨ ȧ ⟩  are in expansion, and exceed . This is best illustrated by Fig. 24 which shows for the 11 RT Vir features matched at all 6 epochs. In fact, is negligible for most sources. IK Tau maser proper motions are clearly dominated by expansion but a small rotational component cannot be ruled out, although this is  < ΔV_th. This is also the case for U Ori although the majority of the proper motion measurements are in scattered directions, as mentioned above.

Figures 25 to 31 show expansion proper motion velocity V_a as a function of mean angular separation from the centre of expansion, $\overline{\mathit{a}}$, for the data in Table 3. We performed a linear error-weighted fit to find the slope ($\mathrm{d}{\mathit{V}}_{\mathrm{a}}\mathit{/}\mathrm{d}\overline{\mathit{a}}$) and found the dispersion in the relationship (the formal errors are very small). This gives an estimate of wind acceleration in the plane of the sky. The RSG, VX Sgr and S Per have $\mathrm{d}{\mathit{V}}_{\mathrm{a}}\mathit{/}\mathrm{d}\overline{\mathit{a}}$ of 0.07 ± 0.02 and 0.11 ± 0.16 km s^-1 AU^-1, respectively. The Miras U Her and IK Tau have $\mathrm{d}{\mathit{V}}_{\mathrm{a}}\mathit{/}\mathrm{d}\overline{\mathit{a}}$ of 0.31 ± 0.28 and 0.31 ± 0.18 km s^-1 AU^-1, respectively. The SRb RT Vir has $\mathrm{d}{\mathit{V}}_{\mathrm{a}}\mathit{/}\mathrm{d}\overline{\mathit{a}}\mathrm{=}\mathrm{0.64}\mathrm{\pm }\mathrm{0.36}$ km s^-1 AU^-1 for the 11 features matched at all epochs; if 2–5-epoch matches are considered, the scatter is greater but the slope is also positive. W Hya has $\mathrm{d}{\mathit{V}}_{\mathrm{a}}\mathit{/}\mathrm{d}\overline{\mathit{a}}\mathrm{=}\mathrm{0.39}\mathrm{\pm }\mathrm{0.35}$ km s^-1 AU^-1 for the 3-epoch matches and additional 2-epoch matches in 2001–2002. The only exceptions to the positive gradients are U Ori, $\mathrm{d}{\mathit{V}}_{\mathrm{a}}\mathit{/}\mathrm{d}\overline{\mathit{a}}\mathrm{=}\mathrm{-}\mathrm{0.2}\mathrm{\pm }\mathrm{0.4}$ km s^-1 AU^-1, and W Hya 1999-2000,  − 0.1 ± 0.1 km s^-1 AU^-1, probably due to a more turbulent maser shell and to greater uncertainties, respectively, as discussed above.

For VX Sgr, S Per, IK Tau, U Her, RT Vir and the better W Hya data, $\mathrm{d}{\mathit{V}}_{\mathrm{a}}\mathit{/}\mathrm{d}\overline{\mathit{a}}$ is positive. It is close to K_grad (Table 2), within the dispersions, for all measurements apart from VX Sgr, where $\mathrm{d}{\mathit{V}}_{\mathrm{a}}\mathit{/}\mathrm{d}\overline{\mathit{a}}$ is slightly greater; this may be due to the axisymmetry of the outflow, modelled by M03. In the other objects the velocity field appears to be approximately spherically symmetric.

These results are consistent with previous proper motion measurements. Yates & Cohen (1994) identified IK Tau features (at lower resolution) surviving for 16 months. Marvel et al. (1998) observed several objects including VX Sgr, S Per, U Her and IK Tau for three epochs during 1995. They deduced that the masers are distributed in multiple thin shells or random clumps. They found that an ellipsoidal model fitted most data better than a spherically symmetric shell, but comparison between their figures and our 1994 plots suggests that fainter, extended emission was resolved out by the VLBA. The proper motions are predominantly directed radially outward, consistent with the LSR maser velocities at distances similar to those given in Table 1. Asaki et al. (2010) performed astrometric monitoring of S Per with VERA to obtain the parallax. The maser proper motions were consistent with a combination of bulk proper motion and a spherically expanding flow. Imai et al. (2003) monitored RT Vir using the VLBA for 5 epochs in 1998. The tendency for the E and W to be dominated by red- and blue-shifted emission is seen in their results as well as in our Figs. 12 and 13 although the brightest masers occur in different locations. They also find expansion (not rotation) proper motions, with at least one feature showing strong acceleration. Note that they use a distance of 220 pc and so obtain proper motion velocities almost double those resulting from the distance we adopt (Table 1).

4. Comparison with single dish data

These targets are among a large number of 22-GHz maser sources which have been monitored by the Pushchino Radio Telescope for several decades. Some results have already been published in Pashchenko & Rudnitskii (1999) (VX Sgr); Lekht et al. (2005) (S Per); Rudnitskij et al. (2000) (U Ori); Mendoza-Torres et al. (1997) and Lekht et al. (1999) (RT Vir). Earlier epochs of W Hya monitoring were published by Rudnitskii et al. (1999).

The sensitivity limit is  ~10 Jy (σ_rms 2–4 Jy) and the spectral resolution is 0.101 km s^-1 (close to the MERLIN resolution of 0.105 km s^-1 at most epochs). We show the MERLIN velocity profiles and the Pushchino spectra taken closest in time, along with a few intermediate epochs, in Figs. 32 to 43. The epochs shown here are superimposed (blue lines) on AAVSO light curves in Figs. 1 to 3. The Pushchino absolute velocities are uncertain at some epochs due to local oscillator problems so we applied a velocity shift V_corr, given in the figure captions, to align the Pushchino spectral peaks with the MERLIN peaks. This correction is only accurate to within a few tenths km s^-1 as some the spectral features show drifts of this magnitude when comparing different epochs taken with the same instrument.

The dotted black lines in Figs. 32 to 43 show the integrated MERLIN velocity profiles. The circles show the individual maser components, outlined in black for features matched at successive epochs. The coloured lines show the Gaussian profiles fitted with  >3σ accuracy to individual features. Comparison with the Puschino data shows that the maser flux densities vary greatly on timescales of a few weeks. Trends in spectral variability and velocity drifts are analysed in detail in the references given at the start of Sect. 4. Here, we consider two questions: does MERLIN recover all the flux, and what does spectral monitoring reveal about cloud survival when compared with imaging?

4.1. Flux detected by MERLIN interferometry

The shortest MERLIN baseline of 20 km means that emission which is smooth on scales  > 100 milli-arcsec cannot be imaged reliably. We tabulate the peak flux densities at the closest epochs in Table 4. The irregular nature of maser variability means that we cannot interpolate reliably between Pushchino measurements bracketing a MERLIN measurement. In some cases, a new feature emerges or disappears in observations made a few weeks apart. We can only investigate whether MERLIN consistently records lower peaks or less integrated flux (I_tot) than Puschino.

RSG. For VX Sgr, the MERLIN:Pushchino (M:P) peak ratio is between 0.42–0.66. However, the MERLIN spectral resolution is much coarser, at 0.42 km s^-1, which can dilute the flux density of narrow peaks. The original 0.21 km s^-1 resolution used in 1994 gave a peak of 375 Jy (Richards 1997) and the integrated flux density at this epoch is very similar for both instruments. For S Per, in 1994 M:P is between 2.8–0.65 and the integrated flux densities are very similar for most adjacent epochs.

Miras. For U Ori, M:P is in the range 0.42–0.89. However, the source is relatively weak and the differences are only significant in 1999. The MERLIN integrated flux densities are slightly closer to the Pushchino values, exceeding it at one epoch. For U Her, M:P lies in the range 0.56–1.1 and for IK Tau M:P is between 0.7–1.7, and the integrated flux densities can be higher or lower at adjacent epochs.

SR. For RT Vir, in 1994 M:P 0.81–2.1 (in the latter case the gap is longer than usual, 15 weeks). In 1996, six MERLIN observations were made over 10 weeks during which the peak flux rose from 389 to 792 Jy but M:P changed from 0.49 to 2.2. The integrated flux densities behave similarly. For W Hya, M:P is in the range 0.30–1.5, but the Pushchino peak and integrated flux densities exceed the MERLIN values for all but the final epoch.

MERLIN clearly recovers most of the flux (within the flux scale uncertainty of 10% for MERLIN and the 10 Jy accuracy limit for Pushchino) for S Per, U Her, IK Tau and RT Vir. This also seems likely for VX Sgr, taking the broader channel width into account, and for U Ori, which is at a very similar distance as U Her and IK Tau, so it is unlikely to be more resolved, but it cannot be ruled out that some of its flux has been missed. It is possible that more than 10% of the emission from W Hya has been resolved out, since it is the closest source, although its low elevation during observations makes the flux scale even less accurate. Therefore, we regard measurements of cloud size or flux density to be lower limits for W Hya.

4.2. Maser features and their parent clouds

Section 2.1 describes how we were able to match some maser features in successive MERLIN epochs (numbers listed in Table 2). Many are spatially distinct in the MERLIN images but spectrally blended in Pushchino observations; however, some of the brightest have unambiguous counterparts. We investigate the spectral variability of these features in the more closely time-sampled Pushchino spectra.

VX Sgr.

M03 described how about half the VX Sgr maser features detected in 1994 were also identified in 1999. The brightest spectral peaks are present in the intervening epochs (Pashchenko & Rudnitskii 1999).

S Per.

Forty S Per MERLIN features were unambiguously matched at both MERLIN epochs. Three peaks dominate the 1994 MERLIN spectrum and the Pushchino spectra at the closest epochs, see Fig. 33. The blue-shifted peak at  − 46 –  − 47 km s^-1 is seen at all Pushchino epochs (an additional, more blue-shifted feature appears after 1997). The red-shifted peak becomes fainter after 1994 and in the 1999 MERLIN spectrum it is under the 10 Jy Pushchino limit. These are seen with a similar internal structure by MERLIN at both 1994 and 1999 (Fig. 22).

The central peak around  − 37.8 km s^-1 exceeds 40 Jy and is spatially compact at both MERLIN epochs. However, it has faded in Pushchino spectra taken at 19940810, a few months after the 1994 MERLIN epoch, and is absent for most of 1995-6 and all of 1997 (Lekht et al. 2005). It re-appears in the 19981228 and later Pushchino spectra as well as in the MERLIN 1999 data.

U Ori.

Nine U Ori MERLIN features were matched at 2000 and 2001. The brightest emission, just under 20 Jy, is close to  − 38 km s^-1 (Fig. 35). This is dominated by the feature at ( − 40,  − 15) milli-arcsec (Fig. 7). The Puschino spectrum shows a corresponding peak of 20–40 Jy at the closest epochs. However, no emission above 10 Jy was detected from U Ori in intermediate spectra, e.g. 20000725.

U Her.

Twenty U Her MERLIN features were matched at 2000 and 2001. The peak, around  − 15.5 km s^-1, was seen at all Pushchino epochs between these dates.

IK Tau.

Twenty-four IK Tau MERLIN features were matched at 2000 and 2001. In 2000, the 75-Jy peak at 25 km s^-1 corresponds to the feature at (–50, 10) milli-arcsec and the 25-Jy peak at 22 km s^-1 corresponds to the feature at (0, 0) milli-arcsec (Fig. 11). Figure 39 shows that the 25 km s^-1 peak declines abruptly after 20000530. The IK Tau H₂O maser shell is well-filled, with several features contributing to the spectral peaks. It is possible that the 25 km s^-1 peak was a flare due to clouds overlapping along the line of sight, similar to that described for W Hya in Sect. 6.2. The 22 km s^-1 feature is present in the 2001 MERLIN spectrum at a similar flux density, 21 Jy, to that seen in 1994. The blue-shifted peaks cannot be distinguished in the 200000819 and 20010409 spectra but are present in 20010529.

RT Vir.

Many features were matched between six epochs of RT Vir MERLIN data taken over 10 weeks in 1996 but there are no intervening Pushchino spectra. Figure 41 shows that the strong blue- and red-shifted RT Vir peaks show systematic increases in  | V_LSR − V_⋆ |  with epoch, totalling 0.43 and 0.17 km s^-1 in 10 weeks, respectively, consistent with the strong acceleration seen in proper motions (Sect. 3; see also Imai et al. 2003).

W Hya.

Between five and nine W Hya MERLIN features were matched between successive epochs in 1999, 2000, 2001 and 2002. The brightest part of the MERLIN spectrum in 1999 is made up of a number of features so their individual behaviour is hard to trace in the Pushchino spectra. However, the 30-Jy feature at 41.1 km s^-1 seen in 1999, located at ( − 95, 120) has a much fainter counterpart in 2000 (Figs. 42 and 14). Faint emission was also detected at the same velocity in 2001 and rises to 160 Jy in 2002, but this is located on the opposite side of the CSE at ( − 70,  − 110) (Figs. 43 and 14). The 2001–2002 flare is discussed in Sect. 6.2.

Figure 1 to 3 show that the intervals over which masers can be matched spans several stellar periods in the case of VX Sgr, S Per and W Hya, and nearly two periods for U Her and IK Tau. U Ori features were matched across one interval close to the stellar period, but not in another, whilst RT Vir features were only matched over intervals less than one stellar period. There is no obvious correlation between maser survival and the number of periods or the phase. In other interferometric studies, Winnberg et al. (2011) found that U Her maser components have an average lifetime of 0.5–1 yr and Winnberg et al. (2008) obtained a lifetime  ~ 1 yr for RX Boo. In general, single dish observations are affected by spatial blending and may overestimate maser lifetimes, but we note that Lekht et al. (1999) analysed 12 yr of Pushchino monitoring of RT Vir and found 22 spectral features persisting for 0.66 yr or more, only 5 of which lasted longer than 2 yr.

The time t taken for the wind to to reach a given distance r in the 22-GHz maser shell is calculated using the velocity gradient K_grad (Table 2), assuming linear, radial acceleration and good coupling of dust and gas velocities for the H₂O maser clouds: $\frac{\mathit{t}}{\mathrm{yr}}\mathrm{=}\mathrm{4.8}{\left(\frac{{\mathit{K}}_{\mathrm{grad}}}{\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{/}\mathrm{AU}}\right)}^{-1}\ln \left[\mathrm{1}\mathrm{+}\frac{\mathrm{(}\mathit{r}\mathrm{-}{\mathit{r}}_{\mathrm{i}}\mathrm{)}{\mathit{K}}_{\mathrm{grad}}}{{\mathit{v}}_{\mathrm{i}}}\right]$(3)and the total shell crossing time from r_i to r_o, is t_cross, given in Table 5. This is about 10–15 yr for SRbs, 20–35 yr for Miras and 35–90 yr for the RSGs. This is the minimum survival time for the clouds as discrete, dense entities, since it is hard to see how they could re-form if disrupted during passage across the water maser shell. On the other hand, most individual maser features become unrecognisable after 2 years or less around AGB stars, and only about half of the RSG features are matched after 5 years, although M03 found that about one fifth survive 16 yr; VY CMa shows similar behaviour (Bowers et al. 1993; Richards et al. 1998).

The comparison with Pushchino data shows that in fact, masers from individual clouds switch off and on again. This is most obvious for S Per, U Ori and IK Tau, where spectral peaks can be clearly identified with spatial features at separate interferometry epochs, but those peaks are missing in the spectra at some intervening single dish monitoring epochs. The smaller the clouds (Table 2), the more short-lived the masers. The sound-crossing times (assuming a sound speed of 3 km s^-1) range from 1.5–5 yr for the AGB clouds and 20–25 yr for the RSG. The maser lifetimes are comparable to the time which would be taken for sub-sonic or mildly supersonic shocks to cross the maser clouds. Figure 22 shows a localised, non-radial disturbance propagating at  ≲ 5 km s^-1. These might temporarily disrupt the maser altogether, or simply cause it to beam in a direction we cannot detect.

5. Acceleration and density of maser clouds

The inner and outer maser shell limit measurements reveal the kinematics (Sect. 5.1). We then use the individual cloud sizes to estimate their contribution to the mass loss (Sect. 5.2).

5.1. Acceleration

Figure 44 shows the shell limits for all sources, taken from Sect. 2.3, Figs. 15 to 21 and Table 2. This demonstrates that the trend for higher velocities at greater distances from the star is similar when comparing different objects which differ in size by an order of magnitude. In all objects, the maser velocity approximately doubles as the wind traverses the maser shells, passing through the gravitational escape velocity V_esc in the process. This implies ongoing acceleration beyond the dust formation zone. Ivezić & Elitzur (2010) model the CSE velocity field using changes in optical depth. However, the consistently strong acceleration over such a wide range of objects seen here is likely to require additional force due to changes in the surface properties of the dust as grains anneal in the wind. This was first suggested by Chapman & Cohen (1986) for VX Sgr and is discussed for others among our maser sources by Richards & Yates (1998). Verhoelst et al. (2009) observed correlations between dust composition and outflow velocities for RSG while Decin et al. (2010) modelled Herschel observerations of IK Tau using gradual acceleration.

Figure 45 shows that the SRs, Miras and RSGs form three distinct groups of increasing 22-GHz shell thickness and this is weakly correlated with shallower velocity gradients, K_grad, which measure the increase in expansion velocity with distance from the star. K_grad is estimated from the V_LSR measurements along the line of sight. The results in Sect. 3.2 demonstrate a similar trend for proper motions in the plane of the sky. The SRb RT Vir has the highest gradient of proper motion velocity with position in the shell, $\mathrm{d}{\mathit{V}}_{\mathrm{a}}\mathit{/}\mathrm{d}\overline{\mathit{a}}$, and the other AGB stars with significant results have higher values of $\mathrm{d}{\mathit{V}}_{\mathrm{a}}\mathit{/}\mathrm{d}\overline{\mathit{a}}$ than the RSG. The thicker circumstellar envelopes show more gradual acceleration, but eventually reach higher velocities.

A theoretical relationship between Ṁ and terminal velocity was established by Cooke & Elitzur (1985), who also showed that Ṁ should determine the 22-GHz shell r_i. The general forms of these relationships (and the acceleration through V_esc) were confirmed observationally by Lane et al. (1987), Yates & Cohen (1994), Bowers & Johnston (1994) and Colomer et al. (2000). However, the observed r_i tend to be much larger than predicted. This implies that the wind is inhomogenous, discussed in Sect. 5.2.

5.2. Water maser clump density and the mass in clouds

We use the relationship between Ṁ and r_i to estimate the density of the H₂O maser clouds. The density of the stellar wind at radius r is ρ = Ṁ/(4πr²v), assuming radial acceleration of a smooth medium. If the average particle mass is 1.37 ×  the mass of molecular hydrogen, ρ = 1.37 × 2m_Hn. The overall number density n of the wind is thus given by: $\frac{\mathit{n}}{{\mathrm{m}}^{-3}}\mathrm{\approx }\mathrm{5}\mathrm{\times }{\mathrm{10}}^{\mathrm{22}}\left(\frac{\mathit{Ṁ}}{{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}{\mathrm{yr}}^{-1}}\right){\left(\frac{\mathit{r}}{\mathrm{AU}}\right)}^{-2}{\left(\frac{\mathit{v}}{{\mathrm{km}\hspace{0.17em}\mathrm{s}}^{-1}}\right)}^{-1}\mathrm{·}$(4)We use this to find the average number density n_i at r_i, where v = v_i. The inner edge of the 22-GHz maser shell is determined by the point at which the collision rate falls below the maser emission rate, at the quenching number density n_q. At 1200 K, n_q ≈ 5 × 10¹⁵ m^-3 (Yates et al. 1997). We assume that this is the density of the clouds at r_i.

The values of Ṁ given in Table 1 are taken from the literature (adjusted to our adopted distances), using high signal-to-noise CO data, supported by interferometry, where available. However, some objects, notably S Per, have faint CO. Other methods of estimating Ṁ, such as IR observations of dust, are more model-dependent; the assumed gas:dust ratios vary between about 80 and several hundred. The variations in values of Ṁ in the literature are a factor of about 5 for most of the objects.

We want to estimate the total number of potentially masing clouds in a shell at any one time, N_tot. However, as described in Sect. 4.2, at any one epoch, some clouds may not be producing masers beamed towards Earth. We attempted to compensate for this using the multi-epoch observations where we should have found reliable matches at neighbouring epochs, i.e. 6 epochs for RT Vir, 4 for W Hya, 3 for U Ori and 2 for the other stars (Sect. 3). This excludes the 1994 data for the AGB stars. N_tot is given by the sum of the number of features NF at the first epoch (e = 1), plus the number of maser features with no previous matches at subsequent epochs (NF − N_prev at 1 < e ≤  end) where e = end at the final epoch. We finally subtract the number of clouds likely to have left the shell during the time interval for all epochs considered, t_all; since t_all ≪ t_cross this is approximated by the uncorrected value of N_tot multiplied by (1 − t_all/t_cross). This approach partially compensates for less sensitive epochs (since a cloud only has to be detected once to be counted) but N_tot is still a lower limit if some clouds are present in the CSE but not detected at any epoch. The degree of underestimation may vary from star to star if some have a persistently lower fraction of detectably masing clouds. ${\mathit{N}}_{\mathrm{tot}}\mathrm{=}\left[\mathit{NF}\mathrm{+}\sum _{\mathit{e}\mathrm{=}\mathrm{2}}^{\mathit{e}\mathrm{=}\mathrm{end}}\mathrm{(}\mathit{N}{\mathit{F}}_{\mathrm{\left(}\mathit{e}\mathrm{\right)}}\mathrm{-}{\mathit{N}}_{\mathrm{prev}\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathrm{)}\right]\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{t}}_{\mathrm{all}}\mathit{/}{\mathit{t}}_{\mathrm{cross}}\mathrm{)}\mathit{.}$(5)The average number of clouds produced per stellar period is given by N_tot/P.

Clouds expanding away from the star under radial acceleration increase in tangential diameter by l(r/r_i) and in the radial diameter by l(r/r_i)^ϵ (Eq. (2)). Thus, the number density in a cloud at r is n_q(r/r_i) ^{− (2 + ϵ)}. We use the brightest part of the maser shell as r_typical (see Sect. 7) to estimate the average mass per cloud in the 22-GHz maser shell: $\frac{{\mathit{M}}_{\mathrm{cloud}}}{{\mathit{M}}_{\mathrm{\odot }}}\mathrm{=}\mathrm{4}\mathrm{\times }{\mathrm{10}}^{-24}{\mathit{n}}_{\mathrm{q}}{\left(\frac{{\mathit{r}}_{\mathrm{typical}}}{{\mathit{r}}_{\mathrm{i}}}\right)}^{\mathrm{-}\mathrm{(}\mathrm{2}\mathrm{+}\mathit{ϵ}\mathrm{)}}{\left(\frac{\mathit{l}}{\mathrm{AU}}\right)}^{\mathrm{3}}$(6)(making the same assumptions as for Eq. (4)).

The mass loss rate in clouds is then given by Ṁ_clouds = M_cloudN_tot/t_cross. The filling factor of the clouds by volume is given by ${\mathit{f}}_{\mathrm{V}}\mathrm{=}{\mathit{N}}_{\mathrm{tot}}\mathrm{(}\mathit{l}\mathit{/}\mathrm{2}{\mathrm{)}}^{\mathrm{3}}\mathit{/}\mathrm{(}{\mathit{r}}_{\mathrm{o}}^{\mathrm{3}}\mathrm{-}{\mathit{r}}_{\mathrm{i}}^{\mathrm{3}}\mathrm{)}$.

We apply Eqs. (4) to (6) to the measurements given in Table 2, averaged for the epochs considered, to derive the parameters tabulated in Table 5. The most striking result is that the H₂O maser cloud density at r_i, n_q, is  ≳ 45n_i for all sources. What are the main causes of uncertainty in this estimate? The collision rate, and hence n_q, is proportional to ^√(temperature). The assumed temperature is likely to be an upper limit, almost double some estimates for the lower-mass stars (Zubko & Elitzur 2000), so our adopted value of n_q is probably a lower limit. The main uncertainty in r_i, 10 milli-arcsec for W Hya and  ≤ 5 milli-arcsec for the other sources, arises from determining the stellar position (Sect. 2.1). At most, this could reduce n_q/n_i to 38 for W Hya, and by a smaller proportion for the other sources. In fact, since one of the main criteria in determining the centre of expansion was maximising r_i, the errors are probably lower. v_i is least well-constrained for U Ori and U Her, but these have lower v_i than any other object considered and any increase would increase n_q/n_i.

Equation (4) is distance-independent, as long as r_i and Ṁ are estimated consistently and both have a quadratic dependence on distance. If the highest values of Ṁ are taken from the literature, this could reduce n_q/n_i about fivefold, but it seems very unlikely that all the errors are in the same sense.

The mass loss rate in clouds, Ṁ_clouds is within  ± 80% of the adopted values of Ṁ. This is a relatively small discrepancy compared to deviations of Ṁ in the literature, as mentioned above. The space between the clouds cannot be completely empty so that if accurate values were known, Ṁ_clouds < Ṁ and thus, for S Per, U Her and U Ori, either Ṁ must be higher or the mass in clouds must be smaller (or both).

In all cases, the cloud filling factor by volume f_V is less than 0.01. Using Eq. (11) of R99, n_c/n must be less than 1/f_V, where n_c/n is the cloud number density at any radius. This constrains n_q/n_i < 105 for U Ori and S Per, or larger values for the other sources, which is not inconsistent with our estimates of n_q/n_i.

There are a number of additional sources of uncertainty in our estimates of M_cloud and derived quantities. t_cross has a linear dependence on distance and M_cloud has a cubic dependence, although Ṁ_cloud/Ṁ is independent. N_tot could be an underestimate if there are clouds which have not beamed masers towards us at any of the epochs used in the analysis, although a much higher value would produce an unrealistically high Ṁ_clouds. For the same reason, although the typical values used for l and r_typical and other assumptions about cloud behaviour are very simplified, these are unlikely to be severe overestimates. R11 used these estimates of l and n to find estimates of parameters such as maser brightness temperatures, beaming and optical depths consistent with existing maser models e.g. Elitzur et al. (1992). Altogether, errors of up to factors of a few are possible.

The inescapable conclusion is that the density of 22-GHz H₂O maser clouds is at least an order of magnitude greater than the average wind density, and probably 50–100 times greater.  ≳ 20% of the mass loss is concentrated in these clouds, with a filling factor  < 1%. Menut et al. (2007) modelled IR interferometry of dust emission from within 300 AU of IRC+10216 using 500 clumps, ten- to a hundred-fold overdense, of the same radius as the star, which is similar to the wind structure deduced for our targets.

6. Time evolution, flares and variability

We model the 3-D changes in maser distribution between epochs in Sect. 6.1. A number of individual maser features brighten or dim by large factors. One, in W Hya, is particularly clearly identified in successive epochs of both MERLIN and Pushchino data and the reason for its isolated flare is explained in Sect. 6.2. We investigate whether the imaging data showed any systematic dependence between flux density variations and location in the maser shell in Sect. 6.3, and look for possible causes of variability in Sect. 6.4.

6.1. Time evolution of maser distributions

The positions of masers (x, y) are measured directly in the plane of the sky but their velocities (u = V_⋆ − V_LSR) are measured in the perpendicular direction. We can calculate the full 3-D coordinates using trigonometry, if we assume spherically symmetric expansion away from the star, obeying v = v_i + K_grad(r − r_i). K_grad is defined in Eq. (1) and r and v are the magnitudes of the radial position and velocity vectors for a maser at any point in the shell, i.e. r = (x,y,z). Thus, r² = a² + z² (where a² = x² + y²) and v² = w² + u² (where w is the velocity component in the plane of the sky, in the direction of a). We use the acceleration along the z-axis to derive an expression which can be solved for the position and velocity components not directly measured: $\mathit{u}\mathrm{=}{\mathit{u}}_{\mathrm{i}}\mathrm{+}{\mathit{K}}_{\mathrm{grad}}\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}_{\mathrm{i}}\mathrm{)}$(7)where u_i is the line-of-sight velocity component at the inner maser shell edge at z_i. By simple trigonometry, z_i = r_i(z/r) and u_i = v_i(z/r), giving $\mathit{u}\mathrm{=}{\mathit{v}}_{\mathrm{i}}\mathit{z}\mathit{/}\mathit{r}\mathrm{+}{\mathit{K}}_{\mathrm{grad}}\mathit{z}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{r}}_{\mathrm{i}}\mathit{/}\mathit{r}\mathrm{)}$(8)Substituting $\mathit{z}\mathrm{=}\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}}$, multiplying by r and then squaring gives ${\mathit{u}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathrm{=}\mathrm{[}{\mathit{v}}_{\mathrm{i}}\mathrm{+}{\mathit{K}}_{\mathrm{grad}}\mathrm{(}\mathit{r}\mathrm{-}{\mathit{r}}_{\mathrm{i}}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}\mathrm{)}$(9)which yields a quartic equation in r$\begin{array}{ccc}\mathrm{0}& \mathrm{=}& \mathrm{[}{{\mathit{K}}_{\mathrm{grad}}^{\mathrm{2}}}^{\mathrm{]}}\hspace{0.17em}{\mathit{r}}^{\mathrm{4}}\mathrm{+}\mathrm{[}\mathrm{2}{\mathit{K}}_{\mathrm{grad}}{\mathit{v}}_{\mathrm{i}}\mathrm{-}\mathrm{2}{\mathit{K}}_{\mathrm{grad}}^{\mathrm{2}}{\mathit{r}}_{\mathrm{i}}\mathrm{]}\hspace{0.17em}{\mathit{r}}^{\mathrm{3}}\\ & & \mathrm{+}\mathrm{[}\mathrm{-}{\mathit{u}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{i}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{K}}_{\mathrm{grad}}{\mathit{v}}_{\mathrm{i}}{\mathit{r}}_{\mathrm{i}}\mathrm{+}{\mathit{K}}_{\mathrm{grad}}^{\mathrm{2}}{\mathit{r}}_{\mathrm{i}}^{\mathrm{2}}\mathrm{-}{\mathit{K}}_{\mathrm{grad}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{2}}\mathrm{]}\hspace{0.17em}{\mathit{r}}^{\mathrm{2}}\\ & & \mathrm{+}\mathrm{[}\mathrm{-}\mathrm{2}{\mathit{v}}_{\mathrm{i}}{\mathit{K}}_{\mathrm{grad}}{\mathit{a}}^{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{K}}_{\mathrm{grad}}{\mathit{r}}_{\mathrm{i}}{\mathit{a}}^{\mathrm{2}}\mathrm{]}\hspace{0.17em}\mathit{r}\\ & & \end{array}$(10)where all the coefficients are composed of measured parameters. We solved this (using python numpy²) for r and hence can calculate any other vector component such as w. A similar process was performed by M03 for VX Sgr.

We calculated the fractional flux density per unit volume in segments of arc and in concentric, spherical sub-shells. In all cases the ranges shown are divided into 25 bins of equal widths in angle θ (as projected against the plane of the sky) or radius r (deduced from 3-D modelling). Since the plots are normalised, the entire shell waxing or waning in response to stellar luminosity changes has no effect. Figures 46 to 52 show the angular and radial distributions of flux as a percentage of the total flux at each epoch.

We used the total expansion velocity v to predict the radial position at the latest epoch, for each component observed at earlier epochs. The corresponding extrapolated radial flux distributions are shown by the faint dotted lines. There are minor discrepancies due to binning different samples into each subshell. No compensation has been made for the effects of changing temperature and density on maser efficiency.

The radial profiles in Figs. 46 to 52 show that the brightest masers occur within the inner half of the shell, frequently peaking around the first quarter (as predicted by e.g. Cooke & Elitzur 1985; Yates et al. 1997). However, there are significant differences between the AGB observations, especially over time intervals of  ≥ 2 years. There are most similarities between epochs with the most matched features, as one would expect. The position angle profiles do not show any systematic shifts which would indicate rotation, nor any obvious preferred directions (note that θ is projected against the plane of the sky and does not provide a full 3-D model of angular distribution, such as was made for VX Sgr by M03).

VX Sgr shows fairly good agreement between the radial locations of the peaks seen in 1999 and the extrapolated 1994 profile. This is consistent with the overall similarities in morphology as well as individual feature matches described by M03. They found good agreement between the trigonometric (w) and proper motion expansion (V_a) velocities in the plane of the sky if the star is at 1.8 ± 0.5 kpc, within the uncertainties of the generally adopted value of 1.7 kpc. The emission at greater r has faded, consistent with Cooke & Elitzur (1985).

When the two brightest peaks in the 1994 S Per radial profile are extrapolated, they appear to overtake the 1999 emission. It is possible that v has been overestimated for some maser features. The mean value of w, the component of motion in the plane of the sky estimated trigonometrically, is 15.2 km s^-1, compared with 4.9 km s^-1V_a but Fig. 26 shows a large scatter due to the large uncertainties in the individual measurements at the large distance of S Per. It is also possible that the velocity field is axisymmetric, with slower expansion in the plane of the sky, where the brightest emission is located. The dashed line in Fig. 47 shows that if we extrapolate the 1994 data at 0.4v, the two highest peaks coincide with the 1999 peaks, with a trough in between. R99 esimated the time when the material producing the masers observed in 1994 left the star, and compared this with AAVSO light curves, finding that the trough corresponded to several years when the optical period amplitude was only about 1.5 mag instead of the more usual  ~ 3.

Among the AGB stars, 1996 RT Vir data peak shifts by the predicted amount (equivalent to a change of about one radial bin from the first to the last epoch) but the 1994 data are significantly different. The mean values of w and V_a are 10 and 6 ± 2 km s^-1, in reasonable agreement given the model uncertainties.

U Ori and W Hya show only limited similarities between epochs, consistent with the fewer matches and greater measurement uncertainties, respectively. IK Tau shows closer apparent similarities, including with 1994 data. However, the IK Tau velocity profiles (Figs. 38 and 39) show that, apart from relatively faint red-shifted emission, the peaks in 1994 were at very different velocities from later epochs.

The brightest U Her emission is located at a similar velocity at all epochs, within the thermal line width plus slight acceleration (Figs. 36 and 37). It is possible that the ring of emission seen in Fig. 8 corresponds to the fainter, outer ring in Fig. 9. The outer ring corresponds to the small peak or knee in the 2000–2001 emission at r ~ 32 AU in Fig. 49, the location of the extrapolated 1994 peak (dotted line). If the 1994 peak was scaled down by a factor of about 20, consistent with maser emission fading at larger r, this would produce the flux densities observed in 2000–2001.

We estimated the date relative to the observation date t_⋆ at which each maser component at position r, velocity v, left the star, based on Eq. (3): $\begin{array}{ccc}& & \frac{{\mathit{t}}_{\mathit{\star }}}{\mathrm{day}}\mathrm{=}\mathrm{-}\mathrm{1736}\\ & & \end{array}$(11)We assume that the material reaches r_i from the star travelling at average velocity v_i and ignore the light travel time from U Her to Earth. Figure 53 shows the flux density summed radially into 25 shells covering this time interval (converted to MJD), (similar to Fig. 49 but with the density profile reversed with respect to the abscissa direction), compared with the AAVSO light curve. The two bright shells of U Her maser emission correspond to MJD around 42 800 and 45 000, but there is nothing obviously unusual about the light curved around these times. The intervening maser dip coincides with a fairly bright U Her maximum but it is not exceptional; inspection of AAVSO data shows that it has exceeded mag 7 twenty times in the past century.

6.2. Cloud overlap causes local flare in W Hya

The 22-GHz peaks from W Hya are generally a few hundred Jy, but exceeded 1000 Jy on several occasions during Pushchino monitoring since 1982 (Rudnitskii et al. 1999). Most recorded emission is confined to the velocity range 36–44 km s^-1 but the MERLIN data for 20000405 show a 100-Jy peak at 45 km s^-1 (Fig. 42). This was not present in the Pushchino spectra just a few weeks before or after.

An even stronger flare is seen in Fig. 43 at 40.4 km s^-1, reaching nearly 1000 Jy in 20010227 and over 2800 Jy in 20011023, was detected by both Pushchino and MERLIN. Figure 54 shows the MERLIN peak in this spectral region in 2000, 2001 and 2002 (located in the NW of the plots shown in Fig. 14). The peak is shown by circles and a fainter, nearby feature, with a distinct, Gaussian spectral profile, is shown by squares. The brighter peak is 7, 136 and 269 Jy and the fainter is 6, 7 and 2 Jy, at the three successive epochs. The relative velocities are  < ΔV_th, the fainter peak being about 1 km s^-1 more redshifted at all epochs, although both show a small drift. Interestingly, the fainter feature is to the N in 2000 but to the S at later epochs, an apparent relative velocity in the plane of the sky of  ~ 5 km s^-1.

The MERLIN data were taken close to stellar minimum whilst the brightest Pushchino peak was seen close to stellar maximum but this cannot be the main reason for the flare since only one feature is affected. It is unlikely that a large-scale shock is involved, since in that case one would expect both adjacent features to be affected. The third possible mechanism for an isolated flare is two overlapping clouds. The brighter cloud is more blueshifted and therefore in the foreground, and could dramatically amplify the background cloud (Kartje et al. 1999). Figure 54 shows that the brigher cloud has a proper motion consistent with outflow; the full extent of the background cloud may be obscured.

6.3. Variability of multi-epoch features

We looked for any systematic trends in the brightness evolution of maser features identified at multiple epochs. The total distance r of each feature from the assumed stellar position was estimated as described in Sect. 6.1. The 3-D modelling is not only vulnerable to measurement errors but also to deviations from the assumptions of spherical symmetry and linear acceleration, but with only two epochs of data (in most cases) these effects cannot be separated. Figures 55 to 61 show the change in flux density as a function of r for matched features.

S Per and W Hya do not show a clear trend, although this may be obscured due to the larger position uncertainties owing to the greater distance and low elevation, respectively. There may be some spurious matches such as the blue-shifted S Per features with very large changes in r, since we did not want to impose an a-priori 3D model but used the same position matching critera at all velocities. However, if spherical expansion dominates, the apparent large proper motions are unlikely for highly blue-shifted features. A possible sign of the effects of asymmetric expansion is described in Sect. 6.1.

Features throughout the maser shell behave similarly in the remaining sources. The majority of VX Sgr, U Her and IK Tau matched features get fainter, generally consistent with Figs. 32,  37 and 39. Most of the U Ori matched features get slightly brighter, consistent with Fig. 35. The tendency for r to decrease in U Ori may be due to a large random component to its proper motions (Sect. 3.2).

RT Vir has the best sampled data. Figure 60 shows that all 11 matched features increase in brightness relative to the first interferometry epoch of 1996. Nine of the features continue to increase, although a few level off or slightly decrease at the final epochs. The other two features reach maximum at intermediate epochs. Figure 24 indicates the fractional increase. Figures 40 and 41 show that the main features of the RT Vir spectrum dropped in flux density between the Pushchino observations in 19960319 and the MERLIN observations in 19960405. The brightness then increased steadily during the interferometric monitoring but dropped sharply again between the last MERLIN epoch in 19960612 and Pushchino observations in 19960619. Figure 3 shows that the MERLIN observations were made during the decline of optical phase.

6.4. Possible causes of variability

In all objects, individual maser features can show idiosyncratic behaviour. Some of this can be attributed to local turbulence (Sect. 4.2). Extreme flares are confined to small spectral and spatial regions in our data, such as the W Hya event (Sect. 6.2). There are several possible types of more general behaviour. Masers switch on, brighten relatively rapidly and then decline slowly as the parent gas enters and eventually leaves the region where maser excitation is possible, over several decades, as predicted by Cooke & Elitzur (1985) and seen in Figs. 46 to 52.

There is no obvious correlation between flux density changes and radial position in the shell of any CSE studied. In those objects which do show systematic behaviour, the whole shell is affected, suggesting a response to stellar illumination (the light crossing time for an entire 22-GHz H₂O shell is a couple of days or less). Figures 1 to 3 show that there is no systematic relationship between maser brighness and the optical phase at the times of our observations. However, this is not surprising, since IR is expected to have more effect on 22-GHz masers. Smith et al. (2006) found 0.1–0.2 P lags of IR variability with respect to the optical light curves for Miras, where P is the stellar period. There was no clear relationship for other sources such as SRbs.

Long-term, well sampled single dish data have shown some correlations between H₂O maser variability and stellar periods. Shintani et al. (2008) monitored a large number of AGB and RSG, including all those in our sample, for 3–4 yr using the VERA Iriki single dish. They fitted periods to the integrated 22-GHz emission and compared these with AAVSO data. Excluding fits marked “bad”, they deduced lags, as fractions of P, for: U Ori 0.47 ± 0.03; U Her 0.34 ± 0.07; IK Tau 0.10 ± 0.03; W Hya 0.19 ± 0.04.  ~ 20 yr of Pushchino monitoring showed similar lags of 0.2–0.4 P for U Ori (Rudnitskij et al. 2000) but W Hya produces semi-regular flares with lags of 0–2.1 P (Rudnitskii et al. 1999). No clear relationship was found for RT Vir. As for the RSG, Pashchenko & Rudnitskii (1999) and Lekht et al. (2005) found phase delays of  ≤ 1 and 0.01–0.5 P for VX Sgr and S Per, respectively. In summary, the clearest relationship is seen for Miras, and tends to align 22-GHz maser and IR periodicity. A relationship between stellar IR and maser brightness is not excluded for the SRbs and RSGs but would not be expected to follow the optical light curve. In fact, OH mainline masers around SRbs show a mixture of cyclic and erruptive behaviour with variability on timescales a short as a month (Etoka et al. 2001), so H₂O masers closer to the star would be expected to be even more variable.

Changing stellar IR will affect the 22-GHz maser in at least 3 ways. Firstly, it is collisionally pumped, so the intensity is sensitive to heating (Cooke & Elitzur 1985). Secondly the pump cycle may be suppressed by a strong IR radiation field (Yates et al. 1997). These phenomena would affect the whole shell at close to light speed, i.e. within a few days. This could explain the behaviour of VX Sgr, U Ori, U Her, IK Tau and RT Vir, where the majority of matched maser features change flux density in concert between our observational epochs. Thirdly, the changing radiation pressure on dust will affect the velocity field and, for example, if acceleration decreases, radial beaming will be enhanced at the expense of tangential beaming (Engels et al. 1997), producing anti-correlated behaviour in the extremes versus the middle part of the spectrum. Our data are not well enough sampled over more than one stellar period, which would be required to test whether there is a correlation between the maser structure and the phase cycle.

Rudnitskii & Chuprikov (1990) suggested that the stellar pulsations shock the maser shell, producing brightness variations which would lag the optical phase by an interval equivalent to the shock crossing time from the stellar surface to the 22-GHz r_i. Shintani et al. (2008) found evidence for lags of  > 0.5P in many objects, taken to support a similar model. Reid & Menten (1997) showed that the pulsation velocity at or above the radio photosphere is unlikely to exceed 5 km s^-1 in AGB stars. The maximum outflow velocities of SiO masers provide another indication of the velocity of possible impacts on the H₂O maser shell. The SiO J = 1 − 0, v = 1 and 2 transitions in all our targets were surveyed by Kim et al. (2010). Although SiO masers are even more variable than H₂O, results from other published observations are generally within a thermal linewidth. The SiO velocities are 9–10 km s^-1 for VX Sgr and S Per, 6 km s^-1 for U Ori, IK Tau and W Hya, 4 km s^-1 for U Her and 5 km s^-1 for RT Vir. This is comparable to v_i and  ≪ v_o for all objects except U Ori and W Hya where it is comparable to v_o.

A shock velocity of 5–10 km s^-1 corresponds to a proper motion of  ≈ 1–2 AU yr^-1 so a single shock would not affect an entire shell; even the smallest (the SRbs) are  ≳ 10 AU thick. It is possible that shocks affect the inner edge of the 22-GHz shell, and possibly penetrate deeper in objects where v_i is smallest, such as U Ori and U Her, which show distinctive “matter bounded” beaming properties for some masers (R11).

Höfner et al. (1995) developed a model in which C-rich dust production is intrinsically episodic. The formation of dust during one time interval leads to back-warming and suppression of further grain growth nearer the star, until radiation pressure has driven the initial dust shell away from the star. The models predict enhanced dust formation episodes at intervals similar to, or longer than the stellar period, but not all episodes are equally efficient. The dust:gas ratio, and thus acceleration, is therefore a function of distance from the star. Jeong et al. (2003) applied this model to an O-rich Mira. The amplitude of velocity variations depend on the extinction coefficients used, giving  ~ 1–7 km s^-1 at 15 R_⋆. Our data could be used to test more precise models for our stars.

7. Relationship between maser cloud size and parent star

The average cloud radius for each epoch ${\mathit{R}}_{\mathrm{c}}\mathrm{=}\overline{\mathit{l}}\mathit{/}\mathrm{2}$. Figure 62 shows that this is an almost linear function of stellar radius R_⋆. The error-weighted fit and dispersion are shown by the solid and dashed lines, where the relationship is $\log {\mathit{R}}_{\mathrm{c}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{0.17}\mathrm{\pm }\mathrm{0.17}\mathrm{)}\mathrm{+}\mathrm{(}\mathrm{1.00}\mathrm{\pm }\mathrm{0.13}\mathrm{)}\log {\mathit{R}}_{\mathit{\star }}\mathit{.}$(12)Note that W Hya cloud sizes are most likely to be underestimated, see Sect. 4.1.

The brightest part of the shell, where the highest number of clouds are detected, is close to the inner radius (see Sect. 6.1), so we use r_typical = r_i + (r_o − r_i)/4 as the typical location in the shell where clouds have radius R_c. If the clouds are formed at the surface of the star and expand uniformly, radially in the wind, the average cloud radius, extrapolated back to the surface of the star, is approximately (0.04–0.09)R_⋆ for the RSG, (0.05–0.12)R_⋆ for the Miras and (0.06–0.07)R_⋆ for the SRbs, i.e.  ~ 5–10% R_⋆ for all stars.

This implies that the size of the clouds is determined by the size of the star. The 22-GHz H₂O masers arise from clouds which are much denser than the average wind density (R99, B03, M03, Sect. 5.2) and which must survive for decades (even if individual masers do not). If these clouds were to be formed in the outflow after it had left the star, the cause would be some local wind phenomenon such as cooling due to dust formation. The length scale of this process would be determined by basic physical properties, implying that all clouds would form on the same size scales. This is not consistent with our data.

Other observations have revealed inhomogeneities and clumps closer to the star. Chiavassa et al. (2009) compared VLTI image reconstruction of VX Sgr at λ 1.6 μm with three-dimensional simulations of surface convection. The best-fitting synthetic image possessed two convective cells of 4–5 milli-arcsec (comparable to R_⋆ at this wavelength), with many 0.5–1 milli-arcsec granules, two being particularly prominent. A third spot appears at  ~10 milli-arcsec separation from the centre, i.e. outside the photosphere. Although the fine details of the models are not uniquely constrained, the stellar surface clearly possesses structure on scales including the extrapolated birth size of the clouds,  ~5–10% R_⋆. Similar models for AGB stars have so far only been published for C-rich objects (Freytag & Höfner 2008), which also produce a few large convection cells. These have a piston effect comparable to traditional pulsations, and thus are potentially capable of enhancing mass loss on sub-surface scales. U Ori was imaged by Pluzhnik et al. (2009) using IOTA. The H₂O band image at 1.78 μm shows three bright clumps just outside the photospheric radius³. Recent VLTI observations of four O-righ Miras (not in our sample) by Wittkowski et al. (2011) show deviations from symmetry in the H₂O band indicating the presence of clumpiness on scales approaching 10% of the stellar size.

At much greater distances, around carbon stars, detatched CO shells have been observed and modelled by Bergman et al. (1993) and Olofsson et al. (1996). The clumps were not individually resolved, but the inhmogeneity was inferred from modelling larger-scale images. The objects studied have mass loss rates in the range (0.4−4) × 10^-5   M_⊙ yr^-1 concentrated in clumps of mass  ~2 × 10^-5   M_⊙, similar to the values for our RSG (Table 5). The model of Bergman et al. (1993) is based on S Sct, using a clump radius R_c ≈ 1350 AU at a distance r ≈ 36   000 AU from the star. Olofsson et al. (1996) model the cloud expansion as R_c ∝ r^0.8 for steady outflow and limited cooling. Extrapolated back to r = 60 AU gives R_c ≈ 8 AU, in agreement with the estimated H₂O cloud radii around S Per and VX Sgr. Recent measurements of dust clumps in two detatched CO shells (Olofsson et al. 2010) yield similar sizes.

8. Summary and further work

We have resolved the detailed structure of the approximately spherical, thick H₂O maser shells around two RSG (VX Sgr, S Per), three Miras (U Ori, U Her, IK Tau) and two SRb stars (RT Vir, W Hya), at multiple epochs. MERLIN detects all the 22-GHz maser emission (or possibly almost all, in the case of the closest object), at high enough resolution to measure the sizes and proper motions of individual clouds. We compare these images with single dish monitoring by the Pushchino radio telescope. Our results provide new insights into the development of the stellar wind in the region where the escape velocity is exceeded.

Survival.

The H₂O maser shells are 10–250 AU thick with crossing times from a few tens of years for AGB stars to nearly a century for RSG, but most individual maser features survive less than 1–2 yr and 1–2 decades, in the AGB stars and RSG, respectively. The disappearance and reappearance of individual masers identified by comparing imaging with single dish monitoring shows that this contradiction is resolved if the masers emanate from long-lived clouds. These undergo subsonic or mildly supersonic turbulence of a few km s^-1, which causes detectable masers to be alternately boosted and supressed without destroying the clouds themselves.

Expansion proper motions.

The proper motions of VX Sgr, S Per, U Her, IK Tau, RT Vir and most epochs of W Hya are unambiguously dominated by spheroidal expansion. There is a hint of rotation in the IK Tau maser proper motions. The small tangential component could be caused by random motions, which is also likely to be a significant factor affecting U Ori. The radial proper motion velocities (with respect to the assumed stellar position) are generally consistent with V_LSR measurements, but the former have larger measurement errors and scatter, especially for the low declination source W Hya. Discrepancies between velocities derived from proper motion and V_LSR measurements could also be due to asymmetry, as seen for VX Sgr (M03).

Acceleration.

The 22-GHz masers lie outside the dust formation zone, where the wind is driven by radiation pressure on grains. The outflow velocity increases twofold or more across the maser shells of all sources. Moreover, comparison of different sources shows that the larger the inner and outer shell radii, the greater the expansion velocities at the outer rims, although the acceleration is more gradual. This is likely to be due to changes in the dust optical properties so a higher proportion of radiation is absorbed, and/or increasing optical depth so that the stellar photons are re-emitted making radiation pressure more efficient. This is consistent with Herschel and other results which show that terminal velocity is not reached until hundreds of R_⋆.

Dense clumps.

The wind at  ~5–50 R_⋆ is concentrated in clumps which are 40–110 times denser than the wind average. They comprise between a fifth and almost all of the mass in the stellar wind in the maser shell, although accurate comparisons are difficult due to uncertainties in measuring the total mass loss rate. The volume filling factor is less than 1%. At least two to six clouds are formed per stellar period.

Flare due to cloud overlap.

Single-dish monitoring of H₂O masers reveals various forms of variability; in some instances a single spectral feature brightens by more than one order of magnitude for between a few weeks and a year or two. One such flare, in W Hya, was imaged at 3 epochs from 2000–2002, bracketing the 2500-Jy maximum captured by Pushchino. The images show that the variable cloud appears to pass in front of another. Two sets of components are spatially differentiated, at slightly different V_LSR, both with distinct Gaussian spectral profiles. The flaring feature is initially to the S of the other cloud, both  < 10 Jy. Two years later it has exceeded 400 Jy and moved to the NW of the fainter feature (in a direction consistent with outflow). This is a dramatic confirmation of the predictions of Kartje et al. (1999) for maser amplification by overlapping clouds.

Variability throughout the maser shell.

The majority of matched features in each of VX Sgr, U Her and IK Tau decrease in flux density between the two imaging epochs regardless of their position in the maser shell; the majority of U Ori matched features increase. The intervals between epochs are of order 1–2 stellar periods, but there is no obvious connection between optical phase and brightening or dimming. Features throughout the RT Vir maser shell, imaged six times during the decline of optical brightness, rise and then slightly decline again in 22-GHz flux density. S Per and W Hya matched features show an apparently random mixture of brightening and dimming. The coordinated behaviour of features in 5 CSEs suggests that radiative effects play a major role in determining maser brightness, probably linked to the IR stellar phase (which lags the optical). The 22-GHz maser shells are 10–200 AU thick, so the light travel time is less than a few days but a shock at 5–10 km s^-1 would take many years to cross even the thinnest shell. Nonetheless, the H₂O masers may be affected by shocks close to the inner rim of the shell.

Cloud size depends on star size.

The average radius of water maser clouds is 0.5–2 AU for the AGB stars and 6–9 AU for the RSG, showing a close dependence on R_⋆. Assuming the clouds expand steadily in the outflow, this corresponds to birth radii  ~ 5–10% R_⋆. Previous CO and dust imaging have hinted at similar clumping scales and overdensity, but only MERLIN H₂O maser observations have well-resolved all the emission from such clumps, close to the star. The dependence on stellar size suggests that stellar surface phenomena such as convection cells determine the scale of clumps.

Future work.

We are in the process of analysing multi-epoch MERLIN and VLBI observations of OH masers which probe different conditions in these objects. A future paper will compare the H₂O morphology with the distribution of OH masers to investigate asymmetry and inhomogeneity. There are very few estimates published for the sizes of SiO or OH maser clouds in our sources, but we will seek available data where the clouds are likely to be resolved but not resolved-out, to test whether the sizes do increase with distance from the star.

The advent of ALMA, e-MERLIN, upgraded VLBI and other interferometers will allow the unanswered questions to be tackled. The current baselines of ALMA, up to 1 km, will locate sub-mm H₂O masers with respect to the 22 GHz transition and the dust formation zone, mapping changes in excitation conditions. Eventual longer baselines and full sensitivity will resolve both dust and molecular species in clumps, showing whether the 22-GHz clouds are indeed dustier and denser. ALMA will also resolve the stars, as will e-MERLIN and the EVLA at its highest frequencies. This will reveal any convective disurbances and timely VLBA SiO maser monitoring will show whether these lead directly to clumpy mass loss, and if the SiO clumps in turn become dusty water maser clouds, resolved by ALMA and e-MERLIN in multiple transitions. Simultaneous imaging of the star and masers at 22 GHz (only previously achieved for a very few objects e.g. Reid & Menten 1990; Reid & Menten 1997) will solve the current astrometric ambiguities.

Online material

Appendix A: Contour plots of H₂O maser emission

These plots show channel-averaged emission (see Sect. 2 for the original channel widths used). The restoring beams are shown in the first or last panel; in most cases these are slightly larger than the beams used in the full-resolution images for analysis (Table 2). The lowest contour is at the average 3σ_rms for the cube but bright channels have been blanked where they are dynamic-range limited. Plots for VX Sgr were published by M03.

The image cubes for these plots were made by averaging 5 to 10 channels (~1 km s^-1) and thus, due to frequency dilution, the peaks are often lower than those measured at full resolution. The same contour levels were used for all epochs of each source, at  ~3–5σ_rms, to aid comparison. Negative contours are not shown to avoid crowding the plots; these generally do not exceed the first positive but blanking has been applied to dynamic-range limited channels where indicated (identified by broken contours).

Appendix B: Measured maser feature parameters

These tables list the maser feature parameters. Lists for VX Sgr were published by M03.

---

Acknowledgments

We recall the tremendous contributions of Jim Cohen (1948–2006) in the early years of this project. We warmly thank the referee, Dieter Engels, for a very thorough reading of the paper which led to great improvements in clarity, Indra Bains for many useful contributions and Robert Laing for trigonometric discussions. We acknowledge, with thanks, our use of data from the AAVSO (American Association of Variable Star Observers), the GCVS (General Catalogue of Variable Stars) and the Vizier and ADS services.

References

All Tables

All Figures