Sub-Parsec-Scale Jet-Driven Water Maser with Possible Gravitational Acceleration in the Radio Galaxy NGC 1052

We report sub-pc-scale observations of the 321-GHz H$_2$O emission line in the radio galaxy NGC 1052. The H$_2$O line emitter size is constrained in $<0.6$ milliarcsec distributed on the continuum core component. The brightness temperature exceeding $10^6$ K and the intensity variation indicate certain evidence for maser emission. The maser spectrum consists of redshifted and blueshifted velocity components spanning $\sim 400$ km s$^{-1}$, separated by a local minimum around the systemic velocity of the galaxy. Spatial distribution of maser components show velocity gradient along the jet direction, implying that the population-inverted gas is driven by the jets interacting with the molecular torus. We identified significant change of the maser spectra between two sessions separated by 14 days. The maser profile showed a radial velocity drift of $127 \pm 13$ km s$^{-1}$ yr$^{-1}$ implying inward gravitational acceleration at 5000 Schwarzschild radii. The results demonstrate feasibility of future VLBI observations to resolve the jet-torus interacting region.


Introduction
The Atacama Large Millimeter/submillimeter Array (ALMA) detected H2O JKa,Kc = 102,9 − 93,6 emission at rest frequency of 321.225677GHz in the radio galaxy, NGC 1052 (Kameno et al. 2023b).This is the first submillimeter maser detection in a radio galaxy and the most luminous 321-GHz H2O maser known to date with the isotropic luminosity of 1090 L⊙.The most plausible interpretation is maser amplification of background synchrotron emis-sion through population-inverted molecular gas in a torus.However, the angular resolution of 0 ′′ .68 × 0 ′′ .56 with a compact array configuration up to 500 m was insufficient to rule out the possibility of thermal emission.A higher angular resolution was required to confirm for maser emission by brightness temperature (TB), which is a criterion to classify the emission mechanism by comparing with the upper energy level of Eu/k b = 1861 K, where k b is the Boltzmann constant.Kameno et al. (2023b) mentioned that baseline length of > 3500 m would provide concrete evidence for non-thermal emission.New ALMA long-baseline observations were deemed necessary to confirm that the emission was of masers, and also allows us to clarify the spatial distribution of population-inverted molecular gas.
NGC 1052 is a unique target emanating well-collimated (Nakahara et al. 2020;Baczko et al. 2022) double-sided sub-relativistic jets with a bulk speed of 0.26c − 0.53c (Vermeulen et al. 2003;Baczko et al. 2019), where c is the speed of light.Kinematic studies clarified that eastern and western jets reside approaching and receding sides, respectively, with the viewing angle > 57 • (Vermeulen et al. 2003) under assumption of intrinsic symmetry.Baczko et al. (2022) questioned the symmetry and confirmed the same sign of the viewing angle.A multi-phase torus in pc-scale vicinity of the core is seen nearly edge-on and has been investigated by various probes such as molecular absorption lines (Omar et al. 2002;Liszt & Lucas 2004;Impellizzeri et al. 2008;Sawada-Satoh et al. 2016;Sawada-Satoh et al. 2019;Kameno et al. 2020;Kameno et al. 2023a), H I absorption line (Shostak et al. 1983;Vermeulen et al. 2003), H2O 22-GHz maser emission (Braatz et al. 1994;Claussen et al. 1998;Braatz et al. 2003;Kameno et al. 2005;Sawada-Satoh et al. 2008), and free-free absorption (FFA) in thermal plasma (Kameno et al. 2001;Kameno et al. 2003;Vermeulen et al. 2003;Kadler et al. 2004).H2O megamaser at 22 GHz is a firm probe for subpc-scale molecular gas in AGNs and its (sub-)Keplerian rotation curve allows us to weigh a central supermassive black hole (SMBH), to study the molecular gas distribution at the highest resolutions around and AGN, and to measure geometrical distance independently from other estimation (Nakai et al. 1993;Miyoshi et al. 1995;Moran et al. 1995;Koekemoer et al. 1995;Herrnstein et al. 2005;Humphreys et al. 2008;Humphreys et al. 2013;Impellizzeri 2022;Gallimore & Impellizzeri 2023).The reach of 22-GHz H2O maser as a cosmological distance ladder is constrained by an angular resolution to clarify the rotation curve.
Since the first detection of the extragalactic 183-GHz and (tentative) 439-GHz H2O masers in NGC 3079 (Humphreys et al. 2005), submillimeter masers have become recognized as new and important probes for AGN.Hagiwara et al. (2013) searched for the 321-GHz H2O maser in five type-2 Seyfert galaxies and detected in the Circinus galaxy.That is followed by the second detection in NGC 4945 (Pesce et al. 2016;Hagiwara et al. 2016;Hagiwara et al. 2021) and the first detection in a radio galaxy NGC 1052 (Kameno et al. 2023b).With a higher angular resolution by 14× with the same baseline length, the 321-GHz H2O maser offers potential to extend the reach of the cosmological distance ladder.
The 321-GHz maser requires a condition with higher temperature and density compared with that of the 22-GHz maser.Neufeld & Melnick (1991) modeled the excitation condition for the 321-GHz H2O maser under a physical temperature of T k = 1000 K and found that the maser emissivity was maximized with a density of hydrogen nuclei ∼ 10 9 cm −3 .Yates et al. (1997) also concluded that the population inversion condition for the 321-GHz H2O maser was found at T k > 1000 K for n(H2) = 8 × 10 8 and 4 × 10 9 cm −3 with n(H2O) = 10 4 cm −3 .Gray et al. (2016) modeled the excitation conditions for possible H2O maser transitions for evolved star environments and clarified that the population-inversion for the 321-GHz maser requires T k > 1000 K and n(H2O) > 10 4 cm −3 .Thus, we expect closer distance to the SMBH, greater gravitational acceleration, wider velocity range, and faster proper motion of disk rotation.
In this paper, we report the results from ALMA observations in the long-baseline array configuration (C-10) for 321-GHz H2O emission toward NGC 1052.This paper is organized as follows: Section 2 describes observation conditions of ALMA and data reduction procedures.Section 3 presents results of continuum images, spectral profiles, size and position of velocity components of the 321-GHz H2O emission.In section 4 we develop discussion about justification of maser emission, velocity drift, velocity gradient, jet-torus interaction, and torus structure.Then we summarize our conclusions in section 5. We employ the systemic velocity of Vsys,LSR = 1492 km s −1 , the luminosity distance of DL = 17.6 Mpc, the angular distance of DA = 17.5 Mpc, and the linear scale of 85 pc arcsec −1 (Kameno et al. 2020).

