$\langle$2023-11-23$\rangle$ $\langle$2024-02-09$\rangle$ $\langle$publication
date$\rangle$ Sub-mm Water Maser in NGC 1052NGC 1052

galaxies: active — galaxies: individual (NGC 1052) — galaxies: nuclei — masers

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

Seiji Kameno 11affiliation: Joint ALMA Observatory, Alonso de Córdova 3107
Vitacura, Santiago 763-0355, Chile 22affiliation: National Astronomical
Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
33affiliation: Department of Astronomy, School of Science, Graduate University
for Advanced Studies (SOKENDAI), Tokyo 181-8588, Japan Yuichi Harikane
44affiliation: Institute for Cosmic Ray Research, The University of Tokyo,
5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan Satoko Sawada-Satoh
55affiliation: Graduate School of Science, Osaka Metropolitan University, 1-1
Gakuen-cho, Naka-ku, Sakai, Osaka, 599-8531, Japan Tsuyoshi Sawada
11affiliationmark:  22affiliationmark:  Toshiki Saito 22affiliationmark:
Kouichiro Nakanishi 22affiliationmark:  33affiliationmark:  Elizabeth
Humphreys 11affiliationmark:  66affiliation: European Southern Observatory,
Karl-Schwarzschild-Strasse 2 D-85748, Garching, Germany and C. M. Violette
Impellizzeri77affiliation: Leiden Observatory, Leiden University, PO Box 9513,
2300 RA, Leiden, The Netherlands<EMAIL_ADDRESS>
###### Abstract

We report sub-pc-scale observations of the 321-GHz H2O emission line in the
radio galaxy NGC 1052. The H2O 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.

##  1 Introduction

The Atacama Large Millimeter/submillimeter Array (ALMA) detected H2O
$J_{Ka,Kc}=10_{2,9}-9_{3,6}$ emission at rest frequency of 321.225677 GHz 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${}_{\Sol}$.
The most plausible interpretation is maser amplification of background
synchrotron emission through population-inverted molecular gas in a torus.
However, the angular resolution of $0^{\prime\prime}.68\times
0^{\prime\prime}.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 ($T_{\rm B}$), which is a criterion to classify the emission
mechanism by comparing with the upper energy level of $E_{u}/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^{\circ}$ (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, 2019; Kameno et al.,
2020, 2023a), HI 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, 2003;
Vermeulen et al., 2003; Kadler et al., 2004).

H2O megamaser at 22 GHz is a firm probe for sub-pc-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, 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, 2021) and the first
detection in a radio galaxy NGC 1052 (Kameno et al., 2023b). With a higher
angular resolution by $14\times$ 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 $\sim 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($H${}_{2})=8\times 10^{8}$ and $4\times
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 $V_{\rm sys,LSR}=1492$ km s-1, the luminosity
distance of $D_{L}=17.6$ Mpc, the angular distance of $D_{A}=17.5$ Mpc, and
the linear scale of $85$ pc arcsec-1 (Kameno et al., 2020).

##  2 Methods

