What Determines the Boundaries of H 2 O Maser Emission in an X-ray Illuminated Gas Disk ?

High precision mapping of H 2 O megamaser emission from active galaxies has revealed more than a dozen Keplerian H 2 O maser disks, which enable a ∼ 4% uncertainty estimate of the Hubble constant as well as providing accurate masses for the central black holes. These disks often have well-defined inner and outer boundaries of maser emission on sub-parsec scales. In order to better understand the physical conditions that determine the inner and outer radii of a maser disk, we examine the distributions of gas density and X-ray heating rate in a warped molecular disk described by a power-law surface density profile. For a suitable choice of the disk mass, we find that the outer radius 𝑅 out of the maser disk predicted from our model can match the observed value, with 𝑅 out mainly determined by the maximum heating rate or the minimum density for efficient maser action, depending on the combination of the Eddington ratio, black hole mass, and disk mass. Our analysis also indicates that the inner radius for maser action is comparable to the dust sublimation radius, suggesting that dust may play a role in determining the inner radius of a maser disk. Finally, our model predicts that H 2 O gigamaser disks could exist at the centers of high-z quasars, with disk sizes of ≳ 10 − 30 pc.


INTRODUCTION
In the nuclear regions of external galaxies hosting active galactic nuclei (AGNs), there exist powerful cosmic masers from the − + = 6 16 − 5 23 transition 1 of the ortho-H 2 O molecule emitting at 22.23508 GHz, which arise either from a sub-parsec circumnuclear disk (e.g.Kuo et al. 2011;Gao et al. 2016) or a nuclear wind or jet (e.g.Claussen et al. 1998;Greenhill et al. 2003;Sawada-Satoh et al. 2008;Kuo et al. 2020a;Surcis et al. 2020).These 22 GHz H 2 O megamasers often display total maser luminosities 10 6 greater than that of typical Galactic maser sources, and their extremely high surface brightnesses permit mapping at sub-milliarcsecond resolution using Very Long Baseline interferometry (VLBI), providing a unique ★ E-mail: cykuo.tara@g-mail.nsysu.edu.tw(NSYSU) 1 Here, is the total angular momentum of the H 2 O molecule, with − and + representing the projections of on two molecular axes (e.g.Cooke & Elitzur 1985).
probe of the gas distribution and kinematics on sub-parsec scales at the centers of distant galaxies (Lo 2005).
In the archetypal H 2 O maser galaxy NGC 4258 (e.g.Haschick et al. 1994;Argon et al. 2007;Humphreys et al. 2008), the masing gas resides between radii of ∼0.13−0.26pc in a thin disk that is viewed in almost edge-on and follows Keplerian rotation.Such disk maser systems enable black hole (BH) mass measurements to percent-level accuracy (e.g.Kuo et al. 2011;Gao et al. 2017) and provide an accurate geometric distance measurement independent of distance ladders and standard candles (e.g.Reid et al. 2009;Gao et al. 2016;Pesce et al. 2020a).
To identify disk megamasers like NGC 4258 for measuring the Hubble constant 0 , the Megamaser Cosmology Project (MCP; Reid et al. 2009;Braatz et al. 2010) has carried out an extensive survey of H 2 O megamaser emission from >4800 AGNs (Kuo et al. 2018b(Kuo et al. , 2020b)), resulting in the detection of 30 candidate disk masers (Pesce et al. 2015).The follow-up imaging of these candidates has increased the number of H 2 O maser disks with high precision VLBI maps by a factor of 4 over the past decade.This progress not only enables an 0 measurement with 4% uncertainty (Pesce et al. 2020b), but also allows us to explore the black holehost galaxy coevolution (e.g.Kormendy & Ho 2013;Greene et al. 2016).In addition, it provides a new tool for constraining the spins of supermassive BHs (Masini et al. 2022).
While progress in the detection and mapping of H 2 O maser disks has been significant in the past decade, the physics of H 2 O maser disks has received less attention and is not well-explored.In particular, it is not yet clear what causes the warps seen in many Keplerian maser disks (e.g.Bregman & Alexander 2009, 2012;Caproni et al. 2007;Kuo et al. 2011;Kamali et al. 2019).In addition, it is still uncertain what determines the inner and outer boundaries of a maser disk and why H 2 O megamaser emissions are typically confined within the narrow range of ∼0.1−1 pc in a circumnuclear disk (e.g.Neufeld & Maloney 1995;Gao et al. 2017;Wardle & Yusef-Zadeh 2012).Addressing questions like these is valuable because it would not only improve our understanding of the maser disks themselves, but would also allow one to explore unrecognized systematics when using the H 2 O megamaser technique for 0 and BH mass measurement.Moreover, a deeper knowledge of maser disk physics can also allow us to investigate whether hyper-luminous H 2 O "gigamaser" disks could exist in the high-redshift universe (Lo 2005).Investigating this possibility is important because such gigamasers, if they exist, would enable the application of the maser technique to high-z galaxies for measuring dynamical BH masses and accurate distances, opening a new avenue to test cosmological models (e.g.King et al. 2014;Lusso et al. 2019) that might address the "Hubble tension problem" (e.g.Di Valentino et al. 2021).
In this paper, we aim to explore whether the physical conditions favorable for population inversion of H 2 O molecules could play the primary role for determining the inner and outer radii of an H 2 O maser disk.We focus primarily on the physics of the masing region, because this is the part of the disk that is relevant for 0 and BH measurements based on the H 2 O megamaser technique.Understanding the mechanisms that determines the physical size of the masing section would allow us to explore whether the megamaser technique could be applied to galaxies in a redshift range beyond the current limit of the MCP (i.e.0.05).In Section 2, we discuss measurements of maser disk radii based on VLBI maps compiled from the literature.In Section 3, we make predictions of the inner and outer radii of a maser disk based on an examination of the physical conditions in a circumnuclear disk, and we compare these predictions with the observed values.A discussion of the possibility of applying the maser technique to high-z H 2 O gigamasers is presented in Section 4, and our present study is summarized in Section 5.