Observations
ALMA observations under C-10 (the most extended) array configuration with the maximum baseline length of 15.2 km have been carried out twice on 2023-07-10 and 2023-07-24 as summarized in table 1. Integration times on the target (NGC 1052) and on a bandpass calibra- tor were 1108 s and 912 s, respectively, for each execution.The first execution, uid://A002/X109d26e/X12c34, was significantly affected by high wind speed that made three antennas shutdown in the beginning and then four antennas later.Absence of antennas in the southern arm degraded the angular resolution in north-south direction poorer than the requested resolution of 0 ′′ .015.The retake, uid://A002/X10a7a20/X65a8, was successful and satisfied the required resolution.Median values of the system noise temperatures were 139 K and 132 K under the condition with the precipitable water vapor (PWV) of 0.58 mm and 0.48 mm, respectively.The distribution of acquired spatial frequencies, also known as (u, v) coverage, is shown in figure 1.
We configured four spectral windows (SPWs) for four 2-GHz basebands (BBs), assigning two in the upper sideband (USB) and two in the lower sideband (LSB).BB4 in USB was tuned to cover the H2O JKa,Kc = 102,9 − 93,6 emission and the SPW was set in the frequency division mode (FDM) with a channel separation of 976.562 kHz corresponding to 0.914 km s −1 .Other three SPWs were set in the time division mode (TDM) to maximize the sensitivity in continuum emission with a effective bandwidth of 5.4 GHz.

