# ALMA Subparsec Resolution Dense Molecular Line Observations of the NGC 1068 Nucleus

Masatoshi Imanishi, Bernd Vollmer, Yoshiaki Hagiwara, Kouichiro Nakanishi, Takuma Izumi, Nozomu Kawakatu

► To cite this version:

Masatoshi Imanishi, Bernd Vollmer, Yoshiaki Hagiwara, Kouichiro Nakanishi, Takuma Izumi, et al. ALMA Subparsec Resolution Dense Molecular Line Observations of the NGC 1068 Nucleus. The Astrophysical Journal, 2025, 995 (2), pp.211. 10.3847/1538-4357/ae1caa. insu-05446284

HAL Id: insu-05446284

https://insu.hal.science/insu-05446284v1

Submitted on 7 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# ALMA Subparsec Resolution Dense Molecular Line Observations of the NGC 1068 Nucleus

Masatoshi Imanishi<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID, Bernd Vollmer<fnref n="4" />, Yoshiaki Hagiwara<fnref n="5" /> ID, Kouichiro Nakanishi<fnref n="1" />,<fnref n="2" /> ID, Takuma Izumi<fnref n="1" />,<fnref n="2" /> ID, and Nozomu Kawakatu<fnref n="6" /> ID

## Abstract

We present the results of our ALMA observations of the dense molecular HCN *J* = 4–3 and HCO$^+$ *J* = 4–3 lines at ≲1 pc (≪14 mas) resolution in the nuclear region of the nearby (~14 Mpc) well-studied AGN NGC 1068. Both emission lines are clearly detected around the AGN along an almost east-west direction, which we ascribe to the dusty molecular torus. The HCN *J* = 4–3 emission is brighter than the HCO$^+$ *J* = 4–3 emission in the compact (≪3–5 pc) torus region. Apparent counterrotation between the inner (≪2 pc) and outer (≫2 pc) parts of the western torus, previously seen in ~1.5 pc-resolution HCN *J* = 3–2 and HCO$^+$ *J* = 3–2 data, is also confirmed in our new ≲1 pc-resolution HCN *J* = 4–3 and HCO$^+$ *J* = 4–3 data. We apply a physically counterrotating torus model, in which a compact dense gas clump collided with the western side of the existing rotating torus from the opposite direction, and we find that this model largely reproduces the observed properties of the combined new ≲1 pc-resolution HCN *J* = 4–3 and HCO$^+$ *J* = 4–3 data, and the previously obtained ≲1.5 pc-resolution HCN *J* = 3–2 and HCO$^+$ *J* = 3–2 data.

*Unified Astronomy Thesaurus concepts: Active galactic nuclei (16); Molecular gas (1073); Submillimeter astronomy (1647); Radio interferometry (1346); Supermassive black holes (1663); Active galaxies (17); Millimeter astronomy (1061); Molecular spectroscopy (2095); Water masers (1790); Interferometry (808)*

## 1. Introduction

