The Astrometric Animation of Water Masers toward the Mira Variable BX Cam

Abstract We report VLBI monitoring observations of the 22 GHz H2O masers toward the Mira variable BX Cam. Data from 37 epochs spanning ∼3 stellar pulsation periods were obtained between May 2018 and June 2021 with a time interval of 3–4 weeks. In particular, the VERA dual-beam system was used to measure the kinematics and parallaxes of the H2O maser features. The obtained parallax, 1.79±0.08 mas, is consistent with Gaia EDR3 and previous VLBI measurements. The position of the central star was estimated relied on Gaia EDR3 data and the center position of the 43 GHz SiO maser ring imaged with KVN. Analysis of the 3D maser kinematics revealed an expanding circumstellar envelope with a velocity of 13±4 km s−1 and significant spatial and velocity asymmetries. The H2O maser animation achieved by our dense monitoring program manifests the propagation of shock waves in the circumstellar envelope of BX Cam.


INTRODUCTION
The long-period variables (LPVs), such as Mira variables, are stars of low to intermediate masses which have reached the late evolutionary stage of Asymptotic Giant Branch (AGB) phase.They are characterised by long-period (>100 d) variations in radius, brightness, and temperature, which are caused by stellar surface instability and radial pulsation (see the latest catalogue of Lebzelter et al. 2022).These stars also have an intense mass-loss phenomenon or a superwind, which leads to formation of a circumstellar envelope (CSE) made of gas and dust (see the review of Höfner & Olofsson 2018).
The stellar mass loss may be driven, primarily, by the pressure from the stellar radiation on dust grains, while pulsation-induced shocks are expected to enhance the mass loss (e.g., Habing & Olofsson 2004;Höfner 2011;Höfner et al. 2016).In addition, the stellar structure is basically spherically symmetric, while the mass loss often exhibits significant asymmetry on very different spatial (1-10 4 AU) and temporal (a few months to a few 10 4 years) scales (Höfner & Olofsson 2018).The mass loss asymmetry observed in the vicinity of the stellar surface (e.g., Khouri et al. 2020) might be linked to the existence of giant gas segments created from convection in the stellar interior.It is interesting to trace how such asymmetry will grow while the gas segments are being accelerated outward.It is likely that the mechanism of dust acceleration in the CSE, the enhancement of the outward flow, and the growth of its asymmetry may be dependent on the property or abundance of the ejected materials (i.e., carbon-rich or oxygen-rich composition) (Höfner & Olofsson 2018).
Circumstellar SiO and H 2 O masers are located at, respectively, ∼2-4 stellar radii (e.g., Diamond et al. 1994) and ∼5-50 stellar radii (e.g., Richards et al. 2012) of a typical oxygen-rich AGB star, while the dust formation layer is co-located between the SiO masers and H 2 O masers (Habing 1996;Wittkowski et al. 2007).The location of SiO masers corresponds to the innermost part of the CSE, where molecules are still rich and outward and inward motions are often observed (e.g.Gonidakis et al. 2013).The location of H 2 O masers corresponds to the zone of the outward acceleration of dust and gas, where molecular gas is cooled and condensed to dust while H 2 O molecules are still enriched.The maser structures are composed of clusters of compact maser features (isolated gas clumps with sizes of ∼1 AU).This enables us to measure the three-dimensional velocity field of the CSE through monitoring observations by Very Long Baseline Interferometry (VLBI) with excellent angular resolution (e.g., Imai et al. 2003).Groups of maser features may trace giant gas/dust clumps, whose origins have been debated in terms of the convection of the stellar interior or just instability on the stellar surface (Höfner & Olofsson 2018).Thus, the maser spatio-kinematics have direct relationship to such gas/dust dynamics and can illuminate important issues on the mechanisms of maser excitation.
The EAVN Synthesis of Stellar Maser Animations (ESTEMA) project has been conducted as one of the KaVA (a combined array of the Korean VLBI Network, KVN, and Japanese VLBI Exploration of Radio Astrometry, VERA) Large Program projects followed by a series of General Proposal projects of the East Asian VLBI Network (EAVN), aiming at intensive VLBI monitoring observations of CSE masers associated with LPVs of different pulsation periods (300-1000 d) over a few stellar pulsation cycles.EAVN is suitable for high cadence VLBI monitoring because more than six telescopes are always available for imaging of CSE masers.The typical angular resolutions are 1.2(K)/0.6(Q)mas with the core array KaVA, and the largest detectable angular scales are 9(K)/5(Q) mas with the short baselines of ∼300 km provided by KVN.The expected image sensitivities are ∼30(K)/50(Q) mJy beam −1 for KaVA and ∼40(K)/60(Q) mJy beam −1 for KVN, with an integration time of 4 hours and narrower bandwidth of 15.625 kHz for maser emission.In particular, it enables us to simultaneously monitor four to eight lines of H 2 O and SiO masers at four frequency (22/43/86/129 GHz or K/Q/W/D) bands with the KVN.Absolute positions, proper motions of maser gas clumps, and the annual parallax of the maser source system can be determined by precise astrometry with the dual-beam system of VERA.The animation of circumstellar masers would demonstrate in detail the behaviors of maser emission, reflecting the dynamics of the maser clumps and variation in regions of maser excitation around the LPVs.
Here we present the first result of ESTEMA, which has been derived from H 2 O masers associated with the Mira variable star BX Camelopardalis (BX Cam).This star is an oxygen-rich Mira variable with similar, but slightly inconsistent, reported pulsating periods of 486 d (AAVSO)1 or 440 d (Matsuno et al. 2020).The spectral type is M9.5 (Solf 1978), and the presence of bright SiO, H 2 O, and OH masers have been confirmed (Ukita & Le Squeren 1984;Matsuno et al. 2020).The trigonometric parallax and the spatial distribution of H 2 O masers with three collimated flows have been reported by Matsuno et al. (2020).This star was selected as one of the ESTEMA targets because of the more flexible allocation of VLBI observation sessions for this Northern-sky star.Similarly to the intensive monitoring observations of SiO (v=1, J=1→0) masers around TX Cam with the Very Long Baseline Array (VLBA) (Diamond & Kemball 2003;Gonidakis et al. 2010Gonidakis et al. , 2013)), the present monitoring project enables us to elucidate the global behaviour of an example of CSE H 2 O masers for the first time.CSE H 2 O masers also should be physically related to the stellar pulsation because of the periodic variation in the maser spectra (Shintani et al. 2008;Kim et al. 2014).However, the true origin of such a maser behavior can be visualized by an animation of the VLBI images.The astrometry on the maser animation is crucial to exactly trace the motions of the individual maser features.After subtracting maser position modulations due to the annual parallax and the secular motion of the maser source system, the maser maps can be registered on the reference frame, which then provides a fixed frame that can be compared with the location of the pulsating star whose astrometric data is now available from the Gaia Early Data Release 3 (Gaia EDR3) (Gaia Collaboration et al. 2021).A maser animation derived based on above astrometric steps enable us to visualize a radial expansion of the CSE masers and some deviations of the maser motions from constant-velocity motions, such as radial accelerations, rotations, etc.

