Evolution of the Outflow in the Water Fountain Source IRAS 18043–2116^*

Single-dish monitoring observations of H₂O and SiO masers toward WFs were conducted in the “Finest Legacy Acquisitions of SiO- and H₂O-maser Ignitions by the Nobeyama Generation” (FLASHING) survey. FLASHING has been described in detail in Imai et al. (2020) and Amada et al. (2022), including the specification and the record of the observations with the NRO 45 m Telescope since 2018 December. The summary paper of FLASHING will be published in a separate paper. Here we focus on the FLASHING observations of H₂O (at the rest frequency of 22.235080 GHz) masers in I18043 conducted during 2019 January–2020 April, including the epochs of the new discovery of the highest-velocity components and those within 1 month of the ATCA observations, described later. Each of the observations lasted for about 1 hr.

We used three spectral windows to obtain a total velocity coverage of ∼1600 km s⁻¹. Each spectral window had 2048 spectral channels, each with a velocity resolution of 0.41 km s⁻¹. Data reduction was carried out using the JavaNewstar package ⁹ in a standard manner, namely, data integration in the time domain followed by spectral baseline fitting in the maser emission-free spectral channels. We adopt a conversion factor of 2.8 Jy K⁻¹ to convert the antenna temperature to the flux density scale.

The ATCA observations presented in this paper were conducted on 2020 April 7 and 30, 2020 June 4, and 2021 March 30 within project C3361, which is a follow-up of the sources monitored in the FLASHING program in order to increase the temporal sampling of spectra, and to obtain information on the evolution of the spatial distribution of the H₂O masers. Two independent spectral windows were observed, centered at 18 and 22 GHz, with a bandwidth of 2 GHz each, obtaining full linear polarization products. We set up the Compact Array Broadband Backend in its 64M-32k mode, which samples each of these spectral windows in 32 broadband channels per polarization with a coarse spectral resolution of 64 MHz, for radio continuum observations. Several of the broadband channels were further zoomed in to observe the H₂O maser line with a finer spectral resolution of 0.42 km s⁻¹. The total velocity coverage for the line was 1727 km s⁻¹ on April 7 and 3453 km s⁻¹ in the later three sessions.

Initial calibration was carried out using standard procedures in Miriad. The absolute flux density scale was calibrated with the source PKS 1934−638. PKS 0537−441 was also used as a bandpass calibrator for the first epoch, while PKS 1613−586 was used for the other three epochs. The sources IERS B1817−254, PMN J1832−2039, IERS B1817−254, and PMN J1755−2232 were used as complex-gain calibrators in each of the four epochs, respectively. Further processing (imaging, deconvolution, and self-calibration) was carried out with the Astronomical Image Processing System (AIPS). For the continuum data, we flagged out the broadband channels containing the maser line, to avoid contamination. In the line data, the spectral channel with strongest maser emission was used for self-calibration as well as reference position for the maps. The phase and amplitude corrections for self-calibration were obtained at 10 s intervals (the integration time of each individual visibility). These solutions were then applied to both the whole spectral data set and the broadband continuum data centered at 22 GHz. Since the line and continuum data share the same calibration corrections, the 1σ relative positional accuracy among individual maser components and between maser and continuum emission is approximately given by the size of the synthesized beam divided by 2 times the signal-to-noise ratio (S/N) of the emission.

We note that the beam is highly elongated along the north–south direction, which strongly worsens the positional accuracy in that direction. The synthesized beams were 253 × 034 at a position angle of −6.78° and 142 × 041 at a position angle of −1.27° on 2020 April 30 and 2021 March 30, respectively. In these two sessions only, the uv coverage was good enough to map the maser emission while the spectra were obtained for all the four sessions previously mentioned. We roughly estimate a positional accuracy of ∼20 mas (1σ) along the minor axis of the beam (almost in the east–west direction) for a component with an S/N ∼ 10. The relative positions of the maser components with respect to the reference one were determined by Gaussian fitting using the AIPS task JMFIT, assuming a single Gaussian component in each velocity channel. This simple assumption looks valid for distinguishing clusters of the maser components described in Section 3.2.