Data Reduction
Reduction scripts and reduced data are available in the GitHub repository. 1  Phase and amplitude calibrations were applied following the standard way, using CASA (CASA Team et al. 2022) 6.4.1.We used bandpass calibrators, listed in table 1, as a flux calibrator assuming 2.104 Jy for J0423 − 0120 and 1.304 Jy for J0238 + 1636.This amplitude calibration yields ∼ 3% uncertainty.We applied smoothed bandpass calibration (Yamaki et al. 2012) with a 9-channel smoothing width to improve the signal-to-noise ratio (SNR) in the 1 https://github.com/kamenoseiji/ALMA-2022.A.00023.S spectra.The target, NGC 1052, was bright and compact enough to allow phase calibration by itself.
We produced continuum-subtracted cross power spectra in BB4.Then, we used the task tclean with the natural weighting to produce image cubes.Spectral profiles of the H2O emission line, shown in figure 2, were sampled at the single phase-center pixel.The continuum-subtracted visibilities were also used to estimate sizes of the emission region and to measure positions of emission, to be presented in section 3.
After subtracting a point CLEAN component at the center position, significant emission remained as an elongated structure spanning ∼ 0 ′′ .001 in PA∼ 70 • .The structure was decomposed by CLEAN components aligned in the same PA.While self calibration with a single central component yields the agreement factors 2 of σ = 3.94 and 4.93 for 2023-07-10 and 2023-07-24, respectively, those with the aligned CLEAN components result in σ = 3.57 and 1.81.Significant improvements in agreement factors justify the reality of the elongated structure.The elongated structure is buried by the central component if the CLEANed images are restored with the main lobe of the synthesized beams.To highlight the structure, we produced super-resolution images with a 1.5-milliarcsec (mas) circular restoring beam.

Results
Table 1 summarizes imaging performance of the observations.Image rms is estimated by statistics in emissionfree area of syhthesized images.While rms in spectral-line channel map was dominated by thermal noise, that in continuum image involves systematic errors such as sidelobe leakages.
The profile mainly consisted of two velocity components 2 The agreement factor in Difmap self calibration is equivalent to χ 2 N , where χ 2 N does not take account of the number of degree of freedom.-the blueshifted narrow component and the redshifted wide component -separated by a local minimum near the systemic velocity.The velocity range of 1280 -1580 km s −1 on 2022-05-12 was shifted redward to 1300 -1700 km s −1 The local minima between red and blue components also redshifted in later observations.The velocity shifts recall the velocity drift of the systemic velocity component of the 22-GHz H2O maser in NGC 4258 (Haschick et al. 1994;Greenhill et al. 1995;Nakai et al. 1995).We will discuss about the velocity drift in section 4.2.
While the blue component showed a simple singlepeaked profile, the red component contained multiple narrow sub-components.The peak intensity of the red component was lower than that of the blue one on 2022-05-12, but it got brighter exceeding the blue one in later epochs.
To characterize the line profiles, we used the non-linear least squares function, nls, in the R statistical language (R Core Team 2023) for spectral decomposition with 7 velocity components with either Gaussian or Lorentzian profile as summarized in table 2 and shown in figure 3. We chose Gaussian or Lorentzian profile for each component to minimize the residuals.Standard errors in table 2 were estimated by nls.See the source code plotMaserSpec.R in the repository for detail.
The decompositions are not unique, especially in the red component consisting of complex sub-components.Thus, the component labels in table 2 do not claim the unique identity between two epochs.Residuals of the fittings results are characterized by χ 2 over the degree of freedoms (d.o.f.) as shown in figure 3.
Identification of the blue component B1 was certain as it was was isolated from others.While the velocity width did not change significantly, the peak flux density increased by 34.5 ± 4.2% and the center velocity shifted by 4.85 ± 0.5 km s −1 in 14 days.
The red component was too complicated to identify every sub component in two epochs.Four spiky sub components with the FWHM < 12 km s −1 (R2, R3, R4, and R5) on 2023-07-10 became unclear behind bumps composed by R3 and R5 on 2023-07-24.We evaluated integrated intensity of the red components by subtracting the Lorentzian blue component model from the observed line profile.Integrated flux densities of the red components increased from 7.69 ± 0.23 Jy km s −1 to 10.12 ± 0.30 Jy km s −1 by 31.6±4.2%.Weighted mean velocity of the red components shifted from 1552.7 ± 1.1 km s −1 to 1559.0 ± 0.8 km s −1 by 6.3 ± 1.4 km s −1 .