OBSERVATIONS AND DATA REDUCTION
The monitoring observations of H 2 O and SiO masers around BX Cam in ESTEMA have been ongoing since 2018 May, except in the EAVN maintenance season (June-August).The EAVN telescopes that have participated in ESTEMA observations include four telescopes of VERA (Mizusawa, Iriki, Ogasawara, and Ishigakijima), three telescopes of the KVN (Yonsei, Ulsan, Tamna), two telescopes of the Chinese VLBI Network (CVN: Tianma, Urumqi in K-band only), and one Japanese telescope (Takahagi in K-band only).The time cadence for a single observational epoch is 3-4 weeks, which is approximately equal to a time resolution of 1/20 stellar pulsation cycles of BX Cam (∼22 d).Each observation epoch corresponds to a pair of K and Q-band experiments over two consecutive days and lasts for 8 hr (2 d×4 hr).We have designed the observations to use all of the unique capabilities of the KVN and VERA in EAVN, using a hybrid setup of simultaneously monitoring H 2 O (22.235080 GHz) and SiO (42.820539, 43.122027, 86.243442, and 129.363358 GHz) maser lines and conducting high-accuracy astrometry by observing the maser and the reference continuum sources simultaneously.
In this paper, we mainly focus on the K-band data as listed in Table 1, spanning observations from 2018 to 2021.The data sample has crossed three light curve maxima of the stellar pulsation from phase φ 0.65 to φ 3.17.Each of the observation sessions had a duration of 4 hr including scans on the target, BX Cam, the fringe finder, NRAO150, and the delay calibrator, J0721+7120 (α J2000 = 07 h 21 m 53.s 4485±0.31mas, δ J2000 = 71 • 20 36.363±0.10mas).We alternated scans between BX Cam for 7 min and the calibrator for 2 min.The angular separation between these two sources is 7.9 • on the sky, which is suitable for source-frequency phase-referencing (SFPR) analysis with KVN data (Dodson et al. 2014).Meanwhile, the phase reference source J0524+7034 (α J2000 = 05 h 24 m 13. s 4333±0.39mas, δ J2000 = 70 • 34 52.906±0.13mas2 ), 2.0 • away from the target, was observed by the second beam of VERA to determine the absolute positions of the maser spots of BX Cam.
The received signals in left circular polarization were recorded in all the EAVN stations with 16 base-band channels (BBCs) each with a bandwidth of 16 MHz, one of which includes the maser emission of BX Cam.In VERA, 15 BBCs were assigned to the signals received by the second beam so as to cover a total width of 496 MHz.The recorded signal data were correlated using three passes with the hardware correlator in Daejeon, Korea (Lee et al. 2015).The first pass was processed for the BBC including H 2 O masers with all EAVN antennas and yielded 2048 spectral channels, corresponding to a velocity spacing of 0.105 km s −1 at 22.235 GHz.The second pass was processed for the phase reference source J0524+7034 observed by the second beam of VERA.The third pass was for simultaneously monitoring of H 2 O and the four SiO maser lines using the KVN.In this pass, all of the 16 BBCs were processed with 1024 spectral channels per BBC.
The data calibration was performed with the Astronomical Image Processing System (AIPS) (Greisen 2003) developed by the National Radio Astronomy Observatory (NRAO).We have developed the pipeline scripts with Parsel-Tongue (Kettenis et al. 2006) for ESTEMA project.A priori amplitude calibration was done using the system noise temperatures and antenna gains logged during observations at KVN and VERA stations.These system noise temperatures were measured with the chopper-wheel method (Ulich & Haas 1976) and were evaluated using the "R-Sky"method (Honma et al. 2004), by observing a reference black body.Since the system temperature measurements were absent for Tianma, Urumqi, and Takahagi stations, we applied the template spectrum method (Reid 1999), in which a template spectrum of the H 2 O maser emission was obtained from the auto-correlation data.The estimated amplitude calibration errors of both methods were less than 10% for EAVN (Cho et al. 2017).For imaging the H 2 O masers with all the EAVN antennas, the velocity channel including a bright, single, and compact maser spot was selected as the reference for fringe fitting and self-calibration.The detection limit was typically 0.4 Jy beam −1 at a 5-σ noise level, in maps without bright maser emission.
For astrometry with VERA, the instrumental relative delays between the dual beams measured by the horn-on-dish method are applied (Honma et al. 2008) and the delay re-calculation tables (Nagayama et al. 2020) were used.In the delay re-calculation tables, the more accurate delay models were updated, with: the tropospheric zenith delay measured by the Global Positioning System (GPS), the ionospheric Total Electron Content (TEC) of the Global Ionosphere Map (GIM) produced by the Center for Orbit Determination in Europe (CODE), the station position measured in the monthly geodetic VLBI observations, and the updated maser positions (Nagayama et al. 2020, and the references therein).Then we used J0524+7034 as the phase reference source to determine the absolute position of masers.5 σ noise levels of the phase referenced map of VERA data are typically 1.0 Jy beam −1 .Therefore, we can measure the absolute positions of the bright masers of BX Cam with the VERA dual-beam bona-fide astrometry, then propagate this to all of the maser spots in the EAVN image by registering the bright masers in the pair of images.We used the SFPR technique described in Dodson et al. (2014) to obtain the astrometric registration of the H 2 O and SiO maser spots from the KVN data.