I18043 was also observed with ALMA on 2019 January 18 (project code: 2018.1.00250.S). The details of these observations are presented by Khouri et al. (2021). The spectral setup covered a frequency range that contains emission of the CO(J = 2 → 1) line at 230.538 GHz. In this paper, we present only the spectrum and image of the CO(J = 2 → 1) in order to estimate the systemic velocity of I18043 and compare the spatio-kinematical distribution of the masers with that of the molecular material. The integration time of ∼11 min yielded an rms noise level of 3 mJy in the spectrum for a velocity resolution of 1.5 km s⁻¹, and a synthesized beam of 124 × 118 with a position angle of −78°. The relative position accuracy is ∼60 mas for a component with an S/N ∼ 10.

Figure 1 shows a time series of the spectra of H₂O masers of I18043 taken with the NRO 45 m telescope and ATCA during 2019 January 2–2021 March 30. Because of the rapid variation in the spectral profile from one observation epoch to another, it is difficult to directly check the consistency of the flux density scales between the spectra of the NRO 45 m telescope and ATCA. However, the maser components around V_LSR ∼82 km s⁻¹ had flux densities of ∼4 Jy in both spectra of the two telescopes during 2020 April 7–16, suggesting that the flux density scales are roughly consistent with each other within ∼10%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The FLASHING session on 2019 January 2 detected the redshifted side of the new velocity components at V_LSR ≃ 376 km s⁻¹, while the ATCA session on 2020 April 30 detected the blueshifted side of the new velocity component at V_LSR ≃ −165 km s⁻¹. These detections break the record of the velocity spread of the maser components in this source (−100 km s⁻¹ > V_LSR > 350 km s⁻¹; Walsh et al. 2009; Pérez-Sánchez et al. 2017); this suggests that the top speed of the jet in I18043 has to be revised. Moreover, this velocity spread (≃540 km s⁻¹) is the largest ever detected in the spectra of H₂O masers of WFs, beating the previously reported one for IRAS 18113−2503 (≃500 km s⁻¹; Gómez et al. 2011).

It is difficult to determine whether these extreme velocity components may be decelerating, as suggested by Orosz et al. (2019) for IRAS 18113−2503, due to their too short lifetimes within a few months. However, their seemingly constant velocities within ∼1 km s⁻¹ on this timescale may rule out such rapid deceleration (up to 10 km s⁻¹ month⁻¹).

Figure 2 shows the maps of H₂O masers taken in I18043 with ATCA on 2020 April 30 and 2021 March 30. Although some of the maser components have relatively large positional uncertainties (up to 500 mas) in the decl. offset, almost parallel to the elongation direction of the synthesized beam, the maser map resolves the structure of the bipolar outflow in the east–west direction. One can see two main clusters of maser components at (A) −75 ≤ V_LSR ≤ 83 km s⁻¹ and (B) 83 ≤ V_LSR ≤ 188 km s⁻¹, as well as another blueshifted cluster (C) −113 ≤ V_LSR ≤ −98 km s⁻¹. An additional redshifted group (D) 206 ≤ V_LSR ≤ 209 km s⁻¹ is also tentatively visible to the west of cluster B. The most redshifted component (V_LSR ∼ 344 km s⁻¹) also may be included in the group D while the most blueshifted component (V_LSR ∼ −165 km s⁻¹) can be included in the group C, taking into account their positional uncertainties mentioned above.

Zoom In Zoom Out Reset image size

Figure 2 also shows the location of the continuum source (X, Y) = (15 ± 23, 165 ± 131) (mas) on 2020 April 30 and (X, Y) = ( −13 ± 4, 270 ± 16) (mas) on 2021 March 30. The continuum source had a flux density of S_ν = 0.75 ± 0.05 mJy on 2020 April 30 and S_ν = 0.97 ± 0.09 mJy on 2021 March 30. The location of the continuum source, detected at S/N ≃ 36, also may have a large positional uncertainty related to the synthesized beam shape. Nevertheless, it is clear that the continuum source is located between the blue- and redshifted H₂O maser clusters.

