Counter-rotation and High-velocity Outflow in the Parsec-scale Molecular Torus of NGC 1068

The unifying model for active galaxies holds that all active galactic nuclei (AGNs) are powered by accretion onto a supermassive black hole (SMBH; Antonucci 1993; Ramos Almeida & Ricci 2017; Hickox & Alexander 2018). The accretion disk and broad-line region (BLR) are surrounded by a dusty, molecular medium, commonly referred to as the "obscuring torus." The torus selectively obscures sight-lines to the central engine, resulting in a dichotomy of AGN spectra. If the torus is viewed more nearly pole-on, we can see the BLR and classify the AGN as type 1; if the torus is viewed more nearly edge-on, the BLR may be obscured and the AGN is classified as type 2 (Lawrence & Elvis 1982). A parsec-scale torus must be geometrically thick to explain the large detection fraction of type 2 AGNs (Osterbrock & Shaw 1988; Tovmassian 2001).

NGC 1068 is one of the nearest, luminous Seyfert type 2 galaxies and has served as an archetype for unifying schemes (Antonucci & Miller 1985). Polarimetry of the reflected nuclear spectrum reveals the hidden BLR and ultraviolet (UV) continuum source (Antonucci & Miller 1985; Capetti et al. 1995; Inglis et al. 1995). The Seyfert nucleus of NGC 1068 is obscured by a dusty, molecular medium (Mason et al. 2006; Gallimore et al. 2016; García-Burillo et al. 2016; Imanishi et al. 2018). A flat spectrum radio source, "S1," marks the location of the central engine (Gallimore et al. 1996; Muxlow et al. 1996). Very Long Baseline Array (VLBA) continuum observations resolve the continuum source into an elongated, pc-scale structure oriented nearly at right angles to the kpc-scale radio jet; the radio source has been interpreted as a plasma disk located inside the obscuring molecular medium (Gallimore et al. 1997, 2004). The cm-wave continuum source associates with the disk of H₂O vapor megamaser emission with systemic velocity 1130 km s⁻¹ (optical convention, LSRK reference frame) and maximum rotation speed 335 km s⁻¹ at a distance of 0.6 pc from the central engine (Greenhill & Gwinn 1997; Gallimore et al. 2001; Gallimore & Impellizzeri 2019).

Recent Atacama Large Millimeter/submillimeter Array (ALMA) observations of the molecular interstellar medium (ISM) in the central few pc have revealed complicated kinematics (Gallimore et al. 2016; García-Burillo et al. 2016; Imanishi et al. 2018). With a ∼ 3 pc (004) beam,¹² ALMA observations of CO (J = 6 → 5) and HCN (J = 3 → 2) resolve an extended, roughly 10 pc by 4 pc distribution of molecular gas oriented nearly at right angles to the radio jet axis but nearly aligned with the major axis of the nuclear plasma disk, S1. We found evidence for high-velocity CO in a bipolar outflow, with radial velocities ranging to ±400 km s⁻¹ relative to systemic, in the central pc (Gallimore et al. 2016). Imanishi et al. observed HCN (J = 3 → 2) at ∼3 pc resolution, and they argued that, unlike CO (J = 6 → 5), the molecular gas traced by HCN better traces the rotational, disk-like region of the torus. However, they also argue that the kinematics are not consistent with pure rotation but show evidence for significant turbulence (i.e., enhanced velocity dispersion revealed by the second spectral moment) and perhaps counter-rotation compared to the molecular gas on hundred parsec scales.

In an effort to untangle the complicated geometry and kinematics of the molecular torus of NGC 1068, we obtained new ALMA observations of HCN (J = 3 → 2) at 1.4 pc (20 mas) resolution. The observations and data reduction are summarized in Section 2. The main results, including moment maps and a kinematic analysis, are presented in Section 3. We summarize our main conclusions in Section 4.