Size of emission regions
We measured sizes of the continuum and line-emitting components by fitting a Gaussian function to visibilities as a function of spatial frequency (also known as u, v distance).Emission-line visibilities of blue and red components were generated by averaging the spectrum in the velocity ranges of 1360 -1500 km s −1 and 1500 -1625 km s −1 , respectively.Figure 4 shows the visibility amplitude as a function of spatial frequencies, together with the confidence bands of the best fits and 99% confidence levels.
The results are summarized in table 3 together with the minimum brightness temperatures.Broader spatial frequency coverage on 2023-07-24 allowed us to constrain the size of emission region in < 0.6 milliarcsec (mas) and the brightness temperature of TB > 10 6 K.

Positioning line emitting components
We performed uvmodelfit task in CASA to estimate the positions of line emission source using continuumsubtracted visibilities binned by 5 spectral channels (4.57km s −1 ).We set an unresolved point source model with a fixed flux density.Threshold for the maser flux density was set to 3× the line image rms listed in table 1.Because self-calibration was applied before continuum subtraction, uvmodelfit determined only relative positions of the line emitters with respect to the continuum core.
Figure 5 shows distribution of the emission line components on 2023-07-24 registered on the continuum map together with the binned spectrum as a color indicator of velocities.The same plot on 2023-07-10 is presented in figure 9 in Appendix.Close-up pictures are shown in figures 6    and 9. Median and interquartile ranges for position errors along major and minor axes were 0.43 (0.33 -0.84) mas and 0.32 (0.24 -0.63) mas on 2023-07-10, 0.23 (0.15 -0.41) mas and 0.16 (0.11 -0.29) mas on 2023-07-24, respectively.Positions on 2023-07-10 were more scattered due to larger position errors compared with those on 2023-07-24.
We identified significant positional differences between blue and red components.While blue components were centered at the continuum peak position, the center of red components offsets along the receding side of the jet by ∼ 1 mas.We applied linear regression of VLSR ∼ jJ +dD, where J and D are offset along PA=68 • .35 and −21 • .65,respectively, to estimate the velocity gradients, j and d along and perpendicular to the jet. Figure 7 shows position-velocity diagram on 2023-07-24 along (J-axis) and perpendicular (D-axis) to the jet.That on 2023-07-10 is presented in figure 10 in Appendix.Best-fit results are summarized in table 4. In both epochs, velocity gradient along the jet is significant while no significant gradient perpendicular to the jet was identified.The sign of velocity gradient along the jet remains the same, though the values are significantly different between two epochs.The smaller gradient in the first epoch can be caused by the larger positional scatter caused by poorer (u, v) coverage.
Comparing with the spatial distribution of the 22-GHz H2O maser consisting of two clusters separated by a 1-mas gap (Sawada-Satoh et al. 2008), the span of the 321-GHz maser fits the size of the 22-GHz gap though the absolute position is uncertain to register these maps.