These ATCA maps taken in 2020–2021 can be directly compared with Figure 3 of Pérez-Sánchez et al. (2017) taken in 2013. Although it is difficult to quantitatively compare the separations of the maser clusters in a time span of ∼7 yr, one can see that the separation between clusters A and B looks consistent. Taking into account the stellar systemic velocity of V_LSR ≃ 87 km s⁻¹ (Deacon et al. 2004), clusters A and B may form a pair with spectral and spatial symmetry with respect to the central star. Regarding the stellar systemic velocity, we revised its value using CO emission detected with ALMA; see Section 3.3.

In the same manner, cluster C and group D will form another symmetric pair. However, cluster C and group D in 2020–2021 have large uncertainties in their locations and distributions. Group D in 2020–2021 may be coincident to cluster B in 2013 (Pérez-Sánchez et al. 2017), while cluster C has been newly visible ∼50 mas east of cluster A.

Figure 3 highlights the maser distributions during 2020–2021 in the east–west direction in which the relative maser positions were better determined. The separation between the two major clusters A and B is discussed in more detail in Section 4.

Zoom In Zoom Out Reset image size

Figure 4 shows the spectrum of the CO(J = 2 → 1) emission observed with ALMA. The line consists of a strong central component and high-velocity wings, which are indicated in the left panel of Figure 4. The high-velocity wings extend over a velocity range >200 km s⁻¹, although the blueshifted wing exhibits contamination of foreground gas that probably belongs to molecular clouds along the line of sight. From the velocity of the line peak of the central component, a systemic velocity of V_LSR,sys = 87 ± 1 km s⁻¹ is derived for I18043. This value is compatible with the systemic velocity obtained by Deacon et al. (2004) using OH maser observations.

Zoom In Zoom Out Reset image size

The spatial distribution of the CO(J = 2 → 1) and 230 GHz continuum emission is shown in Figure 5. There exists an offset of the brightness centroids between the blue- and redshifted components (<02), and the spatially unresolved brightness distributions of these components well support the far distance to this source and the compactness of the jet structure (≲11 kpc × 02 ∼2200 au).

Zoom In Zoom Out Reset image size

The locus of the CO(J = 2 → 1) emission peak in each bin with a velocity width of 6 km s⁻¹ is presented in Figure 6. There is a clear velocity gradient in R.A. offset with respect to the 230 GHz continuum peak, as shown in the right panel of Figure 6.

Zoom In Zoom Out Reset image size

We measured the separation between the two main clusters of H₂O masers, blue- and redshifted with respect to the systemic velocity ≃87 km s⁻¹. This separation of the two clusters was ∼45 ± 5 mas on 2008 July 12 (Walsh et al. 2009). On the other hand, it has changed to ∼74.1 ± 0.4 mas and 87.6 ± 0.5 mas on 2020 April 30 and 2021 March 30, respectively, as shown in Figure 3. Note that the uncertainties in these separations are statistical. This result may indicate a possible expansion of the cluster separation at a rate of ∼2.9 ± 0.5 mas yr⁻¹, namely a growth of a bipolar outflow hosting the two clusters.

We also note that the angular separations of the red- and blueshifted groups of H₂O masers in I18043, ∼74 and 88 mas, correspond to a linear length of ∼800 and 950 au, respectively. They are quite similar to the separation of the red- and blueshifted groups of masers in IRAS 16552−3050 (∼700 ± 100 au; Suárez et al. 2008; Vickers et al. 2015). On the other hand, the velocity spreads are rather different between these two WFs (∼540 and ∼170 km s⁻¹; Suárez et al. 2008). The origin of the difference in the velocity spreads (inclination of the jet major axis, stellar luminosity, evolution phase) will be further investigated.

