Radio Continuum and Water Maser Observations of the High-mass Protostar IRAS 19035+0641 A

We present Very Large Array 1.3 cm continuum and 22.2 GHz H2O maser observations of the high-mass protostellar object IRAS 19035+0641 A. Our observations unveil an elongated bipolar 1.3 cm continuum structure at scales ≲500 au, which, together with a rising in-band spectral index, strongly suggests that the radio emission toward IRAS 19035+0641 A arises from an ionized jet. In addition, eight individual water maser spots well aligned with the jet axis were identified. The Stokes V spectrum of the brightest H2O maser line (∼100 Jy) shows a possible Zeeman splitting and is well represented by the derivatives of two Gaussian components fitted to the Stokes I profile. The measured B los are 123 (±27) and 156 (±8) mG, translating to a preshock magnetic field of ≈7 mG. Subsequent observations to confirm the Zeeman splitting showed intense variability in all the water maser spots, with the brightest maser completely disappearing. The observed variability in a 1 yr timescale could be the result of an accretion event. These findings strengthen our interpretation of IRAS 19035+0641 A as a high-mass protostar in an early accretion/outflow evolutionary phase.


INTRODUCTION
The typical formation time-scale of high-mass stars (M * > 8 M ⊙ ) is much shorter than their lower mass counterparts (∼ 10 5 versus ∼ 10 7 yr, e.g., André 1994;Sabatini et al. 2021), and they are still deeply embedded in their parental cloud by the time hydrogen burning in the core ignites.This makes it extremely challenging to observationally study the early stages of high-mass star formation.According to the turbulent core accretion model (e.g., McKee & Tan 2003), as the young stellar object (YSO) accretes material, jets of ionized gas are ejected from very near the protostellar surface.Even though these ionized jets play a key role in the high-mass star formation process, their nature is still poorly understood, mainly due to a lack of sub-arcsecond resolution radio continuum data.Recently, Rosero et al. (2016Rosero et al. ( , 2019) ) carried out a deep (rms noise ∼ 5 µJy beam −1 ), high angular resolution (θ ∼ 0 ′′ .33)continuum survey at 1.3 and 6 cm toward 58 high-mass star forming regions utilizing the Karl G. Jansky Very Large Array (VLA).They report the detection of 70 individual radio sources, including H ii regions, ionized jets, and ionized jet candidates.The latter were classified as such based on their rising spectral index (expected from ionized jets as per Reynolds 1986) and association with usual jet tracers, such as molecular outflows, mostly from singledish data.However, they were unresolved or slightly resolved in their observations, which made it impossible to completely discard ionization by a zero age main sequence star at the center of a hyper-compact (HC) H ii region as the origin of the radio continuum emission.Nonetheless, their statistical analysis of ionized jets and jet candidates provided significant further evidence supporting the common origin scenario for radio jets driven by YSOs of all luminosities (see Figure 8 in Rosero et al. 2019).This, in turn, supports the idea that high-and low-mass stars are formed in a similar manner, as predicted by the turbulent core accretion model.Ionized jets from such high-mass YSOs are expected to present an elongated morphology at scales of a few hundred au from the driving protostellar object, while H ii regions are more spherical.In order to investigate the shock ionization origin of the radio continuum sources classified as jet candidates on these scales, higher angular resolution observations are essential.Therefore, we carried out a 3× higher resolution (θ ∼ 0 ′′ .1)VLA survey at 1.3 cm toward all jet candidates from Rosero et al. (2016Rosero et al. ( , 2019)).The angular resolution of our data translates to linear resolutions between ∼ 150 and 1200 au (depending on the distance to each source), which will enable us to clearly discern whether the jet candidates are indeed jets or H ii regions.
As a byproduct of our radio continuum survey, we observed the 22.2 GHz (6 1,6 − 5 2,3 ) H 2 O (water) maser emission toward the jet candidates.Water masers are collisionally pumped in high density regions (n ∼ 10 8−10 cm −3 , Elitzur 1992) and are usually associated with jets and outflows from early stage high-mass YSOs (e.g., Moscadelli et al. 2019;Ladeyschikov et al. 2022).The detection of water masers in close proximity to the jet candidates would therefore strongly support the jet interpretation for the radio continuum sources.Moreover, water masers are highly compact and can be extremely bright, which makes them an ideal probe to study shock dynamics at sub-arcsecond resolution (e.g., Moscadelli et al. 2020).In addition, the Zeeman effect detected using H 2 O masers offers direct measurements of the line-of-sight magnetic field strength in star forming regions (e.g., Sarma et al. 2008).However, because the water molecule is weakly sensitive to the Zeeman effect (the splitting factor is 0.0021 Hz µG −1 , versus 2.8 Hz µG −1 for H i), this kind of study has only been carried out in a handful of regions hosting the strongest masers (see the review by Crutcher & Kemball 2019).In general, magnetic field measurements in outflows, are extremely difficult to carry out, thus their role in the launching, collimation, and propagation of ionized jets has been scarcely studied from an observational point of view.In this sense, strong water masers are excellent tools to not only investigate the kinematics of the shocked gas, but also to get key information of the magnetic fields in the post-shock regions.
In this manuscript, we present observations toward one of the sources in our sample, IRAS 19035+0641 A (I19035A hereafter), as a case study that demonstrates the potential of the full survey.I19035A is located at a distance of 2.2 kpc (Osterloh et al. 1997) in a chemically rich environment (e.g., Fuller et al. 2005;Araya et al. 2005Araya et al. , 2007;;López-Sepulcre et al. 2010, 2011;Lu et al. 2014;Taniguchi et al. 2019).Based on HiGal data, Rosero et al. (2019) report a bolometric luminosity of 6 × 10 3 L ⊙ for this region.The ultra-compact (UC) H II region IRAS 19035+0641 B (Rosero et al. 2016) is located at about 2 ′′ to the South-East of I19035A.
The goal of this work is to investigate the jet nature of I19035A, hence we will exclude IRAS 19035+0641 B from the following analysis and discussion, but direct the reader to see Section 2.2 from Rosero et al. (2019) for more information.This paper is organized as follows: in Section 2, the observations, data reduction, and imaging details are given.Our results are presented in Section 3, followed by a discussion in Section 4. Our findings are summarized in Section 5, and complementary figures are given in the Appendix.