Periodic variation of the H 2 O maser spectra
Frequent observations, as well as the time-span of ESTEMA, allow us to study the variability and periodicity of maser intensity on timescales of a few weeks to a few years in terms of the possible correlation with stellar pulsation (Gonidakis et al. 2010).Figure 1 demonstrates the changes in the spectrum with time.Since some of imaging observations of the masers failed and some of the maser images are made by different array configurations, we use the maser intensity on the KVN stations for spectral analysis, e.g., making cross-power spectra on a short baseline in Figure 1, making the total-power spectra on each KVN stations, and analysing the period and phase of the spectra.
The excitation of H 2 O masers is sensitive to both the dynamics of a stellar wind and the stellar radiation with periodic variation.The previous monitoring observations of H 2 O masers around BX Cam with the VERA Iriki telescope (Shintani et al. 2008;Matsuno et al. 2020) did not determine the period and phase lag of the H 2 O maser variability.
Multiple estimates of the optical period of BX Cam have been reported by the American Association of Variable Star Observers (AAVSO) and Matsuno et al. (2020).AAVSO reported a period of 486 d3 without the details of the specific time range of data or a specific method.Matsuno et al. (2020) determined a period of 440 d using the 40 nights in AAVSO database between 2016 February and 2019 December.We also reanalysed 54 nights from the AAVSO database at V-band from 2016 February to 2021 October using two methods.One is the Date Compensated Discrete Fourier Transform (DC DFT) method with the popular analysis tool VStar (Benn 2012), and another is the Lomb-Scargle Periodogram (Lomb 1976;Scargle 1982;VanderPlas 2018), which is equivalent to least-squares fitting of sine waves.We obtained an averaged optical period of 440±5 d using the two methods as shown in Figure 2, which is consistent with Matsuno et al. (2020).
The observed properties of the spectral and flux variability for the H 2 O masers are as follows.
(1) The H 2 O masers are highly variable so that no maser components were detected in some of the maser flux minimum.(2) The maser flux and the optical magnitude exhibit a similar sinusoidal pattern, showing that the maser flux follows the pulsation of the star (same as the SiO masers Gonidakis et al. 2010).(3) The maxima of the periodic maser flux variation do not have a constant magnitude.(4) Time lags of the maser flux maxima with respect to the optical light maxima are also expected, but they are not necessarily constant as seen in the data (same as for the SiO masers Gonidakis et al. 2010).( 5) The maser flux variations do not follow a specific pattern.The number of spectral features, as well as their fluxes, change not only during a single cycle but also from one cycle to another.(6) Systematic radial velocity drifts in the maser features are also found as shown on the spectrum.