The Radius Measurement
H 2 O megamasers in a disk configuration (maser disks hereafter) often display three distinct spatial and velocity groups, including the systemic, redshifted, and blueshifted maser features (Kuo et al. 2011;Pesce et al. 2015).In nearly all cases, the VLBI maps of such maser disks show well-defined inner and outer boundaries of maser emissions on sub-parsec scales, indicated by the positions of the highest and lowest velocity components of either blueshifted or redshifted masers.When the position of the dynamical center of the disk is determined either from disk modeling (e.g.Reid et al. 2013;Kuo et al. 2013Kuo et al. , 2015) ) or rotation curve fitting for the BH mass (e.g.Kuo et al. 2011;Gao et al. 2017), the physical inner ( in ) and outer ( out ) radii of a maser disk can be measured directly from its VLBI map given the distance.Alternatively, if one assumes that the maser emission originates only from a disk, without any components associated with a jet or outflow, one could also infer the physical disk radii indirectly from a single-dish spectrum of a maser disk based on the velocities of the high-velocity maser features and assuming Keplerian rotation2 about a known central mass (see details in Gao et al. 2017).
For both direct and indirect methods, measuring the physical radii of the inner and outer boundaries of a maser disk from its projected angular radii requires the knowledge of the galaxy distance, which can be obtained either from modeling of the maser disk itself or from published values (e.g.NED; https://ned.ipac.caltech.edu/).Since a galaxy distance from disk modeling is available only for a limited number (∼6) of disk maser systems (Pesce et al. 2020b), we simply adopt the Hubble flow distances from NED3 to infer the physical sizes for all disk maser systems discussed in this paper for consistency.
To identify the physical mechanism that determines the physical size of an H 2 O maser disk, we first compile in Table 1 all H 2 O megamasers from literature that display geometrically thin maser disks in their VLBI maps (e.g.Kuo et al. 2011;Reid et al. 2013;Gao et al. 2017;Pesce et al. 2020a).The majority of the radius measurements are drawn from Table 5 in Gao et al. (2017), which provides reliable estimates of disk radii for all "clean"4 H 2 O maser disks that follow Keplerian rotation.We exclude NGC 4388 listed in Gao et al. (2017) from our analysis because the lack of systemic maser features in this system plus the small number of high-velocity maser components (∼4) makes it difficult to be confident that the maser emission originates from a thin, rotating disk.For completeness, we also include the thin maser disk in NGC 1068 in our analysis.This source was first reported to be a complicated maser disk that follows sub-Keplerian rotation (Greenhill et al. 1996), with a disk mass comparable to the BH mass (Lodato & Bertin 2003).However, the latest VLBI observations and analysis (Gallimore & Impellizzeri 2023) suggest an opposite interpretation for NGC 1068, showing that the disk kinematics are consistent with Keplerian rotation and low turbulence speeds, with a disk mass ∼100 times smaller than BH .Therefore, although NGC 1068 could be more complicated than other Keplerian maser systems, considering that the mechanism determining the maser disk size may not be dependent on gas kinematics alone, we still include it.In total, our sample thus consists of 16 maser disks.The BH masses and the inner/outer disk radii of these systems are shown in Column (2) through (6) in Table 1.
In the left panel of Figure 1, we plot disk radii in units of parsec against the BH masses of our maser sample.The top and bottom edges of each black vertical line marked by the filled black circles represent out and in for a maser disk, respectively, revealing that the maser emission is mostly confined within the radial range of ∼ 0.1 − 1 pc for the majority of the sources.The red and blue error bars in the plot show the uncertainties of in and out , respectively.In addition to the uncertainties from the galaxy distance and maser position measurements, the error bars also include the discrepancy between the disk radii directly measured from VLBI maps and indirectly inferred from high-sensitivity single-dish maser spectra (Pesce et al. 2015), which could indicate a systematic uncertainty arising from limitations in the spectral coverage or sensitivity in the VLBI observation.
In the MCP, the single-dish monitoring of the maser galaxies usually covers a ∼3 times wider spectral range than the VLBI observation.In the case of NGC 1194, for example, the high-sensitivity single-dish spectrum (Pesce et al. 2015) shows that the (weak) redshifted maser feature having the highest velocity lies well outside the spectral coverage of the VLBI observations (Kuo et al. 2011), suggesting that the true inner radius of the disk may be smaller than the value derived from the VLBI map.In addition, the spectra of some sources (e.g.NGC 3393) sometimes reveal faint, isolated maser features lying between the systemic and the high-velocity maser complexes in the spectra.Without sensitive VLBI mapping, it is difficult to discern whether such features are systemic masers, outflow components, or high-velocity maser features residing close to the outer boundary of the maser disk, leading to uncertainties in out .For a more detailed discussion on analysis of the uncertainties due to spectral coverage, we refer the readers to Gao et al. (2017).

The Sensitivity Effect on Radius Measurement
It is well known that the flux densities of the maser emissions usually drop quickly at the wings of the maser line complexes in a Keplerian maser disk (e.g.Kuo et al. 2011;Pesce et al. 2015;Gao et al. 2017).Because of limited sensitivity, one may not necessarily be able to detect the weakest maser features that define the intrinsic size of a maser disk.As a consequence, a VLBI observation will inevitably reveal sharp edges of the maser disk in the maser map, and the observed values of the inner and outer radii of the maser disk can vary depending on the sensitivity of the observation.
To estimate the magnitude of the variation in the radius measurement caused by the sensitivity limit, one can consider the most well-known maser disk in NGC 4258, which has a distance of 7.6 Mpc (Humphreys et al. 2013), the closest system in our sample.Its proximity allows for the detection of the weakest maser features among all maser systems studied in this work.We note that if NGC 4258 were placed at a distance of 65 Mpc, the median distance of our maser sample5 , the flux densities of the masers would drop by a factor of ∼70.Given the same VLBI sensitivity (Argon et al. 2007), the observed inner and outer radii of NGC 4258 would change by ∼45% and ∼12%, respectively.We expect that these variations give the magnitude of the uncertainty caused by the sensitivity limit in the radius measurement for our sample if one compares the disk sizes at the same median distance.Since this uncertainty is comparable to the radius errors shown in Table 1, which include uncertainties inferred from sensitive single-dish spectra, the conclusions for this work are not significantly affected by the sensitivity-dependent effect on disk size measurement.

The Characteristic Radii of H 2 O maser Disks
As one can see in the left panel of Figure 1, the mean radius of a maser disk appears to increase with the BH mass.This trend was first reported by Wardle & Yusef-Zadeh (2012) based on the data of 8 maser disks, showing that the outer radii can be approximately described by out = 0.3( BH /10 7 ⊙ ) pc.With the inclusion of 6 more maser disks in their analysis, Gao et al. (2017) obtained a significantly different scaling ( out ∝ 0.57±0.16BH ) and show that both in and out are well correlated with the BH mass6 , suggesting that BH plays an important role in determining the inner and outer radii of a maser disk.Because of the significant correlation between out and BH , we conjecture that the maser disks may reveal an interesting, characteristic scale if the disk size is expressed in units of the Schwarzschild radius S = 2 BH / 2 , where G and c are the gravitational constant and the speed of light, respectively.in , respectively.The horizontal axis denotes the BH mass for each maser system.The black dashed line and dotted line in the plot indicates the radii of 0.12 pc and 0.42 pc, the medians of in and out for all sources, respectively; Right panel: the maser disk size in units of 10 5 Schwarzschild radii ( S ).The majority of the disk megamasers in our sample have inner and outer radii within a factor of ∼2 from in ∼ 1 × 10 5 S and out ∼ 3 × 10 5 S , respectively.The horizontal dashed line and dotted line indicate the medians of in (0.87 × 10 5 S ) and out (2.88 × 10 5 S ) for all sources.
In the right panel of Figure 1, we plot the inner and outer radii of the maser disks in units of 10 5 S .As one can see in this plot, the majority of our maser disks have inner and outer radii within a factor of ∼2 from in ∼ 1× 10 5 S and out ∼ 3× 10 5 S , which are shown by the horizontal dashed and dotted lines, respectively.It is likely that the presence of these characteristic disk radii may be associated with the fine-tuning nature of the H 2 O megamaser phenomenon, and we will explore this in more detail in Section 4.1 based on the disk model presented in the subsequent section.

The Physical Conditions for Maser Pumping
In the interstellar medium, the − + = 6 16 − 5 23 water maser transition can occur naturally through collisional pumping in a warm molecular cloud.In the limit where the medium is optically thick to the far-infrared continuum, the water level populations for a saturated maser source 7 are mainly determined by the local conditions including the gas temperature ( H 2 ), the dust temperature d , the number density ( H 2 ) of molecular gas, and the water abundance (Neufeld 2000).Among these local conditions, the most important ones for maser pumping are H 2 and H 2 , which need to fall within the favored ranges of 400 H 2 1500 K and 10 7 H 2 10 11 cm −3 , respectively (e.g.Neufeld et al. 1994;   7 In the on-going analysis of the H 2 O maser spectra (Kuo et al. in prep.) for disk maser systems listed in Table 1, we find no concrete evidence for line narrowing (e.g.Baudry et al. 2023) which would indicate unsaturated masers in the majority of the systems.Therefore, we make the simple assumption that all maser sources are saturated for the discussion in the rest of the paer.Neufeld 2000;Herrnstein et al. 2005;Lo 2005;Gray et al. 2016).The minimum gas temperature of min ∼ 400 K is required for efficient collisional pumping to the − + = 6 16 level lying at E/k = 643 K above the ground level.Moreover, this temperature threshold is also essential to make a significant enhancement of the water abundance in the cloud ( H 2 O 10 −4 ; e.g.Neufeld 2000) with the reaction O + H 2 → OH + H followed by OH + H 2 → H 2 O + H (e.g.Elitzur & de Jong 1978;Neufeld et al. 1994).Provided that the gas density is below the critical value ( crit 10 11 cm −3 ) for collisional de-excitation, a sufficiently high density ( H 2 10 7 cm −3 ) is also crucial to make maser pumping efficient.
In addition to the suitable physical conditions of the gas, the presence of cold dust in the gas cloud is also important for producing the large maser luminosities (Kuo et al. 2018b) observed in water megamasers (Lo 2005).Given that the gas temperature and density fall in the preferred ranges for population inversion, Collison & Watson (1995) demonstrated that the maser emission can be significantly enhanced if cold dust grains exist in the cloud with a dust temperature dust = H 2 − Δ , where Δ ∼ 50 − 100 K.The presence of such cold dust can absorb nonmasing far-infrared water lines trapped in the cloud, enabling a much larger extent of H 2 O molecules maintaining population inversion without being quenched (Goldreich & Kwan 1974;Collison & Watson 1995;Lo 2005).
To maintain a sufficient temperature for efficient maser pumping, it has been proposed that X-rays from the active nucleus could be the primary heating source for the masing gas (Neufeld et al. 1994;Neufeld & Maloney 1995;Herrnstein et al. 2005).In addition, spiral shock waves travelling through a circumnuclear disk (Maoz & McKee 1998) have also been proposed to pump the masers, but the velocity drifts among high-velocity features predicted by this model have not been detected in sensitive spectral-line observations of eleven H 2 O maser disks (Pesce et al. 2015), making this model less favorable (see also Bragg et al. 2000;Humphreys et al. 2008).
Finally, recent observations of the maser emission from NGC 4258 with space VLBI (Baan et al. 2022) suggest that viscous heating due to the magnetorotation instability (MRI; e.g.Armitage 2022) could also be one of the pumping mechanisms for H 2 O megamaser emission and may explain some spectral variations seen in NGC 4258.Nevertheless, including the MRI-based viscous heating in the modeling of gas disks that could produce maser emission is significantly beyond the scope of this paper and will be deferred to future work.As a consequence, we will focus in the following discussion on the scenario in which H 2 O maser emission arises in the X-ray-dominated region (XDR; Maloney et al. 1996;Lo 2005) in a gas disk, within which the heating of the gas and the chemical composition are dominated by X-rays.

The Role of X-ray Heating
Considering a circumnuclear disk illuminated by the central X-ray source, it is expected that the disk can receive X-ray heating for maser excitation most efficiently if the disk is warped, allowing one side of the disk plane to be irradiated by X-rays directly.In the wellknown picture of maser excitation in NGC 4258 (Neufeld & Maloney 1995, hereafter NM95), it was suggested that the outer edge of the warped maser disk is determined by the critical radius beyond which the disk is more directly illuminated and X-ray irradiation becomes strong enough to dissociate all molecules into atoms.By assuming the viscous gas disk to be in a steady state of accretion, NM95 show that this critical radius can be expressed as cr = 0.040 −0.426 41 where 10 41 41 ergs s −1 is the 2-10 keV X-ray luminosity of the central source, 10 −5 −5 / ⊙ yr −1 is the mass accretion rate normalized by the conventional ( 1) viscosity parameter, and 10 8 8 ⊙ is the BH mass of the maser system. in the above equation is the obliquity parameter defined as = cos , where is the angle at which the disk is illuminated obliquely with respect to the normal direction of the disk.
To explain the presence of the inner edge of the maser emissions in NGC 4258, NM95 observed that the maser disk appears to flatten out close to the inner radius (Miyoshi et al. 1995), suggesting that the obliquity parameter falls to zero and the disk is no longer directly illuminated by the X-ray source, making the gas too cold to mase.As a result, they speculated that the inner edge of the maser disk is determined by the nature of the warp.While the NM95 model seems to provide a plausible explanation for the inner and outer boundaries of the maser emissions in NGC 4258 (e.g.Martin 2008), it is not yet well explored whether this model can be applied to H 2 O maser disks discovered in the past two decades, whose intrinsic AGN luminosities are 2−3 orders of magnitude higher than that of NGC 4258.
Given the possibility that the H 2 O maser disks do not necessarily follow steady-state accretion (e.g.Gammie et al. 1999), Gao et al. (2017) proposed a simple alternative model in which the outer edge of a maser disk is determined by the radius beyond which H 2 400 K. Assuming the gas and dust are well-coupled in a maser disk, Gao et al. (2017) uses the dust temperature d in the opticallythin limit as a proxy for H 2 , finding that the gas temperature is a decreasing function of radius, with the outer radius described by out ∝ 1/2 bol , where bol is the bolometric AGN luminosity.While the derived scaling between out and bol is broadly consistent with the observations, we find that the simple approximation of out in Gao et al. ( 2017) can be further improved by considering the effect of X-ray heating on the gas.It is well-known that the collisions between gas and dust particles do not guarantee d ≈ H 2 in an X-ray irradiated gas.The gas and dust can coexist at substantially different temperatures, with the difference varying depending on the the X-ray heating rate (e.g.Desch et al. 1998;Lo 2005;Namekata & Umemura 2016).In addition, in the region where the gas density is within the favored range for maser pumping, the obscuring column density is required to be large enough (e.g.H 10 23 cm −2 ; Kartje et al. 1999) to avoid molecular dissociation (Neufeld et al. 1994;Neufeld 2000).As shown in Namekata & Umemura (2016), photoheating and thermal emission from dust grains in such a regime is unimportant due to large optical depth.The dust temperature is mainly determined by collisional energy transfer between gas and dust, with d 200 − 300 K (see Section 2 in Namekata & Umemura 2016), suggesting that using d in the optically-thin limit to estimate H 2 may lead to non-negligible systematic errors.To obtain a more reliable estimate of out for a maser disk based on the physical conditions of the interstellar medium, it would be beneficial to explore the gas temperature distribution directly with the X-ray heating rate.

Are Disk Radii Solely Determined by X-ray Irradiation ?
To examine whether the outer edge of H 2 O maser emissions in the maser disks in our sample lie at the molecular-to-atomic transition radius cr as suggested by the NM95 model, we apply Equation 1 to all sixteen maser disks listed in Table 1 and compare the predictions with the observed values.When evaluating cr for our sample, we first estimate the 2-10 keV X-ray luminosity 2−10 X = 10 41 41 ergs s −1 by assuming 2−10 X = 0.1 bol .For all sources except for ESO 558-G009 and CGCG 074-064, we inferred bol from the reddening-corrected [OIII] luminosities, [OIII] , of the maser galaxies using bolometric corrections from Heckman et al. (2004) and Heckman & Best (2014) (see Kuo et al. 2020a).For the two exceptions which do not have [OIII] available, we estimate the bol based on the mid-infrared AGN luminosities derived from SED-fitting (Kuo et al. 2018b), with the bolometric correction given by Gandhi et al. (2009).Since NM95 assumes a steady-state accretion disk, suggesting that the accretion rate is constant with radius, we calculate the mass accretion rate with = bol / 2 , where is the accretion efficiency for a Kerr black hole (e.g.Caproni et al. 2007), which ranges from 5.7% to 42% for the BH spin between 0 ≤ * ≤ 1.To obtain the obliquity parameter , we performed 3-dimensional Bayesian modeling for all maser disks except NGC 5495 and NGC 1068 using the MCP modeling code described in Reid et al. (2013), Humphreys et al. (2013), Kuo et al. (2015), and Gao et al. (2016).These modelings8 are described by 15 global fitting parameters that include the galaxy distance , BH mass BH , the position of the dynamical center, and parameters characterizing the inclination warp and position angle warp (see Equations ( 1) & (2) in Kuo et al. 2013).For all cases, the warped disk models can well fit the observational data of the maser systems with reduced 2 ∼ 1 (e.g.Kuo et al. 2018a).The best-fit warp parameters enable us to calculate the obliquity parameter with = ˆ • ˆ , where ˆ and ˆ are the unit vectors of the disk normal and the unit radial vector for a high-velocity maser spot located at a radius , respectively (Herrnstein et al. 2005).For NGC 5495 and NGC 1068, we simply make a crude estimation of based on their maser maps (e.g.Gao et al. 2017;Greenhill & Gwinn 1997;Gallimore & Impellizzeri 2023), because the quality of the VLBI map for NGC 5495 is not good enough for a reliable 3-D modeling whereas NGC 1068 may have gas distribution and kinematics that deviate from the simple assumptions made in our code.
In Columns ( 8), (9), & (10) in Table 1, we list the obliquity parameters for the inner ( in ) and outer ( out ) edges of each maser disk as well as the ratio = in / out .Based on the values of in / out , we see no evidence of disk flattening toward the inner edges of all maser disks (i.e. in / out → 0) except for UGC 3789.The ratio is of order unity for the majority of the H 2 O maser disks, implying that the X-ray heating at the inner edge of the disk would not be significantly less than that at the outer edge.Therefore, under the conjecture proposed by NM95, it is difficult to account for the presence of the inner edges of most maser disks in our sample.
In Figure 2, we compare the critical radius cr with the observed out for each maser disk, with cr calculated based on the assumption of = 0.25 and = 0.42 (i.e. the maximum accretion efficiency corresponding to BH spin a * = 1; Caproni et al. 2007).The dashed line shows the locations where out = cr .The error bar in cr indicates the possible range of the critical radius given the preferred range of ∼ 0.1 − 0.4 suggested from observations (e.g.King et al. 2007;Kotko & Lasota 2012).This comparison shows that except for three sources including NGC 4258 (the red filled square in Figure 2), the calculated molecular-to-atomic transition radius is substantially larger than the observed outer radius for the majority of our H 2 O maser disks.
We note that it is not unexpected for maser emission to lie well The plot on the top of the figure shows a representative warped maser disk, with the black circle and the red/blue lines indicating the black hole and the redshifted/blueshifted maser features in the maser disk, respectively.To evaluate the X-ray heating rate in the gas disk, we choose the coordinates such that the mid-plane of an approximately flat section of the warped disk lies in the = 0 plane (see the plot in the inset).The X-rays from the black hole vicinity (the green arrow) are incident on the disk at an angle with respect to the normal direction of the disk, where = ( ) is a function of radius due to the warped nature of the disk.In this configuration, the gas clumps (represented by the grey circles) on the < 0 side of the disk are expected to see more X-ray obscuring material compared to those on the > 0 side of the disk.
within the transition radius because the presence of molecules is of course a necessary condition for maser action.Nevertheless, the systematic difference between cr and out for the majority of the maser disks suggests that it is not the boundary of molecular dissociation in these systems, but rather some other effects that are more relevant for setting the outer radius of the maser disk.It is likely that physical conditions of the gas, such as H 2 and H 2 , may have become unfavorable for supporting maser action before reaching the radius of molecular-to-atomic transition.Alternatively, it is also possible that some of the basic assumptions in the NM95 model, such as steady-state accretion, may not hold (e.g.Gammie et al. 1999) in some maser disks, requiring models that relax these assumptions for a better explanation of the observations.

The Outer Radius Determined by the Physical Conditions of the Gas
Considering a dusty, warm medium where the molecular gas is mainly heated by X-rays, instead of using d to approximate H 2 as in Gao et al. (2017), we use the X-ray heating rate to probe the gas temperature in a gas disk.This allows us to examine whether the observed outer edge of a maser disk can be explained by the minimum temperature requirement for efficient maser action as suggested by Gao et al. (2017).
As discussed in the beginning of the section, the water level populations for a saturated maser are primarily determined by H 2 , H 2 , d , and H 2 O in a medium optically thick to far-infrared continuum emission.If one considers an X-ray irradiated medium within which water maser emission can occur, assuming a typical value of the water abundance (e.g.H 2 O ∼ 10 −4 ) that would arise at H 2 400 K, it has been shown by Neufeld (2000) that one can use the equilibrium X-ray heating rate per hydrogen nucleus X / H to replace the gas temper-  Kuo et al. (2020a), which are estimated from the [OIII] luminosity using the bolometric corrections from Heckman et al. (2004) and Heckman & Best (2014).For ESO558-G009 and CGCG074-064, we estimate their bol based on the mid-infrared AGN luminosity derived from SED-fitting (?) with bolometric correction provided by Gandhi et al. (2009); Column (4)−( 6): the mass of the disk in units of 10 4 ⊙ within the outer radius of the maser disk for three representative power-law indices = −1.5, = −1.0,and = −0.5, respectively; Column (7): the disk mass in units of 10 4 ⊙ within 1 pc of the disk assuming the = −1 model; Column (8): the mass ratio between the disk within 1 pc and the black hole; Column (9): the dust sublimation radius prescribed by Nenkova et al. (2008a) ature H 2 as one of the four key variables that determine the level populations (e.g.Wallin & Watson 1997;Babkovskaia & Poutanen 2004).
In their work (i.e.Neufeld 2000), it is demonstrated that the maximum heating rate that allows for efficient maser action is ( X / H ) max ∼ 1.2 × 10 −28 ergs cm 3 s −1 , beyond which the gas will be subject to molecular dissociation.In addition, the minimum heating rate ( X / H ) min that ensures H 2 400 K ranges from ∼5.0 × 10 −31 ergs cm 3 s −1 to ∼6.3 × 10 −30 ergs cm 3 s −1 , depending on d .If the heating rate of a gas is between ( X / H ) min and ( X / H ) max , it is expected that the gas temperature would fall within H 2 ∼ 400 − 1500 K, the favored range for maser excitation.In addition, H 2 would be an increasing function of X / H if the dust temperature and water abundance are roughly constant in the region.In the following analysis, we will assume d ∼ 300 K and H 2 O ∼ 10 −4 in the masing medium, suggesting ( X / H ) min ∼ 3.2 × 10 −30 ergs cm 3 s −1 (see Figure 1 in Neufeld 2000).For the purpose of determining the outer radius of a maser disk, we will simply identify the region within a disk where the X-ray heating rate and gas density fall within the allowed ranges that permit efficient maser action.We will not evaluate the exact gas temperature in the masing region because it will not affect our conclusion.
Following Maloney et al. (1996), we compute the X-ray heating rate per hydrogen nucleus for a gas cloud in the X-ray dominated region as where H is the density of the hydrogen nuclei and eff is the effective ionization parameter defined as eff = 1.26 × 10 −4 X 5 0.9 22 . (3) Here, X = X /4 2 is the unattenuated X-ray flux received by a gas clump at a radius from the central black hole with an intrinsic 1−100 keV X-ray luminosity of X .The two terms in the denominator show the normalized values for the total densities of hydrogen nuclei 5 = H /10 5 cm −3 and the X-ray attenuating column for the gas clump 22 = H /10 22 cm −2 .
To explore the region in the disk where the physical conditions of the gas are suitable for maser excitation, we consider the scenario in which an initially flat, cold molecular disk is warped by a certain mechanism such as resonant relaxation (e.g.Bregman & Alexander 2009, 2012), which can warp a sub-parsec scale disk efficiently in a timescale of ∼10 7 years.After the disk gets warped, one side of it will be subject to direct X-ray illumination that would deposit thermal energy in the disk.To evaluate the density distribution, we assume that the gas motion is dominated by turbulence (Wallin et al. 1998(Wallin et al. , 1999)), with the turbulence velocity g ∼ 2 km s −1 , the typical width of the water maser lines seen in Keplerian maser disks (e.g.Reid et al. 2013;Kuo et al. 2013;Gao et al. 2016).
To facilitate a simpler calculation, we describe the gas distribution in cylindrical polar coordinates ì = ( , , ) and adopt an orientation of the coordinates such that the approximately flat section of the warped disk that could produce maser emissions is aligned with the = 0 plane of the coordinate system (see Figure 3 for an illustration).Following Neufeld & Maloney (1995) and Herrnstein et al. (2005), we take into account the effect of the warp on the X-ray illumination by having the central radiation incident on the disk plane at an angle with respect to the normal direction of the disk, with ≡ cos −1 where is the obliquity parameter determined by the disk warp parameters as described in Section 3.3.In this geometrical setup, the BH lies slightly above the = 0 plane, and the effect of the variation of the warping angle as a function of radius can be accounted for by allowing the X-ray incident angle = ( ) to change with radius.
Since most maser disks are slightly warped and the ratio between in and out is of order unity for the majority of our maser disks, we simply assume = out in our calculation of the X-ray heating rate in the disk.The variations of between inner and outer radii of the maser disks are ignored in the calculation because the gradient term / (or / ) only leads to a mild change of the heating rate distribution within the outer radius of the maser disk and does not result in any change of the conclusion of this paper.
Because of hydrostatic equilibrium, the gas density ( , ) at a radius from the central black hole and an elevation above the mid-plane of the slightly warped, geometrically-thin disk can be expressed as where mid ( ) is the gas density at the mid-plane and ( ) is the scale height of the disk, given by ( ) = g ( BH ) 1/2 3/2 (Armitage 2022).The mid-plane density can be calculated with mid ( ) = Σ( )/[(2 ) 1/2 ] where Σ( ) is the surface density of the disk.Here, we assume that the surface density takes the form of where the power-law index is allowed to vary in the interval of −2 < < 0 and Σ 0 is the surface density at the arbitrarily chosen reference radius 0 of the gas disk (Caproni et al. 2007;Huré et al. 2011;Kuo et al. 2018a).If D represents the disk mass within 0 , one finds that Σ 0 = ( + 2) D /2 2 0 .In the subsequent analysis, we choose the reference radius 0 to be 0 = out , and D thus indicates the disk mass within the outer boundary of the maser disk.
In the context of this power-law surface density profile, the steadystate accretion disk adopted by NM95 corresponds to = −1.5 (Herrnstein et al. 2005;Kuo et al. 2018a), with D expressed as where is the mass accretion rate that can be inferred with = bol / 2 .Note that the disk mass D in the steady-state model can be fully determined by BH , bol , g , and whereas D in the power-law model described by Equation 5 is treated as a free parameter in our modeling.As a result, the analysis based on the power-law model would not necessarily give the same disk mass as derived from the steady-state scenario even if one adopts = −1.5.If one finds that the power-law model assuming = −1.5 suggests a significantly different disk mass from the value predicted by Equation 6, it may imply a deviation from the steady-state accretion.
Given the above expressions, one can evaluate the number density of the molecular gas in the disk with where H 2 = 2.36 is the mean molecular weight per hydrogen molecule (Fischera & Dopita 2008) and H is the mass of the hydrogen atom.Assuming a parcel of gas directly irradiated by the X-ray photons from the > 0 side of the disk is located at the position ì = ( , , ), we calculate the shielding column density of hydrogen nuclei in between the gas and the X-ray source as where the obliquity parameter accounts for the increase in the obscuring column density due to the fact that the disk is illuminated obliquely (e.g.Neufeld & Maloney 1995;Herrnstein et al. 2005), Erf(X) is the standard error function and H ( , ) = 2 H 2 ( , ) is the number of H nuclei.For the region where Thomson scattering becomes important and enhances X-ray obscuration (i.e. the Comptonthick regime; H ( , ) ≥ 1.5×10 24 cm −2 ), we adopt the effective column density eff H ( , ) = T H for evaluating the attenuated Xray heating rate, where the boosting factor is T ∼ 6.65 × 10 −25 H (Neufeld 2000).
Finally, in our calculation of X / H with Equations 2 and 3, we assume that the X-rays that originate from the disk corona in the vicinity of the central BH were emitted isotropically (Namekata & Umemura 2016) and the 1 − 100 keV X-rays account for ∼20% of the bolometric flux (e.g.Neufeld et al. 1994;Netzer 2019), suggesting that X = 0.2 bol /4 2 .In addition to estimating the heating rate, we also compute the density distribution of the molecular gas with Equation 7 to find out the region where H 2 ( , ) falls in the favored range for maser excitation (i.e.H 2 = 10 7 − 10 11 cm −3 ).By comparing this region with the locations where the X-ray heating rate is sufficient to maintain H 2 ∼ 400 − 1500 K, we identify the boundaries of the region within which the level population of water molecules could be inverted.The outer radius out of a maser disk in our model is then defined as the outermost radius of the region within which both gas temperature and density fall within the favored ranges for population inversion.We call such a region the masing region.

The Location and Outer Boundary of the Masing Region
Given the possibility that not all maser disks necessarily follow steady-state accretion, we adopt the more general power-law disk model prescribed by Equation 5 for our analysis, with D and treated as free parameters.We model every maser disk by choosing a set of ( , D ), with varying between −2.0 < < 0.0 in steps of 0.1.For a given value of , our model suggests that the outer radius of the masing region is an increasing function of D .By using the observed outer radius as the constraint and fixing at a specific value, we solve for the disk mass that can predict out consistent with the observation.
Our analysis shows that solutions of D can be found for all values of between −2.0 and 0. In Column (4), ( 5), & (6) of Table 2, we list the solutions for D in units of 10 4 ⊙ for three representative models with = −1.5, = −1.0,and = −0.5, respectively, which are selected to show the typical variations of D when the power-law index varies by Δ ∼ 0.5 − 1.To compare the disk mass within the same reference radius, we also show in Column (7) of Table 2 the total disk mass ˜ D within = 1 pc for each maser system.In this comparison, we adopt the model of = −1 (i.e. the Mestel disk; Mestel 1963), which is often adopted in the study of accretion disks (e.g.Huré et al. 2011;Kuo et al. 2018a;Gallimore & Impellizzeri 2023).Given this model, the disk mass within the reference radius of 0 = 1 pc is evaluated as ˜ D = D (1 pc/ out ).
One can see that the disk masses within 1 pc are comparable for the majority of the maser systems and the disk-to-BH-mass ratios ˜ D / BH shown in Column ( 8) are all significantly smaller than unity, with a mean value of ˜ D / BH ∼ 0.0018.This is consistent with the fact that most maser disks in our sample follow nearly perfect Keplerian rotation.Finally, we also note that the solutions for D Figure 4.The distribution of X-ray heating rate per hydrogen nucleus ( X / H ) and the number density of molecular gas ( H 2 ) for each of the H 2 O megamaser disks in our sample.The horizontal axis gives the disk radius and the vertical axis designates the disk elevation normalized to the scale height .The green dotted and dot-dashed lines indicate the locations in the disk where H 2 = 1 × 10 7 cm −3 and 1 × 10 11 cm −3 , the minimum and maximum density for population inversion, respectively.The black solid line represent the points where the heating rate reaches the critical value for molecular dissociation.The dashed line shows the curve of the minimum heating rate that ensures H 2 400 K.The blue shaded area in each plot shows the region in which both the gas density and temperature fall into the favored ranges for population inversion.The observed inner and outer radii of the maser disk are shown by the two vertical grey bars in each plot.Figure 5.The representative X-ray heating rate per hydrogen ( X ≡ X / H ; the blue curve) and obscuring column density ( H ; the red dashed line) as functions of the normalized disk elevation /H calculated with the disk parameters of NGC 2273 at its outer radius.On the y-axis, the heating rate X is normalized by X (max), which represents the maximum heating rate that enables the presence of molecules.The column density H increases toward the − direction at a given radius because the < 0 side of the disk sees more obscuring materials than the > 0 side in a disk illuminated by X-rays obliquely (see Figure 3).At 0, X decreases rapidly with − as H rises quickly (see Equations 2 & 3).The heating rate reaches the minimum when the exponential drop of H in the < 0 side of the disk starts to dominate over the increase in H .It is the interplay between H ( ) and H ( ) at a given radius that leads to the negative peak position and slight asymmetry of the heating rate curve.are substantially smaller than the predictions from the steady-state accretion model, suggesting that the steady-state assumption may be in question.Given bol ∼ 10 44 ergs s −1 and BH ∼ 10 7 ⊙ for most maser disks (see Tables 1 and 2), Equation 6predicts that ˜ D would be ∼1.0 × 10 6 ⊙ if the accretion is in the steady state assuming = 0.42 and = 0.25.This is ∼1 − 2 orders of magnitudes greater than our solutions for D , suggesting that the steady-state accretion disk is in general too massive to yield an outer radius small enough to be consistent with those found in the observations.
To illustrate how the physical conditions of the gas define the boundaries of the masing region, we show in Figure 4 the X-ray heating rate and gas density distributions for each maser disk based on the Mestel model (i.e.= −1).In each panel, the two vertical grey bars indicate the observed inner and outer radii of the maser disk.The green dotted and dot-dashed lines delineate the locations for the minimum and maximum density for maser excitation, respectively.The black solid line represents the points where the heating rate reaches ( X / H ) max while the black dashed line shows the curve for ( X / H ) min .It is expected that the gas lying beyond the maximum heating curve will become atomic, and the gas bound within the minimum heating curve will be too cold to mase.To produce luminous maser emissions, the gas needs to lie within the masing region marked by the blue shaded area, in which the gas density and temperature would both fall into the favored ranges for population inversion simultaneously.Outside the masing region, either the gas temperature or density is expected to fall outside the preferred range.
To see why the maximum and minimum heating rate curves are roughly symmetric functions which peak at a negative value ∼ −0.7, one can solve for the analytical solutions for the heating rate curves9 by using Equations 2 through 8.This allows one to express the radius of the maximum/minimum heating rate as a function of the disk elevation as H ( ) = ( X ) 0.44 (D) −0.43( / √ 2 ) 2 Erfc 0.78 ( / √ 2 ) pc, ( 9) where Erfc( / √ 2H) is the complementary error function, ( X ) refers to the maximum (or minimum) heating rate per nucleus ( X / H ) m normalized by 3.19 × 10 −28 ergs cm 3 s −1 , and (D) is a function of the accretion disk parameters expressed as : Given a set of disk parameters, H ( ) is dominated by the Gaussian term −0.43( / √ 2 ) 2 , suggesting that it is an approximately symmetric function of .By solving for the equation H out / = 0, one obtains the peak position of the function, which occurs at / = −0.76,displaced from the = 0 plane by the same amount regardless of the disk parameters.The negative displacement in the disk elevation does not change with the magnitude of the disk warp because the obliquity parameter is independent of in the thin disk approximation assumed in our calculation.To see more intuitively the physical origin of the negative displacement and the slight asymmetry of heating rate curves, we refer the readers to the illustration in Figure 5.
In Figure 4, one can also see that the masing region typically lies close to the mid-plane of the disk except for NGC 1194, and the thickness of the masing region between in and out is typically ∼1−4 .In addition, one can also infer that the X-ray heating rate per nucleus always increases with radius in every maser disk, suggesting that the gas temperature would increase with radius, while assuming d is roughly constant.As a result, we argue that the minimum temperature min ∼ 400 K could not be the primary factor that determines out as speculated in Gao et al. (2017).As suggested by Figure 4, the outer radius of a maser disk is primarily determined either by the maximum heating rate (e.g.NGC 4258, UGC 3789, etc.) or the minimum gas density min for maser pumping (e.g.NGC 2960, NGC 5765b, etc.), depending on the combination of bol , BH , D , and the obliquity parameter .

The Critical Black Hole Mass
To identify the most important factor that determines whether the outer radius of a maser disk is primarily limited by the maximum heating rate ( X / H ) max or the minimum gas density min , it is helpful to look for the analytical expressions for the outermost radii of the curves for ( X / H ) max and min shown in consistent with the empirical relationship out ∝ 0.57±0.16BH found by Gao et al. (2017).By equating H out with D out , it can be seen that the outer radii given by the above two equations would be approximately the same ( This critical BH mass separates the BH mass of the maser disk into two regimes.In the regime where BH << crit BH , the outer radius of the disk is primarily limited by the maximum heating rate, with the outer radius determined by H out .When the BH mass of the system is BH >> crit BH , the maser disk falls into the minimum density limited regime where the outer radius is determined by D out .So, given a set of disk parameters, the critical BH mass plays a decisive role in determining whether the outer boundary of the maser disk is mainly confined by ( X / H ) max or min .To identify the regime for each maser system, we label our sample sources with either "H" or "D" in Column (10) in Table 2 to indicate whether they fall into the heating rate limited regime (H) or minimum density limited regime (D).We note that for some heating rate limited sources, their D out are comparable to H out ( D out H out ), suggesting that their BH are close to the critical values.We label these sources with 'H(D)' to indicate that D out is almost as limiting as H out , and they could appear to follow both scaling relations mentioned above because they lie at the transition region between the two regimes.For the special case (NGC 3393) in which D out is nearly the same as H out , we label it with 'H/D' to show that both limiting factors play equally important roles for determining out .
The existence of the two regimes separated by crit BH gives an interesting implication that the seemingly disparate results from Wardle & Yusef-Zadeh (2012) and from Gao et al. (2017) can be unified once one understands that these two regimes of BH mass exist.Indeed, when we examine the eight maser disks studied in Wardle & Yusef-Zadeh (2012), except for NGC 4388 which is not included in our modeling, we find that 5 of the 7 remaining sources fall in the heating rate limited regime (i.e.NGC 4258, NGC 2273, UGC 3789, NGC 6323, and NGC 6264), explaining the out ∝ BH scaling revealed in their work.On the other hand, since Gao et al. (2017) include several additional sources in the density limited regime (e.g.ESO 558-G009, NGC 5765b, CGCG 074-064, and NGC 1194) as well as sources for which BH ≈ crit BH (e.g.IC 2560 and NGC 5495), it is understandable that Gao et al. (2017) would reach a different scaling between out and BH .
Finally, the readers should be aware that the outer radius of a maser disk predicted by our model should be seen as an approximation.For the maser action to take place, the H 2 O molecules are first pumped to a shortlived state at Level 3, followed by spontaneous transitions into a long-lived state at Level 2. The E 3 − E 2 transition that releases far-infrared photons needs to be fast enough to maintain the maser population inversion between E 2 and E 1 .To avoid having the far-infrared photons re-absorbed from Level 2 to Level 3, making the population moving into Level 2 effectively slow, it is useful if the emitted E 3 − E 2 photons are removed from the maser pump cycle e.g., through absorption by cold dust.If the dust has been eliminated by sublimation, this possibility is removed and the maser population inversion can be affected or even destroyed due to the thermalization of the population (e.g.Lo 2005).
In our simple modeling presented above, we ignore the effect of density perturbation of the molecular gas as a result of X-ray heating.Given the turbulence-dominated disk in our analysis, we implicitly assume that the equilibrium density of the molecular gas after X-ray irradiation is comparable to the gas density before any injection of thermal energy.Taking the density perturbation into account would require a more rigorous modeling that involves solving the energy balance equation for the gas.This is beyond the scope of this paper, and we defer this to future work.

The Inner Edge and the Dust Sublimation Radius
One can see in Figure 4 that the density and heating rate requirements for efficient maser pumping do not impose a clear inner bound for the maser emissions in the maser disk.This situation does not change no matter how we vary the model parameters, suggesting that the physical conditions of the gas alone may not be able to define the inner boundary of a maser disk.Other factors might be involved.
Among the additional factors that could affect the production of water maser emissions, dust properties are the most important.As explained in Section 3.1, the presence of cold dust is essential for maintaining population inversion by absorbing nonmasing far-infrared water lines trapped in the cloud.If the amount of dust is substantially reduced at some radius, for example because of dust sublimation, it is possible that the trapping of the far-infrared photons from nonmasing water transitions may become more significant (see Figure 6 for an illustration).This process could substantially diminish or even destroy the population inversion, leading to a boundary of the maser emissions in the disk (e.g.Greenhill et al. 1996;Kartje et al. 1999;Gallimore et al. 2001).
To explore this possibility, we estimate the dust sublimation radius for each maser disk (see Column (9) of Table 2) based on the prescription provided by Nenkova et al. (2008a): indicates where in = sub .Except for NGC 4258, which is marked by the filled red circle, the in values are comparable to sub for our maser disks, given the measurement uncertainties.The yellow square denotes NGC 1068, whose in appears to be an outlier in the right panel of Figure 1, but its inner radius is well consistent with the dust sublimation radius as shown in this figure .which is derived for the clumpy torus models assuming the grain mix has standard interstellar properties (Nenkova et al. 2008b).In our calculation, we use bol from Table 2, with the assumption that dust sublimates at the temperature sub ∼ 1500 K.We estimate the uncertainty sub,Nenkova by adopting bol ∼ 0.54 dex, obtained from the comparison between bol estimated from X-ray spectroscopy and the [OIII] 5007 line for our sample (Kuo et al. 2020a).As noted in Nenkova et al. (2008a), sub,Nenkova is not a sharp boundary within which no dust exists.Instead, it is an approximation that marks the radius across which the environment transitions from being dusty to dust-free as the individual components of the dust mixture gradually sublimate at different radii.It is expected that the largest grains can survive down to the innermost radius of the dusty torus probed by reverberation mapping observations (e.g.Suganuma et al. 2006;Kishimoto et al. 2007), which is ∼3 times smaller than sub,Nenkova .
In Figure 7, we compare the observed inner radius of the maser disk with the dust sublimation radius for all sources in our sample.This figure shows that the majority of the maser disks have their inner radii consistent with sub,Nenkova within ∼ 1 (∼75% of the sources) or ∼ 2 (∼19% of the sources) error.The only significant outlier in this comparison is NGC 4258 (the red square in the figure), for which in is considerably greater than sub,Nenkova given the measurement uncertainties.This result suggests that, except for NGC 4258 (see Section 4.2 for further discussion), dust sublimation may play a role in determining the inner radius of a maser disk.It is likely that dust residing near the sublimation radius would become warmer and its efficiency as a heat sink (Babkovskaia & Poutanen 2004) for absorb-ing far-infrared photons from nonmasing water transitions would get reduced.As a result, the efficiency of maser amplification could get significantly reduced as the dust temperature gradually approaches the sublimation temperature at ∼ sub,Nenkova (see Figure 1 in Gray et al. 2022), making the maser emissions substantially weaker.
In addition to the effect of dust temperature, the reduction of the dust mass could also affect the maser intensity significantly.As dust gradually sublimates away at ∼ sub,Nenkova , maintaining the population inversion becomes more difficult because the trapping of non-masing far-infrared water line emissions by the masing clouds becomes more significant.As shown in Gray et al. (2022), a halving of the dust mass fraction can lead to an extraordinary reduction in the efficiency of maser amplification, suggesting that even before the dust temperature is high enough for dust to totally sublimate away, maser emissions could have dimmed below the detection threshold.This supports the conjecture that dust properties near the dust sublimation radius may become unfavorable for efficient maser action, making dust the culprit for setting the observed inner boundary of a maser disk.

The Primary Physical Factors that Affect the Characteristic Disk Radii
The results of our calculations in the last section suggest that the outer radius of a maser disk could be determined either by the minimum gas density or the maximum heating rate that enables efficient maser action.Moreover, the inner edge of the maser emissions in the disk may result from the quenching of population inversion near the dust sublimation radius.While the physical conditions of the gas and the dust appear to play an important role in defining the inner and outer boundaries of a maser disk, these conditions do not seem to directly explain why the inner and outer radii of the majority of the maser disks are around in ∼ 1.0×10 5 S and out ∼ 3.0×10 5 S as shown in Section 2.
To explain the characteristic inner radius of in ∼ 1 × 10 5 S in light of the dust sublimation radius, it is helpful to replace bol in Equation 14 with bol = Edd Edd , where Edd is the Eddington ratio and Edd = 1.26×10 38 ( BH / ⊙ ) is the Eddington luminosity.This replacement allows one to express the dust sublimation radius in units of S as sub,Nenkova = 1.05 × 10 5 Edd 0.05 The above equation directly implies that the dust sublimation radius in AGN in general does not equal ∼1.0 × 10 5 S since its value depends substantially upon Edd and BH .However, it is well known that the BH masses for the Keplerian megamaser disks are typically ∼10 7 ⊙ , likely resulting from the fact that disk megamasers are preferentially detected in local Seyfert 2 galaxies (Kuo et al. 2011) whose BH mass function peaks at BH ≈ 3×10 7 ⊙ (Heckman et al. 2004).Moreover, the Eddington ratios of the majority of the Keplerian disk maser systems fall in the narrow range of Edd ∼ 0.01 − 0.1, with the median value ≈ 0.04 (see Table 5 and Figure 6 in Kuo et al. 2020a).Given that disk maser systems typically have BH ≈ 10 7 ⊙ and Edd ≈ 0.04, one can infer from Equation 15 that the dust sublimation radii of most Keplerian H 2 O maser disks would be ∼1.0 × 10 5 S , the characteristic inner radius we see in Figure 1.
To explain the characteristic outer radius, one can follow a similar approach to express the outer radii given by Equations 11 & 12 in units of S .For the heating rate limited regime, this gives where the normalization factor 0.0011 is the median value of ˜ D / BH for all sources.These equations suggest that the outer radii of the maser disks would be ∼3×10 5 S , given the typical values of Edd , ˜ D / BH , BH , and for the majority of the Keplerian maser disks.Although the dependence of out on the disk parameters is more complicated, given that these four parameters all fall within relatively narrow ranges for the majority of our maser sources, it is not unexpected that the outer radii of the maser disks have similar values when expressed in units of 10 5 S .
It is likely that the characteristic size of maser disks is deeply connected with the fine-tuning nature of the disk megamaser phenomenon and it reflects the physical state of a gas disk in a certain phase of AGN evolution, with BH masses following the population for low redshift Seyfert 2 galaxies.As suggested in Constantin (2012), the disk megamaser phenomena may only occur in a certain (short) phase in the galaxy-AGN coevolution during which the configuration of the gas accretion transitions from an optically thin, geometrically thick accretion flow to an optically thick, geometrically thin accretion disk, which typically have Eddington ratios Edd 0.01 (e.g.Giustini & Proga 2016).Given this speculation, the general Eddington ratio distributions for local Seyfert galaxies (e.g.Jones et al. 2016) would imply that the disk megamasers living in geometrically thin accretion disks with Edd ∼ 0.01 − 0.1 would outpopulate the ones with Edd 0.1, explaining the typical range of Edd for most Keplerian maser systems.
It can be expected that if H 2 O maser disks also exist at the centers of high redshift AGNs (e.g.quasars at > 2) or local Seyfert galaxies with high Eddington ratios, the combination of their Eddington ratio and BH mass distribution would be considerably different (e.g.Edd 0.1 and B 10 9 BH for high-z quasars), leading to distinctly different characteristic inner and outer radii for these sources.This could explain why NGC 1068 has the largest inner radius (i.e. in = 3.57 × 10 5 S ) among maser disks in our sample, which is ∼4 times greater than the median value for all sources (i.e. in = 0.87 × 10 5 S ; see Figure 1).Given BH = 1.7 × 10 7 ⊙ (see Table 1), the Eddington luminosity of NGC 1068 is Edd = 2.1×10 45 ergs s −1 .The corresponding Eddington ratios are 1.8 and 0.3 if its AGN bolomotric luminosities are inferred from [OIII] luminosity and absorption-corrected X-ray luminosities, respectively (Kuo et al. 2020a).Because of the substantially higher Eddington ratio in comparison with the Keplerian maser disks in our sample, NGC 1068 would have an inner radius substantially greater than the mean inner radius of most maser disks according to Equation 15.Nevertheless, similar to most other maser disks, NGC 1068 has an in well consistent with sub,Nenkova , supporting the conjecture that its inner boundary of maser emissions is constrained by the dust sublimation process.

The Inner Radius of NGC 4258
As shown in Section 3.7 and Figure 7, NGC 4258 is the most prominent outlier ( in >> sub,Nenkova ) in our comparison between in and sub,Nenkova for maser disks in our sample.While one can infer from Equation 15 that the substantially smaller dust sublimation radius in NGC 4258 ( sub,Nenkova = 0.011 pc) is caused by the combination of the extremely low Eddington ratio (i.e.Edd ∼ 0.0001; Kuo et al. 2020a) plus the slightly higher BH mass (i.e.BH = 4.0 × 10 7 ⊙ ; Humphreys et al. 2013), it is unclear why the inner boundary of the maser disk cannot reach down to sub,Nenkova based on what we learn in Section 3.This discrepancy hints that the inner edge of this well-known maser source depends on factors other than the physical conditions of the gas and dust in the disk.
We note that if the conditions that enable the population inversion are not the determining factors that define the inner boundary for NGC 4258, it is possible that its inner radius is characterized by the lower limit in the velocity coherence length that allows for observable maser emissions.By comparing the orientations of all maser disks in our sample, it can be seen that while most maser disks are within ∼ 1 • − 2 • from being edge-on (e.g.Kuo et al. 2011;Gao et al. 2016;Pesce et al. 2020a), NGC 4258 displays the most significant inclination warp.Based on the disk modeling by Humphreys et al. (2013), the disk inclination at a radius ≤ in = 0.11 pc is ≤79.2 • , suggesting that the deviation from the edge-on configuration would be greater than 10 • if there are masers residing inside the observed in .Given such large inclinations, it is possible that the effective coherent path length would be reduced substantially, leading to weak maser emissions below the detection limit.
In addition to the effect of disk inclination, we also note that the coherence length c reduces linearly with the disk radius for a maser disk (i.e.c ∝ ; Lo 2005).Since the maser flux density scales with c as ∝ 3 c , suggesting ∝ 3 (see Section 4.3 for a detailed discussion), it can be expected that the maser flux density would drop quickly when the disk radius is smaller than a certain value at which the physical conditions of the gas no longer support the maximum maser pumping rate (e.g.Gray et al. 2016).As a result, in a maser system with a low Eddington ratio (e.g.Edd << 0.01), the inner radius of the maser disk may not necessarily lie close to the dust sublimation radius because the maser coherence path length at sub,Nenkova may have become too small to produce observable maser emissions before reaching down to the dust sublimation radius.For NGC 4258, it is likely that both the relatively large disk inclination and the reduced coherence length at a smaller radius play a substantial role in causing the discrepancy between in and sub,Nenkova .To confirm this hypothesis requires future work with more rigorous modeling of the maser production rate as a function of disk radius based on the physical conditions of NGC 4258.

The Effect of Increasing Gain Length
Observations of the early universe have revealed strong evidence that galaxies once went through an extremely rapid and luminous phase of evolution, leading to intensive star formation and quasar activity that peak at ∼ 2−3 (e.g.Kotilainen et al. 2009;Lapi et al. 2017).During the most active phase of evolution, it is believed that a high fraction of luminous quasars may harbor 10 9 ⊙ supermassive BHs at their cores (e.g.Vestergaard & Osmer 2009;Li et al. 2021).Lo (2005) speculated that H 2 O gigamasers could be triggered in high-z active galaxies as a result of the immense energy injection into the surrounding gas, and the total maser luminosities of these gigamasers could be 10 3 higher than those of low-z H 2 O megamasers (see also Neufeld & Maloney 1995;Koekemoer et al. 1995; Falcke et al.Over the past two decades, there have been a few maser surveys aiming to detect masers from high-z galaxies (e.g.Impellizzeri et al. 2008;McKean et al. 2011).One of the surveys (Impellizzeri et al. 2008) detected the first 22 GHz H 2 O gigamaser system (MG J0414+0534) at > 1, confirming the existence of H 2 O gigamasers in the high-z universe.However, the analysis of the maser spectra for this source (Castangia et al. 2011) implies that the maser emissions are more likely associated with jet-cloud interaction rather than originating from a nearly edge-on disk, which is the configuration required for precise BH mass and distance measurement.To use the maser technique for studying cosmology, one still needs to look for luminous H 2 O gigamaser systems in disk configurations at high redshifts.
To assess the possibility of discovering high-z H 2 O gigamaser disks quantitatively, we estimate the flux densities of high-z maser disks based on the model presented in the last section.Our estimation suggests that existing radio telescopes, such as the Very Large Array (VLA) and the Green Bank Telescope (GBT), would have sufficient sensitivities to detect these gigamaser disks if their typical radii were 20 − 30 times greater than those of the low redshift maser systems (i.e.∼ 0.3 − 0.8 pc).
As one can infer from Lo (2005), the flux density of a saturated maser source can be expressed as where u is the density of H 2 O molecules in the upper excited state, Δ is the rate of maser transition from the upper to the lower state, ℎ is the energy of a maser photon, g is the gain length for maser am-plification, L is the luminosity distance to the source, and Δ is the observing bandwidth at the observer's frame.The cubic dependence on g shown in the equation suggests that the maser flux density is highly sensitive to the gain length, which is expected to be g c , where c is the velocity coherent path length in a maser disk.In the following discussion, we will make the simple assumption that g ≈ c when trying to provide an optimistic estimation of the detectability of high-z maser disks.The readers should be aware of the caveat that g could be substantially smaller than c in cases where the gas disks are highly clumpy.In such disks, the gain length could be limited to the size of each individual cloud, and g would become long enough to produce luminous maser emissions detectable at high redshifts only when multiple clumpy masing clouds happen to align along the line-of-sight.
Since maser action only takes place when the gas density and temperature fall within the narrow ranges shown in Section 3.1, one could assume that, on average, u Δ in the high-z environment would be comparable to that in local maser sources.Given this assumption and considering the optimal gain length g ≈ c , the key parameters that determine the maser flux density would be the coherence path length c of the disk and the obvious inverse square dependence of the luminosity distance L .
For the high-velocity maser components in a maser disk, the velocity coherent path length c at the tangent points at radius is c = 2( / ) 1/2 , where and are the gas velocity dispersion and the orbital velocity, respectively, suggesting that the gain length would increase linearly with (Lo 2005).If the radius of a maser disk were increased by a factor of 20 − 30, the isotropic luminosity density of the source ≡ 4 2 L /(1 + ) would increase by a factor of 8000 − 27000 due to the increase in c , making the source an H 2 O gigamaser.
Assuming such H 2 O gigamaser disks exist at ∼ 2 − 3, their luminosity distances would be ∼ 150−250 times greater than a maser disk at the distance of ∼100 Mpc assuming standard cosmology.The inverse square dependence of L would then lead to a decrease in the flux density by a factor of ∼ 20000 − 60000.By considering the change in the flux density due to the factor (1 + ) and the increases in L and c based on Equation 18, one would expect that the flux densities of a 22 GHz H 2 O gigamasers at ∼ 2 − 3 could be comparable to that of a local H 2 O megamaser at ∼100 Mpc.Given that the flux densities of known H 2 O maser disks at L ∼ 100 Mpc are typically ∼20 − 40 mJy (e.g.Kuo et al. 2013Kuo et al. , 2015;;Gao et al. 2016), the expected flux densities of the strongest maser features in a high-z gigamaser could range from a few mJy up to 40 mJy, suggesting that a few detections of the strongest lines are possible with a few hours of on-source integration time using the VLA, the GBT, and the High Sensitivity Array (HSA).Note that before applying the above estimation to high-z sources, one needs to explore whether the high-z gigamaser disks could have significantly larger sizes than the ones in the local universe.

The Dependence of Maser Disk Size on the Black Hole Mass
Based on the observations of quasars in Vestergaard & Osmer (2009), it has been shown that the BH mass function of luminous quasars at redshifts 1.5 3 peaks at BH ∼ 1.5 × 10 9 ⊙ , and the Eddington ratios of these quasars tend to be Edd 0.1.To see whether maser disks could exist around these 10 9 ⊙ supermassive black holes, with their sizes significantly larger than the low redshift maser systems, we explore the dependence of maser disk size on black hole mass based on the disk model presented in Section 3.4.In this exploration, we adopt the Mestel disk profile Σ( ) = Σ out ( / 0 ) −1 with g = 2 km s −1 , 0 = 1 pc, and fix the obliquity parameter at = 0.15, a typical value for local H 2 O maser disks.In addition, we set the bolometric luminosity in our model as bol = Edd Edd , with Edd varied between 0.03 to 0.6.Finally, we consider five representative values for the disk-mass-to-BH-mass ratio ˜ D / BH that ranges from 0.001 to 0.03.Given a combination of Edd and ˜ D / BH , we model the heating rate and density distribution in the disk and calculate the outer radius of the masing region for BH mass between 1.0 × 10 7 ⊙ and 1.5 × 10 9 ⊙ .
In the left panel of Figure 8, we show out as a function of BH for Edd = 0.03, 0.06, 0.1, 0.3, and 0.6, with ˜ D / BH = 0.005, which is ∼3 times greater than the mean disk-to-BH mass ratio for our sample of low-z maser disks.The right panel of Figure 8 shows the prediction of the outer radius for ˜ D / BH ranging from 0.001 to 0.03, assuming Edd = 0.1.It can be seen in both plots that the outer radius out for the cases with D / BH 0.005 first increases with BH mass with a steeper slope, and the slope appears to drop distinctly when the BH mass is greater than a critical value that varies depending on the combination of BH , Edd , and ˜ D .We note that this critical value is simply the critical BH mass crit BH discussed in Section 3.6.The slope changes at BH ∼ crit because the maser disk transitions from the heating rate limited regime (i.e.

BH << crit
BH and H out ∝ BH ) to the minimum density limited regime (i.e.BH >> crit BH and D out ∝ 0.6 BH ).As indicated by Equation 13, the critical BH mass would be crit BH 10 8 ⊙ if 0.1 Edd 1 and ˜ D / BH 0.005, suggesting that the outer radius of the maser disk around a 10 9 ⊙ BH in a high-z quasar would be determined by D out if the gas disk is massive enough.As one can see in Figure 8, such high-z maser disks would have out ∼ 10 − 30 pc given ˜ D / BH ∼ 0.005 − 0.03, about ∼ 20 − 60 times greater than the average outer radius of the local maser disks (i.e.out = 0.52 pc).Based on the discussion in Section 4.3.1,we speculate that such large disk sizes could lead to H 2 O gigamasers in high redshift galaxies.
Considering the typical angular-diameter distance A for a galaxy at ∼ 2 − 3 (i.e.A ∼ 1620 − 1760 Mpc), the angular radius of a high-z maser disk with the physical size of ∼ 10 − 30 pc would be ∼ 1.2 − 3.8 milliarcseconds, comparable to those of the local H 2 O maser disks, suggesting that it is possible to apply the H 2 O maser technique to high-z quasars with existing centimeter VLBI facilities.If one could further detect submillimeter water maser emissions (e.g.Pesce et al. 2016Pesce et al. , 2023) ) from these high-z H 2 O gigamaser disks, it is possible that future observations with the Event Horizon Telescope (EHT) at submillimeter wavelengths would provide highly accurate maser imaging with ∼ 20 − 40 microarcsecond resolution at submillimeter wavelengths, leading to maser maps with fractional position uncertainties comparable to NGC 4258 (e.g.Argon et al. 2007;Humphreys et al. 2013).Depending on the redshift and the specific transition, these high-redshift submillimeter masers may even be observable with the Global Millimeter VLBI Array (GMVA) at millimeter wavelengths.For example, a 321 GHz, 325 GHz, and 380 GHz maser (Gray et al. 2016) at a redshift of = 3 would have an observed frequency of ∼80 GHz, ∼81 GHz,and ∼95 GHz, respectively, allowing for observations of these maser transitions with the 3mm receiver of the GMVA.

CONCLUSION
In this work, we combine ideas from multiple previous studies and develop a warped molecular disk model to examine the distributions of gas density and X-ray heating rate in the X-ray illuminated disk, with the gas disk described by a power-law surface density profile Σ( ) ∝ .This allows us to identify the boundaries of the masing region within which both gas density and temperature fall in the favored ranges for maser excitation.We compare the observed radii of sixteen maser disks with the predictions from our model, with the main conclusions summarized as follows : 1.With a suitable choice of disk mass, the predictions from our model agree reasonably well with the observed outer maser disk sizes from observations for all sixteen maser disks in our sample, with the agreement substantially better than previous models.
2. The outer edges of the maser disks are set by either the maximum X-ray heating or by the minimum gas density required for population inversion, depending on the critical black hole mass crit BH for a given system.
3. The critical BH mass crit BH is determined by the combination of Edd , BH , ˜ D , and of a maser system, and it separates the maser systems into the heating rate limited regime and the minimum density limited regime.This understanding explains two otherwise disparate previous results from Wardle & Yusef-Zadeh (2012, out ∝ BH ) and from Gao et al. (2017, out ∝ 0.57±0.16BH ), because it predicts two different scaling laws depending on which regime a maser system falls into.
4. The observed inner radii of the majority of the maser disks are roughly consistent with the dust sublimation radius, suggesting that in of a maser disk is probably set by dust sublimation rather than by disk warping effects as proposed by NM95, though NGC 4258 appears to be a significant outlier, presenting a potential counterexample that requires more investigation.
5. By applying our model to high-z quasars, our work predicts that H 2 O gigamasers should exist in the early universe, with their maser flux densities detectable with existing radio facilities.

Figure 1 .
Figure 1.Left panel: The maser disk size in units of parsecs.The top and bottom of each black vertical bar marked by the filled triangle and black circle indicate the outer ( out ) and inner ( in ) radii of a maser disk, respectively, with the blue and red error bars showing the uncertainties in the measurements of out and

Figure 2 .
Figure 2. The comparison between the observed outer radius out of each maser disk and the corresponding critical radius cr predicted from the NM95 model.The critical radius is calculated with = 0.25 and = 0.42, where and are the viscosity parameter and accretion efficiency, respectively (see Section 3.2 & 3.3).The black dashed line indicates the locations where out = cr .Excepting three sources including NGC 4258 (the filled red square), the maser disks have cr substantially greater than out .

Figure 3 .
Figure 3.The orientation of the warped disk in the chosen coordinate system.The plot on the top of the figure shows a representative warped maser disk, with the black circle and the red/blue lines indicating the black hole and the redshifted/blueshifted maser features in the maser disk, respectively.To evaluate the X-ray heating rate in the gas disk, we choose the coordinates such that the mid-plane of an approximately flat section of the warped disk lies in the = 0 plane (see the plot in the inset).The X-rays from the black hole vicinity (the green arrow) are incident on the disk at an angle with respect to the normal direction of the disk, where = ( ) is a function of radius due to the warped nature of the disk.In this configuration, the gas clumps (represented by the grey circles) on the < 0 side of the disk are expected to see more X-ray obscuring material compared to those on the > 0 side of the disk.
Figure 4. Based on Equations 9 and 10, one can easily show that the outer radius confined by the maximum heating can be expressed as -Zadeh (2012).By combining Equations 4, 5, and 7, one can find that the outer radius of the maser disk in the minimum density limited regime is

Figure 6 .
Figure6.The schematic diagram of the maser excitation process.For the maser action to take place, the H 2 O molecules are first pumped to a shortlived state at Level 3, followed by spontaneous transitions into a long-lived state at Level 2. The E 3 − E 2 transition that releases far-infrared photons needs to be fast enough to maintain the maser population inversion between E 2 and E 1 .To avoid having the far-infrared photons re-absorbed from Level 2 to Level 3, making the population moving into Level 2 effectively slow, it is useful if the emitted E 3 − E 2 photons are removed from the maser pump cycle e.g., through absorption by cold dust.If the dust has been eliminated by sublimation, this possibility is removed and the maser population inversion can be affected or even destroyed due to the thermalization of the population (e.g.Lo 2005).

Figure 7 .
Figure7.A comparison between the observed inner radius in and the dust sublimation radius sub of the H 2 O maser disks.The dashed line in the plot indicates where in = sub .Except for NGC 4258, which is marked by the filled red circle, the in values are comparable to sub for our maser disks, given the measurement uncertainties.The yellow square denotes NGC 1068, whose in appears to be an outlier in the right panel of Figure1, but its inner radius is well consistent with the dust sublimation radius as shown in this figure.

Figure 8 .
Figure 8. Left panel: the prediction of the outer radius of a megamaser disk with the BH mass ranging from 1.0 × 10 7 ⊙ to 1.5 × 10 9 ⊙ , made with the assumption of ˜ D / BH = 0.005, where ˜ D is the disk mass within one pc.The blue, cyan, green, yellow, and red lines represent the predictions assuming Eddington ratios of 0.03, 0.06, 0.1, 0.3, and 0.6 for the accreting BH, respectively; Right panel: The outer radius of a megamaser disk given an Eddington ratio of 0.1.The blue, cyan, green, yellow, and red lines indicate the predictions calculated based on ˜ D / BH = 0.001, 0.003, 0.005, 0.01, 0.03, respectively.

Table 2 .
Properties of the H 2 O Megamaser Systems ; Column (10): The regime into which a maser system falls.The letters H and D indicate that the outer boundaries of the maser disks are primarily limited by H out and D out , respectively; H(D) indicate the cases in which H out is the primary limiting factor for the outer radius but D out is just slightly larger than H out .H/D shows the special case in which H out and D out are approximately the same, implying that both heating rate and density requirements play important roles for limiting out .