Justification for maser emission
The ALMA long-baseline observations revealed the brightness temperature exceeding the excitation energy level and clearly rule out the possibility of thermal emission.Strong intensity variation indicates that the emission is a maser.The timescale of ∆T = 14 days constrains the size of the emitter < c∆T = 0.012 pc.The most plausible explanation for the compact, high-brightness non-thermal emission line on the continuum core, is maser amplification of the background continuum emission via population-inverted molecular gas.
If B1 is identified as the 1333 km s −1 component (Kameno et al. 2023b) on 2022-05-12, the velocity drift would be 35.3 ± 0.9 km s −1 yr −1 .Since the drift rate in 424 days is significantly smaller than that in 14 days, it is unlikely to identify the 1333 km s −1 component as B1.
The velocity drift in red components is more complicated.The change of weighted mean velocity corresponds to 164 ± 37 km s −1 yr −1 or 5.2 × 10 −3 m s −2 .The change is ascribed to either velocity drift or growth of R3 and R5 components.Comparing with the 1467 km s −1 component (Kameno et al. 2023b) on 2022-05-12, we have 73.8 ± 4.4 km s −1 yr −1 in 424 days, which is again inconsistent with the 14-day change.Spectral monitoring with a higher cadence than 14 days is desired to trace the velocity drift of every component.
Nevertheless, all of velocity changes show the same direction, increase of radial velocity, are similar to the velocity drift of the systemic-velocity component of the 22-GHz H2O maser disk in NGC 4258 (Greenhill et al. 1995;Nakai et al. 1995), indicating gravitational acceleration toward the nucleus as the masers locate in front of the core.
We then discuss about the velocity drift of B1.If the velocity drift is ascribed to gravitational acceleration of the SMBH that weighs MBH = 1.5 × 10 8 M⊙ (Woo & Urry 2002), we have the distance of the maser from the SMBH, r, as r = GMBH/g = 0.072 pc or 5000 Rs where Rs = 2GMBH/c 2 is the Schwarzschild radius.While the radial velocity of B1 is blushifted with respect to the systemic velocity, the relative velocity of 118 km s −1 is expected to vanish within 1 yr due to the deceleration of 127 ± 13 km s −1 .Indeed, the maser velocity does not exceed the escape velocity, vesc = 2GMBH/r = 4.2 × 10 3 km s −1 .
Compared with the inner boundary of 0.12 pc or 3.3 × 10 4 Schwarzschild radii for the 22-GHz H2O maser disk in NGC 4258 (Modjaz et al. 2005), the submillimeter maser distance in NGC 1052, scaled by the Schwarzschild radius, is 6.6 times as close as that in NGC 4258.The ratio of kinetic temperatures for 321-GHz to 22-GHz maser excitation conditions follows T k ∝ r −1/2 , as 2000 K : 800 K ∼ √ 6.6 : 1 As we expected in section 1, the submillimeter H2O maser is likely to arise at a closer distance to the central engine.
Alternatively, we admit a chance of misidentifying appearance and disappearance of independent maser components as the velocity drift.To ensure identification of maser components, high-cadence monitoring for the maser spectrum is desired.4) and perpendicular to the jet (right, D-axis in Table 4).The diagram on 2023-07-10 is presented in figure 10 in Appendix.Color and size of markers represents the velocity and the flux density as well as Fig. 5.