The mapped H 2 O maser features
Figure 3 shows the comparison between the maser brightness distributions taken with KVN, VERA and EAVN (KVN+VERA+T6) on 2020Jan10 (φ=2.00).KVN is helpful to detect extended structures, while EAVN can detect a larger number of isolated compact maser features.VERA can measure the bona-fide positional information of bright masers using the dual-beam astrometric system.The H 2 O masers were spatially well resolved into individual maser features with the EAVN synthesized beam.The image including Urumqi baselines shows no obvious difference, since only the brightest and compact maser features located at the map origin were detected on > 3000 km long baselines.  .Integrated flux curves of 22 GHz H2O masers from the cross (red circles) and total (green circles) power spectra of the KVN and V-band optical light curve (blue stars) from the AAVSO monitoring data.The total power spectra is using the weighted average of three KVN antennas (Kt,Ky,Ku), while the cross power spectra is using the short baselines as shown in Figure 1.The V LSR velocity range for estimating the integrated flux densities is from −16 to −7 km s −1 and from 7 to 10 km s −1 .At K band, the beam size of the KVN single-dish telescope is ∼130 as and the synthesized beam size of the projected KYS-KUS baseline is ∼11 mas for BX Cam.The black solid line shows the sinusoidal fit to the optical light curve with a period of 440 d.The AIPS task SAD was used to extract maser spots information, e.g., positions and brightness, from the EAVN image data cube.For the component identification, a non-constant Gaussian was used and features with flux higher than six times the noise level were accepted.Maser feature candidates were identified as components if found in more than three consecutive channels within the cube (Gonidakis et al. 2010).Table 1 gives the possible numbers of maser features detected at individual epochs.We traced the relative positions and radial velocities of maser features at Xu et al.
different epochs, which had been stable from one epoch to another within 1 km s −1 in LSR velocity and within 6 mas yr −1 ×∆t yr in position.The ∆t is the time gap between two epochs and the limit corresponds to the maximum movement velocity in Table 2.In practice it is difficult to automatically cross-identify components just relying on the above velocity and position window, thus we identified the components manually.Figure 4 shows the example of internal structures of the maser features.Among them, 17 maser features with the lifespan longer than one Cycle are labelled with their identifier (MF#).Each maser feature is composed of a compact bright part and an extended structure.There are "spoke-like" (linearly distributed) (Zhang et al. 2012) maser features at all epochs.These maser features appear to be composed of gas flowing outward from the central star and mostly have decreasing (blue-shifted) or increasing (red-shifted) radial velocities.As marked in Table 1, 14 epochs of data yielded successful phase-referencing astrometry with VERA dual-beam system.The dual-beam astrometry failed at other epochs due to various reasons, such as the masers in the trough of the stellar pulsation cycles being too weak, key stations missing, and low sensitivity due to bad weather.The compact calibrator J0524+7034 shown in Figure 6 was used as the phase reference source, thus we expect no or little effects from its source structure.
As shown in Figure 3, few maser features can be detected with VERA alone.Among them, MF1 is the brightest maser feature, which can be detected over one year.The maser feature position can not be well determined using only one maser spot, because the maser spots are blended in position and frequency and the intensity of the blended components can vary with time (as shown in Figure 15).Instead, we used positions in about 5 spectral channels (velocity spacing of ∼0.5 km s −1 ) where the maser spots were at the peak brightness, with compact emission and minimal blending.The positions of the maser spots relative to the calibrator J0524+7034 were measured as a function of time, and the positions were then modeled by the parallax, proper motion, position at a reference epoch, and the possible acceleration.We then assess the astrometric quality based on the magnitude of the post-fit residuals.The formal position uncertainties were usually smaller than the actual error, since the unknown systematic error (i.e., the tropospheric delay residual, maser structures, etc.) often dominate over random noise (Reid et al. 2009;Nagayama et al. 2020).The systematic error was added in quadrature to the priori error so that the reduced chi-square value became nearly unity.Figure 7 and Table 3 show the outcomes of the astrometric fitting for MF1 positions measured with the VERA astrometric technique (Nagayama et al. 2020).In this astrometric fitting, we compared the different parameters, such as, using the data at different timescales and including the acceleration or not.We found different deviations from a common constant proper motion, or an acceleration, in declination between 2020 and 2021.This could come from either the variation in the maser structure (Figure 16) or a non-uniform acceleration in the maser feature.Taking into account the ambiguity in the results from these different estimation procedures, we conservatively obtained a parallax value of 1.79±0.08mas.This is consistent with 1.76±0.10mas from Gaia EDR3 (Gaia Collaboration et al. 2021) and 1.73±0.03mas from the previous VERA result (Matsuno et al. 2020) within 1 sigma level.The proper motions of the individual maser features usually are different from the central star, due to the extra internal motions of the masers.However, the proper motion of the maser feature MF1 (14.37±0.20,-35.46± 0.44 mas yr −1 ) is consistent with Gaia EDR3 (14.29±0.07,-34.63± 0.10 mas yr −1 ) within 1 sigma in east direction and 2 sigma in north direction, due to a small projection distance (less than 7 mas as discussed in Section 3.3.3)between MF1 and central star.

Astrometric fitting for maser features with the EAVN data
We can obtain the absolute positions of the bright and compact maser spots from the 14 epochs with the VERA data.The bright and compact maser feature (MF1 in φ 1.89 ∼ 2.29 and φ 2.93 ∼ 3.03 , MF4 in φ 1.04 ∼ 1.19) was then used to register the VERA and EAVN maps.The registration accuracy can be expressed as σ registration = σ thermalVERA 2 + σ structureVERA 2 + σ thermalEAVN 2 + σ structureEAVN 2 , where σ thermal is the thermal (random) noise determined by an elliptical Gaussian fitting with the AIPS task "SAD" and σ structure is the error effected by maser structure estimated by the blended peak components as shown in Figure 15 in the VERA/EAVN images.Because MF1 and MF4 are bright and compact, the thermal noises are < 0.02 mas and structure errors are ∼ 0.05 mas for both VERA and EAVN maps.We believe the registration errors between VERA and EAVN maps are typically ∼ 0.1 mas.These registration errors are used as the prior errors in the astrometric fitting.Note that this uncertainty may be propagated to the errors in deriving the parallaxes of maser features with the astrometrically registered EAVN maps, as described in Section 4.1.
Table 4 shows the independent parallax fits for the residuals of maser feature positions derived after subtracting the modulation of the annual parallax derived from the VERA data.The results were based on the same fitting procedures as mentioned in Section 3.3.1 and with/without acceleration parameters were compared.Averaging all the fitting results, we obtained a combined parallax residual of 0.07±0.06mas, which is consistent with the estimated error of the parallax (π = 1.79±0.08mas).

The stellar position in the H2O maser distribution
Determining the stellar position with respect to the H 2 O maser distribution is important to estimate the center of 3-D kinematics.Gaia EDR3 provided the stellar position with optical astrometry, however there may exist uncertainty for this AGB star (Xu et al. 2019).On the other hand, the center of the SiO maser ring should correspond to the position of the star.Figure 8 shows the registered map of the H 2 O and SiO masers with the KVN data, to which the SFPR calibration technique was applied (Dodson et al. 2014).The positional reference is the H 2 O maser at V LSR = −11.05km s −1 (MF1).We obtained the position of the central star (−5.1 ± 0.3 mas, −0.8 ± 0.2 mas) as indicated as the center of the solid circle in the least-squares circle method (Yoon et al. 2018).The estimated position of the central star is consistent with Gaia EDR3 result (−6.5 ± 0.3 mas, +1.5 ± 0.4 mas) with an uncertainty of 2.7 mas on 2020 January 10.The differences between the estimated position and Gaia EDR3 result (−6.5 ± 0.3 mas, +1.5 ± 0.4 mas) are 1.4 mas (∼3 sigma) in east direction and 2.3 mas (∼5 sigma) in north direction on 2020 January 10.Because of variable intervals of the observation epochs, including large blanks of the observations, the final version of the animation will be edited with some epoch-interpolation technique (Gonidakis et al. 2013).Even the current version of the animation, one can recognize radial outward motions of the H 2 O maser features, helping us for noting interesting properties of the maser motions described as follows.

The expansion of the BX Cam flow traced by H 2 O masers
Figure 10 and 11 present the combined cube of the maser features synthesized with the cubes from the "High-Res.Imaging" in Table 1.The constant velocity proper motions and spoke-like maser features in BX Cam "point back" towards the central star.However, the point-back directions of the maser features in different groups (e.g, group MF[7,8,9,10,11,12] and group MF[13,14,15]) may converge into different areas, which may also have offsets from the location of the central star and the ring of SiO masers, as shown by auxiliary dashed lines.Thus the H 2 O maser kinematics may have origins different from the central star, if we assume that masers have the constant proper motions.
Figure 12 shows the angular distance-LSR velocity diagram of the H 2 O maser features.To examine the variation of maser flows on long timescales, we compared the maser maps made by VERA from 2012 to 2014 (Matsuno et al. 2020) and made by ESTEMA from 2018 to 2021 (Figure 10). Figure 13 shows the registered map of the ESTEMA maser flows with that of Matsuno et al. (2020).Adopting a radial velocity of BX Cam as 0 km s −1 (Matsuno et al. 2020), we derived the three dimensional expanding velocities of H 2 O maser features using V LSR and internal proper motions.The maser map of Matsuno et al. (2020) showed only the blueshift dominant features with respect to the central star.They obtained a three-dimensional velocity of 14.8±1.4km s −1 with three collimated flows, and concluded that the H 2 O masers trace an outflow with a quite uniform velocity.However, our result shows a wider range of 9 km s −1 (inner radii) to 19 km s −1 (outer radii) at different directions.These differences may be caused by the different observational epochs (different pulsation phases) originated from characteristics of the Mira variable star BX Cam, which traces the different expanding shells relative to Matsuno et al. (2020).b The residual parallax after subtracting the parallax 1.82±0.03mas for with acceleration(Yes) and 1.90±0.02mas for without acceleration(No).
d The possible accelerations/decelerations.
e The root mean square (RMS) of the residuals in RA and Dec.
f The unweighted average and standard deviation.
Base on Figure 9-13, there are three main groups of H 2 O maser features.The first one (Group A: MF [1,2,3,13,14,15,17]) is associated with the inner blue-shifted maser feature closer to the central star, the second (Group B: MF[5,6,7,8,9,10,11,12]) is associated with the outer blue-shifted lobe of the outflow, and the third (Group C: MF4) is associated with the red-shifted maser feature of the outflow.The red-shifted maser feature of the outflow (from behind the central star) is very time variable as shown in Figure 1.The three-dimensional velocities are estimated to be 12.8±3.5 km s −1 for all maser features, 10.5±2.0 km s −1 for Group A, 15.5±2.8km s −1 for Group B, and 10.1 km s −1 for Group C (MF4). ) masers observed on 2020 January 10 using the KVN data, which were taken from velocity-channel integration with weights based on the root-mean-square noise level.Black star indicates the estimated position of the star at the center of the SiO maser ring using the least-squares circle method (Yoon et al. 2018).Red star indicates the position of the star derived from Gaia EDR3.The contour map is drawn with contour levels of 3,6,12,24,...×σ, where 1σ corresponds to 0.025 Jy beam −1 (H2O) and 0.120 Jy beam −1 (SiO), respectively.
We also made a least-squares model-fitting analysis assuming a model of an outflow that is radially expanding at a constant velocity.The detail of the modeling has already been described in Zhang et al. (2012).We adopted weights proportional to the square of the accuracy of a measured proper motion.The expansion center was assumed to be at the ring center of the SiO masers whose map registration was yielded in KVN SFPR (Section 3.3).The systemic velocity is adopted to be ∼0 km s −1 (Matsuno et al. 2020).Thus the expansion velocity of H 2 O masers is estimated to be 13.0 ± 3.7 km s −1 , which is consistent with the result (14.8 ± 1.4 km s −1 ) that found in the VERA observations during 2012-2014 (Matsuno et al. 2020).Figure 14 shows the spatial distribution and velocity vectors of the maser features with respect to the outflow origin projected on the R.A.-Distance (XZ) and Distance-Dec (ZY) planes, which was obtained after the model fit assuming a radially expanding flow (Imai et al. 2000).The outer masers roughly exhibit radial expansions and are located at similar distances from the central star at roughly a constant expansion velocity.Thus the masers seem to be associated with a spherically expanding outflow although the directions of the maser excitation are significantly biased and changed from time to time on the timescale of several years.However, there exist some big deviations from the spherical flow for the inner masers projected closely to the central star in the sky, which may be caused by the result of enhancement of large uncertainties in the maser positions on Z-axis and some random maser motions in the modeling.Our work has evaluated the accuracy of the annual parallax measurement with CSE H 2 O masers, majority of whom are located in the stellar neighborhood and exhibit extended brightness distributions.The importance of astrometry is also demonstrated for the maser animation, which is a key to register the maser images at different epochs and extract the intrinsic motions of maser features associated with the CSE.The astrometric accuracy of the VERA was evaluated in detail by Nagayama et al. (2020), which achieves the parallax accuracy of 10 µas in the best cases.Our targets were observed with a background calibrator 2 • away, and  2.
at elevations of ∼ 50 • .The single-epoch position accuracy for a given baseline can be estimated to be ∼80 µas using the Table 1 of Nagayama et al. (2020), indicating an expected parallax uncertainty better than 30 µas with 10 epochs extended over more than one year.However, our derived parallax (1.79±0.08 mas) has an uncertainty more than twice Xu et al. the expected one, probably due to the variation in the maser structure (MF1 in Figure 16) or an irregular motion due to shocks or radiation-coupled flow as shown in the Figure 7. Nevertheless it is consistent with that previously derived with VERA (1.73±0.03mas) (Matsuno et al. 2020) and Gaia EDR3 (1.76±0.10mas) (Gaia Collaboration et al. 2021) within 1 sigma.These parallaxes also rule out a larger parallax (4.13±0.25 mas) reported in Gaia DR2 (Gaia Collaboration et al. 2018).There is evidence that there may exist uncertainty in the Gaia DR2 parallaxes for the AGB stars from their large sizes, temporal and spatial variation in the surface brightness, and/or extinction by the circumstellar dust (Xu et al. 2019).These factors may prevent the improvement of the parallax accuracy even in future Gaia results using more measurements and better calibration of systematics.As already shown in Section 3.3, the individual maser features also exhibit some deviations from linear motions.The smaller residuals in Figure 7 support the parallax fitting including such an acceleration.However, it is obscure whether the observed acceleration really reflects a real one rather than just an illusion caused by a drift of a maser emission area in the maser gas clump.This should be tested with further analysis of the observed accelerations in the combined maser features around this star.The residual parallaxes of all maser features found in the EAVN maps are probably affected by their extended structures, the errors in the astrometric registration of the EAVN and VERA maps, and/or the possible non-linear motions.Note that an additional error contribution as mentioned above may have possibilities to increase and decrease the parallax value.Therefore, we adopted a parallax value using MF1 only (Section 3.3.1)and give a conservative uncertainty that would cover all of the possible errors (Section 3.3.2).
The stellar position in the H 2 O maser distribution was estimated to be at the ring center of the SiO masers whose locations are determined in the KVN astrometry with the SFPR technique at a single epoch.The SiO maser ring structure will be further clarified with multiple maser lines (Yoon et al. 2018).Even the present results, one can clearly see an offset of the expansion center of the H 2 O maser flow from the center of the SiO maser ring as shown in Figure 10.This may be caused by either the uncertainty of the stellar proper motion measured in Gaia EDR3 and adopted in the maser map registration, or a possible deviation of the outflow from a radial expansion.

