Discovery of a Photoionized Bipolar Outflow toward the Massive Protostar G45.47+0.05

Massive protostars generate strong radiation feedback, which may help set the mass that they achieve by the end of the accretion process. Studying such feedback is therefore crucial for understanding the formation of massive stars. We report the discovery of a photoionized bipolar outflow toward the massive protostar G45.47+0.05 using high-resolution observations at 1.3 mm with the Atacama Large Millimeter/Submillimeter Array (ALMA) and at 7 mm with the Karl G. Jansky Very Large Array (VLA). By modeling the free–free continuum, the ionized outflow is found to be a photoevaporation flow with an electron temperature of 10,000 K and an electron number density of ∼1.5 × 107 cm−3 at the center, launched from a disk of radius of 110 au. H30α hydrogen recombination line emission shows strong maser amplification, with G45 being one of very few sources to show such millimeter recombination line masers. The mass of the driving source is estimated to be 30–50 M⊙ based on the derived ionizing photon rate, or 30–40 M⊙ based on the H30α kinematics. The kinematics of the photoevaporated material is dominated by rotation close to the disk plane, while accelerated to outflowing motion above the disk plane. The mass loss rate of the photoevaporation outflow is estimated to be ∼(2–3.5) × 10−5 M⊙ yr−1. We also found hints of a possible jet embedded inside the wide-angle ionized outflow with nonthermal emissions. The possible coexistence of a jet and a massive photoevaporation outflow suggests that, in spite of the strong photoionization feedback, accretion is still ongoing.