We used ALMA to observe NGC 1068 in Band 6 during ALMA Cycle 5 (project code 2017.1.01666.S). We used the longest baseline configuration available for Band 6 during Cycle 5, C43-10, with baselines ranging from 41.4 m to 16.2 km. The largest recoverable structure is about 032, or 22 pc at the distance of NGC 1068. The source was observed in nine execution blocks between 2017 October 10 and 2017 October 24. The average exposure time per execution block was 1.6 hr, and the total integration time on source was 14.5 hr. Source observations were interleaved with short scans of nearby bright calibrator sources. The phase calibrator was J0239−0234, located 27 from NGC 1068, and the flux calibrator was J0432−0120, with flux density S_ν = 0.97 ± 0.04 Jy.

Spectral windows were tuned to observe three continuum windows, 1.875 GHz total bandwidth each and one tuned to observe HCN (J = 3 → 2; ν₀ = 265.88618 GHz) at recessional velocity 1150 km s⁻¹ (optical convention, LSRK reference frame). The HCN line observation has channel-width ∼1.0 km s⁻¹ and spans ∼1000 km s⁻¹ total bandwidth.

The data were calibrated using CASA pipeline version 40896 (Pipeline-CASA51-P2-B) and imaged using CASA 5.1.1 (McMullin et al. 2007). We performed four cycles of phase-only self-calibration on the continuum data (Cornwell & Fomalont 1989), and the solutions were applied to the spectral line data. The visibility data were Briggs-weighted during Fourier inversion (Briggs 1995). The spectral line data were averaged in frequency to provide approximately 10 km s⁻¹ wide channels. Continuum subtraction of the line data was performed in the (u, v) plane using an average of channels with radial velocities $| {\rm{RV}}| \gt 600$ km s⁻¹ relative to systemic. The common restoring beam is 19.9 × 19.1 mas, PA 798, and the median line channel sensitivity is 1.6 mJy beam⁻¹.

The continuum image was produced by combining all continuum channels into a single image at effective ν = 256.5 GHz (multi-frequency synthesis; Sault & Conway 1999). The restoring beam of the continuum image is 20.2 × 19.6 mas, PA 783, and the rms sensitivity is 6.2 μJy beam⁻¹.

Figure 1 presents images of the HCN (J = 3 → 2) integrated line flux and 256 GHz continuum. The line is detected in absorption at the position of S1, and the brightest integrated emission extends about 7 pc from S1 both southeast and northwest of S1. ALMA resolves the integrated HCN emission, with surface brightness (uncorrected for absorption) peaking roughly 2 pc to either side of S1 and fading with distance from S1. The morphology is irregular, suggesting substructure on scales at least as small as the ∼1.4 pc beam. The orientation of the HCN major axis is PA ∼114°, estimated by drawing a line through the continuum peak and between the outer contours of the integrated HCN flux map. Image moment analysis agrees with this estimate of the major axis orientation to within ∼5°, which corresponds to a difference of less than half of a beam at the extent of the major axis. The major axis is more closely aligned with the milliarcsecond-scale radio continuum structure of S1 but nearly perpendicular to the radio jet on 100 pc scales.

Zoom In Zoom Out Reset image size

At lower surface brightness levels, the ALMA 256 GHz continuum image resolves into an X-shaped structure that resembles a bicone extending perpendicular both to the molecular gas distribution and the major axis of the VLBA continuum source. We have interpreted the VLBA continuum source as arising from a parsec-scale plasma disk (Gallimore et al. 2004). That the major axis of the fainter ALMA continuum emission aligns nearly perpendicular to the VLBA major axis is consistent with a disk-bicone scenario. On the sub-pc scale, the VLBA continuum contours also show X-shaped extensions away from the plasma disk plane that appear to line up roughly with the ALMA continuum extensions. Although a detailed analysis of the continuum is outside the scope of this Letter, we note that the 256 GHz continuum is probably dominated by free–free emission. The integrated flux density of the S1 region in a 013 (9.1 pc) square aperture is 12.7 ± 0.1 mJy. Based on the Very Large Array (VLA) 43 GHz measurement of Cotton et al. (2008) and assuming a 5% uncertainty in the flux scale, the continuum spectral index is α = −0.06 ± 0.03 (S_ν ∝ ν^α; 1σ  uncertainty). Accepting additional uncertainty owing to the difference in resolution (80 mas for the VLA), this measurement is consistent with the spectral index expected for optically thin free–free emission, α = −0.1.