An active galactic nucleus (AGN) shines brightly in the UV-optical because of strong emission from the accretion disk around a mass-accreting supermassive black hole (SMBH). Some AGNs show broad (≫1000 km s⁻<fnref n="1" />) emission lines in the UV-optical spectra (classified as type 1), while others do not (type 2). These AGN properties can naturally be explained by the presence of toroidally distributed, optically thick dust and gas on a ≲10 pc physical scale, the so-called dusty molecular torus (e.g., R. R. J. Antonucci 1993; S. F. Honig 2019). Fast-rotating gas clouds in the vicinity of the SMBH, dominated by its gravity and photoionized by UV radiation from the accretion disk, can display broad (≫1000 km s⁻<fnref n="1" />) UV-optical emission lines. They are detectable if seen from a direction not blocked by the torus (type 1) but become undetectable if obscured by the torus along our line of sight (type 2). This torus-based AGN unification picture was first proposed based on observations of a very nearby (~14 Mpc, *z* = 0.0038, 1" = 70 pc) type 2 AGN, NGC 1068 (R. R. J. Antonucci & J. S. Miller 1985). Detailed observations of the torus in NGC 1068 are thus very important to better understand the torus-based AGN unification paradigm and update it if necessary.

A spatially resolved observational study of the torus is obviously one of the best means to investigate its properties in detail. For this purpose, high-angular-resolution observations

with ≲0.05 are desirable because the physical size of the putative torus is ≲10 pc, which corresponds to ≲0.15 even for very nearby (~15 Mpc) AGNs. ALMA high-angular-resolution (≪0.05) molecular line observations are a very powerful tool to scrutinize the morphological, physical, chemical, and dynamical properties of the torus because (1) very high angular resolution (≪0.02) and high sensitivity are achievable and (2) multiple rotational (J) transition line observations of multiple molecules enable us to investigate the density, temperature, abundance, and kinematics of molecular gas in the torus. Previously conducted ALMA high-angular-resolution (~0.02–0.04 [20–40 mas] or 1.5–3 pc) molecular line observations of NGC 1068, mainly using the *J* = 3–2 lines of CO, HCN, and HCO$^+$, have revealed the presence of compact (≪3–5 pc) dense and warm molecular gas along an almost east-west direction (e.g., S. Garcia-Burillo et al. 2016; M. Imanishi et al. 2018b; S. Garcia-Burillo et al. 2019; C. M. V. Impellizzeri et al. 2019; M. Imanishi et al. 2020), which can naturally be interpreted as originating in the putative dusty molecular torus. However, the observed torus molecular gas properties of NGC 1068 are not in perfect agreement with the simple theory that predicts that molecular gas in the torus should rotate with an almost Keplerian motion dominated by the gravity of the central SMBH (e.g., K. Wada et al. 2009, 2016). It was found that the inner (≪2 pc) and outer (≫2 pc) parts of the NGC 1068 torus do not rotate in the same way but appear to be counterrotating relative to each other (e.g., S. Garcia-Burillo et al. 2019; C. M. V. Impellizzeri et al. 2019; M. Imanishi et al. 2020). The first proposed

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ National Astronomical Observatory of Japan, National Institutes of Natural Sciences (NINS), 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan; masa.imanishi@nao.ac.jp, m.imanishi.astro@gmail.com</footnote>

<footnote n="2">² Department of Astronomy, School of Science, Graduate University for Advanced Studies (SOKENDAI), Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="3">³ Toyo University, 5-28-20, Hakusan, Bunkyo-ku, Tokyo 112-8606, Japan</footnote>

<footnote n="4">⁴ Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 67000, Strasbourg, France</footnote>

<footnote n="5">⁵ Natural Science Laboratory, Toyo University, 5-28-20, Hakusan, Bunkyo-ku, Tokyo 112-8606, Japan</footnote>

<footnote n="6">⁶ Faculty of Natural Sciences, National Institute of Technology, Kure College, 2-2-11 Agaminami, Kure, Hiroshima 737-8506, Japan
Received 2025 September 3; revised 2025 October 22; accepted 2025 November 4; published 2025 December 18</footnote>
<table><caption>**Table 1**
Log of ALMA Cycle 9 Observations of NGC 1068</caption><thead><tr><th rowspan="2">Line<br>(1)</th><th rowspan="2">Date<br>(UT)<br>(2)</th><th rowspan="2">Antenna Number<br>(3)</th><th rowspan="2">Baseline<br>(m)<br>(4)</th><th rowspan="2">Integration<br>(minute)<br>(5)</th><th colspan="3">Calibrator</th></tr><tr><th>Bandpass<br>(6)</th><th>Flux<br>(7)</th><th>Phase<br>(8)</th></tr></thead><tbody><tr><td rowspan="10">HCN/HCO<sup>+</sup> J = 4–3</td><td>2023 Jul 28</td><td>48</td><td>230–16196</td><td>44</td><td>J0238+1636</td><td>J0238+1636</td><td>J0239–0234</td></tr><tr><td>2023 Jul 29</td><td>44</td><td>230–16196</td><td>44</td><td>J0238+1636</td><td>J0238+1636</td><td>J0239–0234</td></tr><tr><td>2023 Jul 31</td><td>46</td><td>230–16196</td><td>44</td><td>J0238+1636</td><td>J0238+1636</td><td>J0239–0234</td></tr><tr><td>2023 Aug 3</td><td>49</td><td>230–16196</td><td>44</td><td>J0238+1636</td><td>J0238+1636</td><td>J0239–0234</td></tr><tr><td>2023 Aug 4</td><td>49</td><td>230–16196</td><td>44</td><td>J0238+1636</td><td>J0238+1636</td><td>J0239–0234</td></tr><tr><td>2023 Aug 9</td><td>44</td><td>83–16196</td><td>44</td><td>J0238+1636</td><td>J0238+1636</td><td>J0239–0234</td></tr><tr><td>2023 Aug 11</td><td>41</td><td>83–16196</td><td>44</td><td>J0238+1636</td><td>J0238+1636</td><td>J0239–0234</td></tr></tbody></table>

**Note.** Column (1): Line. Column (2): Observation date in UT. Column (3): Number of antennas used for observations. Column (4): Baseline length in meters. Minimum and maximum baseline lengths are shown. Column (5): Net on source integration time in minutes. Columns (6)–(8): Bandpass, flux, and phase calibrators for the target source (NGC 1068), respectively.

scenario to explain this peculiar observational result speculated that a collision of a very compact gas clump with the western side of the existing rotating torus from the opposite direction could create a bona fide counterrotating molecular torus (e.g., C. M. V. Impellizzeri et al. 2019; M. Imanishi et al. 2020; B. Vollmer et al. 2022). A second scenario hypothesized that outflow activity could reproduce the observed apparent counterrotation patterns (e.g., S. Garcia-Burillo et al. 2019; D. Williamson et al. 2020; E. Y. Bannikova et al. 2023).

The achievable angular resolution can be even smaller for the *J* = 4–3 lines of HCN and HCO⁺ at ~0.85 mm (~350 GHz) than their *J* = 3–2 lines at ~1.2 mm (~260 GHz) because of the shorter wavelength (higher frequency) of the *J* = 4–3 transitions. Using the longest-baseline (~16 km) observations of ALMA, we can achieve ~0.014 (~14 mas) angular resolution for the *J* = 4–3 lines, which corresponds to subparsec (~1 pc) physical resolution for NGC 1068 (1" = 70 pc). Very recently, V. Gamez Rosas et al. (2025) presented ~0.016 (~1.1 pc) resolution HCO⁺ *J* = 4–3 and CO *J* = 3–2 line data of NGC 1068, and confirmed the apparent counter-rotation between the inner (~2 pc) and outer (≥2 pc) torus regions, which was particularly clear in the western torus. These authors proposed a new third scenario in which an infalling compact gas cloud passes from the west, in front of the ~0.85 mm (~350 GHz) continuum-emitting central engine of the AGN (=mass-accreting SMBH), without strongly interacting with the dusty molecular torus. Because (1) HCN emission is brighter than HCO⁺ emission in the compact (~3–5 pc) molecular torus of NGC 1068, most likely because of HCN-abundance enhancement (M. Imanishi et al. 2018b; M. Imanishi et al. 2020; J. Butterworth et al. 2022; B. Vollmer et al. 2022), and (2) HCN and HCO⁺ have different critical densities (Y. L. Shirley 2015), the highest-angular-resolution HCN *J* = 4–3 line data will enable us to better elucidate the dense molecular gas properties of the NGC 1068 torus. The addition of high-angular-resolution *J* = 4–3 data of HCN and HCO⁺ to the existing ~1.5 pc resolution *J* = 3–2 data of HCN and HCO⁺ (e.g., C. M. V. Impellizzeri et al. 2019; M. Imanishi et al. 2020) can also provide important information on molecular gas excitation conditions inside the ~3–5 pc extended NGC 1068 torus in a spatially resolved manner. In this paper, we present ALMA ~1 pc resolution HCN *J* = 4–3 and HCO⁺ *J* = 4–3 line observations of the NGC 1068 nucleus to further scrutinize various properties of the dense molecular torus around the AGN in NGC 1068.

## 2. Observations and Data Analysis

Our HCN *J* = 4–3 and HCO⁺ *J* = 4–3 line observations of the NGC 1068 nucleus were conducted in our ALMA Cycle 9 program 2022.1.00005.S (PI, M. Imanishi). To achieve the highest angular resolution, the longest-baseline (~16 km) configuration was used. We employed the widest 1.875 GHz mode with 1920 channels for each spectral window. The HCN *J* = 4–3 (rest frame frequency $\nu_{\text{rest}}$ = 354.505 GHz) and HCO⁺ *J* = 4–3 ($\nu_{\text{rest}}$ = 356.734 GHz) lines were simultaneously observed in one sideband (upper sideband; USB). The vibrationally excited HCN $\nu_2$ = 1, *l* = 1f (HCN-VIB) *J* = 4–3 line ($\nu_{\text{rest}}$ = 356.256 GHz) is also covered in the USB, while the vibrationally excited HCO⁺ $\nu_2$ = 1, *l* = 1f (HCO⁺-VIB) *J* = 4–3 line ($\nu_{\text{rest}}$ = 358.242 GHz) is not covered. In lower sideband, the CS *J* = 7–6 line ($\nu_{\text{rest}}$ = 342.883 GHz) is covered. Table 1 summarizes our ALMA Cycle 9 observation log.

We started our analysis with pipeline-calibrated data provided by ALMA, using CASA 6.5.4-9 (The CASA Team et al. 2022; T. R. Hunter et al. 2023). We determined the continuum level by selecting channels that did not show obvious emission or absorption lines and subtracted it using the CASA task “uvcontsub.” Then, we applied the “tclean” task (Briggs-weighting; robust = 0.5, gain = 0.1) to the continuum-subtracted molecular line and continuum data, with a pixel scale of 0.0015 pixel⁻¹ (1.5 mas pixel⁻¹) and a velocity resolution of 10 km s⁻¹. The derived rms noise level for the 10 km s⁻¹ velocity resolution is ~0.3 mJy beam⁻¹ for both HCN *J* = 4–3 and HCO⁺ *J* = 4–3. The maximum recoverable scale of our ALMA Cycle 9 data is ~0.14 (ALMA Cycle 9 Proposer’s Guide), which corresponds to ~10 pc at the distance of NGC 1068, and is thus large enough to properly investigate the morphology and kinematics of the compact (~3–5 pc) dense molecular torus around the AGN in NGC 1068 (e.g., S. Garcia-Burillo et al. 2016; M. Imanishi et al. 2016a, 2018b, 2020; S. Garcia-Burillo et al. 2019; C. M. V. Impellizzeri et al. 2019; B. Vollmer et al. 2022; V. Gamez Rosas et al. 2025). The absolute flux calibration uncertainty is expected to be ~10% at the observed frequency (~350 GHz) of the ALMA Cycle 9 data.

M. Imanishi et al. (2020) presented cleaned maps of ~1.3 pc (~0.019) resolution HCN *J* = 3–2 ($\nu_{\text{rest}}$ = 265.886 GHz) and HCO⁺ *J* = 3–2 ($\nu_{\text{rest}}$ = 267.558 GHz) line data taken in ALMA Cycle 6 (2018.1.00037.S), with pixel scale of 0.003 pixel⁻¹ (3 mas pixel⁻¹). We reanalyzed the HCN *J* = 3–2 and HCO⁺ *J* = 3–2 data with pixel scale of 1.5 mas pixel⁻¹ and
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Integrated-intensity (moment 0) map of (a) HCN *J* = 4–3, (b) HCO⁺ *J* = 4–3, (c) HCN *J* = 3–2, and (d) HCO⁺ *J* = 3–2, created with a pixel scale of 1.5 mas pixel⁻¹. The gray contours indicate simultaneously obtained continuum emission with 50σ, 70σ, 90σ, and 110σ ((a)–(b) at ~350 GHz) and 40σ, 60σ, and 80σ ((c)–(d) at ~260 GHz), where the 1σ value is shown in Table 2. The horizontal thick black bar at the bottom center of each panel corresponds to 1 pc at the distance of NGC 1068. The filled white circle at the lower left of each panel indicates the synthesized beam size, which is 0′′014 × 0′′013 for the *J* = 4–3 lines of HCN and HCO⁺ (a–b), and 0′′020 × 0′′018 for the reanalyzed images of the *J* = 3–2 lines (c–d).</figcaption></figure>

<table><caption>**Table 2**

Continuum Emission</caption><thead><tr><th>Frequency<br>(GHz)<br>(1)</th><th>Peak flux<br>(mJy beam<sup>-1</sup>)<br>(2)</th><th>Peak coordinate<br>(R.A., Decl.) ICRS<br>(3)</th><th>rms<br>(mJy beam<sup>-1</sup>)<br>(4)</th><th>Synthesized beam<br>("×") (deg)<br>(5)</th></tr></thead><tbody><tr><td>340.4–344.2, 352.4–356.2</td><td>6.8 (116σ)</td><td>(02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>709<sup>s</sup>, −00°00′47′′946)</td><td>0.059</td><td>0.014 × 0.013 (−86°)</td></tr><tr><td>263.6–269.1</td><td>6.7 (89σ)</td><td>(02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>709, −00°00′47′′946)</td><td>0.075</td><td>0.019 × 0.017 (56°)</td></tr></tbody></table>

**Note.** Column (1): Frequency range (in GHz) used to extract continuum emission. Frequencies of obvious emission and absorption lines are excluded. Column (2): Flux (in mJy beam⁻¹) at the emission peak. The value at the highest flux pixel (1.5 mas pixel⁻¹) is adopted. The detection significance relative to the rms noise is shown in parentheses. Possible systematic uncertainties coming from absolute flux calibration ambiguity in individual ALMA observations and the choice of frequency range to determine the continuum level are not taken into account. Column (3): Coordinate of the continuum emission peak in ICRS. Column (4): The rms noise level (1σ) (in mJy beam⁻¹), derived from the standard deviation of sky signals. Column (5): Synthesized beam (in arcsec × arcsec) and position angle (in degrees). The position angle is 0° along the north–south direction and increases in the counterclockwise direction.

velocity resolution of 10 km s⁻¹ to compare them with the newly taken, ≲1 pc resolution HCN *J* = 4–3 and HCO⁺ *J* = 4–3 line data in a straightforward manner. The derived rms noise level for the 10 km s⁻¹ velocity resolution is ~0.2 mJy beam⁻¹ for both HCN *J* = 3–2 and HCO⁺ *J* = 3–2. The absolute flux calibration uncertainty of the ALMA Cycle 6 data at ~260 GHz is also expected to be ≲10% (ALMA Cycle 6 Proposer's Guide).

### 3. Results

Table 2 summarizes the continuum emission properties at ~350 GHz and the synthesized beam size (0′′014 × 0′′013), as well as the reanalyzed ~260 GHz data. Figures 1(a) and 1(b) show ~350 GHz continuum (contours) and integrated-intensity (moment 0) maps of HCN *J* = 4–3 and HCO⁺ *J* = 4–3 lines, respectively. Almost east–west oriented, compact (≈3–5 pc) HCN *J* = 4–3 and HCO⁺ *J* = 4–3 emission is clearly detected around the continuum peak (defined as “C-peak”), which we regard as the position of the mass-accreting SMBH. We find that (1) the HCN *J* = 4–3 emission is brighter than the HCO⁺ *J* = 4–3 emission in both the western and eastern torus, and (2) the western torus is brighter than the eastern torus for both the HCN *J* = 4–3 and HCO⁺ *J* = 4–3 lines. Continuum emission at ~260 GHz (contours) and reanalyzed moment 0 maps of HCN *J* = 3–2 and HCO⁺ *J* = 3–2 lines are displayed in Figures 1(c) and

1(d), respectively. In Figures 1(a)–(d), the *J* = 4–3 emission is spatially more compact than the *J* = 3–2 emission for both HCN and HCO⁺. Although the smaller synthesized beam size at *J* = 4–3 could contribute to the observed more compact emission, we interpret that the *J* = 4–3 lines probe warmer and denser molecular gas closer to the central AGN engine than the *J* = 3–2 lines because of the higher excitation energy and critical density of the former (see the last paragraph of Section 3).

We define the locations of the HCN *J* = 4–3 emission peaks in the western and eastern torus in Figure 1(a) as “W-peak” and “E-peak,” respectively, and use these peaks for the following discussion for two reasons: (1) the HCN emission is brighter than the HCO⁺ emission, and (2) the newly taken *J* = 4–3 data have a smaller synthesized beam size than the reanalyzed *J* = 3–2 data, both of which make the HCN *J* = 4–3 data the best for determining the peak position in the most accurate manner. The coordinates of the peak positions of the HCN *J* = 4–3 and other emission lines are summarized in Table 3.

Figure 2 displays beam-sized spectra at the continuum emission peak (C-peak) and at the HCN *J* = 4–3 emission peaks in the western and eastern torus (W-peak and E-peak, respectively). At the C-peak, the HCN *J* = 4–3 and HCO⁺ *J* = 4–3 lines are observed in absorption (Figure 2(a)). At the W-peak and E-peak, the HCN *J* = 4–3 and HCO⁺ *J* = 4–3 emission lines are clearly detected (Figure 2(b), (c)).
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Observed beam-sized spectra at (a) ~350 GHz continuum emission peak (02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>709, −00<sup>o</sup>00′47″946) ICRS (defined as “C-peak”), (b) HCN *J* = 4–3 emission peak in the western torus (02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>7084, −00<sup>o</sup>00′47″940) ICRS (defined as “W-peak”), and (c) HCN *J* = 4–3 emission peak in the eastern torus (02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>710, −00<sup>o</sup>00′47″946) ICRS (defined as “E-peak”). The abscissa is the observed frequency (in GHz) and the ordinate is the flux density (in mJy beam<sup>−1</sup>). The beam size is 0′014 × 0′013 (0.98 pc × 0.91 pc). Downward arrows indicate the expected frequency of HCN *J* = 4–3 (ν<sub>rest</sub> = 354.505 GHz), HCO<sup>+</sup> *J* = 4–3 (ν<sub>rest</sub> = 356.734 GHz), and HCN-VIB *J* = 4–3 (ν<sub>rest</sub> = 356.256 GHz) lines for the systemic optical LSR velocity of NGC 1068, V<sub>sys</sub> = 1130 km s<sup>−1</sup> (e.g., M. Imanishi et al. 2016a, 2018b; C. M. V. Impellizzeri et al. 2019; M. Imanishi et al. 2020; J. F. Gallimore et al. 2024). The horizontal dotted straight line indicates the zero-flux level.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Velocity profile of (a) HCN *J* = 4–3, (b) HCO<sup>+</sup> *J* = 4–3, and (c) HCN-VIB *J* = 4–3 line in the beam-sized spectrum toward the ~350 GHz continuum emission peak ("C-peak"). The abscissa is the optical LSR velocity (in km s<sup>−1</sup>) and the ordinate is the flux density (in mJy beam<sup>−1</sup>). A downward arrow is added at the systemic velocity of NGC 1068 (V<sub>sys</sub> = 1130 km s<sup>−1</sup>), labeled as "V<sub>sys</sub>." The horizontal thin dotted straight line indicates the zero-flux level. In (b), an upward arrow is added at the velocity of V = 1530 km s<sup>−1</sup> (V<sub>sys</sub> + 400 km s<sup>−1</sup>) to indicate an absorption feature that is 400 km s<sup>−1</sup> redshifted relative to the systemic velocity. In (c), a red upward arrow is added at the velocity of V = 1470 km s<sup>−1</sup> (V<sub>sys</sub> + 340 km s<sup>−1</sup>) to indicate redshifted emission at the innermost (≈1 pc) western torus (see Section 3).</figcaption></figure>

<table><caption>**Table 3**
Molecular Emission Line Peak Position</caption><thead><tr><th>Line<br>(1)</th><th>Torus Position<br>(2)</th><th>Peak Flux<br>(Jy km s<sup>-1</sup> beam<sup>-1</sup>)<br>(3)</th><th>Peak Coordinate<br>(R.A., Decl.) ICRS<br>(4)</th><th>rms<br>(Jy km s<sup>-1</sup> beam<sup>-1</sup>)<br>(5)</th><th>Synthesized Beam<br>("&times;") (deg)<br>(6)</th></tr></thead><tbody><tr><td rowspan="3">HCN J = 4-3</td><td>W-torus</td><td>0.77 (13.9&sigma;)</td><td>(02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>7084, -00&deg;00&prime;47&Prime;940)</td><td rowspan="3">0.056</td><td rowspan="3">0.014 &times; 0.013 (-83&deg;)</td></tr><tr><td>E-torus</td><td>0.63 (11.4&sigma;)</td><td>(02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>710, -00&deg;00&prime;47&Prime;946)</td></tr><tr><td>C-abs</td><td>-0.14 (-2.6&sigma;)</td><td>(02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>709, -00&deg;00&prime;47&Prime;949)</td></tr><tr><td rowspan="3">HCO<sup>+</sup> J = 4-3</td><td>W-torus</td><td>0.68 (11.8&sigma;)</td><td>(02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>709, -00&deg;00&prime;47&Prime;939)</td><td rowspan="3">0.058</td><td rowspan="3">0.014 &times; 0.013 (-83&deg;)</td></tr><tr><td>E-torus</td><td>0.41 (7.1&sigma;)</td><td>(02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>710, -00&deg;00&prime;47&Prime;946)</td></tr><tr><td>C-abs</td><td>-0.25 (-4.3&sigma;)</td><td>(02<sup>h</sup>42<sup>m</sup>-40<sup>s</sup>709, -00&deg;00&prime;47&Prime;951)</td></tr><tr><td rowspan="3">HCN J = 3-2</td><td>W-torus</td><td>0.44 (11.7&sigma;)</td><td>(02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>708, -00&deg;00&prime;47&Prime;940)</td><td rowspan="3">0.037</td><td rowspan="3">0.020 &times; 0.018 (54&deg;)</td></tr><tr><td>E-torus</td><td>0.25 (6.8&sigma;)</td><td>(02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>710, -00&deg;00&prime;47&Prime;950)</td></tr><tr><td>C-abs</td><td>-0.18 (-4.9&sigma;)</td><td>(02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>709, -00&deg;00&prime;47&Prime;949)</td></tr><tr><td rowspan="3">HCO<sup>+</sup> J = 3-2</td><td>W-torus</td><td>0.34 (10.8&sigma;)</td><td>(02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>708, -00&deg;00&prime;47&Prime;940)</td><td rowspan="3">0.032</td><td rowspan="3">0.020 &times; 0.018 (55&deg;)</td></tr><tr><td>E-torus</td><td>0.30 (9.4&sigma;)</td><td>(02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>710, -00&deg;00&prime;47&Prime;945)</td></tr><tr><td>C-abs</td><td>-0.22 (-6.9&sigma;)</td><td>(02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>709, -00&deg;00&prime;47&Prime;947)</td></tr><tr><td>CS J = 7-6</td><td>W-torus</td><td>0.072 (5.8&sigma;)</td><td>(02<sup>h</sup>42<sup>m</sup>40<sup>s</sup>709, -00&deg;00&prime;47&Prime;937)</td><td>0.012</td><td>0.014 &times; 0.013 (-82&deg;)</td></tr></tbody></table>

**Note.** Column (1): Molecular line. Column (2): Western torus (W-torus), eastern torus (E-torus), or absorption peak near the continuum emission peak position (C-abs). Column (3): Flux (in Jy km s<sup>−1</sup> beam<sup>−1</sup>) at the emission or absorption peak. The value at the highest or most negative flux pixel (1.5 mas pixel<sup>−1</sup>) in the integrated-intensity (moment 0) map in Figures 1 or 13(a) is adopted. The detection significance relative to the rms noise is shown in parentheses. Possible systematic uncertainties coming from absolute flux calibration ambiguity in individual ALMA observations and the choice of frequency range used to determine the continuum level are not taken into account. Column (4): Coordinate of the emission or absorption peak in ICRS. Column (5): The rms noise level (1&sigma;) (in Jy km s<sup>−1</sup> beam<sup>−1</sup>), derived from the standard deviation of sky signals. Column (6): Synthesized beam (in arcsec &times; arcsec) and position angle (in degrees). The definition of the position angle is the same as that in Table 2.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Upper panels: Observed area-integrated spectra of HCN *J* = 4–3 and HCO⁺ *J* = 4–3, extracted with (a) a 1.7 pc (24 mas) diameter circular aperture centered at the W-peak, (b) a 1.7 pc diameter circular aperture centered at the E-peak, (c) a 3.4 pc (48 mas) diameter circular aperture centered at the C-peak, and (d) a 6.8 pc (96 mas) diameter circular aperture centered at the C-peak. Lower panels: Observed area-integrated spectra of HCN *J* = 3–2 and HCO⁺ *J* = 3–2, extracted with (e) a 1.7 pc (24 mas) diameter circular aperture centered at the W-peak, (f) a 1.7 pc diameter circular aperture centered at the E-peak, (g) a 3.4 pc (48 mas) diameter circular aperture centered at the C-peak, and (h) a 6.8 pc (96 mas) diameter circular aperture centered at the C-peak. The abscissa is the observed frequency (in GHz) and the ordinate is the flux density (in mJy). Downward arrows are added at the expected frequencies of the HCN, HCO⁺, and HCN-VIB lines at *J* = 4–3 and *J* = 3–2, corresponding to the systemic velocity of *V*<sub>sys</sub> = 1130 km s<sup>-1</sup>. The horizontal dotted straight line indicates the zero-flux level.</figcaption></figure>

<table><caption>**Table 4**
Molecular Emission Line Flux</caption><thead><tr><th>Torus Position<br>(1)</th><th>HCN <i>J</i> = 4–3<br>(Jy km s<sup>-1</sup>)<br>(2)</th><th>HCO<sup>+</sup> <i>J</i> = 4–3<br>(Jy km s<sup>-1</sup>)<br>(3)</th><th>HCN <i>J</i> = 3–2<br>(Jy km s<sup>-1</sup>)<br>(4)</th><th>HCO<sup>+</sup> <i>J</i> = 3–2<br>(Jy km s<sup>-1</sup>)<br>(5)</th></tr></thead><tbody><tr><td>W-torus (1.7 pc diameter)</td><td>0.96 ± 0.10</td><td>0.60 ± 0.10</td><td>0.31 ± 0.07</td><td>0.18 ± 0.02</td></tr><tr><td>E-torus (1.7 pc diameter)</td><td>0.73 ± 0.05</td><td>0.47 ± 0.06</td><td>0.17 ± 0.02</td><td>0.14 ± 0.02</td></tr><tr><td>Torus (3.4 pc diameter)</td><td>2.5 ± 0.2</td><td>1.5 ± 0.1</td><td>0.76 ± 0.15</td><td>0.54 ± 0.05</td></tr><tr><td>Torus (6.8 pc diameter)</td><td>4.3 ± 0.3</td><td>2.5 ± 0.2</td><td>2.0 ± 0.1</td><td>1.4 ± 0.1</td></tr></tbody></table>

**Note.** Column (1): Torus region. Columns (2)–(5): Gaussian-fit velocity-integrated emission line fluxes of (2) HCN *J* = 4–3, (3) HCO⁺ *J* = 4–3, (4) HCN *J* = 3–2, and (5) HCO⁺ *J* = 3–2, in units of Jy km s⁻¹.

The velocity profiles of the HCN *J* = 4–3 and HCO⁺ *J* = 4–3 lines toward the C-peak in the beam-sized spectrum are shown in Figure 3. Absorption features at the systemic velocity (*V*<sub>sys</sub> = 1130 km s<sup>-1</sup>) are clearly detected for both lines<sup>7</sup> at the C-peak (Figures 2(a) and 3(a), (b)), as previously seen for the *J* = 3–2 lines of HCN and HCO⁺ at the ~260 GHz continuum emission peak in ~1.5 pc (~0.02) beam-sized spectra (e.g., C. M. V. Impellizzeri et al. 2019; M. Imanishi et al. 2020).

For HCO⁺ *J* = 4–3, an absorption feature is seen at ~400 km s⁻¹ redshifted with respect to the systemic velocity (Figure 3(b)). The ~400 km s⁻¹ redshifted HCO⁺ *J* = 4–3 absorption feature was also reported by V. Gamez Rosas et al. (2025), who interpreted it as due to a highly redshifted (~400 km s⁻¹) infalling gas clump toward the continuum-emitting mass-accreting SMBH, without strong interaction with the torus. However, we note that the observed frequency

spectrally coincides with the frequency at the systemic velocity (*V*<sub>sys</sub> = 1130 km s<sup>-1</sup>), of the HCN-VIB *J* = 4–3 line (*ν*<sub>rest</sub> = 356.256 GHz), which is ~400 km s<sup>-1</sup> redshifted relative to the HCO⁺ *J* = 4–3 line. We discuss these two scenarios in Appendix A.

For HCN *J* = 4–3, a broad absorption tail on the blueshifted side (*V* = 700–1100 km s<sup>-1</sup> or *V*<sub>sys</sub> – (30–430) km s<sup>-1</sup>) is seen (Figure 3(a)). The same absorption tail was previously detected for HCN *J* = 3–2 in the ~1.5 pc (~0.02) beam-sized spectrum toward the ~260 GHz continuum emission peak and was interpreted as originating from a nuclear outflow (C. M. V. Impellizzeri et al. 2019; M. Imanishi et al. 2020). The corresponding absorption tail is less clear in HCO⁺ *J* = 4–3 (Figure 3(b)) and HCO⁺ *J* = 3–2 (M. Imanishi et al. 2020).

Area-integrated spectra of the western and eastern torus, obtained with a 1.7 pc (24 mas) diameter circle around the W-peak and E-peak, and spectra that include both sides of the tori, obtained with 3.4 pc (48 mas) and 6.8 pc (96 mas) diameter circles around the C-peak, are shown in Figure 4. For the HCN emission, an excess emission tail is more significantly recognizable on the lower-frequency side of the main emission component

<footnote n="7"><sup>7</sup> V. Gamez Rosas et al. (2025) mentioned that this absorption feature at the systemic velocity was not detected for HCO⁺ *J* = 4–3 in their analysis using their data with comparable beam size (~0.016 or ~1.1 pc) to ours. We detect this absorption feature for HCO⁺ *J* = 4–3 in our data using the analysis method described in Section 2, as in the case for HCO⁺ *J* = 3–2 (M. Imanishi et al. 2020).</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Intensity-weighted mean velocity (moment 1) map of (a) HCN *J* = 4–3, (b) HCO⁺ *J* = 4–3, (c) HCN *J* = 3–2, and (d) HCO⁺ *J* = 3–2, created with 1.5 mas pixel⁻¹. The displayed velocity range is 1130±320 km s⁻¹ (960–1450 km s⁻¹) for (a)–(d), to clearly display the innermost (≈1 pc) high-velocity redshifted HCN emission in the western torus (e.g., C. M. V. Impellizzeri et al. 2019; M. Imanishi et al. 2020) (see text in Section 3). Moment 1 maps of (e) HCN *J* = 4–3 and (f) HCN *J* = 3–2 are shown with the velocity ranges of 1130±110 km s⁻¹ (= 1020–1240 km s⁻¹) and 1130 ± 60 km s⁻¹ (= 1070–1190 km s⁻¹), respectively, to display overall redshifted and blueshifted motion of the torus dense molecular gas emission with respect to the systemic velocity (*V*<sub>sys</sub> = 1130 km s⁻¹). The plus mark is the ~350 GHz continuum emission peak position (02h 42m 40s709, −00°00′ 47′′946)ICRS defined as C-peak. The open diamond mark is the emission peak in the western torus in the HCN *J* = 4–3 moment 0 map (02h 42m 40s7084, −00°00′ 47′′940) ICRS defined as W-peak. The horizontal thick black bar at the bottom center of each panel corresponds to 1 pc at the distance of NGC 1068. The open circle at the lower-left part of each panel indicates the synthesized beam size, which is 0′014 × 0′013 (≈1 pc) for the *J* = 4–3 lines of HCN and HCO⁺, and 0′020 × 0′018 for the *J* = 3–2 lines. An appropriate cutoff (~4.5σ) is applied to all the moment 1 maps to prevent them from being dominated by noise.</figcaption></figure>

at the systemic velocity in the western torus (Figures 4(a), (e))
than in the eastern torus (Figures 4(b), (f)). Gaussian fits
are applied to the detected HCN and HCO⁺ emission lines in
the area-integrated spectra (Appendix B). Table 4 summarizes the
estimated emission line fluxes.

The CS *J* = 7–6 (ν<sub>rest</sub> = 342.883 GHz) emission line is also
marginally detected in the NGC 1068 compact (≈3–5 pc) torus
region. The details are presented in Appendix C.

Figure 5 displays intensity-weighted mean velocity (moment 1) maps of the *J* = 4–3 lines of HCN and HCO+, together with those of the reanalyzed *J* = 3–2 lines of HCN and HCO+. Figures 5(a)–(d) highlight the redshifted (*V* ≳ 1330 km s⁻¹ or Δ*V*<sub>red</sub> ≳ +200 km s⁻¹) emission component in the innermost (≈1 pc) western torus relative to the systemic velocity (*V*<sub>sys</sub> = 1130 km s⁻¹), which was previously seen for HCN *J* = 3–2 (e.g., C. M. V. Impellizzeri et al. 2019; M. Imanishi et al. 2020) (see also Figure 5(c)). This redshifted (*V* ≳ 1330 km s⁻¹) component in the innermost (≈1 pc) western torus is clearly detected in the new HCN *J* = 4–3 data as well (Figure 5(a)).

The redshifted (ΔV<sub>red</sub> ≳ +200 km s<sup>-1</sup>) HCN emission
component at the innermost (≈1 pc) western torus can contribute
to the excess HCN *J* = 4–3 and HCN *J* = 3–2 emission tails on
the lower-frequency (redshifted) side of the main emission
components at the systemic velocity, which are clearly seen in
the area-integrated spectra at the western torus (Figure 4(a), (e)).
This redshifted emission component is not clearly seen in the
moment 1 maps of the HCO<sup>+</sup> *J* = 4–3 and HCO<sup>+</sup> *J* = 3–2 lines at
the innermost western torus (Figure 5(b), (d)). We will use the

HCN *J* = 4–3 line emission detected in our ≲1 pc resolution data to investigate the innermost dense molecular gas kinematics in the western torus.

VLBI observations at centimeter wavelengths with very high
angular resolution (≈2 mas) show that NGC 1068 has a rotating
22 GHz (1.4 cm) H₂O maser-emitting disk (e.g., L. J. Greenhill
et al. 1996; L. J. Greenhill & C. R. Gwinn 1997; J. F. Gallimore
et al. 2004; J. F. Gallimore & C. M. V. Impellizzeri 2023;
Y. Morishima et al. 2023; J. F. Gallimore et al. 2024). A larger
number of bright, redshifted 22 GHz (1.4 cm) H₂O maser spots
are detected at ~0.5–1.2 pc (~7–17 mas) on the northwestern
(NW) side of the central mass-accreting SMBH (the VLBI radio
22 GHz continuum emission peak, named “S1”) than on the
blueshifted southeastern (SE) side (e.g., L. J. Greenhill &
C. R. Gwinn 1997; J. F. Gallimore et al. 2004; J. F. Gallimore &
C. M. V. Impellizzeri 2023; Y. Morishima et al. 2023;
J. F. Gallimore et al. 2024). The most redshifted H₂O maser
emission has an intensity-weighted velocity of *V* ~ 1470 km s⁻¹
(Δ*V*<sub>red</sub>~ +340 km s⁻¹) at ~0.5 pc NW of the SMBH
(J. F. Gallimore & C. M. V. Impellizzeri 2023). Allowing for a
slight (~9 mas) positional shift between the ALMA ~350 GHz
and VLBI ~22 GHz data, possibly caused by astrometric
uncertainties in both observations (V. Gamez Rosas et al. 2025)<fnref n="8" />,
Figure 6 compares the velocity profiles of the VLBI-detected

<footnote n="8">⁸ According to the ALMA Cycle 9 Proposer’s Guide, the best absolute astrometric accuracy is ~10% of the synthesized beam, namely ≲2 mas for our HCN *J* = 4–3 data. V. Gamez Rosas et al. (2025) adopted an even larger positional shift of ~25 mas between their ALMA and VLBI data than ours (~9 mas).</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Intensity-weighted velocity of the HCN *J* = 4–3 (open circles) and 22 GHz (1.4 cm) H₂O maser emission (red filled triangles) (ordinate) as a function of distance from the SMBH position (continuum emission peak) along the NW direction (PA = 135°) (abscissa). We shift our ALMA HCN *J* = 4–3 data by ~9 mas (~0.65 pc) toward the NW direction (see also V. Gamez Rosas et al. 2025), assuming that the most highly redshifted (*V* ~ 1470 km s⁻¹ or Δ*V*<sub>red</sub> ~ +340 km s⁻¹) HCN *J* = 4–3 emission component is at ~0.5 pc from the VLBI 22 GHz continuum emission peak, where the 22 GHz H₂O maser emission also displays the most redshifted components with *V* ~ 1470 km s⁻¹ (Δ*V*<sub>red</sub> ~ +340 km s⁻¹) (J. F. Gallimore & C. M. V. Impellizzeri 2023). This is because the continuum emission peak positions often do not perfectly agree between ALMA at ~260–350 GHz and VLBI at ~22 GHz, possibly due to nonperfect astrometric accuracy. Red filled triangles: 22 GHz H₂O maser emission. Data in the R1, R2, and R3 regions, and those with >1450 km s⁻¹ in the R4 region (J. F. Gallimore & C. M. V. Impellizzeri 2023) are plotted. Distance in the abscissa is calculated from the offset values from the S1 in Table 1 of J. F. Gallimore & C. M. V. Impellizzeri (2023). Open circles: HCN *J* = 4–3 emission with 2 pixel width. The horizontal dotted straight line indicates the systemic velocity of *V*<sub>sys</sub> = 1130 km s⁻¹. The counterrotation signature, with a velocity change from redshifted at ≲2 pc to blueshifted at ≳2 pc, is seen along the NW direction.</figcaption></figure>

22 GHz H₂O maser-emitting spots and the ALMA-detected HCN *J* = 4–3 emission. In the HCN *J* = 4–3 data, the most redshifted component on the NW side has an intensity-weighted velocity as high as *V* ~ 1470 km s⁻¹ (Δ*V*<sub>red</sub> ~ +340 km s⁻¹) (Figure 6). The rough agreement of the velocity profiles of the redshifted 22 GHz H₂O maser and the HCN *J* = 4–3 emission at ~0.5–1.2 pc along the NW direction (after the slight positional shift) in Figure 6 suggests that the redshifted HCN *J* = 4–3 line-emitting dense molecular gas in the innermost (≪1 pc) western torus is physically connected to the VLBI-detected redshifted H₂O maser-emitting rotating disk.

Figures 5(e) and 5(f) present intensity-weighted mean velocity (moment 1) maps of the HCN *J* = 4–3 and HCN *J* = 3–2 emission using the same redshifted and blueshifted velocity ranges relative to the systemic velocity (*V*<sub>sys</sub> ± 110 km s<sup>-1</sup> and *V*<sub>sys</sub> ± 60 km s<sup>-1</sup>, respectively) to better visualize the overall rotation patterns of the compact (≪3–5 pc) HCN emission. In Figure 5(f), the HCN *J* = 3–2 emission in the western torus is redshifted in the inner part (≪2 pc) and blueshifted in the outer (≫2 pc) part, while in the eastern torus it is blueshifted in the inner part (≪2 pc) and redshifted in the outer (≫2 pc) part, i.e., they appear to be counterrotating (C. M. V. Impellizzeri et al. 2019; M. Imanishi et al. 2020). For the HCN *J* = 4–3 line, a similar apparent counterrotation is seen in the western torus (Figure 5(e)), but in the eastern torus only the inner (≪2 pc) blueshifted emission component is clearly visible, without a

significant outer (≫2 pc) redshifted emission component. This suggests that molecular gas at the outer (≫2 pc) eastern torus is not sufficiently warm and dense because the upper energy level (*E*<sub>upper</sub>) and critical density (*n*<sub>crit</sub>) for HCN *J* = 4–3 (*E*<sub>upper</sub> ~ 43 K and *n*<sub>crit</sub> ~ 1 × 10<sup>7</sup> cm<sup>-3</sup> at 50–100 K) are higher than those for HCN *J* = 3–2 (*E*<sub>upper</sub> ~ 26 K and *n*<sub>crit</sub> ~ 4–6 × 10<sup>6</sup> cm<sup>-3</sup>).

In Figure 6, the velocity changes from an inner redshift (≪2 pc) to an outer blueshift (≫2 pc) along the NW direction for HCN *J* = 4–3, but the outer blueshifted velocity is only Δ*V*<sub>blue</sub> ~ −40 km s<sup>-1</sup> relative to the systemic velocity of NGC 1068 (*V*<sub>sys</sub> = 1130 km s<sup>-1</sup>). This is much slower than the Keplerian rotation velocity of ~130–150 km s<sup>-1</sup> at ~2–2.5 pc for the central SMBH mass of NGC 1068 with *M*<sub>SMBH</sub> ~ 1 × 10<sup>7</sup> M<sub>☉</sub> (e.g., L. J. Greenhill et al. 1996; J.-M. Hure 2002; G. Lodato & G. Bertin 2003; J. F. Gallimore & C. M. V. Impellizzeri 2023; J. F. Gallimore et al. 2024), as has already been noted in previous studies (e.g., M. Imanishi et al. 2020; B. Vollmer et al. 2022; V. Gamez Rosas et al. 2025).

Figure 7 displays the position–velocity diagrams of the HCN *J* = 4–3 and HCO<sup>+</sup> *J* = 4–3 lines, as well as the reanalyzed HCN *J* = 3–2 and HCO<sup>+</sup> *J* = 3–2 lines, along the position angle PA = 114° east of north, following V. Gamez Rosas et al. (2025). Significant emission is seen in the redshifted component at *V* ~ 1330–1530 km s<sup>-1</sup> (Δ*V*<sub>red</sub> = +200–400 km s<sup>-1</sup>) in the western torus, as previously found (e.g., M. Imanishi et al. 2020; B. Vollmer et al. 2022; V. Gamez Rosas et al. 2025). Absorption components at *V* ~ 1430–1580 km s<sup>-1</sup> (Δ*V*<sub>red</sub> = +300–450 km s<sup>-1</sup>) toward the continuum emission peak (=SMBH position) are visible for the HCO<sup>+</sup> *J* = 4–3 and HCO<sup>+</sup> *J* = 3–2 lines (Figure 7(b), (d)), but not for the HCN lines at *J* = 4–3 or *J* = 3–2 (Figure 7(a), (c)).

Intensity-weighted velocity dispersion (moment 2) maps of the *J* = 4–3 lines of HCN and HCO<sup>+</sup>, together with the reanalyzed *J* = 3–2 lines of HCN and HCO<sup>+</sup>, are presented in Figure 8. The *J* = 4–3 data for HCN and HCO<sup>+</sup> clearly reveal high velocity dispersion regions in the western torus compared to the eastern torus, as previously found in the *J* = 3–2 data for HCN and HCO<sup>+</sup> (M. Imanishi et al. 2020). The HCN *J* = 4–3 emission peak “W-peak” spatially coincides with the high velocity dispersion regions in the western torus (Figure 8(a)–(d)), as previously seen for HCN *J* = 3–2 (M. Imanishi et al. 2020).

Figures 9(a) and 9(b) present the maps of (a) the HCN *J* = 4–3 to *J* = 3–2 and (b) the HCO<sup>+</sup> *J* = 4–3 to *J* = 3–2 flux ratios, respectively, after making their synthesized beams identical. The flux ratios tend to be higher at the inner part of the torus than at the outer part, confirming the argument in Section 3 (paragraph 1) that the bulk of the *J* = 4–3 emission comes from more compact regions than the *J* = 3–2 emission for both HCN and HCO<sup>+</sup>. If the emission is thermalized and optically thick, the expected flux ratio is 16/9 (=1.8) when calculated in units of Jy beam<sup>-1</sup> km s<sup>-1</sup>. The innermost (≪1 pc) regions of both the western and eastern tori and regions close to the W-peak show flux ratios with ≳2.5. Such high values cannot be explained solely by the possible absolute flux calibration uncertainty of individual ALMA observations (≪10%), which could produce the ratio with maximum ~2.2 (=1.1/0.9 × 1.8) for the thermalized, optically thick emission. This suggests that HCN and HCO<sup>+</sup> in these regions are highly excited to
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 7.** Position–velocity diagram of (a) HCN *J* = 4–3, (b) HCO⁺ *J* = 4–3, (c) HCN *J* = 3–2, and (d) HCO⁺ *J* = 3–2 lines, along the position angle of PA = 114° east of north. The abscissa is the offset from the ~260–350 GHz continuum emission peak (in pc). West is to the right and east is to the left. The ordinate is the optical LSR velocity (in km s⁻¹).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 8.** Intensity-weighted velocity dispersion (moment 2) map of (a) HCN *J* = 4–3, (b) HCO⁺ *J* = 4–3, (c) HCN *J* = 3–2, and (d) HCO⁺ *J* = 3–2, created with 1.5 mas pixel⁻¹. The displayed velocity range is 0–240 km s⁻¹ for *J* = 4–3 and 0–130 km s⁻¹ for *J* = 3–2. Note that the peak velocity dispersion value is significantly higher for the *J* = 4–3 lines than for the *J* = 3–2 lines, possibly because of the smaller synthesized beam size for *J* = 4–3, which can probe velocity structures at a finer spatial resolution. The horizontal thick bar, open circle, and plus and open diamond marks are the same as in Figure 5.</figcaption></figure>

*J* = 4 and the emission is optically thin (and/or affected by population inversion), likely influenced by the central AGN. The (c) HCN *J* = 4–3 to HCO⁺ *J* = 4–3 and (d) HCN *J* = 3–2 to HCO⁺ *J* = 3–2 flux ratios, also calculated in units of Jy

beam⁻¹ km s⁻¹, are displayed in Figures 9(c) and (d), respectively. HCN-to-HCO⁺ flux ratios greater than unity at both *J* = 4–3 and *J* = 3–2 are confirmed across most of the compact (≤3–5 pc) region of the NGC 1068 torus in a spatially resolved manner.
<figure><img src="image_15.png" /><figcaption>Figure 9. Map of flux ratios calculated in units of Jy beam⁻¹ km s⁻¹ for (a) HCN J = 4–3 to HCN J = 3–2, (b) HCO⁺ J = 4–3 to HCO⁺ J = 3–2, (c) HCN J = 4–3 to HCO⁺ J = 4–3, and (d) HCN J = 3–2 to HCO⁺ J = 3–2. Flux ratios are calculated after modifying the synthesized beam to a 0%021 circular one for (a), (b), and (d), and a 0%015 circular one for (c). The plus and open diamond marks are the same as in Figure 5.</figcaption></figure>

## 4. Discussion

### 4.1. Summary of Multiple Scenarios to Explain the Observed Dense Molecular Gas Properties in the NGC 1068 Torus

As briefly described in Section 1, multiple scenarios have been proposed to account for the apparent counterrotation in the compact (≤3–5 pc) NGC 1068 torus. The first is that the inner (≤2 pc) and outer (≥2 pc) parts of the NGC 1068 torus are indeed counterrotating with respect to each other (e.g., C. M. V. Impellizzeri et al. 2019; M. Imanishi et al. 2020), possibly because of a past collision with a very compact gas clump that had sufficient angular momentum and approached the western side of the originally rotating torus from the opposite direction (e.g., M. Imanishi et al. 2020; B. Vollmer et al. 2022). This configuration is dynamically unstable and short-lived (<10⁵ yr) (B. Vollmer et al. 2022). Thus, the probability of detecting this physically counterrotating phase is very small (e.g., S. Garcia-Burillo et al. 2019), but it may still be detectable (B. Vollmer et al. 2022). In this physically counterrotating gas scenario, angular momentum of the gas in the torus can be efficiently removed, enabling a high mass accretion rate onto the central SMBH (e.g., O. A. Kuznetsov et al. 1999; S. Dyda et al. 2015; D. Quach et al. 2015; E. I. Vorobyov et al. 2015), naturally explaining the fact that NGC 1068 is observed as a luminous AGN ($L_{AGN} \gtrsim 10^{45}$ erg s⁻¹) (e.g., J. J. Bock et al. 2000; S. F. Honig et al. 2008; D. Raban et al. 2009).

The second scenario is that the observed apparent counter-rotation signatures could be produced by outflow activity (e.g., S. Garcia-Burillo et al. 2019; D. Williamson et al. 2020;

difficult to reproduce solely with this outflow-origin model (V. Gamez Rosas et al. 2025).

V. Gamez Rosas et al. (2025) recently proposed a third scenario involving an infalling gas cloud from the western side to the immediate vicinity of the mass-accreting SMBH, without strong interaction with the torus. However, these authors focused only on (1) the observed asymmetry of molecular gas emission between the western and eastern tori and (2) absorption features detected for the CO J = 3–2 and HCO⁺ J = 4–3 lines at ~400 km s⁻¹ redshifted relative to the systemic velocity of NGC 1068 ($V_{sys} = 1130$ km s⁻¹), toward the ~350 GHz continuum emission peak (V. Gamez Rosas et al. 2025). The overall properties of the torus itself were not sufficiently discussed in the context of this third scenario.

We attempt to place further constraints on the dense molecular gas emission and dynamical properties of the NGC 1068 torus, particularly the western torus, by using our highest-spatial-resolution (≤1 pc) bright HCN J = 4–3 line data, combined with the ≤1 pc resolution HCO⁺ J = 4–3 and ~1.5 pc resolution HCN J = 3–2 and HCO⁺ J = 3–2 line data.

### 4.2. Highly Turbulent Dense Molecular Gas at the HCN J = 4–3 Emission Peak of the Torus

We estimate the emission line luminosities of HCN J = 4–3, HCO⁺ J = 4–3, HCN J = 3–2, and HCO⁺ J = 3–2 in the western torus (1.7 pc diameter around the W-peak), the eastern torus (1.7 pc diameter around the E-peak), and the combined torus region (6.8 pc diameter around the C-peak) from their fluxes (Table 4), using the following equation (P. M. Solomon & P. A. Vanden Bout 2005),

$$ \left( \frac{L'_{\text{line}}}{\text{K km s}^{-1} \text{ pc}^2} \right) = 3.25 \times 10^7 \left( \frac{\nu_{\text{rest}}}{\text{GHz}} \right)^{-2} (1+z)^{-1} \left( \frac{D_L}{\text{Mpc}} \right)^2 \left( \frac{S\Delta V}{\text{Jy km s}^{-1}} \right), \quad (1) $$

E. Y. Bannikova et al. 2023). This outflow-origin scenario is longer-lived and therefore more likely to be observed than the first physically counterrotating scenario in terms of timescale. However, it remains unclear what fundamental physical mechanism allows NGC 1068 to sustain its luminous AGN activity by removing the angular momentum of gas in the torus and enabling a high mass accretion rate onto the central SMBH. Furthermore, the significant asymmetry of molecular gas emission observed between the western and eastern tori is

where SΔV is the Gaussian-fit, velocity-integrated emission line flux (in Jy km s⁻¹) and D_L is the luminosity distance (in Mpc). Table 5 summarizes the derived molecular emission line luminosities.

Because the critical densities of the HCN and HCO⁺ lines at J = 4–3 and J = 3–2 are $n_{crit} \gtrsim 10^6$ cm⁻³ at a gas kinetic temperature of 50–100 K (Y. L. Shirley 2015), only dense and warm molecular gas can sufficiently excite HCN and HCO⁺ to J = 4 and J = 3. For most molecular gas in luminous starburst
<table><caption>**Table 5**
Molecular Emission Line Luminosity</caption><thead><tr><th>Torus Position<br>(1)</th><th>HCN <i>J</i> = 4–3<br>10<sup>4</sup> (K km s<sup>-1</sup> pc<sup>2</sup>)<br>(2)</th><th>HCO<sup>+</sup> <i>J</i> = 4–3<br>10<sup>4</sup> (K km s<sup>-1</sup> pc<sup>2</sup>)<br>(3)</th><th>HCN <i>J</i> = 3–2<br>10<sup>4</sup> (K km s<sup>-1</sup> pc<sup>2</sup>)<br>(4)</th><th>HCO<sup>+</sup> <i>J</i> = 3–2<br>10<sup>4</sup> (K km s<sup>-1</sup> pc<sup>2</sup>)<br>(5)</th></tr></thead><tbody><tr><td>W-torus (1.7 pc diameter)</td><td>4.8 ± 0.5</td><td>3.0 ± 0.5</td><td>2.8 ± 0.6</td><td>1.6 ± 0.2</td></tr><tr><td>E-torus (1.7 pc diameter)</td><td>3.7 ± 0.3</td><td>2.3 ± 0.3</td><td>1.5 ± 0.2</td><td>1.2 ± 0.2</td></tr><tr><td>Torus (6.8 pc diameter)</td><td>21.7 ± 1.5</td><td>12.5 ± 1.0</td><td>18.0 ± 0.9</td><td>12.4 ± 0.9</td></tr></tbody></table>

**Note.** Column (1): Torus region. Columns (2)–(5): Emission line luminosities of (2) HCN *J* = 4–3, (3) HCO⁺ *J* = 4–3, (4) HCN *J* = 3–2, and (5) HCO⁺ *J* = 3–2, in units of 10⁴ K km s⁻¹ pc².

galaxies and AGNs, HCN and HCO⁺ luminosities in units of K km s⁻¹ pc² at *J* = 4–3 are usually smaller than those at *J* = 3–2 if these molecules are subthermally excited and the emission is optically thick (e.g., K. K. Knudsen et al. 2007; M. Krips et al. 2008; T. R. Greve et al. 2009; P. P. Papadopoulos et al. 2014; F. P. Israel 2023; M. Imanishi et al. 2023a, 2023b). The luminosities in these units can be comparable between *J* = 4–3 and *J* = 3–2 if HCN and HCO⁺ are thermally excited and the emission is optically thick. In Table 5, the *J* = 4–3 and *J* = 3–2 luminosities in the 6.8 pc diameter circular region centered on the mass-accreting SMBH position are roughly comparable, within ≲20%, suggesting that HCN and HCO⁺ are nearly thermally excited to *J* = 4 and *J* = 3, and that the emission is optically thick when integrating molecular gas emission over the entire NGC 1068 compact (≪3–5 pc) torus region.

However, in both the western and eastern tori with 1.7 pc diameter circular regions around the HCN *J* = 4–3 emission peaks, the luminosities of both HCN and HCO⁺ are a factor of ≳1.7 higher<fnref n="9" /> at *J* = 4–3 than at *J* = 3–2 (Table 5), suggesting that HCN and HCO⁺ are nearly thermally excited and the emission is optically thin in the small areas around the W-peak and E-peak. Large velocity dispersion caused by high turbulence can reduce molecular line opacity and reproduce the observed HCN and HCO⁺ emission line properties in these regions. The first physically counterrotating torus scenario, caused by a physical collision of a compact dense gas clump with the western side of the rotating dense molecular torus from the opposite direction (e.g., M. Imanishi et al. 2020; B. Vollmer et al. 2022), can naturally generate high turbulence in the torus, particularly in small regions around the collision interface, originating from the Kelvin-Helmholtz instabilities (e.g., S. Dyda et al. 2015; D. Quach et al. 2015). The spatial coincidence between the observed large velocity dispersion regions and the W-peak (Figure 8(a)–(d)) can be explained if the putative collision interface in the western torus is close to the W-peak, where the observed velocity dispersion values can increase through the combination of (i) enhancement of intrinsic velocity dispersion produced by the turbulence and (ii) beam smearing of two counterrotating gas components. However, the third scenario of an infalling gas cloud from the west to the SMBH vicinity, without physical collision with the torus (V. Gamez Rosas et al. 2025), has no direct way to reduce the opacity of molecular line emission in the small areas around the W-peak (and E-peak), which is suggested by our new observations.

The physical connection between the redshifted 22 GHz H₂O maser and the HCN *J* = 4–3 dense molecular line emission at the innermost western torus, suggested by their comparable velocity profiles (Figure 6), also supports the first physically counterrotating dense molecular gas scenario (particularly in the western torus) (M. Imanishi et al. 2020), but cannot be naturally explained by the third scenario, in which the redshifted molecular line emission from the innermost western torus is produced by an infalling gas cloud physically unrelated to the H₂O maser-emitting rotating disk (V. Gamez Rosas et al. 2025), as well as outflow-origin models that do not invoke a rotating dense molecular torus (D. Williamson et al. 2020; E. Y. Bannikova et al. 2023). We thus further explore the first collision-induced physically counterrotating torus molecular gas scenario.

### 4.3. Comparison of the Physically Counterrotating Torus Scenario of B. Vollmer et al. (2022) with the Highest-spatial-resolution Dense Molecular Gas Data

In the collision-induced physically counterrotating torus scenario, two more observational characteristics may be produced at the interface between the two counterrotating gas components, particularly in the western torus, in addition to the observed velocity dispersion increase (Section 4.2). First, a possible 0.5–1 pc wide gap may be created (e.g., B. Vollmer et al. 2022). However, we see no clear signature of such a gap, even in our highest-spatial-resolution (≪1 pc) integrated-intensity (moment 0) map of the bright HCN *J* = 4–3 emission (Figure 1(a)). Second, excitation of dense molecular gas can be enhanced by turbulence-induced shocks. We do not see any clear sign of elevated *J* = 4–3 to *J* = 3–2 flux ratios of HCN and HCO⁺ emission, driven by the putative shocks, at the W-peak (Figure 9).

In summary, even using our new ≲1 pc resolution bright HCN *J* = 4–3 line data, we cannot provide definitive observational evidence to conclusively support the first collision-induced physically counterrotating torus scenario, possibly because the achieved ≲1 pc resolution is still insufficient to investigate the predicted collision interface between two counterrotating gas components in detail. Even higher-spatial-resolution (≪1 pc) data may help better scrutinize the putative interface regions.

Even though we cannot clearly detect the putative collision interface, our new ≲1 pc (≪14 mas)-resolution HCN *J* = 4–3 and HCO⁺ *J* = 4–3 line data can be used to update the collision-induced physically counterrotating torus model of B. Vollmer et al. (2022). In particular, B. Vollmer et al. (2022) compared the predictions of their analytical model with the

<footnote n="9">⁹ The possible absolute flux calibration uncertainty of individual ALMA observations (≪10%) again could increase the *J* = 4–3 to *J* = 3–2 flux ratio by a factor of only ≲1.2 (=1.1/0.9).</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 10.** Comparison of position–velocity diagrams between observational data and the model by B. Vollmer et al. (2022). (a): HCN *J* = 4–3. (b): HCO<sup>+</sup> *J* = 4–3. In these plots, the major axis is PA = 110°, consistent with those shown in B. Vollmer et al. (2022). For the observations (left), the images are rotated by 20° to make the elongated structure horizontal (E–W). For the model (right), the data are projected such that the disk plane is horizontal (E–W). In each panel, the curved dotted line indicates the expected Keplerian rotation for the radial distribution of the total enclosed mass, *M*(*R*) = *M*<sub>SMBH</sub> + *M*<sub>0</sub>*R*<sup>5/4</sup>, where *M*<sub>SMBH</sub> is the central SMBH mass (*M*<sub>SMBH</sub> = 1 × 10<sup>7</sup>*M*<sub>☉</sub>) and *M*<sub>0</sub> is the stellar mass distribution (*M*<sub>0</sub> = 0.8 × 10<sup>6</sup>*M*<sub>☉</sub> (pc<sup>-5/4</sup>)), in the same way as presented in B. Vollmer et al. (2022).</figcaption></figure>

<figure><img src="image_17.png" /></figure>

ALMA molecular line data available at that time, including multiple CO J-transition line data and ~1.5 pc (~20 mas) resolution HCN *J* = 3–2 and HCO⁺ *J* = 3–2 line data, and concluded that the observed overall molecular line properties can be quantitatively reproduced by the model within a factor of ~2. However, only low-spatial-resolution HCN *J* = 4–3 (≥7 pc or ≥100 mas) and HCO⁺ *J* = 4–3 (≥2 pc or ≥30 mas) data were used. For HCN *J* = 4–3, no detailed comparison was possible. For HCO⁺ *J* = 4–3, an absorption feature toward the continuum emission peak at the systemic velocity (*V*<sub>sys</sub> = 1130 km s<sup>-1</sup>) was predicted by the model but was not seen in the observed data (Figure 20 of B. Vollmer et al. 2022), most likely because it was diluted by emission in the low-spatial-resolution data. We now have new higher-spatial-resolution (≤1 pc or 14 mas) HCN *J* = 4–3 and HCO⁺ *J* = 4–3 line data, enabling an updated comparison with the model. Indeed, the predicted absorption feature at the systemic velocity (*V*<sub>sys</sub> = 1130 km s<sup>-1</sup>) toward the continuum emission peak was detected in our beam-sized spectra (Figure 3) and position-velocity diagrams (Figure 7) of the HCN *J* = 4–3 and HCO⁺ *J* = 4–3 lines.

Figure 10 shows a detailed comparison of the position–velocity diagrams of the HCN *J* = 4–3 and HCO⁺ *J* = 4–3 lines between our new ≲1 pc resolution data and the best-fit model of B. Vollmer et al. (2022). The main features of the model’s prediction (i.e., absorption at the systemic velocity) are confirmed by our new ≲1 pc resolution HCN *J* = 4–3 and HCO⁺ *J* = 4–3 line data. Observed and predicted spectra at multiple positions along the major axis (PA = 110°) are compared in Figure 11. The overall properties largely agree between the model and observations, as in the case for other molecular lines (B. Vollmer et al. 2022). The observed east–west asymmetry of dense molecular line emission appears not to be adequately reproduced by the best-fit model of B. Vollmer et al. (2022), but this discrepancy largely comes from initial conditions of the modeling. Further model refinement is important because this collision model, consisting of infalling gas clump to the western torus should be able to create significant east–west asymmetry, which is very difficult to naturally explain by the outflow-induced model (e.g., S. Garcia-Burillo et al. 2019; D. Williamson et al. 2020; E. Y. Bannikova et al. 2023).

## 4.4. Further Confirmation of Enhanced HCN Abundance in the Compact Nuclear Region of NGC 1068

### 4.4.1. From Nuclear Outflow Activity

Nuclear dense molecular outflow activity has been proposed based on the detection of a blueshifted HCN *J* = 3–2 absorption wing with up to Δ*V*<sub>blue</sub> ~ −400 km s<sup>-1</sup> in the spectrum toward the ~260 GHz continuum emission peak (the putative location of the mass-accreting SMBH) (C. M. V. Impellizzeri et al. 2019; M. Imanishi et al. 2020). This blueshifted absorption wing was less clear for HCO⁺ *J* = 3–2 (M. Imanishi et al. 2020) and was reported to be undetected for HCO⁺ *J* = 4–3 (V. Gamez Rosas et al. 2025). The blueshifted absorption wing reaching Δ*V*<sub>blue</sub>~−400 km s<sup>-1</sup> (*V* ~ 730 km s<sup>-1</sup>) is recognizable for HCN *J* = 4–3 (Figure 3(a)) but appears weaker for HCO⁺ *J* = 4–3 (Figure 3(b)). The stronger observed outflow-origin absorption wings for HCN than those for HCO⁺ can be explained by the enhanced HCN abundance relative to HCO⁺ found in the compact (≈3–5 pc) nuclear region around the luminous AGN in NGC 1068 (e.g., M. Imanishi et al. 2018b; M. Imanishi et al. 2020; J. Butterworth et al. 2022; B. Vollmer et al. 2022).

### 4.4.2. From Torus Dense Molecular H₂ Gas Mass Estimate

We estimate the HCN *J* = 4–3 emission line flux within the 6.8 pc diameter aperture around the central continuum-emitting AGN to be 4.3 ± 0.3 Jy km s⁻¹ (Table 4). This recovers ≳80% of the integrated HCN *J* = 4–3 flux from the torus (5.3 Jy km s⁻¹) reported in Table 4 of B. Vollmer et al. (2022). The HCO⁺ *J* = 4–3 emission line flux within the same 6.8 pc aperture is estimated to be 2.5 ± 0.2 Jy km s⁻¹ (Table 4), again recovering ≳80% of the integrated HCO⁺ *J* = 4–3 flux (3.0 Jy km s⁻¹) in Table 4 of B. Vollmer et al. (2022). Thus, most of the *J* = 4–3 emission of HCN and HCO⁺ from the entire NGC 1068 torus originates in the compact (≈3–5 pc) region.

Dense molecular hydrogen (H₂) mass is often estimated from (optically thick) HCN *J* = 1–0 line luminosity using the formula *M*<sub>dense-H₂</sub>(HCN) = 6–12 × L(HCN *J* = 1–0) (*M*<sub>☉</sub> [K km s<sup>-1</sup> pc<sup>2</sup>]<sup>-1</sup>) (e.g., Y. Gao & P. M. Solomon 2004; M. Krips et al. 2008; A. K. Leroy et al. 2017). On the compact torus scale (6.8 pc diameter aperture), *J* = 4–3 and *J* = 3–2 emission line luminosities in units of (K km s<sup>-1</sup> pc<sup>2</sup>) are
<figure><img src="image_18.png" /><figcaption>**Figure 11.** Comparison between observational (left) and model (right) spectra by B. Vollmer et al. (2022), at multiple torus positions. (a): HCN *J* = 4–3. (b): HCO<sup>+</sup> *J* = 4–3. In these plots, NW and SE are along PA = 110°, to be consistent with those shown in B. Vollmer et al. (2022). Spectra are extracted along the NW and SE directions, at distances of 0 pc, ±1.4 pc, ±2.8 pc, and ±4.2 pc, with 1.4 pc diameter circular apertures.</figcaption></figure>

roughly comparable for both HCN and HCO⁺ (Table 5), suggesting that they are thermalized and their emission are optically thick. Assuming that the HCN *J* = 1–0 luminosity in these units is also comparable to that of HCN *J* = 3–2 and *J* = 4–3 (~20 × 10⁴; Table 5), we obtain HCN-derived torus dense molecular H₂ gas mass of *M*<sub>dense-H₂</sub> (HCN) = (1.2–2.4) × 10⁶ *M*<sub>☉</sub>.

Dense H₂ mass is also estimated from HCO⁺ *J* = 1–0 line luminosity using the formula *M*<sub>dense-H₂</sub>(HCO⁺) = 2–5 × *L*(HCO⁺ *J* = 1–0) (*M*<sub>☉</sub> [K km s⁻¹ pc²]⁻¹) in the optically

thick regime (e.g., A. K. Leroy et al. 2017). Assuming that the HCO⁺ *J* = 1–0 luminosity in these units is comparable to that of HCO⁺ *J* = 4–3 and *J* = 3–2 (~12 × 10⁴; Table 5), we obtain HCO⁺-derived torus dense H₂ mass of *M*<sub>dense-H₂</sub> (HCO⁺) = (2.4–6.0) × 10⁵ *M*<sub>☉</sub>. This is a factor of ~4–5 smaller than the HCN-derived dense H₂ mass and can be explained by the enhanced HCN-abundance scenario in the compact torus region of NGC 1068 (e.g., M. Imanishi et al. 2018b; M. Imanishi et al. 2020; J. Butterworth et al. 2022), because standard HCN and HCO⁺ abundances are assumed for
the above conversion from line luminosity to dense H₂ mass. The HCO⁺-derived torus dense H₂ mass is roughly comparable to the value adopted (3.6 × 10⁵ M☉) in the modeling of B. Vollmer et al. (2022) (their Table A.1).

## 5. Summary

We presented the results of our new ALMA ≲1 pc (≪14 mas) resolution HCN *J* = 4–3 and HCO⁺ *J* = 4–3 line observations of NGC 1068. The achieved spatial resolution at *J* = 4–3 is ≳30% better than that at *J* = 3–2 because of its shorter wavelength (higher frequency). By combining these new ≲1 pc resolution *J* = 4–3 data with existing and reanalyzed ~1.5 pc resolution HCN *J* = 3–2 and HCO⁺ *J* = 3–2 data, we obtained the following main results.

1. HCN *J* = 4–3 and HCO⁺ *J* = 4–3 dense molecular line emission is clearly detected in the compact (≪3–5 pc) nuclear region along the almost east–west direction from the ~350 GHz continuum emission peak. We interpret this molecular emission as originating from the putative dense molecular torus around a mass-accreting SMBH.

2. In both the western and eastern tori, HCN *J* = 4–3 emission is brighter than HCO⁺ *J* = 4–3 emission, as previously seen at *J* = 3–2. The western torus is brighter than the eastern torus in both HCN *J* = 4–3 and HCO⁺ *J* = 4–3 emission lines.

3. In the western torus, the inner (≪2 pc) part is redshifted and the outer (≫2 pc) part is blueshifted, with respect to the systemic velocity of NGC 1068 ($V_{\text{sys}}$ = 1130 km s⁻¹), for both HCN *J* = 4–3 and HCO⁺ *J* = 4–3 emission, confirming apparent counterrotation as previously found for HCN *J* = 3–2 and HCO⁺ *J* = 3–2 lines.

4. In the eastern torus, both HCN *J* = 4–3 and HCO⁺ *J* = 4–3 emission are blueshifted at the inner (≪2 pc) part, as previously identified at *J* = 3–2 of HCN and HCO⁺. However, a redshifted emission component at the outer (≫2 pc) part, previously detected in HCN *J* = 3–2 and CO *J* = 3–2 lines, is not clearly seen in either HCN *J* = 4–3 or HCO⁺ *J* = 4–3, suggesting that the outer redshifted molecular gas in the eastern torus is not sufficiently warm and dense.

5. The spatial and dynamical properties of the innermost (≪1 pc) redshifted HCN *J* = 4–3 emission in the western torus largely agree with those of the even inner (0.5–1.2 pc) 22 GHz H₂O maser emission spots, previously identified with VLBI very-high-angular-resolution (≪2 mas) observations. This suggests that in the innermost western torus, the highly redshifted HCN *J* = 4–3 dense molecular emission is physically related to the redshifted H₂O maser-emitting rotating disk.

6. Comparison of *J* = 4–3 and *J* = 3–2 fluxes of HCN and HCO⁺ emission extracted with a 1.7 pc diameter circular aperture around the HCN *J* = 4–3 emission peaks in both the western and eastern tori (“W-peak” and “E-peak”, respectively) suggests that dense molecular line emission is optically thin in these small areas. High turbulence of molecular gas can reduce line opacity and reproduce the observed molecular emission line properties there. A physically counterrotating dense molecular gas scenario in the NGC 1068 torus, driven by the collision of a

compact gas clump with the western part of the existing rotating torus from the opposite direction, is a natural explanation if the interface of two counterrotating gas components is located close to these emission peaks.

7. We applied the analytical model based on this assumption presented by B. Vollmer et al. (2022) and found that the overall properties of our new HCN *J* = 4–3 and HCO⁺ *J* = 4–3 data are largely reproduced by this model, as previously confirmed for other molecular lines, including HCN *J* = 3–2 and HCO⁺ *J* = 3–2. However, the observed significant east–west asymmetry of dense molecular line emission properties needs to be better reproduced with this model. Furthermore, even using our highest spatial resolution (≪1 pc) molecular line data so far, (i) a predicted possible ~0.5–1.0 pc wide gap at the interface between two counterrotating gas components and (ii) high molecular gas excitation caused by possible turbulence-induced shocks there are not clearly detected. Even smaller beam-sized data in the future may help further test the physically counterrotating torus scenario in NGC 1068.

8. Enhanced HCN abundance in the vicinity (≪3–5 pc) of the luminous AGN in NGC 1068 is further supported by two additional observational facts: (i) the blueshifted absorption wing toward the mass-accreting SMBH, likely originating from nuclear compact outflow activity, is stronger for HCN *J* = 4–3 than for HCO⁺ *J* = 4–3, and (ii) the HCN-derived torus dense molecular H₂ gas mass is a factor of ~4–5 higher than the HCO⁺-derived one, under the assumption of standard HCN and HCO⁺ abundances for converting their luminosities to dense molecular H₂ gas mass.

Our study demonstrated that the combination of high-spatial-resolution, multiple J-transition molecular line data is a very powerful tool to scrutinize various properties of the compact dense molecular torus in the nearby archetypal AGN NGC 1068.

## Acknowledgments

We thank the anonymous referee for their valuable comments, which helped improve the clarity of this manuscript. This paper made use of the following ALMA data: ADS/JAO.ALMA#2022.1.00005.S and #2018.1.00037.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC 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. M.I. is supported by JP21K03632 and JP25K07359. Y. H. is supported by JP25K07357. Data analysis was in part carried out on the open use data analysis computer system at the Astronomy Data Center, ADC, of the National Astronomical Observatory of Japan. This research has made use of NASA’s Astrophysics Data System and the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

*Facility:* ALMA.

*Software:* CASA (The CASA Team et al. 2022).
## Appendix A

### Physical Origin of the ~400 km s⁻¹ Redshifted Absorption Feature Seen in the HCO⁺ *J* = 4–3 Spectrum

An absorption feature, at ~400 km s⁻¹ redshifted relative to the systemic velocity (V<sub>sys</sub> = 1130 km s⁻¹), is detected in the spectra of HCO⁺ *J* = 4–3 (Figure 3(b) and V. Gamez Rosas et al. 2025), CO *J* = 3–2 (V. Gamez Rosas et al. 2025), and HCO⁺ *J* = 3–2 (M. Imanishi et al. 2020), toward the ~260–350 GHz continuum emission peak position. It is also seen in the position–velocity diagrams of HCO⁺ *J* = 4–3, HCO⁺ *J* = 3–2, and CO *J* = 3–2, at the continuum emission peak (Figure 7(b),(d) and V. Gamez Rosas et al. (2025)), but not seen for HCN *J* = 4–3 and HCN *J* = 3–2 (Figure 7(a),(c)).

The frequency of the ~400 km s⁻¹ redshifted absorption feature seen in the spectra of HCO⁺ *J* = 4–3 and *J* = 3–2 agrees with that of HCN-VIB *J* = 4–3 and *J* = 3–2 lines, respectively, because the rest frequency of the HCN-VIB line is ~400 km s⁻¹ redshifted relative to that of HCO⁺ (vibrational ground level, *v* = 0), both at *J* = 4–3 and *J* = 3–2 (Figure 2(a) and M. Imanishi et al. 2020). If foreground absorbing dense molecular gas at the systemic velocity (V<sub>sys</sub> = 1130 km s⁻¹) is located close to the central continuum-emitting AGN, HCN-VIB *J* = 4–3 and *J* = 3–2 lines can be detected in absorption, given that in the close vicinity (~a few pc) of the AGN in NGC 1068, (1) HCN abundance is derived to be high (e.g., M. Imanishi et al. 2018b; M. Imanishi et al. 2020; J. Butterworth et al. 2022; B. Vollmer et al. 2022) and (2) a sufficient amount of HCN can be vibrationally excited (upper energy level *E*<sub>upper</sub> = 1024 K) by AGN-origin infrared radiation pumping (K. Sakamoto et al. 2010). In fact, the presence of a sufficient amount of vibrationally excited HCN molecules (HCN-VIB) has been found in the vicinity of luminous AGNs (e.g., K. Sakamoto et al. 2010; M. Imanishi & K. Nakanishi 2013; S. Aalto et al. 2015a, 2015b; F. Costagliola et al. 2015; M. Imanishi et al. 2016b, 2016c, 2018a; S. Martin et al. 2016; N. Falstad et al. 2019, 2021; K. Sakamoto et al. 2021). The possible absorption dip by the HCN-VIB *J* = 4–3 line is ~0.5 mJy beam⁻¹ (Figures 2(a) and 3(b)) or ~7% (~0.08) for the ~350 GHz continuum flux of ~6.8 mJy beam⁻¹ (Table 2). A comparable level of absorption dip (~0.08) by the HCN-VIB *J* = 4–3 line has been detected toward the ~350 GHz continuum-emitting mass-accreting SMBH in the luminous obscured AGN NGC 1052 (S. Kameno et al. 2020). The possible absorption dip by the HCN-VIB *J* = 3–2 line toward the ~260 GHz continuum-emitting mass-accreting SMBH in NGC 1068 is also ~9% (~0.09) (M. Imanishi et al. 2020). This HCN-VIB absorption scenario can naturally explain the nondetection of ~400 km s⁻¹ redshifted absorption features for HCN both at *J* = 4–3 and *J* = 3–2 (Figure 7(a),(c)) because there are no corresponding HCN-VIB lines there. We thus argue that the absorption features centered at *V* ~ 1530 km s⁻¹ (Δ*V*<sub>red</sub> ~ +400 km s⁻¹) detected for HCO⁺ *J* = 4–3 and *J* = 3–2 in the beam-sized spectra toward the ~260–350 GHz continuum emission peak (Figure 3(b) and M. Imanishi et al. 2020) and position–velocity diagrams (Figure 7(b), (d)) can be explained by HCN-VIB *J* = 4–3 and *J* = 3–2 lines at the systemic velocity, respectively.

The spatial extent of the ~260–350 GHz continuum emission is ~1–2 pc (Figure 1 and V. Gamez Rosas et al. 2025) and the Keplerian rotation velocity at ~1 pc from the SMBH is ~200 km s⁻¹ for the SMBH mass of *M*<sub>SMBH</sub> ~ 1 × 10⁷ *M*<sub>☉</sub>,

estimated for NGC 1068 (e.g., L. J. Greenhill et al. 1996; J.-M. Hure 2002; G. Lodato & G. Bertin 2003; J. F. Gallimore & C. M. V. Impellizzeri 2023; J. F. Gallimore et al. 2024). Assuming that HCN molecules at ~1 pc from the central AGN are sufficiently vibrationally excited to *v* = 1 level (HCN-VIB), then the observed line width of ~300–400 km s⁻¹ for the ~400 km s⁻¹ redshifted HCO⁺ *J* = 4–3 and *J* = 3–2 absorption features (Figure 3(c) and M. Imanishi et al. 2020) could be explained by the HCN-VIB absorption scenario.

V. Gamez Rosas et al. (2025) interpreted that the ~400 km s⁻¹ redshifted HCO⁺ *J* = 4–3 absorption feature originates in redshifted (~400 km s⁻¹) infalling gas in the close vicinity of the central mass-accreting SMBH because the ~400 km s⁻¹ redshifted absorption feature is seen also for CO *J* = 3–2. In this scenario, both the (i) detection of the absorption feature for HCO⁺ *J* = 3–2 (Figure 7(d)) and (ii) nondetection for HCN *J* = 4–3 and HCN *J* = 3–2 (Figure 7(a), (c)) need to be explained. At the same *J* transitions, the upper and lower energy levels are comparable between HCN and HCO⁺, but the critical density (*n*<sub>crit</sub>) for HCN is a factor of ~5 higher than that for HCO⁺ (Y. L. Shirley 2015). Thus, as long as we compare the same *J* = 4–3 or *J* = 3–2, stronger absorption optical depth for HCO⁺ than for HCN could be explained by molecular gas whose density is high enough to excite HCO⁺ to *J* = 3 and *J* = 2 but not high enough to sufficiently excite HCN to *J* = 3 and *J* = 2. However, it is not immediately obvious whether this scenario can explain the observed significantly stronger absorption optical depth for HCO⁺ *J* = 4–3 (lower energy level *E*<sub>lower</sub> ~ 26 K and *n*<sub>crit</sub> = 2–3 × 10⁶ cm⁻³ at 50–100 K) than for HCN *J* = 3–2 (*E*<sub>lower</sub> ~ 13 K and *n*<sub>crit</sub> = 4–6 × 10⁶ cm⁻³) (Y. L. Shirley 2015).

We ran RADEX calculations (F. F. S. van der Tak et al. 2007) for a foreground absorbing gas cloud with a volume density of 10⁵⁻⁶ cm⁻³, kinetic temperature of 50–150 K, and column density of 3 × 10¹⁵ cm⁻² for HCO⁺ (V. Gamez Rosas et al. 2025) and HCN, and line width of 300 km s⁻¹ (Figure 3(b)). If we adopt a background temperature of ~3 K (i.e., the cosmic microwave background temperature in the nearby universe), we find no solution to reproduce a factor of ~2 stronger absorption optical depth for HCO⁺ *J* = 4–3 than for HCN *J* = 3–2.

The continuum flux density at ~350 GHz and ~260 GHz is 6.8 mJy beam⁻¹ (for a synthesized beam of 0′′014 × 0′′013) and 6.7 mJy beam⁻¹ (0′′019 × 0′′017), respectively (Table 2). We obtain a brightness temperature of ~350 K from these continuum data. By inputting a background temperature of 350 K for our RADEX calculations, we find that the absorption optical depth for HCO⁺ *J* = 4–3 can be a factor of ~2 stronger than that for HCN *J* = 3–2, for a volume density of 10⁵⁻⁶ cm⁻³, kinetic temperature of 50–150 K, and column density of 3 × 10¹⁵ cm⁻² both for HCO⁺ and HCN, and a line width of 300 km s⁻¹. For the kinetic temperature of 50–80 K, the absorption optical depths of both HCO⁺ *J* = 4–3 and *J* = 3–2 can be a factor of ~2 stronger than those of both HCN *J* = 4–3 and *J* = 3–2. Thus, the observed detection of HCO⁺ *J* = 4–3 and *J* = 3–2 absorption features, and nondetection of HCN *J* = 4–3 and *J* = 3–2 absorption features, can be quantitatively reproduced by an infalling gas cloud with a ~400 km s⁻¹ redshifted velocity (V. Gamez Rosas et al. 2025).

The above comparison implicitly assumes that the HCN and HCO⁺ abundances are comparable. In the close vicinity of the AGN in NGC 1068, the HCN abundance is derived to be
significantly enhanced relative to HCO$^+$ (e.g., M. Imanishi et al. 2018b; M. Imanishi et al. 2020; J. Butterworth et al. 2022; B. Vollmer et al. 2022), largely affected by the central AGN. If a similarly enhanced HCN abundance is assumed for the infalling gas cloud, then the stronger HCO$^+$ *J* = 4–3 absorption compared to HCN *J* = 3–2 absorption becomes more difficult to explain. However, it may be possible that the infalling gas cloud is not HCN-abundance-enhanced, despite being located very close to the central AGN, if the infalling cloud comes from a region sufficiently far from the NGC 1068 nucleus and stays in the AGN vicinity only for a short time period (V. Gamez Rosas et al. 2025). In this case, possible AGN effects on the abundance of the infalling gas cloud are limited. However, the elevated HCN abundance in NGC 1068 is reported not only in the compact ($\lesssim$3–5 pc) torus region but also in a more extended ($\sim$100 pc) region (e.g., A. Sternberg et al. 1994; T. Nakajima et al. 2018; S. Takano et al. 2019; J. Butterworth et al. 2022; T. Nakajima et al. 2023). If the $\sim$400 km s$^{-1}$ redshifted HCO$^+$ *J* = 4–3 and *J* = 3–2 absorption features are due to an infalling HCO$^+$ gas clump, then it is likely to come from a non-HCN-abundance-enhanced region, far away ($\gtrsim$100 pc) from the central HCN-abundance-enhanced region of NGC 1068.

In summary, both the HCN-VIB and infalling gas cloud scenarios could reproduce the $\sim$400 km s$^{-1}$ redshifted absorption features detected only for HCO$^+$ *J* = 4–3 and *J* = 3–2, but not for HCN *J* = 4–3 and *J* = 3–2. The former scenario cannot immediately explain the $\sim$400 km s$^{-1}$ redshifted CO *J* = 3–2 absorption feature (V. Gamez Rosas et al. 2025) because there is no HCN-VIB absorption line there. It may come from a low-density ($\lesssim$10$^4$ cm$^{-3}$) infalling gas cloud. The latter scenario is feasible as long as the infalling gas cloud is not HCN-abundance-enhanced.

If HCN-VIB lines are responsible for the $\sim$400 km s$^{-1}$ redshifted absorption features detected for HCO$^+$ *J* = 4–3 and

*J* = 3–2, then HCN-VIB lines from dense molecular gas, which is located in the close vicinity of the central AGN in a direction perpendicular to our line of sight and is strongly illuminated by the AGN radiation, are expected to be detected in emission. However, the detection of the HCN-VIB emission line is usually not easy because it is swamped by the much brighter HCO$^+$ emission line at the vibrational ground (ν = 0) level, which is separated only by $\sim$400 km s$^{-1}$. This is the case in the beam-sized and area-integrated spectra of the W-torus (Figures 2(b) and 4(a)). The HCO$^+$ *J* = 4–3 emission line profile is narrower in the E-torus (Figures 2(c) and 4(b)), but the presence of the HCN-VIB *J* = 4–3 emission line is not clear.

In Figure 3(c), an emission-like feature is recognizable at $\sim$200–400 km s$^{-1}$ redshifted velocity ($V_{\text{sys}}$ + 200–400 km s$^{-1}$ or *V* = 1330–1530 km s$^{-1}$), in the beam-sized HCN-VIB *J* = 4–3 spectrum toward the $\sim$350 GHz continuum emission peak (=position of the putative mass-accreting SMBH). This spectrum is likely to include emission from the innermost ($\lesssim$1 pc) western torus where redshifted $\sim$200–400 km s$^{-1}$ emission components are detected in the intensity-weighted mean velocity (moment 1) maps of HCN *J* = 4–3 (Figure 5(a)) and HCN *J* = 3–2 (Figure 5(c) and M. Imanishi et al. 2020). The redshifted emission-like feature in Figure 3(c) may originate from redshifted HCN-VIB *J* = 4–3 emission components at the innermost ($\lesssim$1 pc) western torus.

## Appendix B: Gaussian Fit of Detected Emission Lines

Figure 12 displays the Gaussian fitting results for the detected HCN *J* = 4–3 and HCO$^+$ *J* = 4–3 emission lines in the area-integrated spectra at multiple torus positions, as well as the HCN *J* = 3–2 and HCO$^+$ *J* = 3–2 emission lines detected in the reanalyzed area-integrated spectra at the same positions.
<figure><img src="image_19.png" /><figcaption>**Figure 12.** Gaussian fit (dashed line) of the detected HCN *J* = 4–3 (top), HCO⁺ *J* = 4–3 (second top), HCN *J* = 3–2 (second bottom), and HCO⁺ *J* = 3–2 (bottom) emission lines at the W-torus (1.7 pc diameter) (left), E-torus (1.7 pc diameter) (second left), torus (3.4 pc diameter centered at the continuum emission peak) (second right), and the torus with an even wider area (6.8 pc diameter) (right). The abscissa is the optical LSR velocity (in km s⁻¹) and the ordinate is the flux density (in mJy). The horizontal dotted straight line indicates the zero-flux level. Lower-frequency (redshifted) emission tails, relative to the molecular line emission at the systemic velocity (*V*<sub>sys</sub> = 1130 km s⁻¹), are discernible in some of the spectra in Figure 4. We apply two Gaussian fits to the HCN *J* = 4–3, HCO⁺ *J* = 4–3, and HCN *J* = 3–2 emission lines at the W-torus and at the torus with a 3.4 pc diameter aperture. For the HCN *J* = 3–2 emission line, a second redshifted emission component is also added at the E-torus (j). For the remaining data, a single Gaussian fit is applied.</figcaption></figure>

## Appendix C
CS *J* = 7–6 Line

Figure 13 presents an integrated-intensity (moment 0) map and area-integrated spectra at multiple torus positions for the CS *J* = 7–6 (*ν*<sub>rest</sub> = 342.883 GHz) line. The CS *J* = 7–6

emission line is stronger (in the moment 0 map in Figure 13(a)) and broader (in the spectra in Figure 13(b), (c)) in the western torus than in the eastern torus, as also seen for HCN and HCO⁺ at both *J* = 4–3 and *J* = 3–2 (Figures 1 and 4).
<figure><img src="image_20.png" /><figcaption>**Figure 13.** (a): Integrated-intensity (moment 0) map of the CS *J* = 7–6 line. In the same way as in Figure 1(a), (i) the gray contours indicate simultaneously obtained ~350 GHz continuum emission with 50σ, 70σ, 90σ, and 110σ, (ii) the horizontal thick black bar corresponds to 1 pc at the distance of NGC 1068, and (iii) the filled white circle indicates the synthesized beam size with 0′′014 × 0′′013. Observed area-integrated spectra around the CS *J* = 7–6 line with (b): a 1.7 pc (24 mas) diameter circle centered at the W-peak, (c): a 1.7 pc diameter circle centered at the E-peak, and (d): a 3.4 pc (48 mas) diameter circle centered at the C-peak. The abscissa is the observed frequency (in GHz) and the ordinate is the flux density (in mJy). Downward arrows are added at the expected frequency of the CS *J* = 7–6 ($ν_{rest}$ = 342.883 GHz) for the systemic velocity of $V_{sys}$ = 1130 km s⁻¹. The horizontal dotted straight line indicates the zero-flux level.</figcaption></figure>

## ORCID iDs

Masatoshi Imanishi https://orcid.org/0000-0001-6186-8792

Yoshiaki Hagiwara https://orcid.org/0000-0002-9043-6048

Kouichiro Nakanishi https://orcid.org/0000-0002-6939-0372

Takuma Izumi https://orcid.org/0000-0001-9452-0813

Nozomu Kawakatu https://orcid.org/0000-0003-2535-5513

## References

* Aalto, S., Costagliola, S., Martin, F., et al. 2015b, A&A, 584, A42

* Aalto, S., Garcia-Burillo, S., Muller, S., et al. 2015a, A&A, 574, A85

* Antonucci, R. R. J. 1993, ARA&A, 31, 473

* Antonucci, R. R. J., & Miller, J. S. 1985, ApJ, 297, 621

* Bannikova, E. Y., Akerman, N. O., Capaccioli, M., et al. 2023, MNRAS, 518, 742

* Bock, J. J., Neugebauer, G., & Matthews, K. 2000, AJ, 120, 2904

* Butterworth, J., Holdship, J., Viti, S., Garcia-Burillo, S., et al. 2022, A&A, 667, A131

* Costagliola, F., Sakamoto, K., Muller, S., et al. 2015, A&A, 582, A91

* Dyda, S., Lovelace, R. V. E., Ustyugova, G. V., Romanova, M. M., & Koldoba, A. V. 2015, MNRAS, 446, 613

* Falstad, N., Aalto, S., Konig, S., et al. 2021, A&A, 649, A105

* Falstad, N., Hallqvist, F., Aalto, S., et al. 2019, A&A, 623, A29

* Gallimore, J. F., Baum, S. A., & O'dea, C. P. 2004, ApJ, 613, 794

* Gallimore, J. F., & Impellizzeri, C. M. V. 2023, ApJ, 951, 109

* Gallimore, J. F., Impellizzeri, C. M. V., Aghelpasand, S., et al. 2024, ApJL, 975, L9

* Gamez Rosas, V., van der Werf, P., Gallimore, J. F., et al. 2025, A&A, 699, A187

* Gao, Y., & Solomon, P. M. 2004, ApJ, 606, 271

* Garcia-Burillo, S., Combes, F., Ramos Almeida, C., et al. 2016, ApJL, 823, L12

* Garcia-Burillo, S., Combes, F., Ramos Almeida, C., et al. 2019, A&A, 632, A61

* Greenhill, L. J., & Gwinn, C. R. 1997, Ap&SS, 248, 261

* Greenhill, L. J., Gwinn, C. R., Antonucci, R., & Barvainis, R. 1996, ApJ, 472, L21

* Greve, T. R., Papadopoulos, P. P., Gao, Y., & Radford, S. J. E. 2009, ApJ, 692, 1432

* Honig, S. F. 2019, ApJ, 884, 171

* Honig, S. F., Prieto, M. A., & Beckert, T. 2008, A&A, 485, 33

* Hunter, T. R., Indebetouw, R., Brogan, C. L., et al. 2023, PASP, 135, 074501

* Hure, J.-M. 2002, A&A, 395, L21

* Imanishi, M., Baba, S., Nakanishi, K., & Izumi, T. 2023a, ApJ, 950, 75

* Imanishi, M., Baba, S., Nakanishi, K., & Izumi, T. 2023b, ApJ, 954, 148

* Imanishi, M., & Nakanishi, K. 2013, AJ, 146, 91

* Imanishi, M., Nakanishi, K., & Izumi, T. 2016a, ApJL, 822, L10

* Imanishi, M., Nakanishi, K., & Izumi, T. 2016b, ApJ, 825, 44

* Imanishi, M., Nakanishi, K., & Izumi, T. 2016c, AJ, 152, 218

* Imanishi, M., Nakanishi, K., & Izumi, T. 2018a, ApJ, 856, 143

* Imanishi, M., Nakanishi, K., Izumi, T., & Wada, K. 2018b, ApJL, 853, L25
Imanishi, M., Nguyen, D. D., Wada, K., et al. 2020, ApJ, 902, 99

Impellizzeri, C. M. V., Gallimore, J. F., Baum, S. A., et al. 2019, ApJL, 884, L28

Israel, F. P. 2023, A&A, 671, A59

Kameno, S., Sawada-Satoh, S., Impellizzeri, C. M. V., et al. 2020, ApJ, 895, 73

Knudsen, K. K., Walter, F., Weiss, A., et al. 2007, ApJ, 666, 156

Krips, M., Neri, R., Garcia-Burillo, S., et al. 2008, ApJ, 677, 262

Kuznetsov, O. A., Lovelace, R. V. E., Romanova, M. M., & Chechetkin, V. M. 1999, ApJ, 514, 691

Leroy, A. K., Usero, A., Schruba, A., et al. 2017, ApJ, 835, 217

Lodato, G., & Bertin, G. 2003, A&A, 398, 517

Martin, S., Aalto, S., Sakamoto, K., et al. 2016, A&A, 590, A25

Morishima, Y., Sudou, H., Yamauchi, A., Taniguchi, Y., & Nakai, N. 2023, PASJ, 75, 71

Nakajima, T., Takano, S., Kohno, K., Harada, N., & Herbst, E. 2018, PASJ, 70, 7

Nakajima, T., Takano, S., Tosaki, T., et al. 2023, ApJ, 955, 27

Papadopoulos, P. P., Zhang, Z.-Y., Xilouris, E. M., et al. 2014, ApJ, 788, 153

Quach, D., Dyda, S., & Lovelace, R. V. E. 2015, MNRAS, 446, 622

Raban, D., Jaffe, W., Rottgering, H., Meisenheimer, K., & Tristram, K. R. W. 2009, MNRAS, 394, 1325

Sakamoto, K., Aalto, S., Evans, A. S., Wiedner, M., & Wilner, D. 2010, ApJ, 725, L228

Sakamoto, K., Martin, S., Wilner, D. J., et al. 2021, ApJ, 923, 240

Shirley, Y. L. 2015, PASP, 127, 299

Solomon, P. M., & Vanden Bout, P. A. 2005, ARA&A, 43, 677

Sternberg, A., Genzel, R., & Tacconi, L. 1994, ApJ, 436, L131

Takano, S., Nakajima, T., & Kohno, K. 2019, PASJ, 71, S20

The CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 114501

van der Tak, F. F. S., Black, J. H., Schoier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627

Vollmer, B., Davies, R. I., Gratier, P., et al. 2022, A&A, 665, A102

Vorobyov, E. I., Lin, D. N. C., & Guedel, M. 2015, A&A, 573, A5

Wada, K., Papadopoulos, P. P., & Spaans, M. 2009, ApJ, 702, 63

Wada, K., Schartmann, M., & Meijerink, R. 2016, ApJL, 828, L19

Williamson, D., Honig, S., & Venanzi, M. 2020, ApJ, 897, 26