Introduction
Massive stars dominate the radiative, mechanical, and chemical feedback to the interstellar medium, thus regulating the evolution of galaxies. However, their formation process is not well understood. One key difference between high-and low-mass star formation is that massive protostars become so luminous that they generate strong radiation feedback, which potentially stops the accretion. For example, radiation pressure was considered a potential barrier for massive star formation (Wolfire & Cassinelli 1987;Yorke & Sonnhalter 2002). Photoionization is another important feedback process, as massive protostars emit large amounts of Lyman continuum photons that ionize the accretion flows to form photoevaporative outflows driven by the thermal pressure of 10 K 4 ionized gas (Hollenbach et al. 1994). However, theoretical calculations and simulations have suggested that, in the core accretion scenario these feedback processes are not strong enough to stop accretion (Krumholz et al. 2009;Kuiper et al. 2010;Peters et al. 2010;Tanaka & Nakamoto 2011). Tanaka et al. (2017) studied the combined effects from various feedback processes in massive star formation, including magnetohydrodynamical (MHD) disk winds, radiation pressure, photoionization, and stellar winds, finding that MHD winds are dominant, while the others processes play relatively minor roles. However, observational confirmation of such a theoretical scenario is still difficult, due to the rarity and typically large distances of massive protostars, especially the most massive type with strong radiation feedbacks.
The protostar G45.47+0.05 (hereafter G45; d=8.4 kpc, Wu et al. 2019) has a luminosity of ∼(2-5)×10 5 L e and a mass of ∼20-50 M e (De Buizer et al. 2017), based on infrared spectral energy distribution (SED) fitting (Zhang & Tan 2018). G45 is associated with an ultracompact (UC) H II region (Wood & Churchwell 1989;Urquhart et al. 2009;Rosero et al. 2019) and OH and H 2 O masers (Forster & Caswell 1989). Molecular outflows are seen in HCO + (1 − 0), with blueshifted emission to the north and redshifted emission to the south, consistent with the elongation of radio continuum emission (Wilner et al. 1996). Infrared observations at 10-40 μm show extended emission offset to the north of the UC H II region (De Buizer et al. 2005, which may come from the near-facing (blueshifted) outflow cavity. The offset between the infrared emission and the UC H II region may be due to the high infrared extinction toward the protostar, suggesting dense molecular gas surrounds the source.
The Very Large Array (VLA) 7 mm observation was performed in 2014 with A configuration and two 4 GHz basebands centered at 41.9 and 45.9 GHz (project ID: 14A-113). Details of the VLA data is summarized by Rosero et al. (2019).

1.3 and 7 mm Continuum
The 1.3 mm continuum emission appears concentrated within 5″ (4.2×10 4 au) from the central source (Figure 1(a)). At low resolution (C3+C6 configuration), the continuum emission is elongated in the north-south direction, with extended structures 53 . 9, panel (b)), ´ 0. 042 0. 028 ( = - P. A. 53 . 6, panel (c), and ´ 0. 043 0. 029 ( = - P. A. 53 . 4, panel (d)). The contour levels are s5 2 n ( = ¼ n 0, 1, ), with s = 1 0.13 K ( -0.33 mJy beam 1 , panel (a)), s = 1 1. . Only regions where both 1.3 and 7 mm emissions are s >5 are included. Here the 1.3 mm image is restored with the same beam size as the 7 mm image. The red rectangle in panel (a) marks the region shown in panels (b)-(f). In panels (b)-(f), the black dashed lines indicate the outflow axis ( =  P.A. 3 ) and projected half-opening angle (40°), and the red dashed line indicates a possible jet ( =  P.A. 15 ). The R.A. and decl. offsets in all the panels are relative to the continuum peak position from the ALMA C9 configuration data, a d = , toward the south. At high resolution (C9 configuration), the continuum shows an hourglass shape aligned in the north-south direction (Figure 1(b)), with a morphology that is highly symmetric with respect to an axis with »  P.A. 3 , which is consistent with the direction of the HCO + outflow (Section 1). The apparent half-opening angle of this hourglass is ∼40°( dashed lines in Figure 1(b)). More extended emissions are recovered by combining with more compact configuration data (panels (c) and (d)). The bipolarity can be still clearly seen in the combined images, and most of the fainter extended emissions are inside of the outflow cavity to the north and south. The 7 mm continuum (panel (e)) has very similar morphology as the highresolution 1.3 mm continuum. At both 1.3 and 7 mm, most of the compact continuum emissions are concentrated within 0 2 (1.7×10 3 au) from the central source. Another emission peak is seen ∼0 4 to the north, which was first identified by Rosero et al. (2019). Its emission appears to be more extended and shows an arc-like shape facing away from the main source in a direction consistent with the main outflow, suggesting that the northern emission structure is part of the main outflow.
Figure 1(f) shows the spectral index map derived using the 1.3 and 7 mm continuum intensities, a nn = n n n I I log log 1 2 1 2 ( ) ( ), where ν 1 =44 GHz and ν 2 =234 GHz. Note that the position of the 1.3 mm emission peak is slightly offset from either the emission peak or the center of symmetry of the 7 mm image (see Figure 1(e)). Therefore, we manually shift the 7 mm image by 7 mas in the R.A. direction and 13 mas in the decl. direction, so that the 1.3 mm emission peak coincides with the 7 mm center of symmetry in making the spectral index map (also see Figure 3(d)). The central region has spectral indices of α ν > 0.5, indicating partially optically thick ionized gas (α ν =2 for completely optically thick free-free emission), while more extended structures have spectral indices of α ν ≈0, consistent with optically thin free-free emission (e.g., Anglada & Rodríguez 2018). Note that, for dust continuum emission, the index is α ν , dust=2 in the optically thick case and 3α ν , dust4 in the optically thin case (assuming a typical dust emissivity spectral index of 1-2). These suggest that in this source, even at 1.3 mm, the compact continuum emission may have a significant free-free contribution, if not dominated by it. This is further supported by the fact that the 1.3 and 7 mm continuum emissions coincide well with the H30α hydrogen recombination line (HRL) emission (Figure 2(a); see below). Furthermore, the 1.3 mm peak brightness temperature is ∼2000 K, higher than that expected from dust continuum (the average dust temperature within 1000 au from a ∼30M e protostar is estimated to be several ×10 2 K from dust continuum radiative transfer (RT) simulations by Zhang & Tan 2018) and the dust sublimation temperature (∼1600 K), but can be naturally explained by thermal emission from ionized gas with a typical temperature of 10 4 K around massive protostars.
However, the extended 1.3 mm continuum emission shown in the low-resolution data should still be dominated by dust emission. The total 1.3 mm continuum flux >3σ within 5″ from the central source measured from the C3+C6 image is 0.77 Jy. The compact structure within 0 9 has a total >3σ flux of 0.34 Jy measured from the C9 image, or 0.61 Jy from the C3+C6+C9 image, which can be considered as an estimate (upper limit) of the free-free flux. Their difference of 0.16 to 0.43 Jy can be used as an estimate (lower limit) for the extended dust continuum flux of the source, which corresponds to a gas mass of (2.1-5.6)×10 2 M e assuming a temperature of 30 K (typically assumed for molecular cores), or (0.95-2.6)×10 2 M e assuming 60 K (suitable for within ∼10,000 au around a 20-50 M e protostar; Zhang & Tan 2018 (Ossenkopf & Henning 1994) and a gas-to-dust mass ratio of 141 (Draine 2011). These mass estimates suggest that there is a large mass reservoir for future growth of the massive protostar. Inside the wide-angle hourglass outflow cavity, a narrower jetlike feature is apparent in the 1.3 mm images with P.A.≈15°(red dashed line in Figure 1). Along this direction, a structure separate from the outflow cavity walls extends from the central source to the northern lobe. To the south, along the same direction, there is an additional emission peak ∼0 1 from the center, and also a stream of extended emission. These are most clearly seen in Figures 1(b) and (c). Along this direction, there are also some emission features in the 7 mm image. However, some of these may arise from strong sidelobe patterns from incomplete uv sampling in the VLA observation. As shown in Figure 1(f), some regions have negative spectral indices of −0.5<α ν <−0.2, inconsistent with pure free-free emission or dust emission, indicating possible contributions from nonthermal synchrotron emission, which is occasionally found associated with massive protostellar jets (e.g., Garay et al. 1996;Carrasco-González et al. 2010;Moscadelli et al. 2013;Beltrán et al. 2016;Rosero et al. 2016;Sanna et al. 2019). Thus, detection of negative spectral indices is also consistent with a jet being embedded inside the wide-angle ionized outflow. Evidence for the coexistence of wide-angle outflows with collimated jets has been found in other massive young stellar sources, such as Cepheus A HW2, where water masers trace a wide outflow and free-free continuum traces a collimated jet (Torrelles et al. 2011).

H30α Line
The position and morphology of the H30α emission are found to coincide very well with the 1.3 and 7 mm continuum emissions (Figures 2(a) and (b)).
The last term results from electrons from He + contributing to free-free emission but not to the HRL, and typically = . Assuming a characteristic ionized gas which is a reference value for the optically thin LTE conditions. As Figure 2(b) shows, most of the H30α emission is stronger than expected under LTE conditions. In some parts, the measured ILTC ratio is >200 km s −1 and reaching ∼300 km s −1 , which is significantly higher than the LTE value of 81 km s −1 , indicating strong maser amplification. Previously, millimeter HRL masers were detected in MWC349A, To our knowledge, G45 is only the third massive young stellar object to show strong millimeter HRL masers, with these reaching a similar level as in MWC349A. Figure 2(b) also shows that the strong maser effect is concentrated toward the east and west of the central source along a direction nearly perpendicular to the outflow axis, i.e., along the disk plane. This behavior is similar to that seen in MWC349A, in which HRL masers are found to be along certain annuli of the disk (e.g., Báez-Rubio et al. 2013;Zhang et al. 2017). In addition, there are two small regions along the outflow cavity walls to the northeast and northwest which also have strong maser levels.
Note that here we have assumed all the compact 1.3 mm continuum is free-free emission, which is reasonably valid, as discussed above. However, if there is a significant dust contribution in the 1.3 mm continuum, the observed ILTC ratio is then underestimated. Furthermore, if the H30α emission is not completely optically thin, the expected ILTC ratio under LTE conditions should be lower than the reference value of 81 km s −1 . Both of these effects would imply even higher levels of maser amplification. In addition, the ionized gas temperature can also affect the LTE value of ILTC ratio, which ranges from 34 km s −1 with T e =2×10 4 K to 156 km s −1 with T e =6× 10 3 K. Even at the low temperature end of this range, maser amplification is needed to explain the observed ILTC ratios.
In addition to the H30α line, the He30α helium recombination line (231.99543 GHz) is also detected (Figure 2(c)). Its spectrum can be fit with a single Gaussian profile at V lsr,H30α =−65 km s −1 (V lsr,He30α =+57 km s −1 ) and = FWHM 83 km s −1 . This line width is consistent with the total H30α line width. He30α integrated emission and ILTC ratio maps are shown in Figures 2(e) and (f). The He30α ILTC ratio increases from~5 km s −1 at the center to ∼20 km s −1 in the outer region. The increasing ILTC ratio toward outside is unexpected, as the outer region should have lower ratio between helium and hydrogen ionizing photon rates. One possible explanation is that, the He30α emission also has maser effects with a similar level to H30α (the expected He30α ILTC ratio in optically thin LTE conditions is 6 km s −1 using similar formula as Equation (1) and T e =10 4 K; note that He ++ is not considered), but at the center, the emission becomes partially optically thick, and the observed ILTC ratio decreases. Note that at V lsr,H30α ≈−50 km s −1 , the dips are due to the absorption of CH 3 OCHO 20 9,12 -20 8,13 line (231.9854 GHz) and CH 3 OCH -13 12 3 0,13 1,12 line (231.9879 GHz) from the cooler surrounding material in the foreground.

Model for Free-Free Emission
We construct a simple model to explain the 1.3 and 7 mm free-free emissions, assuming isothermal ionized gas filling bipolar outflow cavities. The details of the model fitting are described in Appendix A. The best-fit models have electron temperatures around T e =1×10 4 K, electron number densities at the center around n 0 =1.5×10 7 cm −3 and outer radii of the outflow launching region around v = 110 au 0 . The model constraints on the true outflow opening angle (which depends on the inclination) are weak. Although the model with smallest χ 2 has an outflow half-opening angle of θ w =22°, any θ w value within the range of 20°-40°explored by the model grid can generate good-fitting results with similar qualities (Appendix A). Infrared SED fitting of this source provided a range of inclinations of i>55°between the line of sight and outflow axis, corresponding to 35°<θ w <40°(given projected half-opening angle of 40°), which we consider to be a more probable range for the opening angle. The comparisons between the observations and the best-fit model (minimum χ 2 ) are shown in Figure 3. The best-fit model successfully reproduces the absolute values and relative distributions of both 1.3 and 7 mm intensities. There are some emission excesses in the observation at ∼0 04 to the south of the center seen in both bands (stronger at 7 mm), which indicates some substructures in the ionized outflow and/or contributions from nonthermal emissions; these are not taken into account in the model. Although detailed modeling of the H30α emission including the maser effects is outside of the scope of this Letter, we note that the electron temperature and density distributions derived from the free-free model are consistent with the observed strong H30α maser (see Appendix A).
The model requires an ionizing photon rate of (1.4-1.8)× 10 49 s −1 (Appendix A). For zero-age main sequence stars, this corresponds to a stellar mass of 43 to 48M e , i.e., spectral type of O4−O5 (Mottram et al. 2011). According to protostellar evolution calculations with various accretion histories from different initial and environmental conditions for massive star formation (Tanaka et al. 2016;Zhang & Tan 2018), this ionizing photon rate corresponds to a stellar mass of 28-48M e . In all of these protostellar evolution calculations, the central source has started hydrogen burning and reached the main-sequence in the current stage. These mass estimates are consistent with those estimated from infrared SED fitting (∼20-50 M e , De Buizer et al. 2017). Figure 4 shows the observed SED from radio to near-infrared. The infrared SED is well explained by dust continuum emission (De Buizer et al. 2017) using the continuum RT model grid by Zhang & Tan (2018). Based on the same physical model, Rosero et al. (2019) extended the model SED to radio wavelengths by combining photoionization and free-free RT calculations by Tanaka et al. (2016; thin lines in Figure 4). However, the observed radio fluxes are ∼100×higher than these model predictions. The new model presented here not only reproduces the fluxes at 1.3 and 7 mm, but also simultaneously reproduces the 1.3 cm fluxes (only low-resolution data is available at this wavelength). Note that the new free-free model does not fit the SED directly, but rather, fit the 1.3 and 7 mm intensity maps. At 1.3 mm, the freefree model can only reproduce the flux of the compact emission (also see Appendix Figure 6(e)), while the extended emission should be dominated by the dust emission.
In the model presented in Rosero et al. (2019), only photoionization of the MHD disk wind was considered, which has a typical density of 10 4 cm −3 in the outflow cavity, while in the new model, a dense outflow with densities of ∼10 7 cm −3 is included. Such high density suggests that this outflow is not an MHD disk wind, but instead a photoevaporative outflow launched from the disk. The derived outflowing rate of this ionized outflow is -´---V M 2 3.5 10 30 km s yr , which is consistent with theoretical expectations for photoevaporative outflow in later stage of massive star formation (e.g., Tanaka et al. 2017). Note that this outflowing rate is similar to or higher than the typical outflow rates of MHD disk winds from massive protostars, indicating that photoevaporation feedback is important at this stage. However, MHD winds are still the dominant feedback mechanism regulating core-to-star efficiency, as they typically are much faster and can sweep up a much larger amount of ambient gas (Tanaka et al. 2017). Furthermore, in spite of the strong photoionization feedback, the hourglass morphology of the free-free emission suggests that the ionized gas is still confined in the outflow cavity, while the disk/envelope remains neutral along the disk mid-plane.

Dynamics of the Ionized Gas
Figure 5(a) shows that the blueshifted emission of H30α is to the north of the continuum peak, and the redshifted emission is to the south, which indicates organized motion of ionized gas. To better demonstrate the kinematics of the ionized gas, we show the H30α emission centroid distribution in Figure 5(b) (see Appendix B for more details). The centroid distribution shows different patterns in different velocity ranges. At velocities lsr km s −1 and +  V 88 lsr km s −1 , the centroids are distributed roughly along the outflow axis direction, with more blue or redshifted centroids further away from the center, which is consistent with outflowing motion with acceleration. At velocities + + -  V 33 km s 88 1 lsr km s −1 , the centroids are distributed roughly perpendicular to the outflow axis, with blueshifted centroids to the east and redshifted centroids to the west, which is more consistent with rotation kinematics (e.g., Zhang et al. 2019).
These different patterns are better seen in panels (c) and (d), in which we plot the centroid offsets projected perpendicular to the outflow direction and along the outflow direction separately. . For 55°<i<90°determined by the infrared SED fitting, we obtain m dyn =28-40 M e . This mass estimate is consistent with those derived from the free-free modeling (28-48 M e ) and infrared SED fitting (∼20-50 M e ).
Along the outflow direction, the centroid velocity linearly increases with distance at velocities ) for the redshifted emissions, which are among the highest seen in ionized flows around forming massive stars (e.g., Moscadelli et al. 2018). Note that the two high-velocity Gaussian components in the H30α spectrum (Figure 2(c)) have central velocities of = + V 25 lsr and +80 km s −1 , which are close to the velocities where the centroid pattern changes from rotation to outflow. FWHM of these two components are 35 and 42 km s −1 , indicating motions with velocity FWHM of 28 and 36 km s −1 in addition to the thermal broadening with a FWHM of ∼20 km s −1 for ∼10 4 K ionized gas. Therefore it is reasonable to consider that the ionized outflow has a typical line-of-sight outflow velocity of 14-18 km s −1 in addition to the rotation. Considering projection effects, we adopt a fiducial value of = V 30 out km s −1 for the velocity of this outflow (Section 4.1). This outflowing velocity is of order of the sound speed (10 km s −1 ) and thus consistent with the scenario of a photoevaporative outflow (Hollenbach et al. 1994). A model including non-LTE effects is needed to fully understand the H30α kinematics including its transition from rotation to outflowing motion, which we defer to a future paper.

A Possible Embedded Jet
As discussed in Section 3.1, there is a jet-like structure inside the wide-angle ionized outflow seen in the 1.3 mm continuum morphology, along which several regions with negative spectral indices are seen, indicating possible contributions from nonthermal synchrotron emission. Synchrotron emission is sometimes found associated with high-velocity jets from massive protostars (e.g., Garay et al. 1996;Carrasco-González et al. 2010;Moscadelli et al. 2013;Beltrán et al. 2016;Rosero et al. 2016;Sanna et al. 2019). Padovani et al. (2015Padovani et al. ( , 2016 performed model calculations of particle acceleration in protostellar jets under various conditions. According to their model, strong particle acceleration can happen under the The red and black squares show the VLA fluxes integrated within 0 3 and 14″ from the center, respectively. Triangles and error bars show the fluxes at ALMA Band 6 with assumed 10% error. The black triangle shows the ALMA flux observed in C3+C6 configuration and integrated within 5″ from the center. The red and blue triangles show the ALMA fluxes within 0 9 from the center, observed in C3+C6+C9 and C9 configurations, respectively. The solid lines show the combined SEDs from free-free models and dust continuum models, and the dashed lines show the free-free model SEDs. The light dashed lines show the free-free SED from the model presented in Rosero et al. (2019), based on the photoionized MHD disk wind model by Tanaka et al. (2016). The dark dashed lines shows the free-free SEDs from the model presented in this work integrated over a region within 20,000 au (2 4) from the center. Best-fit models for different opening angles (θ w ) are shown but the differences are small. The dust continuum models are from SED fitting by De Buizer et al. (2017) using the dust continuum radiative transfer model grid by Zhang & Tan (2018). appropriate conditions for the G45 photoevaporative outflow (10 6 cm −3 with high ionization fraction), if the magnetic field strength is 1 mG and shock velocities are > 100 km s −1 . This high velocity needed for the particle acceleration suggests that the shock is caused by a fast jet rather than the photoevaporation outflow itself (Ṽ 30 out km s −1 ). Fast jets are commonly considered as a strong indicator of ongoing accretion. Therefore, this massive young star may be still accreting, in spite of strong feedback from its own photoionizing radiation. This provides a direct observational confirmation that photoionization feedback is not stopping accretion and limiting the final mass even for protostars with masses ∼50M e .
For the regions with negative spectral indices <-0.3 at ∼0.1″ to the north and south of the central source, we further estimate minimum-energy magnetic field strengths of B min =15 and 17 mG for the synchrotron sources (see Appendix C). Such values are much larger than those estimated previously in synchrotron emitting regions around massive protostars (e.g., ∼1 mG in G035.02+0.35 Sanna et al. 2019 and0.2 mG in HH80-81 Carrasco-González et al. 2010). However, previously the synchrotron emission was detected further away from the protostar (i.e., 10 4 au in G035.02+0.35 and 10 5 au in HH 80-81), while the synchrotron emission is detected <1000 au from the central source in G45. This magnetic field strength is consistent with some observations. For example, ∼20 mG was measured in Cep A HW2 within 1000 au (Vlemmings et al. 2010). It is also consistent with the magnetic field predicted for the base of the outflow by simulations of collapse of magnetically supported massive cores (e.g., Matsushita

Summary
With high-resolution ALMA and VLA observations at 1.3 and 7 mm, we have discovered a bipolar wide-angle ionized outflow from the massive protostar G45.47+0.05. The H30α recombination line shows strong maser amplification, with this source being one of the very few massive protostars so far known to show such characteristics. By modeling the 1.3 and 7 mm free-free continuum, the ionized outflow is found to be a photoevaporation flow launched from a disk of radius of 110 au with an electron temperature of ∼10 4 K and an electron number density of ∼1.5×10 7 cm −1 at the center. The mass of the protostar is estimated to be ∼30-50M e based on the required ionizing photon rate, or ∼30-40M e based on the H30α kinematics. The kinematics of the photoevaporated gas is dominated by rotation close to the disk plane and outflow motion away from this plane. The derived photoevaporative outflow rate is ∼(2-3.5)×10 −5 M e yr −1 . With robust detection of a resolved bipolar photoevaporative outflow, G45 provides a prototype of photoevaporation outflow in the later stage of massive star formation. Previously, clearly resolved bipolar ionized structure was detected in MWC 349A; however, its evolutionary stage or mass are still uncertain (e.g., Hofmann et al. 2002;Báez-Rubio et al. 2013;Zhang et al. 2017). G45 may also contain an inner, accretion-driven jet. The possible coexistence of a jet and a massive photoevaporative outflow suggests that, in spite of strong photoionization feedback, accretion is still proceeding to this massive young star. This confirms theoretical expectation that radiation feedback plays a relatively minor role in terminating accretion and determining the final mass of a forming massive star. The observed highly symmetric and ordered outflow and rotational motions also suggest that this massive star is forming via core accretion, i.e., in a similar way as low-mass stars. . Note that, in the observed data, the 1.3 mm emission peak is slightly offset from either the emission peak or the center of symmetry of the 7 mm image. Therefore, we manually shift the 7 mm image by 7 mas in R.A. and 13 mas in decl. so that the center of symmetry of the 7 mm image coincides with the 1.3 mm emission peak (see Figure 3(d)). We also adopt a position angle of 3°for the outflow axis.
The Appendix Figure 6 shows the distributions of the χ 2 values with various parameters. As shown in panels (a)-(c), the parameters T e , n 0 , and ϖ 0 are well constrained by the model fitting. The best-fit model (minimum c tot 2 ) has T e ≈10,000 K, »´n 1.5 10 cm the source inclination to be i>55°, corresponding to 35°<θ w <40°, which we consider the more probable range for the opening angle. As the constraint on opening angle (and inclination) is weak, the best-fit model (minimum c tot 2 , marked by a star in Figure 6) should be considered only as an example model among a group of models with good fits to the observations.
In addition to constraining the four free parameters, we also calculate the total fluxes, ionizing photon rates, and mass outflow rates from the models, which are shown in panels (e)-(h). Panels (e) and (f) show the comparison between the model and observation fluxes at 1.3 and 7 mm integrated in regions with observed intensities >0.1I obs,center (the same as the region used in the 2D model fitting). The model fluxes are slightly lower (within 5%) than the observed values. Panel g shows the ionizing photon rate, which is calculated following where α B is the case B recombination rate for which we adopt the fitting formulas from Draine (2011), and the integration is along the axis from 0 to ¥. Note that this is only a lower limit for the ionizing photon luminosity in the outflow, which are just enough to ionize the outflow and ignore the possible effects of dust absorption. The uncertainty of the EUV photon rate is relatively large. Models with low χ 2 values have EUV photon rates spanning from ∼1.4×10 49 s −1 to ∼6×10 49 s −1 . The uncertainty is mostly from the weak constraint on the opening angle (and inclination), as shown in panel d. The best-fit model (minimum c tot 2 ) has an EUV photon rate of 4.5×10 49 s −1 , while the most probable inclination range (55°<i<90°, i.e., 35°<θ w <40°) gives a range 1.4×10 49 s −1 <S < 1.8× 10 49 s −1 . Panel (h) shows the mass outflow rates in the model (assuming a velocity of 30 km s −1 ; see Section 4.2). The best-fit models have rates from 2×10 −5 (v/30 km s −1 ) M e yr −1 to 3.5×10 −5 (v/30 km s −1 ) M e yr −1 . These values are obtained assuming all the ions are H + and the total mass is 1.4 times of the hydrogen mass. If we consider the He + contributions to the freefree emission, the mass outflowing rates are ∼7% lower than the above values, assuming a constant ratio of 0.08 between He + and H + throughout the ionized region.
The comparisons between the observations and best-fit model are shown in Figures 3 and 4, and the Appendix Figure 7, including the images in both bands and spectral index maps, the mm to radio SEDs, as well as intensity profiles along and perpendicular to the outflow axis. The best-fit model successfully reproduces the absolute values and relative distributions of emission intensities in both bands. The profiles perpendicular to the outflow axis are typically well fit, while some asymmetries or excesses along the outflow axis still remain. In particular, at 0 02 to the south, the observed 7 mm image shows an additional peak. At the same location, the observed 1.3 mm image also has a slight excess (see the residual curves in panels (a) and (c) of Appendix Figure 7). This may be caused by some substructures in the ionized outflow that are not taken into account in our modelʼs simple geometry and density distribution. It is also possible that this excess, stronger in the lower frequency, contains contributions of nonthermal synchrotron emission.
Here we briefly discuss whether the observed H30α maser can happen under the temperature and density conditions derived from the free-free model. Strong maser amplification requires the total continuum and HRL optical depth t t + -1 ff HRL  (Báez-Rubio et al. 2013). The non-LTE HRL absorption coefficient can be written as k b k = b n mn HRL HRL,LTE , with k HRL,LTE being the absorption coefficient under LTE conditions, and b n and β mn being the departure coefficients for the recombination line from electronic level m-n (e.g., m=31, n=30 for H30α; Dupree & Goldberg 1970). We adopt the values of b n and β mn derived by Walmsley (1990), which are dependent on the electron temperature T e and density n e . The non-LTE HRL optical depth is then t k b k = = l b l n mn HRL HRL HRL,LTE , with l being the physical length scale. Toward the center, with = = n n e 0 -1.5 10 cm 7 3 and T e =10,000 K from the free-free model, and assuming the length scale l=ϖ 0 =110 au for a first-order estimation, we estimate τ HRL ≈−20. The total continuum and HRL optical depth is then t t + »-19 ff HRL , which is consistent with the observed strong maser amplification (Appendix Figure 8). We further estimate the distribution of non-LTE HRL optical depth over the ionized region, by adopting a simple relation n e ∝r −2 ∼l −2 . The dependence of τ ff +τ HRL on the distance to the center r or the density distribution n e is shown in the Appendix Figure 8. Strong maser amplification appears constrained in the region r400 au for H30α line, which is consistent with the observation that the high maser amplification is concentrated close to the disk mid-plane. It also shows that similar maser effects may be seen in HRLs from about H20α to about H42α, which may be tested by future observations. However, accurate model prediction for HRL maser intensities is difficult as full consideration of the radiation field is need, which is outside of the scope of this Letter.

Appendix B Determining the H30α Emission Positional Centroids
The centroid positions are determined by fitting Gaussian ellipses to the H30α emissions observed in the C9 configuration at channels with peak intensities > 10σ. The accuracy of the centroid positions are affected by the signal-to-noise ratio (S/N) of the data, following Δθ fit =θ beam /(2 S/N) (Condon 1997), where θ beam is the resolution beam size, for which we adopt the major axis of the resolution beam θ beam =38 mas (320 au). The phase noise in the bandpass calibrator also introduces an additional error to the centroid positions through passband calibrations ). The phase noise in

Appendix C Estimating the Magnetic Field Strength from Synchrotron Emission
In the region with negative spectral indices indicating synchrotron emission, we estimate the minimum-energy magnetic filed strength, B min , which minimizes the total energy of the synchrotron source by assuming equipartition between the magnetic field energy and the particle energy. Following the classical formula for estimating minimum-energy magnetic field strength (e.g. where L R is the synchrotron luminosity, R is the source radius, k is the energy ratio between the ions (carrying most of the energy but not radiating significantly) and electrons, and c 12 is a coefficient dependent on the spectral index and the minimum and maximum frequencies in the integration of the spectrum.
Here B min is in Gauss, and L R , R, k, and c 12 are in cgs units. We define the synchrotron emission region as the region with intensity spectral indices of α<−0.3 and both 1.3 and 7 mm emissions >5σ.Here we only consider the negative spectral index regions within ∼0 2 from the center, as the 7 mm image is affected by strong sidelobe patterns further out from the central source. The areas of the regions to the north and south of the central source are 0.00021 arcsec 2 and 0.0013 arcsec 2 , which correspond to effective radii of R=0 008 and 0 020, which we use as the radii of the synchrotron emitting regions. Assuming all the 1.3 and 7 mm emissions from these regions are due to synchrotron emission, we obtain synchrotron fluxes of 0.041 and 0.070 mJy at 1.3 and 7 mm in the northern region, and 1.0 and 1.7 mJy at 1.3 and 7 mm in the southern region. Note that these are only upper limits for the synchrotron fluxes, as these regions only have negative spectral indices ∼−0.3 while typical synchrotron emission has a spectral index of ∼−0.75, suggesting there are both free-free and synchrotron contributions in these regions. However, the magnetic field strength is not very sensitive to the synchrotron flux (see Equation (7)). We then calculate synchrotron luminosities of L R =8.1×10 29 erg s −1 and L R =2.0×10 31 erg s −1 in the northern and southern regions, by integrating the synchrotron fluxes from 1.3 to 7 mm assuming a power-law distribution. For this frequency range, we obtain c 12 =3.2×10 6 following Govoni & Feretti (2004), We also adopt k = 40 following Carrasco-González et al. (2010) and Sanna et al. (2019). These values give B min =15 and 17 mG in the northern and southern regions, respectively.