We measured the radial velocity field of the spectral line cube using Gauss–Hermite polynomial fits to each spectrum in the HCN datacube (see van der Marel & Franx 1993 for a description of this technique); the result is shown in Figure 2. Outside of about 1 pc projected radius, the broad trend is redshifts relative to systemic on the eastern side and blueshifts on the western side. The kinematics are consistent with the results presented by Imanishi et al. (2018) but opposite the sense of motion of the H₂O masers (Greenhill & Gwinn 1997). Inside the central parsec, however, the kinematics more closely resemble the H₂O megamaser disk, with high-velocity redshifted HCN northwest of the nucleus and high-velocity blueshifted HCN southeast of S1.

Zoom In Zoom Out Reset image size

The rotating disk(s) interpretation is complicated by the velocity field northeast and southwest of the nucleus, which does not cleanly resolve from the near-nuclear gas. This axis coincides with the bipolar outflow seen in CO (J = 6 → 5; Gallimore et al. 2016), and so we speculate the complex radial velocity field northeast and southwest of the nucleus results from a beam-convolved blend of clouds in the rotating disk and clouds participating in the outflow. Supporting this picture, we find that the high-velocity HCN (J = 3 → 2) appears to bracket the CO ($J=6\to 5$) outflow. The simplest interpretation would have the HCN trace molecular gas in a rotating disk that collimates the molecular outflow detected in CO ($J=6\to 5$) emission and evidenced by the disturbed HCN radial velocities along the minor axis.

The nuclear absorption profile, provided in Figure 3, shows additional evidence for an outflow component in HCN (J = 3 → 2). We extracted the nuclear spectrum at the position of the nuclear radio source S1. To estimate a correction for line emission in the beam, we averaged the spectra from positions 1.4 pc (one beam) northeast and southwest of the nucleus and subtracted the result from the nuclear spectrum. The resulting corrected absorption profile shows a minimum centered near the systemic velocity, as expected for an edge-on, rotating disk surrounding the nuclear continuum source. However, the line profile is asymmetric: there is a pronounced blue wing on the absorption profile that extends to roughly −450 km s⁻¹ relative to systemic. The blue wing appears on both the raw and corrected spectra, and so it is likely not an artifact of the emission correction. For comparison, the maximum rotation speed of the inner disk is ∼380 km s⁻¹ (see below). The excess blueshifted absorption is consistent with outflow velocities approaching 450 km s⁻¹, exceeding the 400 km s⁻¹, maximum outflow velocities observed in CO ($J=6\to 5$) line emission (Gallimore et al. 2016).

Zoom In Zoom Out Reset image size

Position–velocity diagrams along the major (PA 114°) and minor (PA 24°) axes are shown in Figure 4. On the major axis diagram, the highest velocity emission agrees with the positional offsets and range of velocities spanned by the redshifted H₂O masers (see also Figure 2) and so likely arises from the inner disk. Within offsets ≲3 pc, the region of the brightest continuum emission, the brightest HCN line emission spans radial velocities −170 to 100 km s⁻¹ relative to systemic. Given the coincidence with the radio continuum and the presence of both redshifted and blueshifted emission, we interpret this region as including a significant contribution from spatially unresolved molecular outflow, which is consistent with our interpretation of the velocity map (Figure 2). The outflow emission is most likely overlapping with some of the higher-velocity rotating gas. The presence of both redshifted and blueshifted emission suggests a biconical outflow with an outflow axis aligned near the plane of the sky (see Das et al. 2006). For offsets ≳4 pc, the emission is preferentially redshifted to the southeast and blueshifted to the northwest, as is evident on the velocity map (Figure 2). This region appears to be well separated from the outflow region and runs roughly orthogonal to the known outflow and jet axes. Although some beam-smeared outflow emission may contribute to this outer region, it seems more likely that rotation dominates the kinematics of the outer part of the HCN emission along the major axis.