CSE spatio-kinematics: asymmetric maser distribution and its evolution
The asymmetry/inhomogeneity of the spatial distribution is evident throughout the three-dimensional maser kinematics of H 2 O masers that we have revealed.The maser outflows are dominated by the blue-shifted H 2 O maser features (dominant blue-shifted outflow of bipolar outflow).Compared to the maser features detected with VERA , ZYplane) of the estimated positions and proper motions with respect to the outflow origin, using the method of Imai et al. (2000).The base of the arrow is pointed to the position of the maser feature.The direction and length of the arrow represent the direction and magnitude of the proper motion, respectively.The red and blue arrows show the maser feature motions found in the EAVN (ESTEMA) and VERA observations (Matsuno et al. 2020), respectively.from 2012 to 2014 (Matsuno et al. 2020), we have detected the central and southern components, and also new redshifted components.The northern components previously mapped by Matsuno et al. (2020) have become so faint that they are detected only on KVN baselines, at a few epochs.
Although the velocity field of the masers is generally consistent with the expanding shell model, the expansion is not strictly isotropic.For BX Cam, the distribution of the maser features is limited in a thick-shell region (Section 3.4), which has inner and outer expansion velocities of 9 km s −1 and 19 km s −1 at radii of <20 mas and 60 mas, respectively.There also exist gradual drifts in radial velocities and possible accelerations/decelerations in the proper motion.Here, we suppose a constant/gradual radial acceleration in a CSE (e.g.Richards 2012) for BX Cam.Adopting the distance to the star D =560 pc, one can calculate a travel time from the inner and outer radii of the expanding maser shell to be t travel 7.6 yr and an acceleration to be a 1.3 km s −1 yr −1 , corresponding to μ 0.5 mas yr −2 .This suggests that observed drifts of LSR velocities can be explained by this kind of acceleration while the accelerations/decelerations seen in the maser proper motions (up to 1.8 mas yr −2 ) cannot.Because the accuracy of the observed accelerations/deceleration is currently insufficient for further astrophysical interpretation, we reserve further statistical analysis for a separate paper.Nevertheless, the existence of such accelerations or deviations from constant velocity motion, may be detectable in the systematic residuals after the parallax and proper motion fitting.
The expansion velocity and the projected distance to the central star of a mass-losing flow are seem to be timedependent.The flows found in the ESTEMA data were compared with that in the VERA observations from 2012 to 2014 in Figure 13.Different shells and distributions were detected at different times and with different arrays.
Temporal variation in the morphology of individual maser features is also determined.The shapes of the maser structures are generally correlated with the direction of expansion.Figure 16 also shows the complex variation of the interior structure, such as interior rotation, expansion, sub-structure, and change of direction.This might be due to any or all of: bulk rotation, internal turbulence, or a pattern speed due to shock wave propagation (Richards et al. 2012(Richards et al. , 2020)).