###  2.1 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 2.1. Integration times on the target
(NGC 1052) and on a bandpass calibrator 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^{\prime\prime}.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 $J_{Ka,Kc}=10_{2,9}-9_{3,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.

Observation logs and performances.  ExecBlock UID Date $N_{\rm ant}$ Bandpass
Synthesized beam Image rms (mJy beam-1) uid://A002/ (year-month-day) major
$\times$ minor, PA continuum line (1) (2) (3) (4) (5) (6) (7) X109d26e/X12c34
2023-07-10 37**footnotemark: $*$ J0423-0120 $0^{\prime\prime}.017\times
0^{\prime\prime}.011$, $-7^{\circ}.4$ 0.57 2.3 X10a7a20/X65a8 2023-07-24 50
J0238+1636 $0^{\prime\prime}.012\times 0^{\prime\prime}.010$, $41^{\circ}.5$
0.24 2.0 (1) Execution block UID. (2) Date of observation. (3) Number of
antennas. (4) Bandpass calibrator selected by the operation software. (5)
Major axis, minor axis, and position angle of the synthesized beam. (6) Image
rms of continuum with the bandwidth of 5.4 GHz. (7) Image rms with a spectral
resolution of 976.562 kHz.

**footnotemark: $*$

Among 44 antennas, three antennas were lost before the first scan of NGC 1052
due to high wind speed. Additionally four antennas were lost in the last 21
minutes.

Figure 1: Distributions of spatial frequencies, also known as $(u,v)$
coverages, obtained on 2023-07-10 (left) and 2023-07-24 (right).

###  2.2 Data Reduction

Reduction scripts and reduced data are available in the GitHub
repository.111https://github.com/kamenoseiji/ALMA-2022.A.00023.S

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 2.1, as a flux calibrator assuming 2.104 Jy for J$0423-0120$
and 1.304 Jy for J$0238+1636$. This amplitude calibration yields $\sim 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 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.

We performed continuum imaging with Difmap (Shepherd et al., 1994), using
channel-averaged visibilities of BB 1–3. We employed uniform weighting that
yielded the synthesized beams of $0^{\prime\prime}.017\times
0^{\prime\prime}.011$ and $0^{\prime\prime}.012\times 0^{\prime\prime}.010$,
respectively, as stated in table 2.1.

After subtracting a point CLEAN component at the center position, significant
emission remained as an elongated structure spanning $\sim
0^{\prime\prime}.001$ in PA$\sim 70^{\circ}$. The structure was decomposed by
CLEAN components aligned in the same PA. While self calibration with a single
central component yields the agreement factors222The agreement factor in
Difmap self calibration is equivalent to $\sqrt{\chi^{2}_{N}}$, where
$\chi^{2}_{N}$ does not take account of the number of degree of freedom. of
$\sigma=3.94$ and $4.93$ for 2023-07-10 and 2023-07-24, respectively, those
with the aligned CLEAN components result in $\sigma=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.

Figure 2: Continuum-subtracted spectra of H2O $J_{Ka,Kc}=10_{2,9}-9_{3,6}$
emission in NGC 1052. The bottom red spectrum was taken from Kameno et al.
(2023b) observed on 2022-05-12 with a coarser spectral resolution of 7.8125
MHz. The middle green and top blue spectra stand for
uid://A002/X109d26e/X12c34 on 2023-07-10 and uid://A002/X10a7a20/X65a8 on
2023-07-24, offset by 0.075 Jy and 0.1 Jy, respectively. The dotted vertical
line indicates the systemic velocity of 1492 km s-1. The gray strip around
1150 km s-1 masks the telluric ozone feature at 319.9 GHz.

##  3 Results

Table 2.1 summarizes imaging performance of the observations. Image rms is
estimated by statistics in emission-free 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.

###  3.1 Line profiles

Continuum-subtracted spectra are shown in figure 2. The weighted mean
velocities were 1502.8 km s-1 and 1507.5 km s-1 on 2023-07-10 and 2023-07-24,
respectively. The integrated flux densities of $10.79\pm 0.33$ Jy km s-1 and
$14.32\pm 0.43$ Jy km s-1 correspond to the isotropic luminosities of $1120\pm
34$ L${}_{\Sol}$ and $1486\pm 45$ L${}_{\Sol}$, respectively.

The profile mainly consisted of two velocity components — 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 single-peaked 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 3.1 and shown in figure 3. We chose Gaussian or
Lorentzian profile for each component to minimize the residuals. Standard
errors in table 3.1 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 3.1 do not
claim the unique identity between two epochs. Residuals of the fittings
results are characterized by $\chi^{2}$ over the degree of freedoms (d.o.f.)
as shown in figure 3.

Spectral decomposition.  Component 2023-07-10 2023-07-24 Type**footnotemark:
$*$ $V_{\rm center}$ Peak FWHM $V_{\rm center}$ Peak FWHM (km s-1) (mJy) (km
s-1) (km s-1) (mJy) (km s-1) B1 $1374.2\pm 0.4$ $41.7\pm 0.7$ $50.9\pm 1.3$
$1379.1\pm 0.3$ $56.2\pm 1.1$ $52.1\pm 1.2$ L R1 $1530.1\pm 2.0$ $42.3\pm 1.1$
$113.9\pm 3.7$ $1538.7\pm 4.8$ $26.9\pm 4.8$ $178.4\pm 15.8$ G R2 $1485.6\pm
1.0$ $7.1\pm 1.4$ $10.9\pm 2.6$ $1490.7\pm 1.7$ $7.1\pm 1.6$ $25.0\pm 5.3$ G
R3 $1525.3\pm 0.7$ $11.4\pm 1.4$ $11.9\pm 1.9$ $1547.1\pm 3.0$ $34.2\pm 4.8$
$68.1\pm 8.0$ G R4 $1552.1\pm 0.6$ $9.8\pm 1.6$ $8.4\pm 1.7$ $1567.0\pm 0.6$
$6.6\pm 1.7$ $4.6\pm 1.5$ G R5 $1595.6\pm 0.6$ $8.3\pm 1.8$ $6.4\pm 1.8$
$1603.1\pm 1.6$ $46.7\pm 3.8$ $50.0\pm 2.4$ G R6 $1607.8\pm 1.7$ $27.1\pm 2.1$
$61.3\pm 3.1$ $1610.4\pm 0.5$ $7.2\pm 2.2$ $3.2\pm 1.7$ L **footnotemark: $*$G
and L stand for Gaussian and Lorentzian functions, respectively.

Figure 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
3.1. The residuals are shown offset by $-0.01$ Jy.

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\pm 4.2$% and the center velocity shifted by
$4.85\pm 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\pm 0.23$ Jy km s-1 to $10.12\pm 0.30$ Jy km s-1 by $31.6\pm 4.2$%.
Weighted mean velocity of the red components shifted from $1552.7\pm 1.1$ km
s-1 to $1559.0\pm 0.8$ km s-1 by $6.3\pm 1.4$ km s-1.

###  3.2 Size of emission regions

Figure 4: Visibility amplitudes as a function of spatial frequency (projected
baseline length / wavelength). The vertical axis is presented in logarithmic
scale. Black, blue, and red markers stand for the continuum and H2O line
emission in velocity ranges of 1360 – 1500 km s-1 and 1500 – 1625 km s-1,
respectively. The confidence bands indicate the best-fits and 99% confidence
levels.

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.2 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 $T_{\rm B}>10^{6}$ K.

Estimated Gaussian models and brightness temperatures of the continuum and
line-emitting components  Component 2023-07-10 2023-07-24 Flux density Size
(best) Size (max) $T_{\rm B}$ Flux density Size (best) Size (max) $T_{\rm B}$
(Jy) (mas) (mas) (K) (Jy) (mas) (mas) (K) Continuum 0.624 0.51 0.65
$>1.9\times 10^{7}$ 0.704 0.37 0.38 $>6.2\times 10^{7}$ Blue 0.036 0.63 1.17
$>3.1\times 10^{5}$ 0.052 0.35 0.60 $>1.7\times 10^{6}$ Red 0.045 0.30 0.63
$>1.3\times 10^{6}$ 0.061 0.36 0.55 $>2.4\times 10^{6}$ A 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 figure 4.

###  3.3 Continuum images

The continuum image on 2023-07-24 is shown in figure 5. The image on
2023-07-10 is presented in figure 9 in Appendix. The flux densities, estimated
from visibility amplitudes (figure 4 and table 3.2), increased from $0.624\pm
0.019$ Jy to $0.704\pm 0.021$ Jy by $12.8\pm 4.2$%. They were dominated by the
unresolved core component with the peak intensities of 0.581 Jy and 0.649 Jy
in the maps (figures 5 and 9). The extended structure (‘jets’) in east
(PA$=68^{\circ}.2$) and west (PA$=-111^{\circ}.5$) spans $\sim
0^{\prime\prime}.006$ to the east (PA$=68^{\circ}.2$) and west
(PA$=-111^{\circ}.5$). The extension is comparable to the size of the
synthesized beam and is distinct in the super-resolution images. The
orientations were consistent with $61^{\circ}$ to $72^{\circ}.8$ and
$-111^{\circ}.9$ to $-120^{\circ}.9$ for east and west jets in VLBI images
(Kameno et al., 2001, 2003; Vermeulen et al., 2003; Kadler et al., 2004;
Baczko et al., 2019, 2022).

Figure 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 $\pm 3\sigma_{\rm
cont}$ multiplied by powers of 2 where $\sigma_{\rm 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.

###  3.4 Positioning line emitting components

We performed uvmodelfit task in CASA to estimate the positions of line
emission source using continuum-subtracted visibilities binned by 5 spectral
channels (4.57 km s-1). We set an unresolved point source model with a fixed
flux density. Threshold for the maser flux density was set to $3\times$ the
line image rms listed in table 2.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 $\sim 1$ mas. We applied linear regression of $V_{\rm LSR}\sim jJ+dD$,
where $J$ and $D$ are offset along PA=$68^{\circ}.35$ and $-21^{\circ}.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 3.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.

Figure 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.

Figure 7: Position–velocity diagrams along the jet (left, $J$-axis in Table
3.4) and perpendicular to the jet (right, $D$-axis in Table 3.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.

Linear regression for velocity gradient  Epoch 2023-07-10 2023-07-24 gradient
$P$ value gradient $P$ value Axis (km s-1 mas-1) (km s-1 mas-1) $J$ $j=-94\pm
19$ $5.2\times 10^{-6}$ $j=-173\pm 19$ $1.0\times 10^{-13}$ $D$ $d=16\pm 15$
$0.29$ $d=-18\pm 36$ $0.602$

##  4 Discussion

###  4.1 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 $\Delta T=14$ days constrains the size of the emitter
$<c\Delta 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.

###  4.2 Velocity drift

B1 maser component showed the velocity drift of $127\pm 13$ km s-1 yr-1 or
acceleration of $g=4.0\times 10^{-3}$ m s-2. The drift rate is an order of
magnitude greater than that of well-known 22-GHz H2O maser in NGC 4258; 7.5 km
s-1 yr-1 (Haschick et al., 1994), $9.5\pm 1.1$ km s-1 yr-1 (Greenhill et al.,
1995), or $9.6\pm 1.0$ km s-1 yr-1 (Nakai et al., 1995).

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\pm 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\pm 37$ km s-1 yr-1 or $5.2\times
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\pm 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 $M_{\rm
BH}=1.5\times 10^{8}$ M${}_{\Sol}$ (Woo & Urry, 2002), we have the distance of
the maser from the SMBH, $r$, as $r=\sqrt{GM_{\rm BH}/g}=0.072$ pc or $5000\
R_{\rm s}$ where $R_{\rm s}=2GM_{\rm BH}/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\pm 13$ km s-1. Indeed, the maser velocity
does not exceed the escape velocity, $v_{\rm esc}=\sqrt{2GM_{\rm
BH}/r}=4.2\times 10^{3}$ km s-1.

Compared with the inner boundary of $0.12$ pc or $3.3\times 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}\propto r^{-1/2}$, as 2000 K : 800 K $\sim\sqrt{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.3 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 $\sim 0.1$ pc. The direction of the velocity gradient is consistent
with the 100-pc-scale bipolar outflow traced by [OIII] (Sugai et al., 2005;
Cazzoli et al., 2022) and by [ArII] and [NeIII] (Goold et al., 2023). As
argued in Kameno 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\pm 43$ K implying excitation by jet-torus
interaction (Kameno et al., 2023a). The observed velocity gradient supports
this scenario.

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
$\nu^{-2.1}$, where $\nu$ is the frequency, the opacity of FFA at 1 GHz,
$\tau_{0}\simeq 1000$, yields $\tau\simeq 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\times 10^{-3}$ m s-2 is ascribed to gravitational acceleration
and the derived distance is the radius of a rotating disk, we expect the
rotating speed of $V_{\rm rot}=c\sqrt{R_{\rm s}/2r}=3000$ km s-1 and the
velocity gradient of $V_{\rm rot}/r=4.2\times 10^{4}$ km s-1 pc-1 or
$3.5\times 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 $\Delta V\sim 400$ km
s-1 would be $r\Delta V/V_{\rm rot}=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\times 10^{3}$ $R_{s}$ inside a distance
of $10^{4}$ Rs from the central engine. The estimated maser emission section
width of $\sim 700$ $R_{s}$ 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 $\sim 0.01$ mas is desired.

A rotating disk with $V_{\rm rot}=3000$ km s-1 would yield high-velocity
components at large impact parameters. Expected frequency offsets of $\pm 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.

###  4.4 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 $\rho_{j}v^{2}_{j}=1.5\rho_{t}v^{2}_{s}$,
where $\rho_{j}$ and $v_{j}$ are density and speed of the jet, and $\rho_{t}$
and $v_{s}$ are the density of the torus and the advancing speed of the shock,
respectively. Taking $v_{j}=0.26c$ (Vermeulen et al., 2003),
$\rho_{t}=1.1\times 10^{-12}$ kg m-3 estimated from the maser excitation
condition of $n_{\rm H2}\gtrsim 3.3\times 10^{8}$ cm-3 (Gray et al., 2016), we
have $\rho_{j}=3.5\times 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\pm 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.

Figure 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^{\circ}$. (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.

###  4.5 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\pm
4.2$% and $31.6\pm 4.2$% in 14 days while the continuum flux density increased
by $12.8\pm 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 $\sim 20$% increase of the maser amplification gain
variation with the timescale of $1.2\times 10^{6}$ s$/\log 1.2=6.6\times
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 $\leq$ 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 population-inverted
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 $\lambda_{J}=c_{s}/\sqrt{G\rho}\sim 5.3\times 10^{14}$
m or 0.017 pc, where $\rho=1.1\times 10^{-12}$ kg m-3 is the density for maser
excitation condition of $n_{\rm H2}>3.3\times 10^{8}$ cm-3 (Gray et al., 2016)
and $c_{s}$ is the sound speed for the kinematic temperature of 2000 K.
Crossing time of $\lambda_{J}/V_{\rm rot}=1.8\times 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\pm 4$ K and $344\pm 43$ K, respectively. Smaller Jeans length for lower-
temperature 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.

##  5 Summary

ALMA follow-up observations with the widest array configuration led findings
summarized below.

  1. 1.

The 321-GHz H2O emission in NGC 1052 is definitely a maser with the brightness
temperature $T_{\rm B}>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. 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\pm 4.2$% and $31.6\pm 4.2$% in 14 days
while the continuum flux density increased by $12.8\pm 4.2$%. The variability
in amplification gain implies inhomogeneity of the population-inversion zone
and/or change of excitation condition.

  3. 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
sub-components, the velocity drift of $127\pm 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 $M_{\rm BH}=1.5\times 10^{8}$
M${}_{\Sol}$, the distance of the maser will be 5000 $R_{s}$ 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 $V_{\rm rot}=3000$ km s-1.

  4. 4.

We identified velocity gradient of $-174\pm 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 $\sim 0.16$ mas. To resolve expected velocity gradient
for a rotating disk, a spatial resolution sharper than 0.1 mas is required.

  5. 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 $\sim\pm 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.

## Acknowledgments

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2021.1.00341.S
and ADS/JAO.ALMA#2023.A.00023.S. ALMA is a partnership of ESO (representing
its member states), NSF (USA) and NINS (Japan), together with NRC (Canada),
MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the
Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and
NAOJ. This work is supported by JSPS KAKENHI 18K03712 and 21H01137.

## Results on 2023-07-10

The main text addresses results on 2023-07-24 for discussion about the maser
distribution and the velocity gradient because the image quality is
significantly better than that on 2023-07-10. They are presented in figures 9
and 10.

Figure 9: Continuum and maser map (left) and close-up of maser components
(right) on 2023-07-10, as well as figures 5 and 6. The contour levels are $\pm
3\sigma_{\rm cont}$ multiplied by powers of 2 where $\sigma_{\rm cont}=0.57$
mJy beam-1.

Figure 10: Position–velocity diagrams along and perpendicular to the jet on
2023-07-10, as well as figure 7.

## References

  * Baczko et al. (2019) Baczko, A.-K., Schulz, R., Kadler, M., et al. 2019, A&A, 623, A27
  * Baczko et al. (2022) Baczko, A.-K., Ros, E., Kadler, M., et al. 2022, A&A, 658, A119
  * Braatz et al. (1994) Braatz, J. A., Wilson, A. S., & Henkel, C. 1994, ApJ, 437, L99
  * Braatz et al. (2003) Braatz, J. A., Wilson, A. S., Henkel, C., et al. 2003, ApJS, 146, 249
  * CASA Team et al. (2022) CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 114501
  * Cazzoli et al. (2022) Cazzoli, S., Hermosa Muñoz, L., Márquez, I., et al. 2022, A&A, 664, A135
  * Claussen et al. (1998) Claussen, M. J., Diamond, P. J., Braatz, J. A., et al. 1998, ApJ, 500, L129
  * Fromm et al. (2019) Fromm, C. M., Younsi, Z., Baczko, A., et al. 2019, A&A, 629, A4
  * Gallimore & Impellizzeri (2023) Gallimore, J. F. & Impellizzeri, C. M. V. 2023, ApJ, 951, 109
  * Goold et al. (2023) Goold, K., Seth, A., Molina, M., et al. 2023, arXiv:2307.01252
  * Gray et al. (2016) Gray, M. D., Baudry, A., Richards, A. M. S., et al. 2016, MNRAS, 456, 374
  * Greenhill et al. (1995) Greenhill, L. J., Henkel, C., Becker, R., et al. 1995, A&A, 304, 21
  * Hagiwara et al. (2013) Hagiwara, Y., Miyoshi, M., Doi, A., et al. 2013, ApJ, 768, L38
  * Hagiwara et al. (2016) Hagiwara, Y., Horiuchi, S., Doi, A., et al. 2016, ApJ, 827, 69
  * Hagiwara et al. (2021) Hagiwara, Y., Horiuchi, S., Imanishi, M., et al. 2021, ApJ, 923, 251
  * Haschick et al. (1994) Haschick, A. D., Baan, W. A., & Peng, E. W. 1994, ApJ, 437, L35
  * Herrnstein et al. (2005) Herrnstein, J. R., Moran, J. M., Greenhill, L. J., et al. 2005, ApJ, 629, 719
  * Humphreys et al. (2005) Humphreys, E. M. L., Greenhill, L. J., Reid, M. J., et al. 2005, ApJ, 634, L133
  * Humphreys et al. (2008) Humphreys, E. M. L., Reid, M. J., Greenhill, L. J., et al. 2008, ApJ, 672, 800
  * Humphreys et al. (2013) Humphreys, E. M. L., Reid, M. J., Moran, J. M., et al. 2013, ApJ, 775, 13
  * Impellizzeri (2022) Impellizzeri, C. M. V. 2022, Nature Astronomy, 6, 885
  * Impellizzeri et al. (2008) Impellizzeri, V., Roy, A. L., & Henkel, C. 2008, The role of VLBI in the Golden Age for Radio Astronomy, 9, 33
  * Kadler et al. (2004) Kadler, M., Ros, E., Lobanov, A. P., et al. 2004, A&A, 426, 481
  * Kameno et al. (2001) Kameno, S., Sawada-Satoh, S., Inoue, M., et al. 2001, PASJ, 53, 169
  * Kameno et al. (2003) Kameno, S., Inoue, M., Wajima, K., et al. 2003, PASA, 20, 134
  * Kameno et al. (2005) Kameno, S., Nakai, N., Sawada-Satoh, S., et al. 2005, ApJ, 620, 145
  * Kameno et al. (2020) Kameno, S., Sawada-Satoh, S., Impellizzeri, C. M. V., et al. 2020, ApJ, 895, 73
  * Kameno et al. (2023a) Kameno, S., Sawada-Satoh, S., Impellizzeri, C. M. V., et al. 2023a, ApJ, 944, 156
  * Kameno et al. (2023b) Kameno, S., Harikane, Y., Sawada-Satoh, S., et al. 2023b, PASJ, 75, L1
  * Koekemoer et al. (1995) Koekemoer, A. M., Henkel, C., Greenhill, L. J., et al. 1995, Nature, 378, 697
  * Liszt & Lucas (2004) Liszt, H. & Lucas, R. 2004, A&A, 428, 445
  * Miyoshi et al. (1995) Miyoshi, M., Moran, J., Herrnstein, J., et al. 1995, Nature, 373, 127
  * Modjaz et al. (2005) Modjaz, M., Moran, J. M., Kondratko, P. T., et al. 2005, ApJ, 626, 104
  * Moran et al. (1995) Moran, J., Greenhill, L., Herrnstein, J., et al. 1995, Proceedings of the National Academy of Science, 92, 11427
  * Nakai et al. (1993) Nakai, N., Inoue, M., & Miyoshi, M. 1993, Nature, 361, 45
  * Nakai et al. (1995) Nakai, N., Inoue, M., Miyazawa, K., et al. 1995, PASJ, 47, 771
  * Nakahara et al. (2020) Nakahara, S., Doi, A., Murata, Y., et al. 2020, AJ, 159, 14
  * Neufeld & Melnick (1991) Neufeld, D. A. & Melnick, G. J. 1991, ApJ, 368, 215
  * Omar et al. (2002) Omar, A., Anantharamaiah, K. R., Rupen, M., et al. 2002, A&A, 381, L29
  * Pesce et al. (2016) Pesce, D. W., Braatz, J. A., & Impellizzeri, C. M. V. 2016, ApJ, 827, 68
  * R Core Team (2023) R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/
  * Sawada-Satoh et al. (2008) Sawada-Satoh, S., Kameno, S., Nakamura, K., et al. 2008, ApJ, 680, 191
  * Sawada-Satoh et al. (2016) Sawada-Satoh, S., Roh, D.-G., Oh, S.-J., et al. 2016, ApJ, 830, L3
  * Sawada-Satoh et al. (2019) Sawada-Satoh, S., Byun, D.-Y., Lee, S.-S., et al. 2019, ApJ, 872
  * Shepherd et al. (1994) Shepherd, M. C., Pearson, T. J., & Taylor, G. B. 1994, BAAS
  * Shostak et al. (1983) Shostak, G. S., van Gorkom, J. H., Ekers, R. D., et al. 1983, A&A, 119, L3
  * Sugai et al. (2005) Sugai, H., Hattori, T., Kawai, A., et al. 2005, ApJ, 629, 131
  * Vermeulen et al. (2003) Vermeulen, R. C., Ros, E., Kellermann, K. I., et al. 2003, A&A, 401, 113
  * Watson & Wallin (1994) Watson, W. D. & Wallin, B. K. 1994, ApJ, 432, L35
  * Woo & Urry (2002) Woo, J.-H. & Urry, C. M. 2002, ApJ, 579, 530
  * Yamaki et al. (2012) Yamaki, H., Kameno, S., Beppu, H., et al. 2012, PASJ, 64, 118
  * Yates et al. (1997) Yates, J. A., Field, D., & Gray, M. D. 1997, MNRAS, 285, 303


# On the Interpretability of Quantum Neural Networks

Lirandë Pira<EMAIL_ADDRESS>University of Technology Sydney,
Centre for Quantum Software and Information, Ultimo NSW 2007, Australia  Chris
Ferrie  University of Technology Sydney, Centre for Quantum Software and
Information, Ultimo NSW 2007, Australia

###### Abstract

Interpretability of artificial intelligence (AI) methods, particularly deep
neural networks, is of great interest. This heightened focus stems from the
widespread use of AI-backed systems. These systems, often relying on intricate
neural architectures, can exhibit behavior that is challenging to explain and
comprehend. The interpretability of such models is a crucial component of
building trusted systems. Many methods exist to approach this problem, but
they do not apply straightforwardly to the quantum setting. Here, we explore
the interpretability of quantum neural networks using local model-agnostic
interpretability measures commonly utilized for classical neural networks.
Following this analysis, we generalize a classical technique called LIME,
introducing Q-LIME, which produces explanations of quantum neural networks. A
feature of our explanations is the delineation of the region in which data
samples have been given a random label, likely subjects of inherently random
quantum measurements. We view this as a step toward understanding how to build
responsible and accountable quantum AI models.

##  I Introduction

Artificial intelligence (AI) has become ubiquitous. Often manifested in
machine learning algorithms, AI systems promise to be evermore present in
everyday high-stakes tasks [1, 2]. This is why building fair, responsible, and
ethical systems is crucial to the design process of AI algorithms. Central to
the topic of trusting AI-generated results is the notion of interpretability,
also known as explainability. The interpretability of AI models is a pivotal
concern in contemporary AI research, particularly with the ubiquity of deep
neural networks. This has given rise to research topics under the umbrella of
interpretable machine learning (or IML) and explainable AI (or XAI), noting
that the terms interpretable and explainable are used synonymously throughout
the corresponding literature. Generically, interpretability is understood as
the extent to which humans comprehend the output of an AI model that leads to
decision-making [3]. Humans strive to understand the “thought process” behind
the decisions of the AI model — otherwise, the system is referred to as a
“black box”. Even though the terms “interpretation” and “explanation” are used
colloquially in the literature with varying degrees, throughout this paper, we
use them as follows. Interpretable entails a model that humans can understand
and comprehend through direct observation of its internal workings or outputs.
It implies that the model is intuitive to the human observer and that no
further tools are required. On the other hand, explanation refers to the
output of a tool used to articulate the behavior of a model. In our case, and
indeed many others, an explanation is a model itself, which is justifiably
called such when it is interpretable. Complex models and their simpler
explanations are often called black-box and white-box models, respectively.

Figure 1: Categorization of interpretability techniques as they apply to
classical and quantum resources. Here, the well-known QML diagram represents
data and an algorithm or device, which can be classical (C) or quantum (Q) in
four different scenarios. We consider a reformulation of interpretable
techniques to be required in the CQ scenario. In the QC and QQ quadrants, the
design of explicitly quantum interpretable methods may be required. The scope
of this paper covers CQ approaches.

The precise definition of a model’s interpretability has been the subject of
much debate [4, 5]. Naturally, there exist learning models which are more
interpretable than others, such as simple decision trees. On the other hand,
the models we prefer best for solving complex tasks, such as deep neural
networks (DNNs), happen to be highly non-interpretable, which is due to their
inherent non-linear layered architecture [6]. We note that DNNs are one of the
most widely used techniques in machine learning. Thus, the interpretability of
neural networks is an essential topic within the field of interpretable ML
research [7, 8]. In this work, we focus on precisely the topic of
interpretability as we consider the quantum side of neural networks.

In parallel, recent years have witnessed a surge of research efforts in
quantum machine learning (QML) [9, 10]. This research area sits at the
intersection of machine learning and quantum computing. The development of QML
has undergone different stages. Initially, the field started with the quest
for speedups or quantum advantages. More recently, the target has morphed into
further pursuits in expressivity and generalization power of quantum models.
Nowadays, rather than “competing” with classical models, quantum models are
further being enhanced on their own, which could, in turn, improve classical
machine learning techniques. One of the key techniques currently used in QML
research is the variational quantum algorithm, which acts as the quantum
equivalent to classical neural networks [11]. To clarify the analogy, we will
refer to such models as quantum neural networks (QNNs) [12].

Given the close conceptual correspondence to classical neural networks, it is
natural to analyze their interpretability, which is important for several
reasons. Firstly, QNNs may complement classical AI algorithm design, making
their interpretability at least as important as classical DNNs. Secondly, the
quantum paradigms embedded into QNNs deserve to be understood and explained in
their own right. The unique non-intuitive characteristics of quantum nature
can make QNNs more complicated to interpret from the point of view of human
understandability. Finally, with the growing interest and capabilities of
quantum technologies, it is crucial to identify and mitigate potential sources
of errors that plague conventional AI due to a lack of transparency.

In this work, we define and establish key notions pertaining to the
interpretability of quantum neural networks. Through this formalization, we
aim to contribute a robust foundation for advancing the discourse on the
interpretability of quantum neural networks within the broader context of
quantum machine learning research. In doing so, we generalize some well-known
interpretable techniques to the quantum domain. Consider the standard
relationship diagram in QML between data and algorithm (or device) type, where
either can be classical (C) or quantum (Q). This entails the following
combinations (CC, CQ, QC, and QQ), shown in Fig. 1. Classical interpretable
techniques are the apparent domain of CC. We will discuss, but not dwell on,
the potential need for entirely new techniques when the data is quantum (QC
and QQ). In CQ, the domain that covers the so-called quantum-enhanced machine
learning techniques, although the data is classical, the output of the quantum
devices is irreversibly probabilistic. CQ can be seen as a specialization of
the quantum data techniques above. Generalizing classical notions of
interpretability to this domain is the subject of our work.

The question of interpretability in quantum machine learning models more
broadly, as well as of QNNs more specifically, has already started to receive
attention [13, 14, 15], particularly involving the concept of Shapley values
[16], which attempt to quantify the importance of features in making
predictions. In [13], interpretability is explored using Shapley values for
quantum models by quantifying the importance of each gate to a given
prediction. The complexity of computing Shapley values for generalized quantum
scenarios is analyzed in [17]. In [15], Shapley values are computed for a
particular class of QNNs. Our work complements these efforts using an
alternative notion of explainability to be discussed in detail next.

##  II Interpretability in AI

###  II.1 Taxonomy

There are several layers to the design of interpretability techniques. To
start, they can be model-specific or model-agnostic. As the name suggests,
model-specific methods are more restrictive in terms of which models they can
be used to explain. As such, they are designed to explain one single type of
model. In contrast, model-agnostic methods allow for more flexibility in usage
as they can be used on a wide range of model types. At large, model-agnostic
methods can have a global or local interpretability dimension. Locality
determines the scope of explanations with respect to the model.
Interpretability at a global level explains the average predictions of the
model as a whole. At the same time, local interpretability gives explanations
at the level of each sample. In another axis, these techniques can be active
(inherently interpretable) or passive (post-hoc). The state of these
interpretable paradigms implies the level of involvement of interpretable
techniques in the outcome of the other parameters. Active techniques change
the structure of the model itself by leaning towards making it more
interpretable. In contrast, passive methods explain the model outcome once the
training has finished. In comparison to model-agnostic methods, which work
with samples at large, there also exist example-based explanations that
explain selected data samples from a dataset. An example of this method is the
$k$-nearest neighbors models, which average the outcome of $k$ nearest
selected points.

Other than the idea of building interpretable techniques, or more precisely,
techniques that interpret various models, there exist models that are
inherently interpretable. Such models include linear regression, logistic
regression, naive Bayes classifiers, decision trees, and more. This feature
makes them good candidates as surrogate models for interpretability. Based on
this paradigm, the concept of surrogate models exists, which uses
interpretable techniques as a building block for designing other interpretable
methods. Such important techniques are, for example, local interpretable
model-agnostic explanations (LIME) [18] and Shapley additive explanations
(known as SHAP) [16].

###  II.2 Interpretability of neural networks

The interpretability of neural networks remains a challenge on its own. This
tends to amplify in complex models with many layers and parameters.
Nevertheless, there is active research in the field and several proposed
interpretable techniques [7]. Such techniques that aim to gain insights into
the decision-making of a neural network include saliency maps [19], feature
visualization [20, 21], perturbation or occlusion-based methods [18, 16], and
layerwise relevance propagation (also known by its acronym LRP) [22].

To expand further on the above-mentioned techniques, saliency maps use
backpropagation and gradient information to identify the most influential
regions contributing to the output result. This technique is also called pixel
attribution [4]. Feature visualization, particularly useful for convolutional
natural networks, is a technique that analyses the importance of particular
features in a dataset by visualizing the patterns that activate the output. In
the same remark, in terms of network visualizations, Ref. [23] goes deeper
into the layers of a convolutional neural network to gain an understanding of
the features. This result, in particular, shows the intricacies and the rather
intuitive process involved in the decision-making procedure of a network as it
goes through deeper layers. Occlusion-based methods aim to perturb or
manipulate certain parts of the data samples to observe how the explanations
change. These methods are important in highlighting deeper issues in neural
networks. Similarly, layerwise relevance propagation techniques re-assign
importance weight to the input data by analyzing the output. This helps the
understanding by providing a hypothesis over the output decision. Finally, the
class of surrogate-based methods mentioned above is certainly applicable in
neural networks as well.

The importance of these techniques is also beyond the interpretability
measures for human understanding. They can also be seen as methods of
debugging and thus improving the result of a neural network as in Ref. [23].
Below, we take a closer look at surrogate model-agnostic local interpretable
techniques, which are applicable to DNNs as well.

###  II.3 Local interpretable methods

Local interpretable methods tend to focus on individual data samples of
interest. One of these methods relies on explaining a black-box model using
inherently interpretable models, also known as surrogate methods. These
methods act as a bridge between the two model types. The prototype of these
techniques is the so-called local interpretable model-agnostic explanations
(LIME), which has received much attention since its invention in 2016 [18].
Local surrogate methods work by training an interpretable surrogate model that
approximates the result of the black-box model to be explained. LIME, for
instance, is categorized as a perturbation-based technique that perturbs the
input dataset. Locality in LIME refers to the fact that the surrogate model is
trained on the data point of interest, as opposed to the whole dataset (which
would be the idea behind global surrogate methods). Eq. (1) represents the
explanation $\xi$ of a sample $x$ via its two main terms, namely the term
$L(f,g,\pi_{x})$ representing the loss, which is the variable to be minimized,
and $\Omega(g)$ which is the complexity measure, which encodes the degree of
interpretability. Here $f$ is the black-box model, $g$ is the surrogate model,
and $\pi_{x}$ defines the region in data space local to $x$ (See Alg. II.3
pseudocode). In broader terms, LIME is a trade-off between interpretability
and accuracy,

$\xi(x)=\operatorname*{argmin}_{g\in G}L(f,g,\pi_{x})+\Omega(g)$ (1)

In the following, we make use of the concept of local surrogacy to understand
the interpretability of quantum models using LIME as a starting point. Much
like LIME, we develop a framework to provide explanations of black-box models
in the quantum domain. The class of surrogate models, the locality measure,
and the complexity measure are free parameters that must be specified and
justified in each application of the framework.

LIME Algorithm [18]

1:function LIME($f,x,D,K$)

2: for $i=1$ to $K$ do

3: Sample $z_{i}$ from $D$

4: $g_{i}=f(z_{i})$ $\triangleright$ $g_{i}$ is the prediction of the
classifier on the perturbed sample

5: end for

6: $\xi(x)=\operatorname*{argmin}_{g\in G}L(f,g,\pi_{x})$ $\triangleright$ $L$
is a loss function, $\pi_{x}$ is a locality measure to the data point, and $G$
is the set of surrogate models

7: return $\xi(x)$

8:end function

9:function Loss($f,g,\pi_{x}$)

10: $L(f,g,\pi_{x})=\sum_{z_{i}\in D}\pi_{x}(z_{i})[f(z_{i})\neq
g(z_{i})]+\Omega(g)$ $\triangleright$ Penalize the difference between
predictions

11: return $L(f,g,\pi_{x})$

12:end function

13:function Locality($\pi_{x},z,x$)

14: $\pi_{x}(z)=\exp\left(-\frac{d(x,z)^{2}}{\sigma^{2}}\right)$
$\triangleright$ $d$ is a distance metric, $\sigma$ is a bandwidth parameter

15: return $\pi_{x}(z)$

16:end function

##  III The case for quantum AI interpretability

As mentioned in Section I, interpretability in the quantum literature in the
context of machine learning can take different directions. We consider the
case when data is classical and encoded into a quantum state, which is
manipulated by a variational quantum circuit before outputting a classical
decision via quantum measurement. Our focus is on interpreting the classical
output of the quantum model.

A quantum machine learning model $f$ takes as input data $x$ and first
produces quantum data $|\psi(x)\rangle$. A trained quantum algorithm, such as
the QNN, subsequently processes the quantum data and generates a
classification decision based on the outcome of a quantum measurement. This is
not conceptually different from a classical neural network beyond the fact
that the weights and biases have been replaced by parameters of quantum
processes, except for one crucial difference — quantum measurements are
unavoidably probabilistic.

Probabilities, or quantities interpreted as such, often arise in conventional
neural networks. However, these numbers are encoded in bits and directly
accessible, so they are typically used to generate deterministic results
(through thresholding, for example). Qubits, on the other hand, are not
directly accessible. While procedures exist to reconstruct qubits from
repeated measurements (generally called tomography), these are inefficient —
defeating any purpose of encoding information into qubits in the first place.
Hence, QML uniquely forces us to deal with uncertainty in interpreting its
decisions.

In the case of probabilistic decisions, the notion of a decision boundary is
undefined. A reasonable alternative might be to define the boundary as those
locations in data space where the classification is purely random (probability
$\frac{1}{2}$). A data point here is randomly assigned a label. For such a
point, any explanation for its label in a particular realization of the random
decision process would be arbitrary and prone to error. It would be more
accurate to admit that the algorithm is indecisive at such a point. This
rationale is equally valid for data points near such locations. Thus, we
define the region of indecision as follows,

$R=\left\\{x^{\prime}\in X:\left\lvert
P(f(x^{\prime})=1)-\frac{1}{2}\right\rvert<\epsilon\right\\},$ (2)

where $\epsilon$ is a small positive constant representing a threshold of
uncertainty tolerated in the classification decision. In a certain context,
points situated within this rationale lack a discernible rationale, hence
explanation, for their specific outcomes, relying instead on random
circumstances.

At present, although certain data points possess labels that were randomly
assigned, one may pose the question as to the underlying rationale for such
assignments. In other words, even data points lying within the region of
indecision demand an explanation. Next, we will show how the ideas of local
interpretability can be extended to apply to the probabilistic or quantum
settings. In doing so, we aim to provide a nuanced and comprehensive
understanding of model predictions.

###  III.1 Quantum LIME: Probabilistic local interpretability

Here, we define Quantum LIME (or Q-LIME) as an extension of the classical
counterpart to quantum models, namely quantum local interpretable model-
agnostic explanations. A specialization of the definition of Q-LIME is the
hybrid version in the $CQ$, such as in Fig. 1.

In the context of the LIME algorithms, the loss function is typically chosen
to compare models and their potential surrogates on a per-sample basis.
However, if the model’s output is random, the loss function will also be a
random variable. An obvious strategy would be to define loss via expectation:

$\xi(x)=\operatorname*{argmin}_{g\in G}\mathbb{E}[L(f,g,\pi_{x})]+\Omega(g).$
(3)

Nevertheless, even under these circumstances, a definitive assertion regarding
$\xi$ as an explanation remains elusive. This stems from the recognition that
its predictions are only capturing the average behavior of the underlying
model’s randomness. In fact, the label provided by $\xi$ may be the opposite
of that assigned to $x$ by the model in any particular instance, which clearly
does not capture a reasonable notion of explanation.

To mitigate this, we generalize the classical definition and call an
explanation the distribution $\Xi$ of trained surrogate models $g$ obtained
through the application of LIME to the random output of any probabilistic
classifier. Note again that $g$ is random, trained on synthetic local data
with random labels assigned by the underlying model. Thus, the explanation
inherits any randomness from the underlying model. It is not the case that the
explanation provides an interpretation of the randomness per se — however, we
can utilize the distribution of surrogate models to simplify the region of
indecision, hence providing an interpretation of it.

The symbol $\Xi$ represents the Q-LIME explanation, emphasizing the
interpretability framework tailored for QNNs. Additionally, Alg. III.1 refers
to a Monte Carlo approximation method employed within the Q-LIME process. This
algorithm plays a pivotal role in facilitating the computation of
interpretable insights by utilizing a Monte Carlo sampling approach, thereby
enhancing the efficiency of the interpretation process.

Quantum LIME (Q-LIME)

1:function QLIME($U,x,D,K$)

2: $\Xi\leftarrow$ Empty list to store quantum surrogate models

3: for $i=1$ to $K$ do

4: $D_{i}\leftarrow$ Generate synthetic quantum data locally around $x$

5: $\gamma_{i}\leftarrow$ LIME($U,x,D_{i},K$) $\triangleright$ Apply LIME to
generate a quantum surrogate model (see Alg. 1)

6: Add $\gamma_{i}$ to $\Xi$

7: end for

8: return $\Xi$

9:end function

###  III.2 Local region of indecision

Figure 2: Depiction of the concept of the local region of indecision. The
space within the dashed lines represents the region in the decision space
where data samples exhibit ambiguous classification due to randomness. The
figure showcases a two-class classification task in a two-dimensional space
with two features represented along the horizontal and vertical axes. Here,
$\epsilon$ is the pre-defined threshold. Likely, data samples, either inside
or close to the band, are interpreted as randomly assigned.

In this section, we define the local region of indecision. In a broad sense,
this is the region of indecision interpreted locally through a distribution of
surrogate models. Suppose a particular data point lies within its own local
region of indecision. The explanation for its label is thus the data was given
a random label, even though its explanation may not be a satisfying one.
Moreover, this is a strong statement because — in principle — all possible
interpretable surrogate models have been considered in the optimization.

The local region of indecision can be defined as the region of the input space
where the classification decision of the quantum model is uncertain. More
formally, we can define the local region of indecision $B$ for a data point
$x$ in a dataset $X$ as,

$B=\left\\{x^{\prime}\in X:\left\lvert
P(g(x^{\prime})=1|f,\Xi,\pi_{x})-\frac{1}{2}\right\rvert<\epsilon\right\\},$
(4)

where $\epsilon$ is again a small positive constant representing a threshold
of uncertainty tolerated in the classification decision — this time with
reference to the explanation rather than the underlying model.

The exact relationship between $B$ and $R$ is complicated, depending on the
details of the loss function. However, for a $0-1$ loss function and $\pi_{x}$
defining a closed boundary around $x$, $B=R$ within $\pi_{x}$. In a more
general case, $B$ approximates $R$ up to the freedom that elements of $\Xi$
are allowed to deviate from $f$. Note that the distribution in Eq. (4) is over
$\Xi$ as each $g$ provides deterministic labels. This has been formalized in
Alg. III.2.

In much the same way that an interpretable model approximates decision
boundaries locally in the classical context, a local region of indecision
approximates the proper region of indecision in the quantum (or probabilistic)
context. The size and shape of this region will depend on several factors,
such as the choice of the interpretability technique, the complexity of the
surrogate model, and the number of features in the dataset. We call the region
a “band” as it describes the simplest schematic presented in Fig. 2.

In the examples presented in this paper, we chose the interquartile range
(IQR) as a summary measure for the set of quantum surrogate models $\Xi$.
However, other options are undoubtedly possible. It is important to note that
in higher-dimensional spaces, the complexity of quantum models and the
diversity of decision boundaries may necessitate alternative summary measures.

Range

1:$B\leftarrow$ Empty list to store quantum decision boundaries

2:$\Xi\leftarrow$ QLIME($U,x,D,K$) $\triangleright$ Generate quantum surrogate
models using Quantum LIME (see Alg. 2)

3:for $i=1$ to $M$ do $\triangleright$ Repeat for $M$ iterations

4: $\xi\leftarrow$ Sample from $\Xi$

5: $B_{i}\leftarrow$ Boundary($\xi$) $\triangleright$ Obtain the quantum
decision boundary from the sampled quantum surrogate model

6: Add $B_{i}$ to $B$

7:end for

8:$IQR\leftarrow$ InterquartileRange($B$) $\triangleright$ Compute the
interquartile range of quantum decision boundaries

9:return $IQR$

10:function InterquartileRange($B$)

11: Sort the quantum decision boundaries in $B$

12: $Q_{1}\leftarrow$ The median of the lower half of $B$

13: $Q_{3}\leftarrow$ The median of the upper half of $B$

14: $IQR\leftarrow Q_{3}-Q_{1}$ $\triangleright$ Interquartile range

15: return $IQR$

16:end function

##  IV Numerical experiments for interpreting QNNs

We use the well-known Iris dataset [24] for our numerical experiments. To
enhance the explicability of our methodology, we opt to simplify the analysis
by framing it as a binary classification problem. This entails the exclusive
utilization of two out of the three available classes within the dataset,
coupled with the restriction to employ only two out of the four available
features. In consideration of our indifference toward classifying the two
categories of flowers, we abstract the names of these classes and labels
below.

The trained quantum model to be explained is a hybrid QNN trained using
simultaneous perturbation stochastic approximation (better known by its
acronym SPSA) [25], which is an optimization algorithm that efficiently
estimates the gradient of an objective function. This model has been built and
simulated using the Qiskit framework [26]. Each data point is encoded into a
quantum state with the angle encoding [27], which acts by mapping numerical
values to specific angles, facilitating their integration into quantum models.
The QNN model is an autoencoder with alternating layers of single qubit
rotations and entangling gates [28]. Since our goal here is to illustrate the
local region of indecision, as in Eq. 4, we do not optimize over the
complexity of surrogate models and instead fix our search to the class of
logistic regression models with two features.

The shaded background in each plot of Figs. 3 and 4 show the decision region
of the trained QNN. Upon inspection, it is clear that these decision regions,
and the implied boundary, change with each execution of the QNN. In other
words, the decision boundary is ill-defined. In Fig. 3, we naively apply the
LIME methodology to two data points — one in the ambiguous region and one deep
within the region corresponding to one of the labels. In the latter case, the
output of the QNN is nearly deterministic in the local neighborhood of the
chosen data point, and LIME works as expected — it produces a consistent
explanation.

Nonetheless, in the first example, the data point receives a random label. It
is clearly within the region of indecision for a reasonably small choice of
$\epsilon$. The “explanation” provided by LIME (summarized by its decision
boundary shown as the solid line in Fig. 3) is random. In other words, each
application of LIME will produce a different explanation. For the chosen data
point, the explanation itself produces opposite interpretations roughly half
the time, and the predictions it makes are counter to the actual label
provided by the QNN model to be explained roughly half the time. Clearly, this
is an inappropriate situation to be applying such interpretability techniques.
Heuristically, if a data point lies near the decision boundary of the
surrogate model for QNN, we should not expect that it provides a satisfactory
explanation for its label. Q-LIME and the local region of indecision rectify
this.

Figure 3: The classification uncertainty for a chosen data sample (yellow).
The two shadings represent the decision boundary of the QNN, which is clearly
randomly defined. The synthetic data and the produced surrogate linear
decision boundary are shown for a single instance of QNN labels. (top) A clear
example of an undefined classification result for the selected data point.
(bottom) Decision boundary unambiguously separates the who classes with
respect to the selected data point.

Figure 4: The approximated local region of indecision. (top) An example of a
marked data point that lies on the local region of indecision. (bottom) A data
point that is outside of this region can be assessed for interpretability as
per the interpretable techniques.

For the same sample data points, the local region of indecision is shown in
Fig. 4. Data points within their own band should be regarded as having been
classified due purely to chance — and hence have no “explanation” in the
conventional sense for the label they have been assigned. Note that the band
itself is unlikely to yield an analytic form. Hence, some numerical
approximation is required to calculate it in practice. Our approach,
conceptually, was to repeat what was done to produce Fig. 3 many times and
summarize the statistics to produce Fig. 4. A more detailed description
follows, and the implementation to reproduce the results presented here can be
found at [29].

Approximating the local region of indecision requires a trained QNN, a sample
data point, a function defining “local”, and some meta-parameters defining the
quality of the numerical approximation, such as the amount of synthetic data
to be used and the number of Monte Carlo samples to take. Once these are
defined, the general procedure is as follows, which is formalized in Alg.
III.1 and Alg. III.2.

  1. 1.

Create synthetic data locally around the chosen data point;

  2. 2.

Apply LIME to generate a surrogate model;

  3. 3.

Repeat 1 and 2 several times to produce an ensemble of decision boundaries.

  4. 4.

Take the interquartile range (for example) of those boundaries.

To make evident, Fig. 3 is produced through a single application of steps 1
and 2 while Fig. 4 is produced after the final step. It is important to
recognize that, in a practical or real-world application, it would be
infeasible to plot the entire decision regions to any reasonable level of
accuracy. The computational demands and complexities involved in plotting such
extensive regions can be prohibitive. This makes the region of indecision both
valuable and convenient in interpreting QNNs.

##  V Discussion & Future Work

In summary, in interpreting the classification output of a QNN, we have
defined Q-LIME and the local region of indecision, a simple description of a
region where the QNN is interpreted to be indecisive. The data samples within
their associated band should not be expected to have an “explanation” in the
deterministic classical sense. While directly useful for hybrid quantum-
classical models, we hope this stimulates further research on the fundamental
differences between quantum and classical model interpretability. In the
remainder of the paper, we discuss possible future research directions.

Our results are pointed squarely at the randomness of quantum measurements,
which might suggest that they are “backward compatible” with classical models.
Indeed, randomness is present in training classical DNNs due to random weight
initialization or the optimization techniques used. However, this type of
randomness can be “controlled” via the concept of a seed. Moreover, the
penultimate output of DNNs (just before label assignment) is often a
probability distribution. Nevertheless, these are produced through arbitrary
choices of the activation function (i.e., the Softmax function), which force
the output of a vector that merely resembles a probability distribution. Each
case of randomness in classical DNNs is very different from the innate and
unavoidable randomness of quantum models. While our techniques could be
applied in a classical setting, the conclusions drawn from them may ironically
be more complicated to justifiably action.

In this work, we have provided concrete definitions for QNNs applied to binary
classification problems. Using a probability of $\frac{1}{2}$ would not be a
suitable reference point in multi-class decision problems. There are many
avenues to generalize our definitions, which would mirror standard
generalizations in the move from binary to multi-class decision problems. One
such example would be defining the region of indecision as that with nearly
maximum entropy.

We took as an act of brevity the omission of the word local in many places
where it may play a pivotal role. For example, the strongest conclusion our
technique can reach is that a QNN is merely locally inexplicable in the
classical sense. In such cases, we could concede that (for some regions of
data space) the behavior of QNN is inherently stochastic, lacking any
discernible patterns or deterministic explanations. Alternatively, we can use
this conclusion to signal that an explanation at a higher level of abstraction
is required. This shift toward higher-level abstraction becomes imperative for
gaining more insights into QNNs, acknowledging their inherent quantum nature
and the limitations imposed by classical interpretability paradigms.

Classically, should a data point inquire, ”Why has this label been assigned to
me?”, within a quantum framework, an initial response may invoke the inherent
stochastic nature of quantum systems. Nonetheless, persistent interrogation
from the data point could necessitate a quantum extension of global
interpretability methodologies to clarify the specific features or attributes
contributing to the assigned label.

Referring back to Fig. 1, we have focused here on CQ quantum machine learning
models. Nevertheless, the core idea behind local surrogate models remains
applicable in the context of quantum data — use interpretable quantum models
as surrogate models to explain black box models producing quantum data. Of
course, one of our assumptions is the parallels between the classical
interpretable models we mentioned above, with their quantum equivalents. This
can be a line for future work. The ideas here encapsulate inherently quantum
models such as matrix product states or tensor network states, which can act
as surrogate models for quantum models as they may be considered more
interpretable. These surrogate models can serve as interpretable proxies,
facilitating a comprehension of the underlying quantum processes, thus paving
the way for enhanced interpretability of quantum models.

Furthermore, the idea behind interpreting or “opening up” black-box models may
be of interest in control theory [30, 31, 32]. In this scenario, the concept
of “grey-box” models — portions of which encode specific physical models —
give insights into how to engineer certain parameters in a system. These grey-
box models can thus be considered partially explainable models. The proposed
algorithm in [33] may also be of interest in terms of creating intrinsically
quantum interpretable models, which would act as surrogates for other more
complex quantum models.

Advanced ideas for future work entail exploring and developing specific
metrics tailored for assessing interpretability in QML models. This may
involve defining measures that capture the degree of coherence or entanglement
contributing to model predictions.

An obvious open question that inspires future research remains to investigate
the difference in computational tractability of interpretability methods in
quantum versus classical. This will lead to understanding whether it is more
difficult to interpret quantum models as opposed to classical models. We hope
such results shed light on more philosophical questions as well, such as is
inexplicability, viz. complexity, necessary for learning?. This exploration
holds the potential to advance our understanding not only of the
interpretability challenges in quantum models, but also of the fundamental
relationship between complexity and the learning process itself.

For completeness, the case for the interpretability of machine learning models
does not go without critique. Certain viewpoints contend that the pursuit of
insights into a model’s decision-making process should not come at the expense
of sacrificing performance, and, in practical terms, this prioritization might
not always be feasible or deemed necessary [34, 35]. It is often posited that
simpler models inherently possess a greater degree of explainability.
Nonetheless, it is the more complex models that require explanations, as they
may be more likely employed in critical applications.

Regardless of the two distinct camps of beliefs, the niche field of
interpretable machine learning keeps growing in volume. An argument is that
having a more complete picture of the model’s performance can help improve the
performance of the model overall. As QML becomes increasingly relevant to AI
research, we expect the need for quantum interpretability will also be in
demand. This research contributes to a broader understanding of explicability
in quantum AI models, paving the way for the development of accountable and
transparent systems in the quantum computing and AI era.

Acknowledgments: LP was supported by the Sydney Quantum Academy, Sydney, NSW,
Australia.

## References

  * Russell and Norvig [2010] S. Russell and P. Norvig,  _Artificial Intelligence: A Modern Approach_ , 3rd ed. (Prentice Hall, 2010).
  * Mitchell [1997] T. M. Mitchell,  _Machine Learning_ , 1st ed. (McGraw-Hill, Inc., USA, 1997).
  * Molnar [2022] C. Molnar,  _Interpretable Machine Learning_, 2nd ed. (2022).
  * Molnar  _et al._ [2020] C. Molnar, G. Casalicchio, and B. Bischl, Interpretable machine learning – a brief history, state-of-the-art and challenges, in  _ECML PKDD 2020 Workshops_ (Springer International Publishing, 2020) pp. 417–431.
  * Lipton [2018] Z. C. Lipton, The mythos of model interpretability, Queue 16, 31 (2018).
  * Goodfellow  _et al._ [2016] I. Goodfellow, Y. Bengio, and A. Courville,  _Deep Learning_ (MIT Press, 2016) http://www.deeplearningbook.org.
  * Zhang  _et al._ [2021] Y. Zhang, P. Tino, A. Leonardis, and K. Tang, A survey on neural network interpretability, IEEE Transactions on Emerging Topics in Computational Intelligence 5, 726 (2021).
  * Guidotti  _et al._ [2018] R. Guidotti, A. Monreale, S. Ruggieri, F. Turini, F. Giannotti, and D. Pedreschi, A survey of methods for explaining black box models, ACM computing surveys (CSUR) 51, 1 (2018).
  * Biamonte  _et al._ [2016] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Quantum Machine Learning, Nature 549, 195 (2016).
  * Schuld and Petruccione [2021] M. Schuld and F. Petruccione,  _Machine Learning with Quantum Computers_ (Springer International Publishing, 2021).
  * Cerezo  _et al._ [2022] M. Cerezo, G. Verdon, H.-Y. Huang, L. Cincio, and P. Coles, Challenges and opportunities in quantum machine learning, Nature Computational Science (2022).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Variational quantum algorithms, Nature Reviews Physics 3, 625 (2021).
  * Heese  _et al._ [2023] R. Heese, T. Gerlach, S. Mücke, S. Müller, M. Jakobs, and N. Piatkowski, Explaining quantum circuits with shapley values: Towards explainable quantum machine learning, arXiv preprint arXiv:2301.09138 10.48550/arXiv.2301.09138 (2023).
  * Mercaldo  _et al._ [2022] F. Mercaldo, G. Ciaramella, G. Iadarola, M. Storto, F. Martinelli, and A. Santone, Towards explainable quantum machine learning for mobile malware detection and classification, Applied Sciences 12, 10.3390/app122312025 (2022).
  * Steinmüller  _et al._ [2022] P. Steinmüller, T. Schulz, F. Graf, and D. Herr, explainable ai for quantum machine learning, arXiv preprint arXiv:2211.01441 10.48550/arXiv.2211.01441 (2022).
  * Lundberg and Lee [2017] S. M. Lundberg and S.-I. Lee, A unified approach to interpreting model predictions, in  _Advances in Neural Information Processing Systems_, Vol. 30 (Curran Associates, Inc., 2017).
  * Burge  _et al._ [2023] I. Burge, M. Barbeau, and J. Garcia-Alfaro, A quantum algorithm for shapley value estimation, arXiv preprint arXiv:2301.04727 10.48550/arXiv.2301.04727 (2023).
  * Ribeiro  _et al._ [2016] M. T. Ribeiro, S. Singh, and C. Guestrin, ”why should i trust you?”: Explaining the predictions of any classifier, in  _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_, KDD ’16 (Association for Computing Machinery, New York, NY, USA, 2016) p. 1135–1144.
  * Simonyan  _et al._ [2013] K. Simonyan, A. Vedaldi, and A. Zisserman, Deep inside convolutional networks: Visualising image classification models and saliency maps, arXiv preprint arXiv:1312.6034 10.48550/arXiv.1312.6034 (2013).
  * Olah  _et al._ [2017] C. Olah, A. Mordvintsev, and M. Tyka, Feature visualization: How neural networks build up their understanding of images, Distill 2 (2017).
  * Yosinski  _et al._ [2015] J. Yosinski, J. Clune, A. Nguyen, T. Fuchs, and H. Lipson, Understanding neural networks through deep visualization, in  _Proceedings of the 32nd International Conference on Machine Learning (ICML)_ (2015) pp. 1582–1591.
  * Bach  _et al._ [2015] S. Bach, A. Binder, G. Montavon, F. Klauschen, K.-R. Müller, and W. Samek, On pixel-wise explanations for non-linear classifier decisions by layer-wise relevance propagation, PLoS ONE 10, 10(7): e0130140 (2015).
  * Zeiler and Fergus [2014] M. D. Zeiler and R. Fergus, Visualizing and understanding convolutional networks, European conference on computer vision , 818 (2014).
  * Fisher [1936] R. A. Fisher, The use of multiple measurements in taxonomic problems (1936).
  * Spall [1998] J. C. Spall, An overview of the simultaneous perturbation method for efficient optimization, Johns Hopkins APL Technical Digest 19, 482 (1998).
  * IBM [2021] IBM, Qiskit: An Open-Source Framework for Quantum Computing (2021), accessed on May 1, 2023.
  * Schuld and Petruccione [2018] M. Schuld and F. Petruccione,  _Supervised Learning with Quantum Computers_, Quantum Science and Technology (Springer International Publishing, 2018).
  * Farhi and Neven [2018] E. Farhi and H. Neven, Classification with quantum neural networks on near term processors, arXiv preprint arXiv:1802.06002 10.48550/arXiv.1802.06002 (2018).
  * Pira and Ferrie [2023] L. Pira and C. Ferrie, Interpret QNN: Explicability and Inexplicability in the Interpretation of Quantum Neural Networks, https://github.com/lirandepira/interpret-qnn (2023), accessed on July 24, 2023\.
  * Youssry  _et al._ [2023a] A. Youssry, Y. Yang, R. J. Chapman, B. Haylock, F. Lenzini, M. Lobino, and A. Peruzzo, Experimental graybox quantum system identification and control, arXiv preprint arXiv:2206.12201 10.48550/arXiv.2206.12201 (2023a).
  * Youssry  _et al._ [2020] A. Youssry, G. A. Paz-Silva, and C. Ferrie, Characterization and control of open quantum systems beyond quantum noise spectroscopy, npj Quantum Information 6, 95 (2020).
  * Youssry  _et al._ [2023b] A. Youssry, G. A. Paz-Silva, and C. Ferrie, Noise detection with spectator qubits and quantum feature engineering, New Journal of Physics 25, 073004 (2023b).
  * Weitz  _et al._ [2023] G. Weitz, L. Pira, C. Ferrie, and J. Combes, Sub-universal variational circuits for combinatorial optimization problems (2023).
  * Sarkar [2022] A. Sarkar, Is explainable AI a race against model complexity?, arXiv preprint arXiv:2205.10119 10.48550/arXiv.2205.10119 (2022).
  * McCoy  _et al._ [2022] L. McCoy, C. Brenna, S. Chen, K. Vold, and S. Das, Believing in Black Boxes: Must Machine Learning in Healthcare Be Explainable to Be Evidence-Based?, Journal of Clinical Epidemiology 10.1016/j.jclinepi.2021.11.001 (2022).