In a previous analysis of Very Long Baseline Array (VLBA) data from 2008 October 26 and 2009 January 25 in I18043, Orosz (2017) was able to distinguish two distinct arc-shaped groups of masers, blue- and redshifted with respect to the systemic velocity of ∼87 km s⁻¹. The blueshifted maser emission was within a velocity range from −111 km s⁻¹ to 78 km s⁻¹ and the redshifted maser emission was within a velocity range from 94 km s⁻¹ to 176 km s⁻¹. Comparing these results with those by Walsh et al. (2009), it seems that the high-velocity components detected previously at velocities higher than 200 km s⁻¹ were not included in the VLBA epochs, with a velocity coverage up to 227 km s⁻¹ (Day et al. 2011). In our recent observations, we were able to detect higher-velocity components than those previously reported, in both the blue- and redshifted sides of the spectra at velocities up to V_LSR ≃ −165 km s⁻¹ and ≃376 km s⁻¹, respectively. This velocity spread is the largest ever detected in a WF ≃540 km s⁻¹. The emergence of the new highest-velocity components at the outer outflow lobes may indicate a rapid growth of the outflow triggered by an increase in the maximum outflow velocity.

From two epochs of the VLBA data, 2008 October 26 to 2009 January 25, Orosz (2017) estimated roughly the expansion rate of the outflow to be ∼1 mas yr⁻¹, with an average separation of the two maser clusters of ∼63 mas. If we extrapolate the expansion rate of the outflow from 2009 to 2020, the separation of the two clusters would be ∼74 mas, which is consistent with the results of the ATCA observations. Thus, the dynamical age of the outflow is estimated to be t_jet ≲ 30 yr, using the expansion rate and the separation of the two clusters from the VLBA data. Orosz (2017) used the full velocity spread of the H₂O masers seen in their VLBA maps (+176 to −111 km s⁻¹) and derived the inclination angle of the jet’s major axis to be i_jet ∼ 75° with respect to the plane of the sky. Using this inclination angle, we estimate a 3D expansion velocity for the redshifted side of V_jet,3D ∼ 170 or 300 km s⁻¹, when assuming mean proper motions of ∼1 and 2.9 mas yr⁻¹, respectively, considering the most redshifted components observed with ATCA of V_LSR ∼ 284 and 344 km s⁻¹, respectively.

The CO(J = 2 → 1) emission consists of a spatially extended low-velocity component, shown as green contours in Figure 5, and a compact high-velocity component, shown as blue and red contours. The green contours represent the velocity-integrated emission in the velocity range +65 to +110 km s⁻¹, and the blue and red contours represent the velocity-integrated emission in the velocity ranges −125 to +65 km s⁻¹ and +110 to +300 km s⁻¹, respectively. The peaks of the blue- and redshifted components have an offset in the east–west direction by approximately 02, with the blueshifted emission located toward the east. Thus, the velocity gradient seen in the CO emission is in agreement with the velocity gradient of the blue- and redshifted masers.

The right panel of Figure 6 shows the locus of the CO(J = 2 → 1) emission peak in each bin with a velocity width of 6 km s⁻¹. The observed dispersion of the locus in the R.A. offset indicates a wide angle of the outflow. We note that there exists a clear velocity gradient in the R.A. offset of the CO emission from −20 to 20 km s⁻¹ with respect to the systemic velocity. This velocity gradient is roughly followed by the blue- and redshifted CO emission, considering the uncertainties and a S/N of 10.

The wanders of the peak positions in the blue- and redshifted CO emission components, V_offset > 70 km s⁻¹, may be affected by the large positional uncertainty, while those in the range V_offset = 10–70 km s⁻¹ are roughly consistent with the ones of H₂O masers. This suggests that the H₂O masers are associated with the CO outflow itself, or possibly shock regions in the outflow. The absence of CO emission at the fastest velocity range ∣V_offset∣ > 100 km s⁻¹, where the H₂O masers are visible, is different from the case of W43A in which the H₂O masers are associated with the slower-moving entrained material rather than the faster jet (Tafoya et al. 2020). The CO velocity wanders also indicate a large opening angle of the outflow, within which the locations of different velocity components are scattered due to an inhomogeneity of the outflow velocity. Taking into account the direction coverage of the wandering-velocity components with respect to the originating point of the outflow, the full opening angle of the outflow is estimated to be ∼60°.

Also looking at the right panel of Figure 2 in more detail, the blueshifted components at V_LSR ∼ −104 to −101 km s⁻¹ present an offset from the two major clusters of masers A and B. This separation also may imply the presence of a wide-angle outflow, indeed supported by the CO(J = 2 → 1) observations.