Variability of the H 2 O maser
A correlation between the 22 GHz H 2 O maser and optical light curves with a time lag has been reported in AGB stars (Benson et al. 1993).BX cam is in also such a case, and the measured time lag between the integrated water maser and optical light curves is about ∼25 d.However, several bright masers often dominate the total integrated intensity, and thus it is hard to depict the variation of individual maser features or spots without spatially resolved maser maps.Our VLBI images allow us to trace individual maser features (shown in Figure 15) and groups of maser features (shown in Figure 5).We find that the light curves of the individual maser features or groups have different time lags with the optical light curve.We note that H 2 O masers are pumped by collision rather than radiation field, and hence other factors should be considered to interpret the correlation between maser and optical light curves, such as shock waves caused by stellar pulsation or giant planets, revolving the central star (Struck-Marcell 1988;Rybicki & Denis 2001;Richards et al. 2012).

SUMMARY
We have measured the trigonometric parallaxes, proper motions, drifts in LSR velocities, and the possible accelerations/decelerations of the H 2 O masers around the Mira variable BX Cam using EAVN observations.Our parallax is 1.79 ± 0.08 mas, which is consistent with Gaia EDR3 and previously measured VLBI parallaxes within their joint uncertainties.The stellar position with respect to the H 2 O masers is registered by a ring of SiO masers using the KVN SFPR astrometry, which is consistent with that of Gaia EDR3 within 3 mas.There are asymmetries/inhomogeneities in both the spatial and velocity distributions of the H 2 O masers although their motions are roughly modeled with a spherically expanding outflow.In fact, the maser locations have been dominant in a blue-shifted outflow with respect to the central star.The 3D H 2 O maser kinematics indicates that the circumstellar envelope is expanding at a velocity of 13 ± 4 km s −1 .A time lag of about 25 d is also measured between the integrated water maser and optical light curves.More detailed maser kinematics for the CSE around BX Cam will appear in separate papers, covering multiple maser-line movies with SiO masers and the further discussion on the possibility of the maser clump acceleration.