Velocity gradient
As presented in section 3.4, we identified velocity gradient along the jet.While the blue components reside on the continuum core, the red components offset by ∼ 0.1 pc.The sub-pc-scale 22-GHz H2O maser (Claussen et al. 1998;Sawada-Satoh et al. 2008) shows opposite direction of the velocity gradient in the western cluster toward the receding jet, while the eastern cluster did not show obvious velocity gradient.Sawada-Satoh et al. (2008) found that the 1-mas (0.1 pc) gap between two clusters of the 22-GHz H2O maser coincides the plasma obscuring torus that caused FFA.As the FFA opacity is proportional to ν −2.1 , where ν is the frequency, the opacity of FFA at 1 GHz, τ0 ≃ 1000, yields τ ≃ 1.5 at 22 GHz.Thus, the 22-GHz H2O emission inside the torus is obscured and the maser must locate outside the plasma torus where dynamical interaction with the jet does not significantly impact the radial velocity.They claimed that the velocity gradient in the western cluster indicate acceleration infalling gas.On the contrary at 321 GHz, the FFA opacity will be 0.005 and the submillimeter maser is visible at a 0.2-pc vicinity of the core inside the plasma where jet-torus interaction plays a role in the velocity field.This explains the opposite directions of velocity gradient between 22-GHz and 321-GHz masers.Gravitational acceleration by the central SMBH is the major factor in both inside and outside the torus to cause the velocity drift.
Alternative model ascribes the velocity gradient to projected velocity of a rotating disk, which successfully modeled for 22-GHz H2O megamaser sources (Miyoshi et al. 1995).However, we conclude that the velocity gradient does not represent such a rotating disk model by following reasons.If the velocity drift of g = 4.0 × 10 −3 m s −2 is ascribed to gravitational acceleration and the derived dis-tance is the radius of a rotating disk, we expect the rotating speed of Vrot = c Rs/2r = 3000 km s −1 and the velocity gradient of Vrot/r = 4.2 × 10 4 km s −1 pc −1 or 3.5 × 10 3 km s −1 mas −1 .The observed velocity gradient is not consistent with the expectation for the rotating disk model.Furthermore, the orientation of the gradient along the jet is unlikely for a simple rotating disk model.
If there exists a rotating disk perpendicular to the jet, the width of the maser emission section to support the velocity width of ∆V ∼ 400 km s −1 would be r∆V /Vrot = 0.01 pc or 0.1 mas.The position accuracy of the maser components in our observations were not sufficient to measure the expected velocity gradient in that width, as shown in figure 7 (right).VLBI observations up to 43 GHz (Nakahara et al. 2020;Baczko et al. 2022) measured the width of the jet and revealed that the jet is collimated in a cylindrical shape with a width of 1.3 × 10 3 Rs inside a distance of 10 4 Rs from the central engine.The estimated maser emission section width of ∼ 700 Rs is comparable to the estimated jet width.The constraint on the maser emission section supports the scenario that the population-inverted gas at the inner side of the torus amplifies the background continuum emission of the jet.To trace the disk rotation, position accuracy of ∼ 0.01 mas is desired.
A rotating disk with Vrot = 3000 km s −1 would yield high-velocity components at large impact parameters.Expected frequency offsets of ±3.2 GHz are outside of the spectral coverage in our observations.An expected impact parameter of r = 0.072 pc or 0.85 mas is too small to be resolved.Future observations with the spectral setup targeting the terminal velocity component would clarify presence or absence of a rotating disk.A proper motion of a rotating disk with 3000 km s −1 would be identified by 8000-km VLBI observations within 240/SNR days if position accuracy would be obtained by a synthesized beam size divided by a SNR.

Jet-torus interaction
Following the former section, the velocity gradient along the jet axis is vindication of jet-torus interaction.Here we argue about physical parameters following discussion in Kameno et al. (2023a).Adopting overpressure factor of 1.5 (Fromm et al. 2019), we have ρjv 2 j = 1.5ρtv 2 s , where ρj and vj are density and speed of the jet, and ρt and vs are the density of the torus and the advancing speed of the shock, respectively.Taking vj = 0.26c (Vermeulen et al. 2003), ρt = 1.1 × 10 −12 kg m −3 estimated from the maser excitation condition of nH2 > ∼ 3.3 × 10 8 cm −3 (Gray et al. 2016), we have ρj = 3.5 × 10 −18 kg m −3 .This is one order of magnitude greater than the estimation by SO absorption feature (Kameno et al. 2023a).Since the required temperature for the submillimeter H2O maser is significantly higher than 344 ± 43 K estimated by submillimeter SO absorption lines, the population-inversion zone is likely to locate upstream of the jet compared with the SO evaporation region.Higher pressure in the population-inversion zone than that in the SO evaporation region may indicate collimation and acceleration of the jet via jet-torus interaction.