Zoom In Zoom Out Reset image size

The minor axis position–velocity diagram places some constraints on the scale of the outflow. Allowing for beam-smearing, HCN (J = 3 → 2) emission extends roughly 1 pc northeast of the nucleus, but the HCN emission is much weaker to the southwest (see Figure 1). The CO (J = 6 → 5) emission similarly spans roughly 1–2 pc both southwest and northeast of the nucleus. Of course, the presence of line absorption, sensitivity limitations, and, potentially, confusion with overlapping emission from the inner disk complicate the measurement of the extent of the molecular outflow in HCN. On the other hand, molecular gas lifted from the shielding of the nuclear disk(s) would be soon exposed to the ionizing continuum from the AGN, and we speculate that the molecular outflow might become ionized and feed into to the larger, ∼400 pc ionized gas outflow (see Cecil et al. 1990). The implication would be acceleration from maximum line-of-sight velocities of roughly 450 km s⁻¹ for outflowing molecular clouds (Figure 3) to 1000 km s⁻¹ for [O iii] emission line clouds within the central 14 pc (Das et al. 2006).

We measured the parsec-scale rotation curve from the HCN (J = 3 → 2) data using the conventional tangential velocity technique (Kerr & Westerhout 1965). To perform this analysis, we used an automated numerical technique to determine the tangential velocity as a function of offset along the major axis of the integrated HCN map. For projected r ≥ 1.2 pc, we extracted spectra at each pixel along the major axis (1 pixel = 0.2 pc). For r < 1.2 pc, we instead extracted two spectra, one northwest of the nucleus and one southeast along the axis defined by the maser spot distribution (PA 131°) and at the locations of the brightest peak line emission in that region. The northeast spectrum was extracted at a position located r = 0.5 pc from the nuclear continuum source. Owing to confusion with the blue wing of the nuclear absorption line, the extraction for the southeast position is farther from the nucleus, at projected radius r = 0.9 pc.

A Gauss–Hermite polynomial was fitted to each extracted HCN spectrum, and Monte Carlo realizations of the spectrum were generated by adding a randomly generated noise spectrum scaled to the observed rms. The procedure next bracketed the region surrounding the spectral peak and the nearest velocity where the simulated spectrum falls below the background noise. The average and standard deviation of 10⁴ trials, respectively, provide an estimate of the tangential velocity and its uncertainty. The median uncertainty is 30 km s⁻¹. Note that this method does not correct for the unknown, intrinsic velocity dispersion of the gas. However, intrinsic velocity dispersions as high as σ_r = 90 km s⁻¹, comparable to the inferred rotation speeds, introduce corrections that are smaller than the measurement uncertainties.

The results are shown in Figure 4, which also includes the radial velocities of the nuclear H₂O megamasers and the Keplerian rotation curve that best fits the H₂O megamaser kinematics (Gallimore & Impellizzeri 2019).

The r > 1.2 pc rotation curve is consistent with an extrapolation of the Keplerian rotation curve. There are bumps and wiggles that deviate slightly from Keplerian rotation perhaps owing to confusion with the molecular outflow, i.e., overlap between the outflow and rotating gas within our beam; however, from inspection of Figure 4, these deviations are small compared to the measurement uncertainties. More specifically, the largest statistical deviation from the predicted Keplerian rotation curve is located 3.4 pc northwest of the nucleus. There, the measured rotational speed deviates from the predicted Keplerian rotation curve by 70 ± 40 km s⁻¹, which is a <2σ deviation.