bFigure 1 .
Figure1.The cross power H2O maser spectra on the Ky-Ku baseline (305 km).Few epochs on Ky-Kt baseline (476 km) when Ku is not available.All of the 37 epochs data were used.

Figure 2
Figure2.Integrated flux curves of 22 GHz H2O masers from the cross (red circles) and total (green circles) power spectra of the KVN and V-band optical light curve (blue stars) from the AAVSO monitoring data.The total power spectra is using the weighted average of three KVN antennas (Kt,Ky,Ku), while the cross power spectra is using the short baselines as shown in Figure1.The V LSR velocity range for estimating the integrated flux densities is from −16 to −7 km s −1 and from 7 to 10 km s −1 .At K band, the beam size of the KVN single-dish telescope is ∼130 as and the synthesized beam size of the projected KYS-KUS baseline is ∼11 mas for BX Cam.The black solid line shows the sinusoidal fit to the optical light curve with a period of 440 d.

Figure 3 .
Figure 3. Left to Right: Moment-zero maps of the water maser emission in BX Cam from KVN, VERA, and EAVN (KVN+VERA+T6) observations at the same epoch, Jan 10, 2020.The maps are from images created by collapsing each spectral line cube to a single plane using the AIPS task SQASH, taking the weighted averages of each pixel.An ellipse in the bottom-left corner of each panel indicates the synthesized beam pattern of the image synthesis.The contour map is drawn with contour levels of 3,6,12,24,...×σ, where 1σ corresponds to 0.03 Jy beam −1 (KVN), 0.04 Jy beam −1 (VERA), and 0.03 Jy beam −1 (EAVN), respectively.