OBSERVATION & DATA REDUCTION
The K-band (18 − 26.5 GHz) observations were carried out with NRAO's 1 VLA on March 27, 2022 in the A configuration utilizing the 3-bit samplers.A total of 57 spectral windows (SPWs) with 128 × 1 MHz channels were used for the continuum.A narrow 16 MHz SPW with 1024 × 15.6 kHz (i.e., 1024 × 0.21 km s −1 ) channels was set up to simultaneously observe the 22235.08MHz H 2 O(6 1,6 − 5 2,3 ) maser transition2 .The target phase center was set to RA(J2000) = 19 h 06 m 01 s .1,Dec(J2000) = 6 • 46 ′ 35 ′′ .0 and the total time on source was ∼ 11 minutes.The primary beam and largest recoverable scale of these observations are ∼ 2 ′ and 2 ′′ .4,respectively.Observations of the quasar 3C286 were used for flux density scale and bandpass calibration, while the complex gain calibration was done with observations of the calibrator source J1851+0035, with a cycle time of 4 minutes.
The line SPW was calibrated using the CASA pipeline version 5.3.0-147,modified for spectral line observations.Self-calibration for both phase and amplitude was performed on the strongest water maser spot (maser spot #4, see below) using NRAO's Astronomical Image Processing System (AIPS).The self-calibration solution ta-bles were then applied to the entire line SPW data.We used the AIPS task IMAGR with a grid weighting intermediate between natural and uniform (ROBUST = 0) to create the Stokes I water maser image cube for total intensity analysis.This resulted in a channel width, synthesized beam size, PA, and typical rms noise of 0.21 km s −1 , 0 ′′ .083× 0 ′′ .076,8.1 • , and 13.2 mJy beam −1 per channel (in line-free channels), respectively.Finally, for the Zeeman effect analysis, we additionally made Hanning smoothed Stokes I and V cubes3 .

Radio Continuum Emission
Our high angular resolution observations unveiled two 1.3 cm continuum peaks, which are indicated as A2 and A3 in the right panel of Figure 1, with a significance of > 5σ toward I19035A.The positions of these peaks are listed in Table 1.Peaks A2 and A3 are aligned in the NE-SW direction, separated by ∼0 ′′ .12(∼ 265 au), and were reported as a single radio continuum source at both 1.3 and 6 cm by Rosero et al. (2016) in their lower angular resolution observations.In addition, our detailed inspection of the images from the Rosero et al. (2016) observations revealed a second, weaker 6 cm continuum peak (> 10σ) located ∼ 0 ′′ .65 (∼ 1400 au) to the NE of A2, as seen in the left panel of Figure 1.We will refer to this second, weaker 6 cm peak as A1.While the lower resolution 1.3 cm continuum emission is dominated by a rather compact source, as reported by Rosero et al. (2016), it appears to be elongated toward A1 (see the lowest red contour in the left panel of Figure 1).However, this weak, extended 1.3 cm continuum emission toward A1 was not recovered in our higher angular resolution observations.
The spectral index (α, with S ν ∝ ν α ) between the 6 and 1.3 cm for the dominant continuum source, as seen in the lower resolution images, has a value of 0.9 (Rosero et al. 2016).Through the re-inspection of their data, as noted above, we find that the spectral index flattens to 0 as we reach A1.In our high angular resolution 1.3 cm continuum data, the measured in-band spectral indices are 0.9 and 0.6 for peaks A2 and A3, respectively.
We fitted 2-D Gaussians to each continuum peak using the highest resolution images available on them, that is, the 6 cm for A1 and our 1.3 cm for A2 and A3.The obtained positions, sizes, PA, peak intensities, and flux densities are presented in Table 1.Peaks A2 and A3 are marginally resolved, with upper limits to their deconvolved minor axes sizes at full width half maximum (FWHM) as reported in the third column of the table.The sum of the flux densities of A2 and A3 is consistent with the value reported by Rosero et al. (2016) for the 1.3 cm continuum source in their image.

22.235 GHz H 2 O Masers
We detected 8 water maser spots that seem to be associated with I19035A.The maser spots are distributed along the same direction as the continuum emission and form a narrow linear structure in the NE-SW direction, as shown in Figure 2. Maser #1 is the farthest from the radio continuum source, located at about 0.2 pc (∼19 ′′ ) to the NE, while the other maser spots are found within ∼2500 au from the continuum peaks.Masers #2 and #3 lay toward peak A1, while maser #4 is found close to the center of A2.Masers #5 − #8 are distributed to the SW of I19035A.
We report the peak position and intensity of each maser spot in columns 2, 3, and 4 of Table 2, respectively.We measured the velocity of the peak, its full width at zero power (FWZP), and the total velocity extent of the emission using the peak pixel spectrum for each maser spot.We list these in columns 5, 6, and 7 of Table 2, respectively, while the spectra can be found in the Appendix.Maser spots #4 and #6 span a wide velocity range.This likely indicates that these spots have substantial substructure that we are not able to spatially resolve.The brightest maser spot, maser #4, has a peak intensity of 102.8 Jy beam −1 at V peak LSR = 20.1 km s −1 .

Zeeman splitting in maser #4
The Stokes I and V profiles of the strongest maser #4 feature are shown with a histogram-like line in the top and bottom panels of Figure 3, respectively.The Stokes I spectrum shows an asymmetric shape, indicative of at least two Gaussian components.We used the AIPS task XGAUS to fit the Stokes I profile.The best fit, shown with a red line in the right panels of Figure 3, corresponds to two Gaussian components, shown separately in the left panels with a blue and a green line.The S-shape appearance observed in the Stokes V profile is usually considered a Zeeman effect detection.We derived the line-of-sight magnetic field B los by fitting the Stokes V profile (e.g., Momjian & Sarma 2017, 2019) to the expression (Troland & Heiles 1982;Sault et al. 1990): The first term on the right-hand side of the equation accounts for small leakage terms and, in this case, the  Values in parentheses are the uncertainties in the measurements.Note that the values for A1 were obtained using the C-band data from Rosero et al. (2016).
scale factor a was on the order of 10 −3 .The second term is the frequency derivative of the Stokes I profile multiplied by a fit parameter b, which contains the magnetic field information.Namely, b = zB los , where z is the Zeeman splitting factor.For the 22.2 GHz water transition, z = 2.1 Hz mG −1 (Nedoluha & Watson 1992).We used the AIPS task ZEMAN to fit equation 1, which derives a b value for each Gaussian component by fitting their derivative to the Stokes V profile.As seen in Figure 3, the Stokes V profile of maser #4 is well represented by the derivatives of the fitted Gaussians, and the measured B los values are 123 (±27) and 156 (±8) mG.These values are at a level of 4σ and 19σ (respectively), hence the measured magnetic field is significant in both cases.We note that, by convention, a positive value for B los means the line-of-sight magnetic field is pointing away from the observer.

Follow-up observations: extreme H2O maser variability
The detection of the ∼ 100 Jy H 2 O maser peak in I19035A led us to investigate the Zeeman splitting of this line and the magnetic field in the post-shock region.However, since the original observations were not designed for Zeeman effect analysis (e.g., the use of 3-bit samplers), we carried out follow-up VLA Bconfiguration observations on March 14, 2023, to confirm the Zeeman splitting.These observations, which utilized the 8-bit samplers, resulted in ∼ 30 minutes of on-source time.We used a single 8 MHz wide SPW with 2048 × 3.91 kHz channels, which translates to a The columns show the maser spot number (#), peak position, peak brightness (Iν ), velocity of the peak brightness (V peak LSR ), FWZP of the peak (∆V peak FWZP ), and the total velocity extent of the maser spot emission (∆V ).Values in parentheses are formal errors from the 2-D Gaussian fits.52 m s −1 channel width.The data were manually calibrated and imaged using CASA.The synthesized beam size and typical rms noise of the final image cube are ∼0 ′′ .3 and 5.5 mJy beam −1 per channel in line-free channels, respectively.The new observations recovered most maser spots seen in the 2022 observations, with varying degrees of variability.In Figure 4, we show a comparison of the two epochs for maser spot #4.While most peaks observed in the original spectrum are present, albeit with different intensities, the bright feature that showed Zeeman splitting was not detected in the new data.Thus, while we were not able to confirm the Zeeman splitting detection, the complete disappearance of the ∼100 Jy peak in maser #4 indicates extreme H 2 O maser variability in I19035A.

Nature of the Radio Continuum Emission
Based on the compact nature of our new 1.3 cm continuum data, one might consider the possibility of peaks A2 and A3 being two individual YSOs in a tight (d ≈ 265 au) protobinary system (e.g., Kraus et al. 2017;Zhang et al. 2019).In this scenario, the radio contin- uum emission would be attributed to early-stage H II regions.If A2 and A3 have a spherical geometry of diameter equal to their geometric mean (calculated using the sizes reported in Table 1), then their radii are approximately 40 and 30 au, or about 1.5 × 10 −4 pc.These sizes are consistent with those expected for HC H II regions (< 0.03 pc, Kurtz 2005).Using the ionizing photon rate (N Ly ) equation from Sánchez-Monge et al. ( 2013), assuming optically thin emission from a homogeneous H II region and an electron temperature of 10 4 K, we calculate N Ly values of 1.5 × 10 44 and 9.5 × 10 43 s −1 for A2 and A3, respectively, corresponding to B2-B3 Zero Age Main Sequence stars.Note that the cm spectrum is still rising in K-band, which for the above calculation was assumed to mark the turnover point between optically thick to optically thin conditions.A rising spectrum can also be explained by a stratified electron density distribution, which natu-rally would occur in an ionized stellar wind.We can use equation 5 in Felli et al. (1982) to estimate the stellar wind mass loss rate.Assuming a terminal wind velocity of 600 km s −1 , we obtain mass loss rates of the order of 10 −7 M ⊙ yr −1 .This number is several magnitudes higher than what is expected from near-main-sequence B2-B3 type stars (e.g., Cohen et al. 1997).Nonetheless, for stars in evolutionary stages prior to the main sequence, such high values are perhaps not unexpected.
While the above analysis shows that our data are consistent with a model where the radio continuum peaks A2 and A3 are caused by a high-mass binary, there is much morphological evidence that an ionized jet model is preferable.First, the observed extended and elongated morphology of the radio continuum emission at 6 cm can hardly be attributed to two compact objects located at the position of A2 and A3.Moreover, this elongated structure is very well aligned with A2 and A3, in the NE-SW direction.Second, the water maser spots #2 − #7 are also very well aligned in the direction of the radio continuum emission, along a structure much larger than the distance between A2 and A3, which is inconsistent with a wind-driven shell scenario (e.g., Beltrán et al. 2007).Ultimately, even though we cannot entirely discard a high mass binary for the origin of the observed 1.3 cm continuum emission, it appears more likely to be tracing an ionized jet.We will now discuss this scenario.Rosero et al. (2019) had classified I19035A as a possible ionized jet based on its association with a bipolar molecular outflow observed in CO(2−1) and HCO + (1 − 0) (Beuther et al. 2002a;López-Sepulcre et al. 2010).The orientation of the outflow at the angular resolution of these studies (11 ′′ and 29 ′′ , respectively) is SE-NW, i.e., perpendicular to the NE-SW orientation of the 6 cm continuum.However, the classification as a high-mass ionized jet was still made based on a number of additional observational results, like the coincidence of the 6 cm continuum with a cone-shaped 2.2 µm structure, presumably tracing scattered light from an outflow cavity, the presence of SiO(2−1) and (3−2) in the region (López-Sepulcre et al. 2011), as well as detections of Class II 6.7 GHz CH 3 OH masers in the vicinity of I19035A (Beuther et al. 2002b;Hu et al. 2016;Ouyang et al. 2019).Furthermore, with a radio luminosity of S 5GHz d 2 = 1.05 mJy kpc 2 , I19035A falls in the high-mass end of the S ν d 2 −L Bol relation shown in figure 8 from Rosero et al. (2019) for ionized jets, and it is well below the radio luminosity values for UC/HC H ii regions (≈ 10 2 mJy kpc 2 ).This further indicates that the I19035A radio luminosity at 5 GHz is likely due to shock ionization rather than photoionization from the central YSO.The fact that the large scale outflow is perpendicular to the jet direction can be explained by precession of the jet axis, for which there are a number of known cases discussed in the literature (e.g., Rodríguez et al. 2021).High angular resolution observations of the bipolar molecular flow are necessary to verify this conjecture and unveil the true outflow orientation.
Recently, Fedriani et al. (2023), using the Stratospheric Observatory for Infrared Astronomy (SOFIA) 7.7 − 37.1 µm observations together with multiwavelength data from the literature, compiled the detailed SED of IRAS 19035+0641 and carried out radiative transfer modeling based on the turbulent core accretion theory (Zhang & Tan 2018, and references therein).Their average model predicts a central protostellar mass of 15 M ⊙ accreting at a rate of 2.1×10 −4 M ⊙ yr −1 , and a bolometric luminosity of 4.9×10 4 L ⊙ .The mass outflow rate of the CO was estimated by Beuther et al. (2002a) as 0.9 × 10 −4 M ⊙ yr −1 , i.e., consistent with the idea that part of the energy released by the accretion process powers the outflow.
What do our new continuum data add to this picture?The radio luminosity at 5 GHz discussed above was derived from a source consisting of a compact core with a NE-SW extension (Rosero et al. 2016).The higher angular resolution 1.3 cm observations reported here resolve the source into two peaks (A2 and A3) aligned in the NE-SW direction.This shows that the elongated bipolar structure persists on a much smaller scale of 500 au, as expected from an ionized jet.Moreover, the rising in-band spectral indices for both A2 and A3 (0.9 and 0.6, respectively) are in good agreement with the ideal jet model from Reynolds (1986).In conclusion, our observations strengthen the interpretation that I19035A is an ionized jet from a high-mass protostar.

Association with the I19035A Jet
Water masers in high-mass star forming regions can be associated with both jets/winds and accretion disks (e.g., Imai et al. 2006;Sarma et al. 2008;Moscadelli et al. 2022).It is therefore pertinent to analyze whether the water masers detected in I19035A probe the ionized jet or a disk.In the latter case, water maser spots arise in different regions of the disk, and sometimes show Keplerian velocity gradients and arc-shape structures of a few au in size that are only resolved using Very Large Baseline Interferometry (VLBI) techniques (e.g., Torrelles et al. 1998;Franco-Hernández et al. 2009;Rodríguez et al. 2012).In contrast, our observations reveal a well defined linear maser spot distribution covering about 2500 au (ex-cluding maser #1, located at ∼ 40 000 au from the continuum source), and the water maser positions are well aligned with the ionized jet axis.If these were associated with a disk, a distribution perpendicular to the jet would be expected.Therefore, the water masers in I19035A appear to be mostly associated with the ionized jet observed in the radio continuum.However, we note that maser spot #4, located close to the radio continuum peak A2, shows several emission peaks of varied intensity at red-and blue-shifted velocities from the adopted V LSR (see Figure 4), thus it deserves additional consideration.We carried out 2-D Gaussian fits to each maser #4 spectral feature to obtain their positions, which we show color-coded by velocity on top of the 1.3 cm continuum emission in Figure 5.The individual components of maser spot #4 are located along a N-S line, and there is an apparent velocity gradient along this structure.If this gradient was associated with a Keplerian accretion disk, the equation M = RV 2 /G can be used to estimate the mass of the central object.Using this expression, the derived mass is ∼ 95 M ⊙ , which appears too large for an early stage high-mass YSO.Therefore, we believe it is unlikely that maser spot #4 traces an accretion disk.More likely, the observed velocity gradient is reflecting the velocities of individual shocks in the outflowing gas, and might be tracing an expanding or rotating flow.VLBI observations and a proper motion study are required to determine the nature of maser spot #4.
In the following, we discuss the association of the maser spots #1 and #8 with the I19035A jet.Maser #1 is found at ∼ 0.2 pc (∼ 40000 au) from the ionized jet, thus their relation is not immediate.Although Moscadelli et al. (2020) found that 84% of water masers are found within 1000 au from the driving protostar, they report detections at distances up to 18000 au (0.09 pc), and strong jet bow-shocks where water masers can be pumped are usually found at pc-scale distances from the high-mass YSOs (e.g., HH 80/81 at ∼1 pc from IRAS 18162−2048, Reipurth & Graham 1988).As seen in Figure 2, the PA of maser #1 is consistent with that of the ionized jet and the other water maser spots.Also the peak velocity of maser #1 is within 10 km s −1 from the systemic velocity of the region, and it is consistent with the velocity range in which the other masers in the jet present emission.Moreover, we detected variability in the intensity of maser #1 in our follow up observations (see Section 3.2.2),which suggests a common origin for this and the other maser spots in the jet.Finally, we did not detect a radio continuum source in the vicinity of maser #1 that could be associated with its emission, although we do not discard the possibility of this source being outside our field of view or fainter than ∼ 40 µJy beam −1 (i.e., our 3σ limit).Meanwhile, maser #8 is found in close proximity (projected) to the UC H ii region IRAS 19035+0641 B. The association of water masers with H ii regions has been investigated in the past.Observational studies found that water masers can be actively pumped in the edges of expanding winddriven, ionized shells (e.g., Beltrán et al. 2007).It has also been suggested that they might coexist with H ii regions and eventually disappear while the YSO remains visible (e.g., Sánchez-Monge et al. 2013).The PA and peak velocity of maser #8, as in the case of maser #1, are consistent with the other water masers in the jet.We did not detect maser #8 in our follow-up observations, which might be indicating extreme variability and a common pumping mechanism with the other masers that arise in the jet.Additionally, maser spot #8 was the only one detected nearby the UC H II region, and we do not see an arc-or ring-like distribution, indicative of an expanding shell, around IRAS 19035+0641 B. This could be due to an anisotropic shell structure, however, there is not enough information in the literature about the kinematics of this UC H II region, such as the expanding shell velocity.In summary, our observations suggest maser #1 is likely associated with the I19035A jet and are insufficient to unambiguously determine the origin of maser #8.A subsequent water maser proper motion study could shed light on the connection between these masers and the jet.Alternatively, observations of other shock probes and molecular outflow tracers could be used to determine whether there is a spatial correlation between masers #1 and #8 and the jet/outflow structure.

Measured Magnetic Field
We can use the more significant line-of-sight magnetic field component of B los ∼ 156 mG, measured from our possible Zeeman splitting detection to estimate some interesting jet and shock parameters.Firstly, we will estimate the number density (n post , n pre ) and magnetic field strength (B post , B pre ) in the pre-and post-shock gas.If we assume that the magnetic field is amplified proportionally with the density in the shocked region (e.g., Sarma et al. 2002), then we can relate these values as: In addition, Crutcher (1999) showed that, statistically, the total post-shock magnetic field (B post ) can be well approximated by B post = 2B los .We can also relate the B pre and the pre-shock gas density (ρ pre ) as: where β ≈ 5.8 × 10 5 in cgs units (calculated empirically by Sarma et al. 2002).Combining these expressions, and taking n post ∼ 10 9 cm −3 (Elitzur 1992), we find that the density in the pre-shock region is ρ pre = n pre m ∼ 10 −16 g cm −3 , where m is the mean molecular weight (i.e., m = 2.8m p , with m p being the proton mass).This yields a pre-shock gas number density of the order of n pre ∼ 2 × 10 7 cm −3 , which in turn translates to a pre-shock magnetic field strength of B pre ∼ 7 mG.We can also investigate what process is more likely to dominate the motion in the post-shock region by comparing the magnetic energy density (M) and the kinetic energy density (K).The magnetic energy density is defined as: Crutcher (1999) showed that, for an ensemble of measurements, a good approximation for the post-shock magnetic field value squared is given by B 2 post = 3B 2 los .For B los = 156 mG, M = 2.9 × 10 −3 erg cm −3 .Meanwhile, the kinetic energy density can be expressed as: where σ is the velocity dispersion of the FWHM linewidth ∆v, i.e., σ = ∆v(8 ln 2) −1/2 .If ∆v is the FWHM of the brightest peak in maser #4, which is 1.16 km s −1 , then K ∼ 1.7 × 10 −5 erg cm −3 in the postshock region.However, the maser emission is amplified along the path of greatest velocity coherence, which translates to a narrowing of the maser line.Therefore, a ∆v value obtained from a maser line is a lower limit of the real dispersion of the gas in the region.If instead we use ∆v = 7.0 km s −1 from single-dish CH 3 CN observations (Araya et al. 2005), the kinetic energy density would be of the order K ∼ 6.2 × 10 −4 erg cm −3 .This is still an order of magnitude smaller than M, indicating that the magnetic energy is likely dominating the motion in the post-shock region.
In addition, the shock velocity v s can be estimated by equating the magnetic energy density and the ram pressure (e.g., Sarma et al. 2008): For B los = 156 mG and the ρ pre estimated above, we obtain a shock velocity of 62 km s −1 .Typically, shocks are classified according to their velocity as C -or J -type shocks, with the former having v s 30 km s −1 and the latter v s > 30 km s −1 .The derived shock velocity values are in good agreement with maser spot # 4 arising in a J -type shock.
As noted previously, the above calculations make use of the B los measurement obtained from our Zeeman splitting detection in the brightest peak from the maser #4 spectrum.Even though we were unable to confirm this detection, as noted in Section 3.2.2,we briefly discuss the implications of the calculated physical parameters on the jet.The estimated magnetic field in the jet is B pre ∼ 7 mG, which is comparable to the higher values obtained by Curran & Chrysostomou (2007) for a set of high-mass YSOs at scales of about 0.1 pc (< 0.1 ∼ 6 mG).We expect the value of the magnetic field along the jet axis to remain fairly constant at scales of ∼ 0.5 pc (e.g., Carrasco-González et al. 2010), therefore it is likely that the energetic parameters derived above represent the magnetic field in the I19035A jet.However, when comparing the magnetic energy density and the kinetic energy density derived from the pre-shock magnetic field (B pre ∼ 7 mG) and density (n pre ∼ 2 × 10 7 cm −3 ) values, we find that the magnetic energy density (M pre ∼ 1.9 × 10 −6 erg cm −3 ) is an order of magnitude smaller than the kinetic energy density in the pre-shock region (K pre ∼ 1.2 × 10 −5 erg cm −3 ).Therefore, while the magnetic field appears to dominate the post-shock gas dynamics, the motions along the jet seem to be dominated by the gas kinematics.

Maser Variability
There has been increasing observational evidence that variations in the underlying protostellar object could be responsible for maser flares and maser variability of several different species.Brogan et al. (2018) discuss the case for water maser variability in an outflow caused by a burst of infrared radiation due to an accretion event propagating along the flow cavity.In this scenario, extreme increase or decrease in water maser intensity has a radiative origin, even though the principal pumping mechanism is collisions (see also Gray et al. 2016).A similar scenario was invoked by Hirota et al. (2021) to explain water maser variability in the highmass protostellar object S255 NIRS 3.Alternatively, Trinidad et al. (2003) demonstrate that random Gaussian fluctuations of an unsaturated maser (i.e., exponential amplification) in the line opacity can account for maser variability in the span of a few months.The results in Andreev et al. (2017) are also consistent with changes in the amplification path length.
In I19035A, we detected strong variability within one year in all water maser spots, although we do not see a global increase or decrease trend in the intensities.In all cases, some spectral features become brighter and others dimmer, new lines appear, and others become completely undetected.No position dependence of the variability along the NE-SW line of H 2 O masers was found either.While we find many very small scale intensity variations, the weakening of a ∼ 100 Jy maser to below 16.5 mJy (3× the rms noise in our follow-up observations) makes a statistical fluctuation unlikely.A possible scenario that would explain the observed variations is that accretion events occurred in I19035A.Monitoring of the H 2 O maser line, VLBI observations, and observations of other probes (such as the radiatively pumped Class II CH 3 OH masers) toward I19035A are needed to explore in detail the origin of the water maser variability in this high-mass protostellar object.

SUMMARY & CONCLUSIONS
We report VLA A-configuration observations of the 1.3 cm continuum and 22.2 GHz H 2 O masers toward the high-mass protostar I19035A.Our discoveries are summarized as follows: 1. We resolved the 1.3 cm continuum emission toward I19035A into a bipolar structure consisting of two peaks aligned in the NE-SW direction.We showed that the elongated, bipolar structure reported by Rosero et al. (2016Rosero et al. ( , 2019) ) persists at scales 500 au.The in-band spectral index calculated is consistent with what is expected from an ionized jet.
2. We detected eight 22.2 GHz H 2 O maser spots well aligned along the ionized jet axis.The intensity of the brightest peak in each maser spot ranges between 0.11 (maser #7) and 102.78 Jy beam −1 (maser #4).
3. We detected possible Zeeman splitting in the brightest line observed.The Stokes I profile was well fitted by two Gaussian components, and the measured line-of-sight magnetic field values are 123 and 156 mG.This yields a pre-shock magnetic field of ∼ 7 mG.
4. We carried out a second set of observations of the water maser emission in I19035A spaced one year from the original to confirm our Zeeman splitting detection.We found intense variability in all maser spots.Particularly, the bright ∼ 100 Jy beam −1 peak was not detected.
In conclusion, our observations of I19035A strengthen the original classification as an ionized jet from a high mass protostar based on the bipolar nature of the 1.3 cm continuum on scales 500 au, the occurrence of H 2 O masers along the jet axis, and the presence of extreme variability which could have been caused by an accretion event.

Figure 2 .
Figure2.Left: The 6 cm continuum emission in IRAS 19035+0641 is shown in color and cyan contours(Rosero et al. 2016).The contours are the same as in the left panel of Figure1.Right: Enlarged version of the area marked with the red rectangle in the left-hand panel.The 6 cm continuum emission is shown in color (note the scale difference between the panels), and our 1.3 cm continuum emission is shown in white contours.The contours are the same as in the right panel of Figure1, and the outlined and filled ellipses in the bottom left represent the 6 cm (∼0 ′′ .33)and our 1.3 cm (∼0 ′′ .08)synthesized beam sizes, respectively.The white rectangle shows the region covered in the right panel of Figure1.In both panels, the red + symbols mark the position of masers #1 through #8, which are also numbered.

Figure 3 .Figure 4 .
Figure 3. Stokes I (top) and V (bottom) profiles of the brightest peak of maser spot #4 are shown with a black histogram-like line.The green and blue curves in the left panel show the two Gaussian components fitted to the Stokes I profile (top) and their derivatives (bottom).The red curves in the right panel show the sum of the two components.

Figure 5 .
Figure 5. Positions obtained from 2-D Gaussian fits to each maser #4 spectral feature are shown with + symbols overlaid to the 1.3 cm continuum in black contours.The data points are color-coded by velocity the contours are the same as in the right panel of Figure 1.The dashed black ellipse in the bottom left represents the synthesized beam size.We note that the difference in the positions obtained are much larger than the formal fit errors provided by the CASA task imfit.

Figure
Figure A.1.Peak pixel spectra of all the detected maser spots.The vertical dashed line marks the adopted systemic velocity (Vsys = 33.8km s −1 ,Araya et al. 2005).The solid red and blue lines represent residual side lobes contamination from masers #4 and #6, respectively.

Table 1 .
Properties of the continuum peaks

Table 2 .
Properties of the Detected Maser Spots