Torus structure
This section addresses clumpy torus structure based on the time variation of the maser profile.
The blue and red components of H2O maser showed amplitude increase by 34.5 ± 4.2% and 31.6 ± 4.2% in 14 days while the continuum flux density increased by 12.8 ± 4.2%.As maser intensity is a product of the brightness of the background emission source and the amplification (or absorption) gain through population-inversion (or thermal) zone, the breakdown of the increase of maser intensity involves ∼ 20% increase of the maser amplification gain variation with the timescale of 1.2 × 10 6 s/ log 1.2 = 6.6 × 10 6 s.
The variation in amplification gain indicates inhomogeneity inside the molecular torus.VLBI images for spatially resolved HCN (Sawada-Satoh et al. 2016) and HCO + (Sawada-Satoh et al. 2019) absorption lines estimated the upper limit of the clump size of ≤ 0.1 pc in the molecular torus of NGC 1052.ALMA millimeter/submillimeter absorption-line studies (Kameno et al. 2020) showed that the molecular gas in the torus is clumpy with an estimated covering factor of 0.17 +0.06  −0.03 .Passage of populationinverted clumps across the line of sight toward the continuum background can cause variation of amplification gain.The Jeans length in the population-inversion zone is estimated to be λJ = cs/ √ Gρ ∼ 5.3 × 10 14 m or 0.017 pc, where ρ = 1.1 × 10 −12 kg m −3 is the density for maser excitation condition of nH2 > 3.3 × 10 8 cm −3 (Gray et al. 2016) and cs is the sound speed for the kinematic temperature of 2000 K. Crossing time of λJ /Vrot = 1.8 × 10 8 s is 27 times as long as the variation timescale.Thus, passage of a population-inverted clump in front of the continuum background is unlikely to explain the maser gain variation.
Another possibility is the change of excitation condition supported by the increase of the continuum flux density.Increase of depth in the population-inversion zone, generated by interaction with more powerful jet, can cause increase of the maser flux density.
Attenuation through thermal molecular gas outside the population-inversion zone would be another origin of the gain variation.The spectral minimum near the systemic velocity can be caused by self absorption through thermal gas (Watson & Wallin 1994).Kameno et al. (2023a) estimated the temperatures of SO molecules showing millimeter and submillimeter absorption lines to be 26 ± 4 K and 344 ± 43 K, respectively.Smaller Jeans length for lowertemperature yields shorter crossing time by one order of magnitude than that of a population-inverted clump and thus comparable to the maser variation timescale.
Figure 8 illustrates the clumpy torus model.Inside of the torus is excited by interaction with the jet to form population-inversion zone.Outer region of the torus consists of colder thermal clumps that attenuate maser emission to bear the spectral minimum near the systemic velocity and absorption lines in thermal transitions.

Summary
ALMA follow-up observations with the widest array configuration led findings summarized below.
1.The 321-GHz H2O emission in NGC 1052 is definitely a maser with the brightness temperature TB > 10 6 K.The intensity variation in 14 days also supports that the emission is a maser.This indicates feasibility for VLBI observations to clarify spatially-resolved structures and proper motions of the maser together with the core and the jets.A 8000-km-baseline array would offer a resolution of 0.002 pc which allows resolving the structure of maser emitting gas as discussed in the previous section.2. The maser profile consists of mainly two velocity components, blue and red, straddling a local minimum near the systemic velocity.Both velocity components increased the flux density compared with the previous observation 424 days ago.They still grew by 34.5 ± 4.2% and 31.6±4.2% in 14 days while the continuum flux density increased by 12.8 ± 4.2%.The variability in amplification gain implies inhomogeneity of the populationinversion zone and/or change of excitation condition.3.Both velocity components showed redward velocity drifts.While the red component was too complex to distinguish the drift from growth of R3 and R5 subcomponents, the velocity drift of 127 ± 13 km s −1 yr −1 in the blue component is more confident.If the drift is real and is ascribed to gravitational acceleration onto the SMBH with MBH = 1.5×10 8 M⊙, the distance of the maser will be 5000 Rs from the SMBH.Since the maser velocity relative to the systemic velocity is smaller than the escape velocity, the maser emitting gas must be bound in the SMBH gravitational potential.For a rotating disk model, corresponding rotation velocity will be Vrot = 3000 km s −1 .4. We identified velocity gradient of −174 ± 19 km s −1 mas −1 along the jet.The gradient indicates that the population-inverted gas is driven by the jet.No significant velocity gradient perpendicular to the jet was identified under positional errors of ∼ 0.16 mas.To resolve expected velocity gradient for a rotating disk, a spatial resolution sharper than 0.1 mas is required.5.The pressure of population-inversion zone is one order of magnitude greater than that of the SO evaporation region.This implies that the maser locates upstream of the jet compared with the SO evaporation region and jet acceleration and collimation process by jet-torus interaction.
For future works, spectral monitoring is desired to ensure the velocity drift possibly caused by gravitational acceleration.Designed spectral setup targeting the terminal velocity emission at ∼ ±3000 km s −1 would allow us to determine the rotation velocity, if the maser associates in a rotating disk.VLBI observations are feasible to resolve the structure of the population-inverted zone and detect a proper motion for disk rotation or outflows.Submillimeter masers would yield a potential to clarify the inner edge of molecular torus and to determine geometric distances.