Figure 4 .
Figure 4. EAVN images of the spatial structures of maser features (top left), with a zoom-in of the colored blocks (top right, bottom left, bottom right), observed on 2020 Jan 10 (φ 2.00).These images are obtained after component fitting from the image data cube with AIPS task SAD.Each of the filled circles shows a velocity component (maser spot), and each maser feature (MF#) has a velocity channel spacing of 0.105 km s −1 .The position reference feature (MF1) is around the map origin.

Figure 5 .
Figure 5. Temporal variation of maser features in different regions.The locations of the different regions are shown in Figure 4, while Region 1 includes MF[1,2,3,17], Region 2 includes MF[4,13,14,15], Region 3 includes MF[5,6], Region 4 includes MF[7,8,9,10,11,12].The velocity gradient is classified by color and the size of maser spots means the flux density.The color and size of maser spots are the same to Fig 15.Optical light curve found in the AAVSO photometric data (blue stars).

Figure 4 .bd
Absolute proper motions are defined as µx = µ α cos δ and µy = µ δ .c The fitted accelerations/decelerations.The root mean square (RMS) of the residuals in RA and Dec. e The unweighted average and standard deviation.

Figure 7 .
Figure 7. Parallax fitting for MF1 with VERA astrometry.Top panel: An example of the fitted parallax curve using the data points during the period 2019.8-2021.3without an acceleration.Bottom-left panel: Residuals of the fitting using the data points during the period 2019.8-2021.3.The data points from the fitting including an acceleration have a constant horizontal offset from those without an acceleration for clarity.Bottom-right panel: Same as the bottom-left panel but the period 2020.2-2021.3.

Figure 9
Figure9displays the astrometric animation of H 2 O masers at the 33 epochs from 2018 May 24 to 2021 Jun 7. Because of variable intervals of the observation epochs, including large blanks of the observations, the final version of the animation will be edited with some epoch-interpolation technique(Gonidakis et al. 2013).Even the current version of the animation, one can recognize radial outward motions of the H 2 O maser features, helping us for noting interesting properties of the maser motions described as follows.Figure10and 11 present the combined cube of the maser features synthesized with the cubes from the "High-Res.Imaging" in Table1.The constant velocity proper motions and spoke-like maser features in BX Cam "point back" towards the central star.However, the point-back directions of the maser features in different groups (e.g, groupMF[7,8,9,10,11,12]  and groupMF[13,14,15]) may converge into different areas, which may also have offsets from the location of the central star and the ring of SiO masers, as shown by auxiliary dashed lines.Thus the H 2 O maser kinematics may have origins different from the central star, if we assume that masers have the constant proper motions.The CSE of BX Cam were reported with the LSR velocity range from -5 to +5 km s −1 for SiO masers(Matsuno et al. 2020), from -16 to -4 km s −1 and from +7 to +10 km s −1 for H 2 O masers (Figure12), from +12 to +17 km s −1 for OH masers(Slootmaker et al. 1985) , and up to 21 km s −1 for CO J = 1 → 0 line(Knapp & Morris 1985), which shows a gradual acceleration from SiO maser to CO J = 1 → 0 line.Figure12shows the angular distance-LSR velocity diagram of the H 2 O maser features.To examine the variation of maser flows on long timescales, we compared the maser maps made by VERA from 2012 to 2014(Matsuno et al. 2020) and made by ESTEMA from 2018 to 2021 (Figure10).Figure13shows the registered map of the ESTEMA maser flows with that ofMatsuno et al. (2020).Adopting a radial velocity of BX Cam as 0 km s −1(Matsuno et al. 2020), we derived the three dimensional expanding velocities of H 2 O maser features using V LSR and internal proper motions.The maser map ofMatsuno et al. (2020) showed only the blueshift dominant features with respect to the central star.They obtained a three-dimensional velocity of 14.8±1.4km s −1 with three collimated flows, and concluded that the H 2 O masers trace an outflow with a quite uniform velocity.However, our result shows a wider range of 9 km s −1 (inner radii) to 19 km s −1 (outer radii) at different directions.These differences may be caused by the different observational epochs (different pulsation phases) originated from characteristics of the Mira variable star BX Cam, which traces the different expanding shells relative toMatsuno et al. (2020).

Figure 8 .
Figure8.Registered map of the H2O and SiO(v=1, J=1→0) masers observed on 2020 January 10 using the KVN data, which were taken from velocity-channel integration with weights based on the root-mean-square noise level.Black star indicates the estimated position of the star at the center of the SiO maser ring using the least-squares circle method(Yoon et al. 2018).Red star indicates the position of the star derived from Gaia EDR3.The contour map is drawn with contour levels of 3,6,12,24,...×σ, where 1σ corresponds to 0.025 Jy beam −1 (H2O) and 0.120 Jy beam −1 (SiO), respectively.
Astrometry for CSE H 2 O masers

Figure 9 .
Figure 9. Astrometric Animation of H2O masers around BX Cam.The coordinate system is astrometrically fixed with the star whose position is determined using Gaia EDR3.Maser feature motions on the sky with modulation due to the annual parallax (1.79±0.08 mas) and the proper motion (14.37±0.20,−35.46 ± 0.44 mas yr −1 ) of the referenced maser feature MF1 subtracted.A red-filled circle denotes the central star whose position (parallax: 1.76±0.10mas, proper motion: 14.29±0.07,−34.63 ± 0.10 mas yr −1 ) has been determined in Gaia EDR3.The lower panel is the optical light curve of AAVSO and the red vertical line follows the phase.An animated version of this figure is available.It covers 33-epoch maser maps from 2018 May 24 to 2021 Jun 7. The video duration is 16 seconds.

Figure 10 .
Figure 10.Expansion of the BX Cam flow traced by H2O masers.The positional reference is Gaia EDR3 (red star).Black circle denotes the fitted ring of SiO masers in Section 3.3.3.The different colors of the maser spots indicate different observation epochs, which are astrometrically fixed using the method in Figure 9.An auxiliary dashed line indicates a possible trajectory of an maser feature from 2012 to 2022 at a constant velocity vector, while the nearby lines in light grey show the position error of the back-extrapolated motion vectors.The lines of motion are drawn with the positions and proper motions of the maser feathers estimated in Table2.

Figure 11 .
Figure 11.Same as Figure 10 but the different colors of the maser spots indicate different LSR velocities.

Figure 12 .
Figure 12.Angular distance-LSR velocity diagram.The position of the star was used as the center of the SiO maser ring.

Figure 13 .
Figure 13.The registered maser flows from EAVN (ESTEMA) and VERA (Table 4 of Matsuno et al. 2020) observed at different pulsation circles.The registration is performed with Gaia EDR3.The numbers around the maser features are the three-dimensional velocities in km s −1 .
Figure 14.The Top view (Left panel: R.A.-Distance, XZ-plane) and the East-side view (Right panel: Distance-Dec., ZYplane) of the estimated positions and proper motions with respect to the outflow origin, using the method ofImai et al. (2000).The base of the arrow is pointed to the position of the maser feature.The direction and length of the arrow represent the direction and magnitude of the proper motion, respectively.The red and blue arrows show the maser feature motions found in the EAVN (ESTEMA) and VERA observations(Matsuno et al. 2020), respectively.

Figure 15 .
Figure 15.Temporal variation of the internal maser features.The size of the circle indicates the logarithm of the flux density of the maser spot (see the fourth sub-panel for the each maser feature as a legend).The color of the circle indicates the LSR velocity of the spot (see the third sub-panel).

Table 1 .
EAVN observations of BX Cam at K band Num.d RMS e High-Res.d VERA d Problem f

Table 2 .
The radial velocity drifts and relative proper motions of the maser features MF a Cycle b Radial Velocity c

Table 3 .
Astrometric fitting for MF1 with VERA

Table 4 .
Astrometric fitting for positions of the maser features with EAVN, registered using VERA astrometry MF a