Figure 7 shows a schematic picture of the WF I18043, including the current observational panorama discussed in the previous sections. The continuum emission, close to the central stellar system, may correspond to an ionized area at the inner part of a biconical outflow. Farther away from the central system, around several hundreds of au, thick biconical cavities surrounding the outflow are expected from the observed clusters of H₂O masers, which are likely excited in the cavities by shock interaction with the outflow. The blue- and redshifted H₂O masers may be more amplified along the front and rear sides of the inner cavities, where a longer amplification path is expected close to the line of sight. Moreover, the CO emission is tracing an envelope of several thousands of au. The velocity wanders in the blue- and redshifted components of the CO emission indicate a large opening angle of the biconical outflow (a full angle of ∼60° in each component). This is also consistent with a large scatter of the maser distribution in the direction perpendicular to the major axis of the outflow.

Zoom In Zoom Out Reset image size

Apart from the H₂O maser emission detected at 22 GHz, this source also presents H₂O maser emission at 321 GHz. These maser components span a velocity range similar to that at 22 GHz, indicating that these H₂O masers at both frequency transitions probably coexist. The intensity of the submillimeter masers is comparable to that of the 22 GHz masers, implying the kinetic temperature of the region where the submillimeter masers originate is T_k > 1000 K (Tafoya et al. 2018).

The origin of the radio continuum emission in this source is also relevant to understand its evolution. Pérez-Sánchez et al. (2017) obtained flux density data at four frequency bands between 1.5 and 22 GHz. They reported spectral indices (α, defined as S_ν ∝ ν^α ) < 2 between 1.5 and 5.5 GHz. Specifically, we estimated a value of 0.49 ± 0.09 from their data between 1.5 and 5.5 GHz, and α = −0.2 ± 0.2 between 5.5 and 22 GHz. These values are consistent with free–free radiation from a plasma, which ranges from α = 2 (optically thick regime) to −0.1 (optically thin), with an optical depth decreasing at higher frequencies. As mentioned by these authors, the spectral indices should be taken with care, since these observations were not simultaneous: the data at 22 GHz were taken in 2015, two years after those at lower frequencies, and the radio continuum emission in evolved stars is known to be variable (Cerrigone et al. 2011). Pérez-Sánchez et al. (2017) argued that if the radio continuum emission happened to fade with time, then the flux density at 22 GHz at the time of the low-frequency observations (2015) would have been lower, leading to a spectral index at high frequencies significantly lower than −0.1, which in its turn would point out to a nonthermal contribution. Actually, a decreasing flux density has been suggested in itself as being a signature of nonthermal radio continuum emission (Suárez et al. 2015; Cerrigone et al. 2017). However, following this line of reasoning, our data at 22 GHz do not show any evidence for a nonthermal component, since the radio continuum emission at 22 GHz seems to be consistently increasing: while a flux density of 0.70 ± 0.10 in 2013 July was obtained by Pérez-Sánchez et al. (2017), we measured S_ν = 0.75 ± 0.05 mJy on 2020 April 30 and S_ν = 0.97 ± 0.09 mJy on 2021 March 30. This is a factor 1.3 increase in one year in our observations, with an apparently higher increase rate in our observations than in the previous ones by Pérez-Sánchez et al. (2017).

Considering, then, that the current radio continuum data are compatible with thermal free–free emission from electrons in a plasma, we note that the spectral index α = 0.49 ± 0.09 in the Pérez-Sánchez et al. (2017) data between 1.5 and 5.5 GHz indicates a partially optically thick emission. It is also close to the standard value of n_e ∝ r⁻² expected in the case of an isotropic distribution of electrons with a radial variation of electron density of n_e ∝ r⁻², where r is the distance to the central star (Olnon 1975; Panagia & Felli 1975). This would be the case of an isotropic ionized wind or a photoionized region at the beginning of the PN phase (since the previously ejected envelope would also have an r⁻² density dependence). Under the assumption of an isotropic wind, Pérez-Sánchez et al. (2017) obtained a mass-loss rate of ≃3 × 10⁻⁵ (D/11 kpc) M_⊙ yr⁻¹.