Fig. 3 .Fig. 4 .
Fig. 3. Decomposition of the line profiles on 2023-07-10 (left) and 2023-07-24 (right).Solid blue lines show the best-fit model to the observed spectra shown in black filled circles with error bars.Dashed lines with labels represent decomposed components whose parameters are listed in Table 2.The residuals are shown offset by −0.01 Jy.

Fig. 5 .
Fig. 5. Continuum and maser maps on 2023-07-24.The contours show continuum image at 310 GHz combining BBs 1-3.The contour levels are ±3σcont multiplied by powers of 2 where σcont = 0.24 mJy beam −1 .The synthesized beam and the restoring beam are shown as the gray ellipse at the map center and the circle at the left bottom corner.Locations of maser components are drawn in filled colored circles with error bars.The marker size is proportional to line flux density and the color represents the velocity as shown in the spectral profiles in the upper right corner.The map on 2023-07-10 is presented in figure 9 in Appendix.

Fig. 6 .
Fig. 6.Close-up of maser components on 2023-07-24.Coordinates are the same with figure 5.The cross at the origin indicates the continuum peak position.The color represents the velocity as shown in the spectral profiles in the bottom left corner.The size of markers is proportional to flux density of at the channel.The plot on 2023-07-10 is presented in figure 9 in Appendix.Color indicator is shown in figure 5.

Fig. 8 .
Fig. 8.A schematic of the clumpy torus and jets in NGC 1052.The black marker at the center indicates the core which emanates double-sided jets presented in white.The torus surrounding the core consists of molecular clumps (filled circles) and plasma (unseen filler).Color of clumps indicates excitation condition: population-inversion zone in red to orange and thermal gas in yellow, green, and blue.While yellow clumps reside in warm zone, green and blue clumps in the equatorial plane and outer area of the torus are colder.The population-inversion zone amplifies the background continuum emission (core and the jets) to produce the 321-GHz H2O maser.The thermal gas casts molecular absorption lines on the background emission.(a):observer's view with a viewing angle of 70 • .(b):pole-on view of the torus.The dashed lines stand for the jet width and the maser amplification area in the cross section with the population-inversion zone.(c):cross cut of the torus and the jet.Lines of sight through the torus stand for blue and red maser components while systemic velocity component is attenuated through the thermal gas on the plane of the equator.
* G and L stand for Gaussian and Lorentzian functions, respectively.

Table 3 .
Estimated Gaussian models and brightness temperatures of the continuum and line-emitting componentsA size stands for FWHM of a Gaussian component in milliarcsec (mas).Best and max sizes are derived by the most likelihood and 99% confidence band shown in figure4.

Table 4 .
Linear regression for velocity gradient (Kameno et al. 2023a)zzoli et al. 2022)is consistent with the 100-pc-scale bipolar outflow traced by [O III](Sugai et al. 2005;Cazzoli et al. 2022)and by [Ar II] and [Ne III](Goold et al. 2023).As argued inKameno et al. (2023b), the velocity range of the submillimeter H2O maser in NGC 1052 coincides with the submillimeter SO absorption lines at the estimated temperatures of 344 ± 43 K implying excitation by jet-torus interaction(Kameno et al. 2023a).The observed velocity gradient supports this scenario.