Again we emphasize that, between r = 1.2 pc and r ∼ 7 pc, the outer HCN rotation curve runs counter to the inner maser disk rotation; that is, the masers recede toward the northwest, but the outer HCN disk recedes to the southeast. However, inside r = 1.2 pc, the HCN rotation follows the sense of the maser disk rotation. The inferred rotation speeds of the inner HCN disk reach 380 ± 30 km s⁻¹, exceeding the maximum rotation speed of the H₂O megamaser disk, 335 km s⁻¹. This result suggests that we have detected molecular gas closer to the central engine than the gas traced by the H₂O megamasers. Note also that the masers align along PA 131° (from the nucleus toward the approaching masers; Gallimore & Impellizzeri 2019), misaligned from the integrated, larger-scale HCN major axis by nearly 20°.

We used ALMA to observe 256 GHz continuum and the transition HCN (J = 3 → 2) at roughly parsec resolution in the molecular torus of NGC 1068. Comparing to VLA 43 GHz measurements, the continuum has a flat spectrum that is consistent with free–free emission. The continuum structure resembles a bicone viewed from the side, perpendicular to the symmetry axis.

Turning to the HCN observations, we identify three kinematically distinct regions: (1) an outflow component apparent in emission on the HCN position–velocity diagram and detected as a blueshifted wing in absorption against the nuclear continuum source, with projected outflow speeds approaching ∼450 km s⁻¹; (2) an inner disk spanning 0.5 ≲ r ≲ 1.2 pc; and (3) an outer disk extending to r ∼ 7 pc. The two disks counter-rotate, and the kinematics of the inner disk agree with the H₂O megamaser disk mapped by the VLBA (Greenhill & Gwinn 1997). The outer disk shows a Keplerian rotation curve consistent with an extrapolation of the rotation curve of the inner disk. We also find that the HCN radial velocity field is more complex along the molecular outflow axis, which suggests that we are also detecting but not fully resolving HCN emission associated with the outflow.

Conservation of angular momentum precludes the formation of counter-rotating disks out of a single epoch of accreting gas, and, to our knowledge, accretion disk warping mechanisms cannot induce counter-rotation (e.g., Tremaine & Davis 2014 and references therein). It seems likely that the counter-rotating outer disk formed from gas that has been more recently introduced to the >1 pc-scale environment, perhaps from an errant cloud falling out of the host galaxy or a captured satellite dwarf galaxy. Imanishi et al. (2018) demonstrated that the outer disk of the parsec-scale molecular torus also counter-rotates with respect to the molecular ring at 100 pc scales (see their Figure 1) and, we note, the sense of rotation of the galactic disk on kiloparsec scales (see, e.g., Helfer & Blitz 1995 and Brinks et al. 1997).

Because the central SMBH appears to dominate the gravitational potential in this region, the counter-rotating disks should be stable until, owing to dissipation, the outer disk starts to collapse onto the inner disk. When that occurs, we would expect the outer disk to collapse on an orbital timescale, t_orb ∼ 3 × 10⁴ yr. The interaction between counter-rotating disks may enhance the accretion rate by one or two orders of magnitude (Kuznetsov et al. 1999; Nixon et al. 2012; Dyda et al. 2015). King & Pringle (2006, 2007) described a chaotic accretion scenario for the cosmological growth of black holes, namely, rapid accretion owing to randomly oriented and episodic accretion events. They cite the random orientation of Seyfert radio jets and ionization cones with respect to host galaxies as evidence for chaotic accretion (see Schmitt et al. 2003; Gallimore et al. 2006). In NGC 1068, we find direct evidence that the molecular obscuring medium consists of counter-rotating and misaligned disks on parsec scales.

We thank Santiago García-Burillo, the referee, and the participants of the TORUS 2018 workshop (Puerto Varas, Chile; 2018 December, https://www.torus2018.org/), where we first presented the discovery of the counter-rotating outer disk, for helpful feedback that improved this work. This Letter makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.01666.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. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Facility: ALMA. -

Software: CASA (McMullin et al. 2007), astropy (Astropy Collaboration et al. 2013).