However, I18043 clearly shows the presence of a collimated jet traced by H₂O maser emission, therefore assuming an isotropic ionized wind may not be appropriate. In this sense, Pérez-Sánchez et al. (2017) mentioned that the mass-loss rate they obtained under this assumption is an upper limit to the actual value. On the other hand, a collimated, ionized jet with n_e ∝ r⁻² would also show a similar spectral index in its radio continuum emission (Reynolds 1986), with a particular value depending on the jet geometry. In this model, a biconical jet (constant opening angle) would give α = 0.6, while lower or higher values correspond to a jet whose width decreases or increases with distance from the star, respectively. Here the jet half-width is given by w ∝ r , where is related to the spectral index as α = 1.3–0.7/ following Reynold’s model. The case of a standard spherical wind yields = 1 and α = 0.6. On the other hand, an isothermal, constant-velocity, fully ionized flow yields 0.5 < ≤ 1 and 0.6 ≥ α > −0.1. For the case of I18043, α = 0.49 and ≃ 0.89, therefore the radio continuum emission is congruous with a fully ionized collimated outflow.

In an ionized jet, the flux density of the radio continuum emission can provide an estimate of the mass-loss rate. Using Equation (8) in Anglada et al. (2018) and the parameters from the jet traced by H₂O masers, we estimate M_⊙ yr⁻¹ and ≃6 × 10⁻⁶ M_⊙ yr⁻¹, considering the flux density at 5.5 GHz and the spectral index α = 0.49 measured by Pérez-Sánchez et al. (2017) in 2013, and a jet velocity of 170 and 300 km s⁻¹, respectively. In this calculation we assumed an injection opening angle of the jet equal to 60°, a distance to the source of 11 kpc, an electron temperature of 10⁴ K, an inclination angle 75° with respect to the plane of the sky, and a turnover frequency of ≃5.5 GHz. On the other hand, we obtain the following values for the mass-loss rates from 2 to 4 × 10⁻⁶ M_⊙ yr⁻¹ when they are derived from the flux density measurements at 22 GHz in 2020−2021, although these are likely to be lower limits, considering that the continuum emission at this frequency may be in the optically thin regime (Pérez-Sánchez et al. 2017). This would suggest that, if the continuum emission arises from a jet, its mass loss is increasing. In any case, better estimation of the physical parameters would require an up-to-date measurement of the continuum spectrum of the source, with (nearly) simultaneous observations at a wide range of frequencies. Even so, these values of the mass-loss rates are not very accurate; they seem to be higher than those expected in models of post-AGB stars (≃10⁻⁷–10⁻⁸ M_⊙ yr⁻¹; Vassiliadis & Wood 1994), although these models follow the evolution of a single star. Lately, observations and models suggest that the presence of collimated jets in WFs could be due to a binary system in the center, leading to the formation of accretion disks from which the collimated outflow can be launched. In particular, these objects seem to have gone through a common-envelope phase (García-Segura et al. 2021; Khouri et al. 2021), which greatly enhances mass-loss rates.

An alternative explanation for the origin of the radio continuum emission is the onset of photoionization at which this source is entering the PN phase. An increase in the flux density at radio wavelengths has also been observed in several young PNe (Kwok et al. 1981; Knapp et al. 1995; Christianto & Seaquist 1998; Gómez et al. 2005; Tafoya et al. 2009; Cala et al. 2022). Such a flux density increase has been interpreted as an expansion of the ionized region. If this is the case, then an increase by a factor of 1.3 in flux density would correspond to an expansion by a factor of 1.14 in radius in one year for an isotropic ionized region. We should consider that the size of the emitting region is unresolved in our observations, therefore we are not able to estimate other physical parameters such as the mass of the ionized region. Furthermore, the scenario of a photoionized region around this source has also been discussed by Pérez-Sánchez et al. (2017), but they concluded that a shock-ionized wind is more likely to explain the observed radio continuum emission. Whether it is a shock- or radiative-ionized region around I18043, further observations would be necessary to identify the actual scenario.