© ESO, 2014

1. Introduction

The study of the content, distribution, and kinematics of interstellar gas is a key to understanding the origin and maintenance of nuclear activity in galaxies. The processes involved in the fueling of active galactic nuclei (AGNs) encompass a wide range of scales, both spatial and temporal (Combes 2003; 2006; Jogee 2006). Current mm-interferometers are instrumental in providing a sharp view of the distribution and kinematics of molecular gas, the dominant gas phase in galaxy nuclei, through extensive CO line mapping. Combined with high-resolution near infrared images, interferometric CO maps are used to derive the angular momentum transfer budget in the circumnuclear disks of AGNs (e.g., García-Burillo et al. 2005; Haan et al. 2009; Meidt et al. 2013). Maps of CO in nearby AGNs unveil a wide range of large-scale and embedded m = 1,2 instabilities in their central 1 kpc circumnuclear disks (CNDs). Results from the NUclei of Galaxies (NUGA) project, a CO interferometric survey of a sample of 25 nearby low-luminosity AGNs (García-Burillo et al. 2003), indicate that molecular gas is frequently stalled in rings, which are the signposts of gravity torque barriers, and that only ~1/3 of galaxies show smoking-gun evidence of AGN fueling (García-Burillo & Combes 2012). In agreement with the picture drawn from observations, the most recent state-of-the-art numerical simulations show that several mechanisms, namely large-scale and nuclear stellar bars as well as slowly precessing lopsided m = 1 instabilities, are expected to cooperate to drain the gas angular momentum at the different spatial scales of galaxy disks (Hopkins et al. 2010b; 2011; 2012).

Furthermore, the use of molecular tracers specific to the dense gas phase can probe the feedback of activity on the chemistry and energy balance/redistribution of the interstellar medium of galaxies. Observations suggest that the excitation and chemistry of the main molecular species in AGNs are different with respect to those found in purely star-forming galaxies (Tacconi et al. 1994; Kohno et al. 2001; Usero et al. 2004; Graciá-Carpio et al. 2008; Krips et al. 2008; 2011; García-Burillo et al. 2010; Imanishi et al. 2007; 2013; Aladro et al. 2013). Although theoretical models have tried to explain these differences, the underlying reasons for the apparent AGN specificity are still debated (Lepp & Dalgarno 1996; Maloney et al. 1996; Meijerink & Spaans 2005; Meijerink et al. 2007; Yamada et al. 2007; Harada et al. 2013). Molecular outflows, which are considered as a footprint of the mechanical feedback of activity, are being discovered in a growing number of nearby active galaxies including ultra luminous infrared galaxies (ULIRGs), radio galaxies, and Seyferts (Feruglio et al. 2010; Sturm et al. 2011; Alatalo et al. 2011, 2012; Chung et al. 2011; Dasyra & Combes 2012; Combes et al. 2013; Morganti et al. 2013; Cicone et al. 2012; 2014). Radiative and mechanical feedback is often invoked as a mechanism of self-regulation in galaxy evolution (Di Matteo et al. 2005; 2008). Observations of nearby AGNs, where the distribution and kinematics of molecular gas can be spatially resolved, are thus instrumental if we are to understand if and how gas accretion can self-regulate in galaxies.

1.1. The prototypical Seyfert 2 galaxy NGC 1068

NGC 1068 is a prototypical nearby (D = 14 Mpc; Bland-Hawthorn et al. 1997) Seyfert 2 galaxy. It has a large-scale oval and a nuclear bar with a pseudo-bulge which is overly massive with respect to its central black hole (e.g., Kormendy & Ho 2013). NGC 1068 has been the subject of numerous campaigns using molecular line observations to study the fueling and the feedback of activity. Schinnerer et al. (2000) used the Plateau de Bure Interferometer (PdBI) to map the emission of molecular gas in the central r ~ 1.5–2 kpc disk using the J = 1−0 and J = 2−1 lines of CO. The CO maps spatially resolved the distribution of molecular gas in the disk, showing a prominent starburst (SB) ring of ~1–1.5 kpc – radius, which contributes significantly to the total CO luminosity of NGC 1068 (see also Planesas et al. 1991; Helfer et al. 1995; and Baker 2000). Furthermore CO emission is detected in a central r ~ 200 pc CND that surrounds the AGN. The CO(J = 3−2) map recently obtained by Tsai et al. (2012) with the Submillimeter Array (SMA) gives a similar picture of the large-scale distribution of molecular gas.

The kinematics of molecular gas in NGC 1068 were first interpreted as being due to the action of two embedded bars. The gas response seen in the CND, characterized by strong non-circular motions observed in CO, SiO, and CN maps, indicates that molecular clouds would be trapped between the Inner Lindblad Resonances (ILRs) of the inner nuclear bar and thus cannot fuel the AGN at present (Schinnerer et al. 2000; Baker 2000; García-Burillo et al. 2010). Alternative scenarios invoke the existence of large-scale outflow motions in the CND (Galliano & Alloin 2002; Davies et al. 2008; García-Burillo et al. 2010; Krips et al. 2011); these models would also suggest that the AGN feeding is presently thwarted.

However, closer to the nucleus, at r ≤ 50 pc, the kinematics of the molecular gas revealed by the 2.12 μm H₂ 1–0 S(1) map of Müller-Sánchez et al. (2009) give a completely different picture. These data, sensitive to hot (T_K ≃ 10³ K) and moderately dense (n(H₂) ≃ 10³ cm^-3) molecular gas, show elliptical streamers that bridge the CND and the central engine. Müller-Sánchez et al. (2009) suggest that these structures correspond to gas feeding the AGN. Although on different spatial scales, inflowing and outflowing gas may therefore coexist in the CND. However, the H₂ map only traces a small fraction of the total molecular gas reservoir of the CND, which is known to be much denser (n(H₂) ≃ 10^{5 − 6} cm^-3; Sternberg et al. 1994; Usero et al. 2004; Krips et al. 2008; 2011; Pérez-Beaupuits et al. 2007; 2009). The use of molecular tracers which are more representative of the total H₂ content is thus essential to derive the mass and energetics associated with the different inflowing/outflowing gas components at the CND.

Interferometric images of NGC 1068, obtained in tracers specific to the dense molecular gas, have also revealed the existence of a strong chemical differentiation in the disk of this Seyfert. Essential diagnostic line ratios are different in the SB ring compared to the CND. Tacconi et al. (1994) derived a high HCN/CO intensity ratio (~1) in the CND, about a factor of 5–10 higher than the ratio measured in the SB ring (Usero et al., in prep.). Radiative transfer calculations showed that the abundance of HCN relative to CO is globally enhanced in the CND: HCN/CO ~ 10^-3 (Sternberg et al. 1994; Usero et al. 2004; Krips et al. 2011; Kamenetzky et al. 2011). The detection of molecular ions like HOC⁺ and H₃O⁺ has been interpreted as the signature of X-ray processing (Usero et al. 2004; Aalto et al. 2011). García-Burillo et al. (2010) analyzed the likely drivers of chemical differentiation inside the CND with high-resolution observations of CN and SiO. The abundances of SiO and CN are enhanced at the extreme velocities of gas associated with non-circular motions/shocks close to the AGN (r< 70 pc). On the other hand, the correlation of CN/CO and SiO/CO ratios with hard X-ray irradiation suggests that the CND is a giant X-ray-dominated region (XDR). Although these results imply a strong radiative and mechanical feedback in the CND, the mechanism that drives the excitation and chemistry of molecular gas is yet to be elucidated.

1.2. This project

We use here the Atacama large millimeter array (ALMA) to map the emission of a set of molecular gas tracers (the J = 3 − 2 and J = 6 − 5 lines of CO, the J = 4 − 3 lines of HCN and HCO⁺, and the J = 7 − 6 line of CS) and their underlying continuum emission in the central r ~ 2 kpc of NGC 1068 with spatial resolutions ~0.3″ − 0.5″ (20–35 pc). These line transitions span a range of critical densities (n_crit [CO(3 − 2)] ~ a few 10⁴ cm^-3, n_crit [CO(6 − 5)] ~ a few 10⁵ cm^-3, n_crit [HCO⁺(4 − 3)] ~ 10⁷ cm^-3, n_crit [CS(7 − 6)] ~ a few 10⁷ cm^-3, and n_crit [HCN(4 − 3)] ~ 10⁸ cm^-3). The actual average densities of molecular gas probed by these lines in NGC 1068 are seen to be ~10⁵⁻⁶ cm^-3 due to subthermal excitation of the higher J-lines and radiative trapping (Krips et al. 2011; Viti et al. 2014, hereafter, Paper II). The line and continuum maps of the dense molecular gas and dust emission greatly improve the sensitivity and spatial resolution of any previous interferometric study of NGC 1068 in the (sub)millimeter range. We analyze the distribution, kinematics, and excitation of the molecular gas and study the processes associated with the fueling and the feedback of activity in this Seyfert. In Paper II we use the line ratio maps derived in this work to model the excitation and chemistry of molecular gas in the different environments of the CND and SB ring. Star formation laws will be analyzed in a future paper (García-Burillo et al., in prep.; Paper III).

We describe in Sect. 2 the ALMA observations and the ancillary data used in this work. Section 3 presents the continuum maps obtained at 349 GHz and 689 GHz. Section 4 discusses dust masses derived from continuum observations and illustrates the use of the CLUMPY torus models (Nenkova et al. 2008a,b) to constrain the parameters of the torus. The distribution of the molecular gas derived from the CO, HCN, HCO⁺, and CS line maps is discussed in Sect. 5. We describe the kinematics of the molecular gas and derive the main properties of the outflow component in Sect. 6. A first description of line ratio maps is presented in Sect. 7. The main conclusions of this work are summarized in Sect. 8.

2. Data

2.1. ALMA data

We observed the CO(J = 3 − 2) emission and the CO(J = 6 − 5) emission in NGC 1068 with ALMA during Cycle 0 using Band 7 and Band 9 receivers (project-ID: #2011.0.00083.S). The data in the two bands were calibrated using the ALMA reduction package CASA¹ while the calibrated uv-tables were subsequently exported to GILDAS² where the mapping and cleaning were performed as detailed below. Hereafter we adopt a distance to NGC 1068 of D ~ 14 Mpc (Bland-Hawthorn et al. 1997); this implies a spatial scale of ~70 pc/″.

2.1.1. Band 7 maps

In order to cover the CND and the SB ring of NGC 1068 we used an eleven-field mosaic in Band 7 with a field-of-view of 17″ per mosaic pointing. In total four tracks were observed between June and August 2012 resulting in initial on-source times (i.e., before flagging) of ~24–39 min per track or a total of 138 min. Between 18 and 27 antennas were available during the observations with projected baselines ranging from 17 m to 400 m. Weather conditions were good with median system temperatures of T_sys = 120 − 220 K and peak values not exceeding 300 K. Four spectral windows with a bandwidth of 1.875 GHz each and a spectral resolution of 488 kHz were placed in Band 7, two in the lower sideband (LSB) and two in the upper sideband (USB); the centers of the two sidebands are separated by 12 GHz. The four windows were centered on the following sky frequencies: 341.955 GHz and 343.830 GHz in the LSB and 353.830 GHz and 354.705 GHz in the USB. This setup allowed us to simultaneously observe CO(J = 3 − 2) (345.796 GHz at rest) and HCO⁺(J = 4 − 3) (356.734 GHz at rest) in the LSB bands, and HCN(J = 4 − 3) (354.505 GHz at rest) and CS(J = 7 − 6) (342.883 GHz at rest) in the USB bands. J0522-362 and 3C454.3 were used as bandpass calibrators and J0217+017 was used to calibrate the amplitudes and phases in time. To set the absolute flux scale, Uranus and Callisto were observed using the Butler-JPL-Horizons 2012 catalogue to model their visibilities as both were resolved at the given resolution of the observations. We estimate that the absolute flux accuracy is about 10–15%.

The angular resolution obtained using natural weighting was ~ at a position angle of ~60° in all the line and continuum data cubes. The conversion factor between Jy beam^-1 and K is 37 K Jy^-1 beam. The line data cubes were binned to a common frequency resolution of 14.65 MHz (equivalent to ~12–13 km s^-1 in Band 7). The point source sensitivities in the line data cubes were derived selecting areas free from emission in all channels. They range from 2 mJy beam^-1 (in the CS line) up to 2.8 mJy beam^-1 (in the CO line) in channels of 12.4–12.8 km s^-1 width. Images of the continuum emission were obtained by averaging in each of the four sub-bands those channels free of line emission. The resulting maps were averaged to produce an image of the continuum emission of the galaxy at 349 GHz. The corresponding point source sensitivity for the continuum is 0.14 mJy beam^-1.

As our observations do not contain short-spacing correction, we expect that a significant amount of flux will be filtered out on scales ≥6″ for continuum emission in Band 7. A comparison with the single-dish fluxes measured by Papadopoulos & Seaquist (1999) with the James Clerk Maxwell Telescope (JCMT) in apertures of ~15″ indicates that we are filtering up to ≃2/3 of the total flux on these scales. This will affect mostly the low level emission that may extend on large-scales in the interarm and spiral arm region. However, fluxes estimated in small (~1 − 4″) apertures centered on the brightest spots of emission of the CND and the SB ring are not expected to be significantly underestimated. The percentage of missing flux in the continuum maps in Band 7 represents the worst case scenario: the same estimate points to significantly lower missing fluxes in the molecular line ALMA maps. In order to estimate the missing flux due to the lack of short-spacing correction in the CO(3–2) ALMA map we used as references the fluxes measured in two studies of the galaxy done with single-dish telescopes: the single-pointed observations done at 22″-angular resolution with the Heinrich Hertz Telescope (HHT) by Mao et al. (2010), and the fully-sampled map of the central 1′ region of NGC 1068 done at 14″-angular resolution with the JCMT by Israel (2009a). Compared on these apertures (14–22″), the ALMA map recovers 80–90% of the CO(3–2) flux at the CND and over the SB ring. This percentage is lowered to 60–70% at the edges of the area mapped by the ALMA mosaic and over the interarm region. As expected, the clumpy distribution of gas aided by the velocity structure of the emission helps us recover the bulk of the flux over the map on these spatial scales. Fluxes estimated on smaller (~1 − 4″) apertures centered on the brightest emission spots of the CND and the SB ring are thus not expected to be significantly affected.

We set the phase tracking center of the central field of the mosaic to α₂₀₀₀ = 02^h42^m40.771^s, δ₂₀₀₀ = − 00°00′47.84″, which is the galaxy’s center according to SIMBAD (taken from the Two Micron All Sky Survey–2MASS survey; Strutskie et al. 2006). By default, (Δα, Δδ)-offsets are relative to this position. Rest frequencies for all the lines were corrected for the recession velocity initially assumed to be v_o(LSR) = 1135 km s^-1 = v_o(HEL) = 1146 km s^-1. Systemic velocity is nevertheless re-determined in this work (Sect. 6.1) as v_sys(LSR) = 1116 km s^-1 = v_sys(HEL) = 1127 km s^-1. Relative velocities throughout the paper refer to v_sys.

2.1.2. Band 9 maps

Given the small field-of-view (~9″) and the difficulty of Band 9 observations, we used only one field to cover the CND in NGC 1068. In total four different tracks were observed during July and August 2012 of which one had to be discarded because of poor quality. Weather conditions during the remaining three tracks were excellent with median system temperatures of T_sys = 600 K for a range of 400–1200 K. This results in an initial (i.e., prior to flagging) on-source time of 7–29 minutes per track or a total of ~52 min. Between 21 and 27 antennas were used during each track with projected baselines ranging between 17 m and 400 m. Three spectral windows with a spectral bandwidth of 1.875 GHz and a spectral resolution of 488 kHz were placed in the USB at sky frequencies of 686.899 GHz, 688.865 GHz, and 690.884 GHz, while only one window with a spectral bandwidth of 1.875 GHz and a spectral resolution of 488 kHz were placed in the LSB at a sky frequency of 674.944 GHz; the centers of the two sidebands are separated by ~16 GHz in Band 9. This setup allowed us to cover the CO(J = 6 − 5) line entirely and having enough line-free channels to determine the continuum emission. Bandpass calibration was performed on either Uranus or Ganymede while phase and amplitude calibration in time was done on J0339-017. The absolute flux scale was set by using models from the Butler-JPL-Horizons 2012 catalogue for Uranus, Ganymede, and/or Callisto. For most tracks, we had at least two solar bodies available so that we could crosscheck the flux calibration. However, we assume a conservative accuracy of ~30% for the absolute flux scale. We performed a self-calibration on the final maps using the strongest peak emission, located in the eastern knot of the CND, which purposely coincides within 0.1″ with the adopted phase tracking center. Self-calibration helped to improve the fidelity and more importantly the dynamic range of the final images.

The synthesized beam obtained using natural weighting is ~ at a position angle of ~50°. The conversion factor between Jy beam^-1 and K is 27 K Jy^-1 beam. The line data cube was binned to a frequency resolution of 29.30 MHz (equivalent to ~12.8 km s^-1). The point source sensitivity in the line data cube, derived by selecting areas free of emission is 23 mJy beam^-1 in channels of 12.8 km s^-1 width. Images of the continuum emission were obtained by averaging in each of the two nearby sub-bands at 686.899 GHz and 690.884 GHz those channels free of line emission. The resulting maps were combined to produce an image of the continuum emission of the galaxy at 689 GHz. The corresponding point source sensitivity for the continuum is 1.95 mJy beam^-1.

As our observations do not include zero-spacings, we expect to filter a non-negligible amount of flux on scales ≥3″ for continuum emission at 689 GHz. The total integrated flux measured inside the 9″ field-of-view amounts to ~815 mJy. A comparison with the fluxes measured with the JCMT by Papadopoulos & Seaquist (1999), and with Herschel by Hailey-Dunsheath et al. (2012), on apertures of ~10″ indicates that we are filtering at most ≃30–45% of the total flux on these scales. However, fluxes estimated on small (~1 − 3″) apertures centered on the brightest spots of emission of the CND are not expected to be significantly affected. As it happens in Band 7 observations, a similar estimate would also suggest that significantly lower missing fluxes are likely to be filtered out in the CO(6–5) maps, due to the high clumpiness and velocity structure of the emission.

The phase tracking center of the central field and the reference velocity are the same as the ones assumed in Band 7 (see Sect. 2.1.1 and Fig. 1).

2.2. Ancillary data

We retrieved the HST NICMOS (NIC3) narrow-band (F187N, F190N) images of NGC 1068 from the Hubble Legacy Archive (HLA). These images were completely reprocessed by the HLA with a fully revamped calibration pipeline which includes temperature-dependent dark frames, improved temperature measurements from bias levels, and other routines to reduce the impact of image artifacts and improve the calibrated data. The pixel size of the HLA images is 0.̋10 square. The redshift of NGC 1068 is sufficiently small that Paα falls within the F187N filter. We checked the alignment and background subtraction, then scaled the F190N and subtracted it from the F187N image. Although the two images are very close together in wavelength, their ratio is not exactly unity. Following Thompson et al.(2001), we determined the scaling factor by examining the ratio of the two images at the nucleus, where the F190N emission peaks. Because the Paα line is not centered at the filter central wavelength, we corrected the flux and converted the data numbers in the image to flux densities using the IRAF task synphot.

We also use the hard X-ray image obtained in the 6–8 keV band by Chandra. The X-ray map was obtained by applying an adaptive smoothing on the Chandra archive data available in this wave-band for NGC 1068, following the same procedure described by Usero et al. (2004). This image is identical to the hard X-ray map published by Ogle et al. (2003) (see also Young et al. 2001).

3. Continuum emission

3.1. Observations at 349 GHz

Figure 1a shows the continuum emission map obtained in Band 7 (at 349 GHz). The emission detected inside the inner r ~ 25″ (1.7 kpc) region imaged by ALMA, with a total integrated flux ~355 mJy, is distributed in three distinct regions:

• 1.

  The CND: the brightest emission stems from this centralcomponent, which appears as an asymmetric elongated ring of – (non-deprojected) size as shown in Fig. 1b. The intrinsic (deprojected) shape of the ring (~370 pc× 200 pc), derived assuming the disk geometry fit of Sect. 6.1 (PA = 289 ± 5° and i = 41 ± 2°), reinforces the idea that the ring is of elliptical shape. The ring, which shows a rich substructure, with two knots of emission located east and west of the nucleus, is remarkably off-centered relative to the location of the AGN. The latter is identified as the strongest emission peak in the continuum map at [Δα, Δδ] ≈ [–0.9″, –0.1″] = [α₂₀₀₀ = 02^h42^m40.71^s, δ₂₀₀₀ = − 00°00′47.94″]. This maximum coincides within the errors with the position of the AGN core (S1) determined from VLBI radiocontinuum maps of the galaxy (Gallimore et al. 1996; 2004). The striking off-centered ring morphology of the CND is echoed by the molecular line emission maps discussed in Sect. 5. This particular geometry questions the simplest scenario that interprets the gas ring as the signature of the ILR region located at r ≤ 5″ (Schinnerer et al. 2000; García-Burillo et al. 2010).

• 2.

  The bar: inside the domain occupied by the stellar bar, which extends along PA ~ 46 ± 2° up to a corotation radius R_COR ~18 ± 2″ (~1.3 kpc) (Scoville et al. 1988; Bland-Hawthorn et al. 1997; Schinnerer et al. 2000; Emsellem et al. 2006), we identified an arc of dust emission on the northeast side of the disk at distances ~4 − 7″ from the AGN (deprojected radii r ~ 350 − 650 pc). The source, hereafter referred to as the bow-shock arc, has a V-shaped morphology on its far side. This feature coincides in position with the northern edge of the AGN nebulosity identified in ionized gas emission, and with the northeast radio lobe. Wilson & Ulvestad (1987) interpreted the limb-brightened shape and the high polarization of the northeast radio lobe as the signature of a bow-shock in the galaxy disk. Leaving aside some isolated clumps, no significant emission is detected anywhere else inside the bar region.

• 3.

  The SB ring: most of the emission in the disk (≥72% of the flux integrated within the ALMA mosaic) is detected from a two-arm spiral structure that starts from the ends of the stellar bar and unfolds in the disk over ~180° in azimuth forming a pseudo-ring at r ~ 18″(~1.3 kpc), i.e., the position of the ILR resonance of the outer oval (Schinnerer et al. 2000). Continuum emission at 349 GHz over the SB ring is highly clumpy.

3.2. Observations at 689 GHz

Figure 1c shows the continuum emission map obtained in Band 9 (at 689 GHz) in the inner ~9″ field-of-view imaged by ALMA, which covers completely the CND region³. Similar to the picture drawn from observations in Band 7, the CND emission in Band 9 stems from a clumpy ring that is off-centered relative to the AGN. The prominent east and west knots identified in Band 7 are spatially resolved into several clumps, due to the higher spatial resolution at 689 GHz. Furthermore, while emission is also detected at the AGN locus, like in Band 7 observations, the strongest peak coincides with the position of the east knot. There is also a secondary peak of emission in Band 9 which is shifted to the northeast forming a spatially resolved ~35 pc – size elongated structure that connects the AGN with the CND along PA ~ 25°. This feature is reminiscent of the ~north-south elongated source present in the mid-infrared (MIR) maps of the central r ~ 1 − 2″ of the galaxy. This component has been attributed to emission of hot dust from narrow-line-region (NLR) clouds (Bock et al. 1998; 2000; Alloin et al. 2000; Tomono et al. 2001; 2006; Galliano et al. 2005; López-Gonzaga et al. 2014). Gallimore et al. (2001) mapped the subarcsecond radio jet and found evidence of a jet-cloud interaction around (20 pc) north and PA ~ 20° relative to the AGN locus, based on the detection of strong H₂O and OH maser emission and the bending of the radio jet at this location. This is very close to the position of a dust cloud identified in our maps, as shown in Fig. 3.

4. Dust masses

We combined observations of continuum emission in Bands 7 and 9 of ALMA with observations done in two PACS bands at 70 μm and 160 μm by Hailey-Dunsheath et al. (2012), and with high-resolution observations done at near- and mid-infrared (NIR and MIR) wavelengths by Alonso-Herrero et al. (2011) to construct the spectral energy distributions (SED) and derive the mass of the dust in two different environments of NGC 1068: the r ≤ 20 pc region centered at the AGN (Sect. 4.1), and the central r ≤ 400 pc of the galaxy, a region that includes the CND and the bow-shock arc (Sect. 4.2).

4.1. The central r ≤ 20 pc

4.1.1. Thermal versus non-thermal emission

Figure 2 shows the map of the 689 GHz–to–349 GHz continuum flux ratio (hereafter, B9 /B7) in the CND derived at the (lower) spatial resolution of observations in Band 7, after convolution of the 689 GHz continuum image with a Gaussian beam with full width at half maximum (FWHM) ~0.4″. The B9 /B7 ratio is seen to change from ≃5 close to the AGN () up to a maximum of ≃20–25 in a ring-like region of radius r ≃ 1.5″. The highest ratios can be explained if we assume a steep dependence of the dust emissivity κ_ν, which scales as ~ν^β, on frequency (β ≥ 2.5), and attribute high temperatures to dust T_dust ≥ 100 K in the ring. On the other hand, the low value of B9 /B7 close to the AGN is indicative of a non-negligible contribution of non-thermal emission, which can be particularly relevant at 349 GHz. Hönig et al. (2008) and Krips et al. (2011) discussed the relevance of three different mechanisms of non-thermal emission at the vicinity of the AGN: electron-scattered synchrotron emission, pure synchrotron emission, and thermal free-free emission. The inclusion of the two ALMA band fluxes (integrated in 1″-apertures) into the AGN SED confirms that pure synchrotron emission is a poor representation of the SED in the submillimeter range, as also argued by Krips et al. (2011): this scenario over-predicts by 50% the total flux measured in Band 7, i.e., well beyond the associated uncertainties (<10%). On the contrary the two alternative models discussed by Krips et al. (2011) account equally well for the AGN SED. In either case, the non-thermal contamination can at best be estimated to range between 30 and 65% in Band 7. However, this percentage is much lower at the higher frequencies of Band 9 (≤5–18%).

4.1.2. CLUMPY torus models: NIR/MIR and ALMA observations

We used the ALMA Band 9 and Band 7 continuum thermal fluxes at 435 μm and 860 μm, respectively, to investigate whether we can set further constraints on the properties of the putative torus of NGC 1068 using CLUMPY torus models (Nenkova et al. 2008a,b). To do so we took into account the NGC 1068 SED and MIR spectroscopy presented in Alonso-Herrero et al. (2011). We added the ALMA Band 9 measurement at the full resolution of , which is similar to that of the MIR data. The ALMA Band 7 continuum measurement has a slightly worse resolution () and a higher contribution from non-thermal emission as discussed in Sect. 4.1.1, and therefore we used the estimated thermal component as an upper limit. For the ALMA Band 9 flux uncertainties we included both the error in the photometric calibration (~30%) and the uncertainties associated with the modeling of the non-thermal component at this wavelength estimated above (~10%). The details of the fitting procedure are explained in Appendix A.

The fitted torus parameters are all similar to those of Alonso-Herrero et al. (2011) except for the torus size now fitted to a factor ~10 larger value, which corresponds to a torus outer radius of $\mathrm{2}{\mathrm{0}}_{-10}^{\mathrm{+}\mathrm{6}}\hspace{0.17em}$pc, derived using the AGN bolometric luminosity of ~$\mathrm{4.}{\mathrm{2}}_{-1.1}^{\mathrm{+}\mathrm{1.4}}\mathrm{\times }{\mathrm{10}}^{\mathrm{44}}$ erg s^-1 estimated from scaling the best fit model to the data. The torus size is not well constrained because the large value of the fitted index of the radial distribution of the clouds, which implies that most of the clouds are located close to the inner radius of the torus. For comparison the modeled 12 μm interferometric half radius of the resolved and unresolved components of NGC 1068 is 1.6 pc (Burtscher et al. 2013). These differences might be explained because the NIR and MIR emission is probing warm dust that is on average closer to the AGN, whereas in the submillimeter range we are more sensitive to cold dust distributed over larger distances from the AGN.

The best CLUMPY model fit and the 1σ uncertainty to the nuclear emission of NGC 1068 is presented in Fig. 3 (upper panel). It is clear from this figure that the measured continuum ALMA fluxes are above the CLUMPY torus fit. This could be explained if cold dust not associated with the torus was included in the ALMA Band 7 and 9 flux measurements. Figure 3 (lower panel) shows a close-up of the nuclear region of NGC 1068 at 435 μm at full resolution. Clearly there is much cold dust in this region including the ionization cone, as discussed in Sect. 3.2, but more importantly there is no point source associated with the position of the AGN. This seems to confirm our suspicion that even at the angular resolution of the ALMA Band 9 observations, there can be a contribution of cold dust not associated with the clumpy torus. We note that the fitted outer torus radius would correspond to an angular diameter of $\mathrm{0.5}{\mathrm{7}}_{-0.43}^{\mathrm{+}\mathrm{0.17}}$″. The lower limit is consistent with not having resolved the torus in Band 9 with the current angular resolution. Incidentally, the NIR data points are also above the CLUMPY torus model fit. The extra flux at these wavelengths might come from the dust at the base of the ionization cone, as can be the case with other Seyfert galaxies (see, e.g., Hoenig et al. 2013).

We can finally estimate the gas mass in the torus of NGC 1068 based on the fit, according to Eq. (A.1). We obtained: M_torus = 2.1( ± 1.2) × 10⁵ M_⊙, where we considered that the main uncertainties come from the relatively unconstrained torus size and from the scatter around the adopted N_H2/A_V scaling ratio taken from Bohlin et al. (1978) (see Appendix A). This mass estimate is comfortingly similar to the estimated molecular gas mass detected inside the central r = 20 pc aperture derived from the CO(3–2) emission, as discussed in Sect. 5.4.

4.2. The central r ≤ 400 pc region: the CND and the bow-shock arc

ALMA observations in Bands 7 and 9 were combined with PACS observations obtained at 70 μm and 160 μm in comparable apertures by Hailey-Dunsheath et al. (2012) to constrain the overall SED of the dust emission in the central r = 400 pc of the galaxy. The estimate of M_gas based on this fit and the implied conversion factor for CO (X_CO), discussed in Sect. 5.4, are used in Sect. 6 to derive the mass load of the outflow identified in this region.

The SED was fit using a modified black-body (MBB) model to derive the dust temperature (T_dust), dust mass (M_dust) and emissivity index (β) in this region. In this approach, the measured fluxes, S_ν, can be expressed as S_ν = M_dust × κ_ν × B_ν(T_dust)/D², where the emissivity of dust, κ_ν, scales as ~κ_{352 GHz} × (ν[GHz]/ 352)^β, with κ_{352 GHz} = 0.09 m² kg^-1 (a value rounded up from κ_{352 GHz} = 0.0865 m² kg^-1 used by Klaas et al. 2001), B_ν(T_dust) is the Planck function, and D is the distance. Prior to the fit, the fluxes in the two ALMA bands were corrected for the contamination by non-thermal emission in the central 1″ aperture at the AGN, as estimated in Sect. 4.1.1. The best MBB fit is found for M_dust = (8 ± 2) × 10⁵ M_⊙, T_dust = 46 ± 3 K and β = 1.7 ± 0.2. The errorbars on the parameters of the fit reflect the uncertainties due to the estimated range of missing flux in the ALMA bands, derived in Sects. 3.1 and 3.2.

The value of M_dust can be used to predict the associated neutral gas mass budget in this region. Applying the linear dust/gas scaling ratio of Draine et al. (2007) (see also Sandstrom et al. 2013) to the gas-phase oxygen abundances measured in the central 2 kpc of NGC 1068 (~12 +log (O/H) ~ 8.8; Pilyugin et al. 2004; 2007), which yields a gas-to-dust mass ratio of ~$\mathrm{6}{\mathrm{0}}_{-30}^{\mathrm{+}\mathrm{30}}$, we estimate that the total neutral gas mass in the central 10″ aperture is M_gas = (5 ± 3) × 10⁷ M_⊙. This number is a good proxy for the molecular gas mass in this region because the HI distribution in the disk of NGC 1068, studied by Brinks et al. (1997), shows a bright ring between 30″ and 80″ with a true central hole.

5. Molecular line emission

In this section we analyze our newly acquired ALMA maps and compare the morphology of the dense-gas tracers.

5.1. CO maps

Figure 4 shows the CO(3–2) and CO(6–5) velocity-integrated intensity maps of NGC 1068 obtained with ALMA. The CO maps reveal the distribution of the dense molecular gas in NGC 1068 with unprecedented high-dynamic range capabilities: ≥600 and ≥200 for the CO(3–2) and CO(6–5) maps, respectively. As for the dust emission, we identify three main regions in the disk of NGC 1068:

• 1.

  The CND: the brightest CO(3–2) line emission comes from theCND, which appears as a closed asymmetric elliptical ring of6″ × 4″–size as shown in Fig. 4b. The substructure of the CO(3–2) ring reveals two strong emission peaks located ~1″ east and ~1.5″ west of the AGN. Similar to the dust continuum ring, the CO(3–2) ring is noticeably off-centered relative to the location of the AGN: the two emission knots are bridged by lower-level emission that completes the ring north and south of the AGN and leaves a gas-emptied region to the southwest. However, unlike for Band 7 continuum emission, CO(3–2) emission does not peak at the AGN. The morphology of the ALMA map of the CND is to a large extent similar to that of the SMA map of Krips et al. (2011). Nevertheless, the order of magnitude higher dynamic range of the ALMA image of the CND helps reveal that the ring closes south of the AGN, i.e., similarly to the molecular ring imaged with the Very Large Telescope (VLT) in the 2.12 μm H₂ line by Müller-Sánchez et al. (2009). The brightest CO(3-2) emission features are also detected and spatially resolved in the CO(6–5) map of Fig. 4c.

• The bar: the CND appears connected to lower-level emission which extends farther out in the disk into the stellar bar region. On these scales the CO(3–2) emission is detected along two lanes, offset by 2 − 4″, which run mostly parallel to the bar’s major axis (PA = 46 ± 2°; Scoville et al. 1988; Bland-Hawthorn et al. 1997; Schinnerer et al. 2000; Emsellem et al. 2006). These are reminiscent of the gas leading edge morphology that is expected to prevail between the corotation of the bar and the ILR region. This is illustrated in Fig. 5, which shows the emission of CO(3–2) at two velocity channels symmetrically located relative to v_sys: v − v_sys = − 50 km s^-1 and +50 km s^-1, with v_sys(HEL) = 1127 ± 3 km s^-1 as determined in Sect. 6.1.1. Gas emission along the bar’s leading edges appears at the velocity range expected for gas falling to the nucleus: at redshifted velocities on the southwestern (near) side and blueshifted velocities on the northeastern (far) side. We nevertheless detect a gas component with anomalously redshifted velocities on the northeastern side of the disk at distances ~4 − 7″ from the AGN. This CO feature, closely associated with the bow-shock arc identified in dust continuum emission, is interpreted in Sect. 6 as the signature of a molecular outflow.

  Fig. 5 Zoom in of the inner disk of NGC 1068 to show the overlay of two CO(3–2) velocity channel maps obtained for v − vsys = − 50 (black contours: 5, 10, 20 to 100% in steps of 20% of the peak value = 1 Jy) and 50 km s-1 (white contours: 5, 10, 20 to 100% in steps of 20% of the peak value = 0.5 Jy), with vsys(HEL) = 1127 ± 3 km s-1, on the NIR K-band image (in color scale and arbitrary units) of NGC 1068 obtained by the Two Micron All Sky Survey (2MASS). We highlight the location of gas emission along the bar’s leading edges as well as the existence of an anomalous component associated with the outflow. The orientation of the stellar bar’s major axis along PA = 46° is marked by a dashed line.

  Fig. 6 a) Overlay of the CO(3–2) ALMA intensity contours (levels as in Fig. 4a) on the Paα emission HST map (color scale as shown in counts s-1pixel-1). b) Same as a) but zooming in on the CND region. c) Overlay of the CO(6–5) ALMA intensity contours (levels as in Fig. 4c) on the HST Paα emission map (color scale as shown). The filled ellipses at the bottom right corners represent the CO beam sizes.

  Fig. 7 a) Overlay of the CO(3–2) intensity contours (levels as in Fig. 4a) on the dust continuum emission at 349 GHz (color scale in Jy beam-1 units as indicated). b) Same as a) but zooming in on the CND region. c) Overlay of the CO(6–5) intensity contours (levels as in Fig. 4c) on the dust continuum emission at 689 GHz (color scale in Jy beam-1 units as indicated). The filled ellipses at the bottom right corners represent the CO beam sizes.

• 3.

  The SB ring: similar to the dust emission, most of the CO(3–2) flux in the disk mapped by ALMA (≃63% of the grand total) is detected in the SB ring. The two spiral arms, which unfold over ~180° in azimuth, are spatially resolved. The (transverse) widths of the arms go typically from 100 pc to 350 pc. Figure 6 shows that the SB ring concentrates most of the ongoing massive star formation in the disk identified by strong Paα emission. The emission of CO(3–2) over the SB ring is highly clumpy and it appears to be organized as coming from molecular cloud associations of ≥50 pc-size which are also identified in the dust continuum map, as shown in Fig. 7.

We also report the detection of CO(3–2) emission at different locations throughout the interarm region. These interarm complexes, which remained undetected in dust emission due to the lower dynamic range of the Band 7 continuum map, are organized into a network of filaments that extends out to the edge of our mapped region.

5.2. HCN, HCO⁺, and CS maps

Figure 8 shows the integrated intensity maps of NGC 1068 obtained with ALMA in the HCO⁺(4–3), HCN(4–3), and CS(7–6) lines. In stark contrast with the CO(3–2) map, most of the emission in these likely denser gas tracers that are characterized by ~a factor 100 comparatively higher critical densities, stems from the CND. Notwithstanding, a few (~12) isolated clumps are detected in the HCO⁺(4–3) and HCN(4–3) lines at significant levels towards the SB ring (Fig. 8a shows the HCO⁺ emission in the SB ring). The different distributions of the CO(3–2), HCN(4–3), and HCO⁺(4–3) line emission in the disk is reflected in the different line ratios measured in the CND and in the SB ring (see Sect. 7 and Paper II). Other molecular tracers also show remarkably different distributions in the SB ring and the CND (e.g., Takano et al. 2014).

The overall morphology of the CND in HCO⁺, and HCN is similar to that of the ALMA CO maps, i.e., an asymmetric closed molecular ring, which is noticeably off-centered relative to the AGN. However, the lower signal-to-noise ratio of the ALMA CS(7–6) map restricts the detected emission in this line mostly to the two prominent knots of the CND. The superior capabilities of ALMA dramatically improve the picture drawn for the distribution of the dense molecular gas traced by the HCN and HCO⁺ lines in the CND compared to the previous SMA data of Krips et al. (2011): while most of the emission detected in the HCO⁺(4–3) SMA map is restricted to the western knot, the ALMA map shows the emission to be widespread over the whole CND. The integrated flux of this line over the CND is a factor ~1.7 higher in the ALMA map (104 ± 10 Jy km s^-1) compared to the SMA map (60 ± 10 Jy km s^-1).

5.3. Molecular emission near the AGN

With the exception of the CS(7–6) line, emission is detected in all the molecular gas tracers probed by ALMA at the position of the AGN. Figure 9 shows the CO(3–2), CO(6–5), HCN(4–3), and HCO⁺(4–3) emission line profiles towards the AGN. As expected, molecular emission is within the errors centered around v = v_sys(HEL) = 1127 ± 3 km s^-1 (determined in this work: see Sect. 6.1.1) and extends in contiguous channels across 200–300 km s^-1 over significant levels (>1σ) in all tracers. Gaussian fits to the line profiles extracted from the same apertures (~40 pc) indicate that the higher density tracers (HCN(4–3) and HCO⁺(4–3)) show significantly wider lines (FWHM ~ 180 ± 10 km s^-1) compared to CO(3–2) (FWHM ~ 106 ± 3 km s^-1). This result indicates that the excitation of molecular gas, estimated from the HCN(4–3)/CO(3–2) and HCO⁺(4–3)/CO(3–2) line ratios, is enhanced at the highest velocities, which in all likelihood correspond to gas lying closer to the central engine. This trend runs in parallel with the observed tendency to find higher velocity-integrated ratios at the smallest radii inside the CND (see discussion in Sect. 7).

As shown in Fig. 9 the half width of the CO(6–5) line derived using the velocity channels that show contiguous emission above a 1σ-level is ~125 km s^-1. Adopting an inclination angle i = 40 − 41° (Bland-Hawthorn et al. 1997; Sect. 6.1), the implied spherical mass enclosed at (10 pc) is M_dyn~8 − 9 × 10⁷ M_⊙. This is a factor ~7–8 higher than the black hole mass (1 − 1.2 × 10⁷ M_⊙) estimated from the H₂O maser kinematics in the inner r ~ 0.7 pc of the galaxy (Greenhill et al. 1996; Gallimore et al. 2001), an indication that most of the dynamical mass at (10 pc) is contained in the central stellar cluster.

5.4. The CO–to–H₂ conversion factor

Using the CO(1–0) map of Schinnerer et al. (2000), we find that the CO(1–0) flux (S_{CO(1 − 0)}) integrated in the central 10″(700 pc)-aperture of the galaxy is ~100 Jy km s^-1. A comparison with the CO(1–0) flux derived using a similar aperture on the BIMA map of Helfer et al. (1995), which, unlike the Schinnerer et al. map, includes zero spacings, suggests that the PdBI map recovers ≃80% of the total flux in this region. The implied conversion factor X_CO = N(H₂) /I_{CO(1 − 0)} that is required to match the value of M_gas derived from the dust SED fitting discussed in Sect. 4.2 is ~$\mathrm{1}\mathit{/}\mathrm{\left(}{\mathrm{4}}_{-1}^{\mathrm{+}\mathrm{6}}\mathrm{\right)}\mathrm{\times }{\mathit{X}}_{\mathrm{CO}}^{\mathrm{MW}}$, where ${\mathit{X}}_{\mathrm{CO}}^{\mathrm{MW}}\mathrm{=}\mathrm{2}\mathrm{\times }{\mathrm{10}}^{\mathrm{20}}$ cm^-2 (K km s^-1)^-1 is the average CO conversion factor assumed to hold in the molecular clouds of the Milky Way disk (Strong et al. 1988; Strong & Mattox 1996; Dame et al. 2001). A similar conclusion, pointing to an X_CO ~a factor of 4–8 lower than ${\mathit{X}}_{\mathrm{CO}}^{\mathrm{MW}}$ in the CND of NGC 1068, was reached by Usero et al. (2004) who used LVG models to fit the CO line ratios observed in this region. This result is also in line with the observational work of Israel (2009a,b) and Sandstrom et al. (2013), who found that X_CO can be up to a factor 4–10 lower than ${\mathit{X}}_{\mathrm{CO}}^{\mathrm{MW}}$ in the central ~1 kpc of a subset of galaxy disks of solar metallicities. Bell et al. (2007) also derived conversion factors which are typically ~a factor of 10 lower than ${\mathit{X}}_{\mathrm{CO}}^{\mathrm{MW}}$ in their models of dense PDRs (n(H₂) ≥ 10⁴ cm^-3) of starburst nuclei. Lower X_CO factors have also been commonly assumed for ULIRGs; however, it is still debated whether this is a common feature of extreme starburst systems (e.g., see discussion in Papadopoulos et al. 2012).

An estimate of the molecular gas mass detected inside the central r = 20 pc aperture (M_gas[AGN]) can be made based on the detected CO(3–2) emission, provided that we assume a conversion factor for this line. Taking the average value derived above for X_CO, which is representative for the 1–0 line in the CND, the conversion factor for the 3–2 line has to be scaled down by an additional factor 3, based on the 3–2/1–0 line ratio ~3 (in T_mb units) derived in the neighborhood of the central engine of NGC 1068 (see Sect. 7). This implies that M_gas[AGN]$\mathrm{=}\mathrm{1.}{\mathrm{6}}_{-0.9}^{\mathrm{+}\mathrm{0.5}}\mathrm{\times }{\mathrm{10}}^{\mathrm{5}}{\mathit{M}}_{\mathrm{\odot }}$, i.e., the typical mass of a GMC and a value compatible within the errors with the gas mass derived from the dust emission in Sect. 4.1. This is equivalent to an average H₂ column density N_H2 = 5 × 10²¹ mol cm^-2.

Similarly, if we use the CO(6–5) integrated intensities in the ALMA aperture for Band 9, i.e., r = 10 pc, and the factor of ~3–4 scaled-down version of X_CO for this line estimated from the 6–5/1–0 ratio ~3–4 (in T_mb units) measured at the AGN, we derive M_gas[AGN]$\mathrm{=}\mathrm{1.}{\mathrm{2}}_{-0.7}^{\mathrm{+}\mathrm{0.4}}\mathrm{\times }{\mathrm{10}}^{\mathrm{5}}{\mathit{M}}_{\mathrm{\odot }}$. Considering the smaller aperture used in the 6–5 line, this value is, not surprisingly, slightly below the estimate derived from 3–2. This is equivalent to an average H₂ column density N_H2 = 1.2 × 10²² mol cm^-2.

6. Gas kinematics: the molecular outflow

Figures 10 and 11 show the mean-velocity field of molecular gas in the disk of NGC 1068 derived from the CO(3–2), CO(6–5), HCN(4–3), and HCO⁺(4–3) lines. By default, isovelocities are derived by integrating the emission above 5σ-levels throughout the disk in all tracers. The clipping is lowered to 3σ–levels when we zoom in on the CND region for HCN(4–3) and HCO⁺(4–3) in Fig. 11.

The CO(3–2) isovelocities of Fig. 10a, which sample the kinematics of molecular gas throughout the central ~40″ (2.8 kpc)-region, show the expected pattern for a rotating disk characterized by an overall east-west orientation of its kinematic major axis (PA = 278° ± 10°; Bland-Hawthorn et al. 1997; Schinnerer et al. 2000⁴). At close sight, the orientation of the major axis is seen to change from the CND region (PA ~ 330°) to the bar region (PA ~ 290°) and farther out to the spiral arm region (PA ~ 260°). The PA of the CND is close to the orientation of the kinematic major axis derived for the much smaller r ~ 0.7 pc rotating disk traced by H₂O maser emission published by Greenhill et al. (1996) (PA = 315°). The PA values for the CO disk, derived from a qualitative fit of isovelocities, are confirmed by the kinematic modeling of Sect. 6.1.1. If we assume that gas orbits are roughly coplanar, the observed trend can be interpreted as the footprint of non-circular motions on the gas flow. As discussed in Sect. 6.1.1, these distortions of the gas kinematics are caused by different mechanisms which are at work on different spatial scales: the spiral arm structure in the outer disk, the stellar bar on intermediate scales and the nuclear outflow in the CND. All molecular tracers show a similar tilt of the CND kinematic major axis relative to the large-scale disk as shown in Fig. 11.

The observed line widths (FWHM), displayed in Fig. 10d, go from 15–20 km s^-1 in the well detected (>5σ) molecular complexes of the interarm region, to an average value of 35–40 km s^-1 in the SB ring. The corresponding velocity dispersion estimates would be σ_v ~ 6 − 8 km s^-1 in the interarm region and σ_v ~ 15 − 17 km s^-1 in the SB ring. The measured widths at the CND, where FWHM values range from 50 to 200 km s^-1, likely overestimate σ_v for molecular gas in this region as they reflect the superposition of different velocity components inside the beam associated with rotation but also with an underlying strong outflow pattern (see Sect. 6.1). Figures 10e and 10f show that the region in the CND where FWHM values are above ~100 km s^-1 has a ring-like morphology. The widespread distribution of high line width values in this region suggests that the outflow, which is found in Sect. 6.1 to be a mode superposed onto rotation, spatially extends over most of the CND.

We analyze below the evidence supporting the existence of a massive molecular outflow in NGC 1068 based on the modeling of the velocity field derived from the CO(3–2) data (Sect. 6.1). We examine in Sect. 6.2 that the evidence for an outflow is also found in the molecular gas traced at higher densities by CO(6–5), HCN(4–3), and HCO⁺(4–3). The possible powering sources of the molecular outflow are discussed in Sect. 6.3. We explore in Sect. 6.4 if other alternative scenarios can explain the observed kinematics. Section 6.5 compares the properties derived for the molecular outflow detected by ALMA with those derived from previous studies of NGC 1068.

6.1. The CO(3–2) molecular outflow

6.1.1. Fourier decomposition of the velocity field

The general description of the two-dimensional line-of-sight velocity field of a galaxy disk can be expressed as ${\mathit{v}}_{\mathrm{los}}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathrm{=}{\mathit{v}}_{\mathrm{sys}}\mathrm{+}{\mathit{v}}_{\mathit{\theta }}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\cos \mathit{\psi }\sin \mathit{i}\mathrm{+}{\mathit{v}}_{\mathit{R}}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\sin \mathit{\psi }\sin \mathit{i,}$(1)where (v_R,v_θ) is the velocity in polar coordinates, ψ is the phase angle measured in the galaxy plane from the receding side of the line of nodes, and i is the inclination angle restricted to the range 0 <i<π/ 2. With this convention v_θ> 0 always, while v_R> 0 (v_R< 0) indicates outflow (inflow) for counterclockwise rotation and inflow (outflow) for clockwise rotation, respectively. We note that in the case of NGC 1068, where the receding side of the major axis is to the west, rotation should be counterclockwise in the likely scenario where spiral arms are trailing (e.g., Contopoulos 1971; Toomre 1981; Binney & Tremaine 1987).

Alternatively, the line-of-sight velocity field can be divided into a number of elliptical ring profiles defined by (PA, i) for a given radius r, and v_los(r,ψ) is decomposed as a Fourier series with harmonic coefficients c_j(r) and s_j(r), where $\begin{array}{ccc}{\mathit{v}}_{\mathrm{los}}\mathrm{\left(}\mathit{r,\psi }\mathrm{\right)}\mathrm{=}{\mathit{c}}_{\mathrm{0}}\mathrm{+}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{n}}\mathrm{\left[}{\mathit{c}}_{\mathit{j}}\mathrm{\right(}\mathit{r}\mathrm{\left)}\cos \mathrm{\right(}\mathit{j\psi }\mathrm{)}\mathrm{+}{\mathit{s}}_{\mathit{j}}\mathrm{(}\mathit{r}\mathrm{\left)}\sin \mathrm{\right(}\mathit{j\psi }\mathrm{\left)}\mathrm{\right]}\mathit{.}& & \end{array}$(2)In the most general case, c₁ reflects the contribution from circular rotation while all remaining terms represent contributions from noncircular motions (Schoenmakers et al. 1997; Schoenmakers 1999). It is known that expanding the series of Eq. (2) out to n = 3 provides a fair description of v_los in most models (Trachternach et al. 2008). In a disk with simple (axisymmetric) circular rotation (v_c), c₀ = v_sys and c₁ = v_csini and all remaining terms can be neglected. In the case of a pure (axisymmetric) radial flow (v_R), c₀ = v_sys and s₁ = v_Rsini, with the rest of coefficients being 0.

We derived the Fourier terms that best describe the v_los extracted from the CO(3–2) mean-velocity field presented in Sect. 5. We adopted an iterative three-step process using the software package kinemetry developed by Krajnovic et al. (2006). We used in the fit 28 ellipses with semi major axes covering the disk from to r = 24″, with a spacing Δr = 0.5″ in the inner 4″ and Δr = 1″ farther out. At each step, kinemetry finds the best fitting ellipses by minimizing the contribution of noncircular motions (v_nonc), evaluated as ${\mathit{v}}_{\mathrm{nonc}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\sqrt{{\mathit{s}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}{\mathit{s}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}{\mathit{c}}_{\mathrm{2}}\mathrm{(}\mathit{r}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{s}}_{\mathrm{3}}^{\mathrm{2}}\mathrm{(}\mathit{r}\mathrm{)}\mathrm{+}{\mathit{c}}_{\mathrm{3}}^{\mathrm{2}}\mathrm{(}\mathit{r}\mathrm{)}}\mathit{.}$(3)In a first step we left the position angle PA(r) and the inclination i(r) as free parameters in the fit and assumed that the dynamical center coincides with the AGN. We then derived the average values of PA(r) and i(r): ⟨ PA ⟩ = 290 ± 5° and ⟨ i ⟩ = 41 ± 2° (excluding outliers). In a second step, we fixed ⟨ i ⟩ = 41 ± 2° and re-determined the PA(r) profile. This profile shows a changing PA value from the CND region (PA ~ 330° ± 10° at r< 5″) to the bar (PA ~ 290° ± 10° at 5″<r< 10″) and the spiral arm region (PA ~ 260° ± 10° at r> 12″). From this second iteration, we derived ⟨ PA ⟩ = 289 ± 5°. These values of ⟨ PA ⟩ and ⟨ i ⟩ are compatible within the errors with previous determinations of the disk orientation available in the literature (PA = 286 ± 5° and i = 40 ± 3°, e.g., see compilation by Bland-Hawthorn et al. 1997). Finally, we derived the Fourier decomposition of the NGC 1068 velocity field fixing PA = 289 ± 5° and i = 41 ± 2° at all radii. Implicit in this assumption is that the gas kinematics at all radii can be satisfactorily modeled by coplanar orbits (alternative scenarios are described in Sect. 6.4). As an output parameter of this final iteration we obtained the best fit for the systemic velocity, v_sys(HEL) = 1127 ± 3 km s^-1. This coincides within the errors with the value of v_sys inferred from the kinematics of the water maser emission detected in the central r ~ 0.7 pc, as derived from the Very Long Base Interferometry (VLBI) observations of Greenhill & Gwinn (1997).

Figures 12a and b compare the magnitude of the (projected) circular component of the CO(3–2) velocity field (v_circ = c₁) obtained in the fit with that of the non-circular motions (v_nonc) defined in Eq. (3). While the value of v_nonc/v_circ stays within the range ~0.2–0.6 over a sizable fraction of the outer disk (3″<r< 20″), reflecting the expected order of magnitude of the contribution from the bar and the spiral structure to v_nonc in this region, the same ratio reaches extremely high values (~0.8–1.30) in the inner r ≤ 3″ of the galaxy. As further illustrated in Fig. 12c, the main contribution to v_nonc comes from the s₁ term. In particular, the sign (>0) and normalized strength of s₁ (s₁/c₁> 0.3) indicates that strong outflow motions prevail at r ≤ 3″. The s₁ term changes sign at r ~ 8″ and stays moderately strong and negative (− 0.5 <s₁/c₁< − 0.1) farther out in the disk; this behavior fairly reflects the expected influence of the bar and the spiral structure on the gas flow provided that we are inside corotation of the relevant perturbation (bar or spiral) at these radii (Schinnerer et al. 2000; Wong et al. 2004; Emsellem et al. 2006).

Figure 13, where we compare the s₁ and s₃ terms of the Fourier decomposition of v_los, also supports that bar-and-spiral induced streaming motions (inside corotation) are the simplest explanation for the s₁ profile in the outer disk: the dominance of the s₁ term over the s₃ term indicates that we remain inside corotation in this region. This supports the conclusion that the molecular gas pseudo-ring at r ~ 18″ (1.3 kpc) is not part of the inner bar pattern, but that it constitutes an independent wave feature characterized by a lower pattern speed (e.g., see Rand & Wallin 2004). By contrast, this mostly excludes the streaming motion scenario to explain the strong outward radial motions identified in the inner disk, which we rather interpret as a molecular outflow: values of s₁ ≫ 0 would be prevalent outside corotation for the assumed geometry of NGC 1068.

Figure 14 shows the residual mean-velocity field (⟨ V_res ⟩) obtained after subtraction of the rotation component derived in the analysis above. The residuals clearly show a pair of approaching-receding regions in the outer disk. This morphology follows the expected 2D-pattern produced by the combined action of the bar and the spiral on the gas flow provided that we are inside corotation of the perturbations on these scales (Canzian 1993; Sempere et al. 1995; Colombo et al. 2014). This is in agreement with the predominance of an inward radial flow that characterizes the fit to non-circular motions in the outer disk.

Closer to the nucleus the comparison between the residual velocity field, shaped by outward radial motions, and the gas/dust distribution in the inner disk suggests that a significant fraction of the gas in this region is associated with the outflow, as shown in Figs. 15a and b. Most remarkably, Figs. 15c and d show a noticeable spatial correlation between the AGN ionized nebulosity, traced by Paα emission, the radio jet plasma, traced by the radio continuum emission at 22 GHz, and the molecular outflow signature identified in the CO velocity field. This close association, which applies to a large range of radii and to a wide angle in the disk, suggests that a sizable fraction of the dense molecular gas traced by CO(3–2) is being entrained due to AGN feedback in the CND, and farther out in the disk, in the bow-shock arc region (located at r ~ 400 pc).

6.1.2. Mass outflow rate

To derive the mass load of the outflow we need to estimate the characteristic mass (M_mol) as well as the projected values of the radial size (R_out) and velocity (V_out) of the outflow, and also assume a certain geometry (i.e., a certain angle α between the outflow and the line of sight). If we assume a multi-conical outflow uniformly filled by the outflowing clouds (Maiolino et al. 2012; Cicone et al. 2014), the mass load rate (dM/ dt) can be estimated from the expression $\frac{\mathrm{d}\mathit{M}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{3}\mathrm{\times }{\mathit{V}}_{\mathrm{out}}\mathrm{\times }{\mathit{M}}_{\mathrm{mol}}\mathit{/}{\mathit{R}}_{\mathrm{out}}\mathrm{\times }\tan \mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathit{.}$(4)As discussed by Maiolino et al. (2012), if instead of a homogeneous multi-conical outflow, we assume a single shell-like geometry of thinness dR_out (where generally dR_out ≪ R_out), the above estimate should be replaced by $\frac{\mathrm{d}\mathit{M}}{\mathrm{d}\mathit{t}}\mathrm{=}{\mathit{V}}_{\mathrm{out}}\mathrm{\times }{\mathit{M}}_{\mathrm{mol}}\mathit{/}\mathit{d}{\mathit{R}}_{\mathrm{out}}\mathrm{\times }\tan \mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathit{.}$(5)In the following we use Eq. (4), as it provides a more conservative lower limit to the outflow rate estimated globally over the CND.

The mass (M_mol) has been calculated from the CO(3–2) data cube, after subtraction of the projected rotation curve (derived in Sect. 6.1), by integrating the emission of the line outside a velocity range ⟨ V_res ⟩ = [− 50, +50] km s^-1, which encompasses most of the expected virial motions around rotation. We determine that ≃50% of the total CO(3–2) flux in the CND stems from the outflow component. A similar percentage is derived for the bow-shock arc region, where we also find the signature of the outflow at larger radii. This translates into a total molecular mass ${\mathit{M}}_{\mathrm{mol}}\mathrm{~}\mathrm{1.}{\mathrm{8}}_{-1.1}^{\mathrm{+}\mathrm{0.6}}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}{\mathit{M}}_{\mathrm{\odot }}$ for the CND, including the mass of helium, if we assume a CO conversion factor ~$\mathrm{1}\mathit{/}\mathrm{\left(}{\mathrm{4}}_{-1}^{\mathrm{+}\mathrm{6}}\mathrm{\right)}$ of the MW value (see discussion in Sect. 5.4). We note that this estimate is rather conservative as the X_CO conversion factor for the fraction of molecular gas that participates in the outflow could be higher if this component consists of an ensemble of optically thick dense clumps embedded in a diffuse medium.

The average projected radial extent of the outflow in the CND is R_out ~ 1.5″ (100 pc) and the projected radial velocities V_out are close to ~100 km s^-1 according to the residual velocity field shown in Figs. 14–16.

The implied outflow rate given by Eq. (4) is $\mathrm{d}\mathit{M}\mathit{/}\mathrm{d}\mathit{t}\mathrm{~}\mathrm{5}{\mathrm{4}}_{-32}^{\mathrm{+}\mathrm{18}}\mathrm{\times }\tan \mathrm{\left(}\mathit{\alpha }\mathrm{\right)}{\mathit{M}}_{\mathrm{\odot }}$ yr^-1 , where α reflects the unknown angle between the outflow and the line-of-sight. If we assume that the outflow is coplanar with the main disk, tan(α) = 1 / tan(i), where i = 41°, then $\mathrm{d}\mathit{M}\mathit{/}\mathrm{d}\mathit{t}\mathrm{~}\mathrm{6}{\mathrm{3}}_{-37}^{\mathrm{+}\mathrm{21}}{\mathit{M}}_{\mathrm{\odot }}$ yr^-1. This is significantly above the mass load rate estimated on similar spatial scales for the ionized gas outflow: ~9 M_⊙ yr^-1 (Müller-Sánchez et al. 2011). The molecular mass load rate implies a very short gas depletion timescale of ≤1 Myr in the CND.

Krips et al. (2011) used a simplified model to fit by eye the CO(3–2) spectra of the CND obtained by the SMA. While based on a different approach, Krips et al. (2011) also concluded that a significant fraction (≥30%) of the molecular gas in the CND could be participating in an outflow, in qualitative agreement with our findings. Similarly, Davies et al. (2008) argued that the kinematics of molecular gas at the CND traced by the 2.12 μm H₂ line are suggestive of an outflow with typical projected velocities V_out ~ 100 km s^-1, i.e., in quantitative agreement with the values derived in this work. While the NIR line traces a minor fraction of the total H₂ content and a distinctly different molecular gas phase to that probed by the ALMA data, the similar picture drawn from these different data sets supports the outflow scenario in the CND.

The same estimate applied to the outflow detected in this work in the bow-shock arc region, with V_out ~ 75 km s^-1, ${\mathit{M}}_{\mathrm{mol}}\mathrm{~}{\mathrm{9}}_{-5.4}^{\mathrm{+}\mathrm{3}}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}{\mathit{M}}_{\mathrm{\odot }}$ and R_out ~ 5″ (~400 pc), gives $\mathrm{d}\mathit{M}\mathit{/}\mathrm{d}\mathit{t}\mathrm{~}{\mathrm{6}}_{-3.6}^{\mathrm{+}\mathrm{2}}{\mathit{M}}_{\mathrm{\odot }}$ yr^-1 for the uniformly filled model. The bow-shock arc component was not detected in the SMA map of Krips et al. (2011), nor is there a signature of the outflow on these scales in the VLT map of Davies et al. (2008).

6.2. The molecular outflow in dense gas tracers (n(H₂) ≥ 10⁵⁻⁶ cm^-3)

In the following we will show that the analysis of the kinematics derived from CO(6–5), HCN(4–3), and HCO⁺(4–3) confirms that the molecular outflow detected in the CO(3–2) line has a higher density (n(H₂) ≥ 10⁵⁻⁶ cm^-3, according to Paper II) counterpart in NGC 1068.

As discussed in Sect. 5, the mean-velocity field of the CND shows in all tracers a similar pattern: the kinematic major axis of the CND is tilted by ≥ 40° relative to that of the large-scale disk traced by CO(3–2), as shown in Fig. 11. This distinct kinematic feature has been modeled as an outflow in CO(3–2) (see Sect. 6.1), and as such can be similarly interpreted as a signature of outflow in the dense gas tracers. Figure 16 shows the position-velocity (p-v) plots taken along the minor axis, determined in Sect. 6.1 (PA = 19°), for various dense molecular tracers. In these diagrams any deviation of the emission from v_sys beyond the virial range, determined by the expected cloud-cloud velocity dispersion, is indicative of radial inflow/outflow motions. An inspection of this figure shows that the outflow signatures identified in CO(3–2) out to radii r ~ 400 pc are echoed in CO(6–5), HCN(4–3), and HCO⁺(4–3) on the scales of the CND (r ~ 100 − 200 pc). A sizable fraction of the emission of these molecular tracers lies outside the expected range of virial motions attributable to rotation and dispersion: on average, emission is ≥50 km s^-1 – redshifted on the northern side of the CND, while it is ≥100 km s^-1 – blueshifted on the southern side. This reflects the sign and the right order of magnitude of the mean-velocity field deviations seen in the CO(3–2) map of Fig. 15 at the CND.

The limited velocity coverage of the CO(3–2) data at the blue velocity end (⟨ V_res ⟩ < − 170 km s^-1) implies that we underestimate the most extreme outflow velocities on the southern side of the CND in this line. We note that outflow velocities up to ⟨ V_res ⟩ ~ + 280 km s^-1 are detected in CO(3–2) on the northern side. The reality of these highly redshifted velocities in the outflow is confirmed by the HCN(4–3) data, which show emission from ⟨ V_res ⟩ = − 260 to +260 km s^-1 in the CND, as shown in Fig. 16. This implies that the outflow mass load rate estimated in Sect. 6.1.2 from CO(3–2) should be taken as a lower limit.

6.3. The powering source of the molecular outflow: star formation or AGN jet driven?

Evidence of a young stellar population in the CND has been recently found by Storchi-Bergmann et al. (2012), who located the young star formation (SF) episode inside the expanding molecular ring at r ~ 100 pc. While the presence of young stars is well established, the total energy possibly injected by SF in the CND is much lower compared to the AGN. Davies et al. (2007) (see also discussion in Hailey-Dunsheath et al. 2012) estimated that only a fraction of the total FIR continuum luminosity (L_FIR) from the CND can be attributed to ongoing/recent SF. Based on the luminosity measured in the K-band in the inner r ~ 35 pc, Davies et al. (2007) conclude that L_FIR due to SF is ~(1.7–3) × 10⁹L_⊙; this implies an integrated star formation rate of SFR ~ 0.4 − 0.7 M_⊙ yr^-1. This estimate is similar to the nuclear SFR measured from the 11.3 μm PAH luminosity in the inner r ~ 12 pc of NGC 1068 by Esquej et al. (2014), who estimated that SFR_nuclear ~ 0.4 M_⊙ yr^-1. The SFR estimated by Esquej et al. (2014) for the circumnuclear region out to a radius r ~ 2″ (140 pc) is about 1 M_⊙ yr^-1. The total SFR in the CND is therefore about an order of magnitude lower than the estimated mass load rate of the molecular outflow. Taken at face value this discrepancy suggests that SF is not able to drive the molecular outflow in the CND of NGC 1068 (see Murray et al. 2005 and Veilleux et al. 2005 for a general discussion).

The kinetic luminosity of the CND outflow can be derived from the expression ${\mathit{L}}_{\mathrm{kin}}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{2}\mathrm{\times }\frac{\mathrm{d}\mathit{M}}{\mathrm{d}\mathit{t}}\mathrm{\times }{\left(\frac{{\mathit{V}}_{\mathrm{out}}}{\cos \mathrm{\left(}\mathit{\alpha }\mathrm{\right)}}\right)}^{\mathrm{2}}\mathrm{·}$(6)Assuming that the gas in the outflow is coplanar, ${\mathit{L}}_{\mathrm{kin}}\mathrm{~}{\mathrm{5}}_{-3}^{\mathrm{+}\mathrm{1.7}}\mathrm{\times }{\mathrm{10}}^{\mathrm{41}}$ erg s^-1. AGN feedback models require that a significant fraction of the radiated luminosity (~5%L_bol; di Matteo et al. 2005) should be coupled to the ISM to produce an outflow. This fraction is lowered to ~0.5% in the two-phase feedback model of Hopkins & Elvis (2010a). The bolometric luminosity of the AGN in NGC 1068 (L_bol), estimated from MIR and X-ray wavelengths, is at least three orders of magnitude larger than L_kin: L_bol ≥ 10^{44 − 45}erg s^-1 (Bock et al. 2000; Laurent et al. 2000; Matt et al. 2000; Raban et al. 2009; Prieto et al. 2010; Alonso-Herrero et al. 2011 and Sect. 4.1.2). This result indicates that the AGN can power the outflow in the CND of NGC 1068.

The momentum flux of the CND outflow can be computed from the expression $\frac{\mathrm{d}{\mathit{P}}_{\mathrm{out}}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathrm{d}\mathit{M}}{\mathrm{d}\mathit{t}}\mathrm{\times }\frac{{\mathit{V}}_{\mathrm{out}}}{\cos \mathrm{\left(}\mathit{\alpha }\mathrm{\right)}}\mathrm{·}$(7)In the case of coplanar gas, Eq. (7) yields $\mathrm{d}{\mathit{P}}_{\mathrm{out}}\mathit{/}\mathrm{d}\mathit{t}\mathrm{~}{\mathrm{6}}_{-3.6}^{\mathrm{+}\mathrm{2}}\mathrm{\times }{\mathrm{10}}^{\mathrm{34}}$g cm s^-2. Compared to the momentum provided by the AGN photons, derived as L_bol/c ~ (0.3 − 3) × 10³⁴g cm s^-2, dP_out/ dt is a moderate factor 1–27 larger. The required boost is only ~1.7–6 if we adopt L_bol = 4.2 × 10⁴⁴erg s^-1 from Sect. 4.1.2. The molecular outflow could thus be in the momentum-conserving regime (where cooling is fast), and radiation pressure could be the driving mechanism. This is also within the range of values of the momentum boost factors predicted by AGN feedback models under the assumption that molecular outflows are energy-conserving (adiabatic Sedov phase): (dP_out/ dt) / (L_bol/c) ~ 10 − 50 (Faucher-Giguère & Quataert 2012).

Alternatively, the radio jet of NGC 1068 could also inject the required outflow power. The jet power (W_jet) can be estimated from the monochromatic luminosity at 1.4 GHz, according to Bîrzan et al. (2008). In the case of NGC 1068, the spatially integrated flux of the jet at 1.4 GHz, including the main components (known in the literature as NE, C, and S), is ~840 mJy from the radio continuum map of Gallimore et al. (1996). This implies W_jet = 1.8 × 10⁴³ erg s^-1. Since W_jet ~ (30 − 100) × L_kin, we conclude that the jet can drive the molecular outflow in the CND even assuming a low coupling efficiency.

6.4. Alternatives to the molecular outflow

The velocity field of the outer disk from r ~ 400 pc out to ~1.8 kpc shows a regularly rotating pattern with conspicuous distortions that reflect the expected perturbation on the gas flow due to the bar and the spiral structure inside corotation of the patterns (Sect. 6.1.1). While the outer disk and, also, the molecular gas component detected at the AGN (see Sect. 5.3) seem to share the latter as the common dynamical center, the CND is a strongly off-centered ring, which is apparently rotating but with a different kinematic axis, as shown in Fig. 10. As an alternative to the outflow scenario described in Sects. 6.1 and 6.2, the abrupt shift of ≥40° degrees in the kinematic PA of the two disks can be interpreted as due to the CND being a non-coplanar disk which is far from being dynamically relaxed. Two types of non-coplanar instabilities can be invoked to account for the decoupled kinematics observed in the CND: a nuclear warp or a non-coplanar lopsided instability.

6.4.1. A nuclear warp

Schinnerer et al. (2000) explored a nuclear warp scenario in an attempt to model the kinematics of molecular gas in the nuclear region of NGC 1068, based on CO(2–1) observations obtained with the PdBI. They concluded that gas motions could be equally fit with either a warp or a nuclear bar in the CND. One of the main predictions of their warp model was that the CO disk should become edge-on at a radius of ~70 pc. However, we do not see the signature of an edge-on disk on these spatial scales in the higher-resolution ALMA maps presented in this work. Furthermore the internal kinematics of the CND do not show the typical S-shaped distortion attributable to a nuclear warp instability.

The diagram shown in Fig. 13, originally introduced by Wong et al. (2004), can also be used as a diagnostic tool to identify nuclear warp signatures in the velocity field of galaxies in the particular case where the hypothesized tilted orbits share a common center. An inspection of Fig. 13 indicates that the continuous line, which represents the least-squares fit to the NGC 1068 data points, shows no correspondence to the expected location of the warp line, which would relate the s₁ and s₃ terms in the case of a warp having the AGN as a common center of the tilted non-coplanar orbits. This disagreement, more evident in the inner disk, leads us to conclude that the simplest nuclear warp scenario described above is not a satisfactory explanation for the velocity residuals observed in this region.

6.4.2. A non-coplanar lopsided instability

If we abandon the restriction of having the AGN as common orbiting center for both the CND and the outer disk, the scenario of non-coplanarity remains nevertheless viable. In this scenario, a non-coplanar CND would be orbiting around a secondary nucleus at (Δα, Δδ) ~ ( − 1.5″, − 1″), i.e., 1″= 70 pc offset to the southwest relative to the AGN.

Two types of mechanisms have been described as potential triggers of lopsided non-coplanar gas instabilities in galactic nuclei:

• 1.

  External trigger: the postulated secondary nucleus could bethe footprint of a recent minor merger with a nucleated satellite. Inthis scenario dynamical friction would make the satellite quicklysink toward the nucleus of the host dragging the accreted gasdisk to the central region (Taniguchi & Wada 1996; Taniguchi1999; 2013). Minor mergers could thus explaina random orientation of circumnuclear gas disks in the NLRof Seyfert galaxies, as the orbital plane of the resulting nucleardisk would be determined by the orbital parameters of the ac-creted satellite (Taniguchi 2013). As the CND is noticeablyoff-centered, this would imply that NGC 1068would be at an early stage in the merging process where thesecondary nucleus has not yet sunk toward the nucleus. How-ever, the absence of any signature of a secondary nucleus in thehigh-resolution HST/NICMOS images of the CND, whichshould be at the position of the apparent guiding center of the non-coplanar disk, makes the minor merger scenario implausible inNGC 1068.

• 2.

  Internal trigger: alternatively, a lopsided instability triggered on the gas could have an internal origin (Jog & Combes 2009). The CND could be the result of a recent episode of gas infall produced by the feedback of a star formation burst in the disk. This feedback action would be able to eject gas perpendicular to the plane, at heights of 100–200 pc, as frequently seen in numerical simulations (e.g., Emsellem et al., in prep.). The gas fountain would fall back while settling initially in an inclined disk, the perpendicular orientation being preferred since differential precession is then canceled. This would explain the offset between the AGN and the guiding center of the instability. While settling back to the disk, the nuclear gas would trigger a slow and long-lasting lopsided perturbation, which could remain during more than 100 dynamical times at this radius, i.e., 400 Myr for the CND (e.g., Jog & Combes 2009). In this scenario it is difficult to predict what evolutionary stage during the settling of the disk best represents the case of NGC 1068. However, at the end of the process, the m = 1 instability and the main molecular disk would become coplanar.

As will be discussed in Sect. 6.5, the orientation of the jet and the geometry of the ionized gas outflow in NGC 1068 imply that these can efficiently interact with molecular gas in a coplanar CND. Overall, the outflow hypothesis is the simplest explanation for the distorted kinematics of molecular gas at the CND, but also for the bow-shock arc region, the signature of the interaction at larger radii. In this scenario the remarkable off-centering of the molecular ring would be triggered before the jet or the ionized gas wind hits the disk and would launch the outflow rather than being the result of it. Furthermore, the lopsided morphology in a coplanar CND, which is characterized by a strong near side/far side asymmetry in NGC 1068, could also be enhanced due to opacity effects in molecular lines, as discussed by Boone et al. (2011).

6.5. The molecular outflow in context: relation to other tracers

NGC 1068 harbors a wide-angle (FWHM ~ 50 − 60°) biconical outflow of ionized gas (Macchetto et al. 1994; Arribas et al. 1996; Crenshaw et al. 2000; Tecza et al. 2001; Cecil et al. 2002; Mueller-Sánchez et al. 2011). The orientation of the outflow is not perpendicular to the galaxy disk: i_outflow ~ 70 − 80°, while i_disk ~ 41°. This particular geometry favors interaction with the molecular disk out to a galactocentric radius of r ~ 400 pc. An interaction between the jet and the ISM was already identified close to the AGN (r ~ 20 − 40 pc) by Gallimore et al. (1996; 2001) (see also Bicknell et al. 1998).

Figure 17 shows a simple kinematic model of the NLR, first proposed by Tecza et al. (2001) (see also: Cecil et al. 2002; Das et al. 2006; 2007), here revised to account for the molecular outflow detected by ALMA. The figure shows a scaled crosscut of the NLR along the projected direction of the radio jet (PA ~ 30°). The observer’s line-of-sight, which makes an angle i_disk ~ 41° relative to the direction orthogonal to the molecular disk, as shown, is at 90° with respect to the viewer of this figure. The northeast side of the ionized gas cone is at the upper right of the figure. The CND, represented by two molecular knots located asymmetrically north and south of the AGN (N and S in figure), is embedded in the molecular disk.

According to Tecza et al. (2001) (see also Cecil et al. 2002 and Mueller-Sánchez et al. 2011), high-ionization lines, like [SiVI] or [OIII], are produced where the radio lobes encounter diffuse material located outside the galaxy plane. The interaction of the radio lobes with this medium generate two bow-shock fronts that are responsible for the simultaneous detection of redshifted and blueshifted emission in both ionization cones, as shown in Fig. 17. The low-ionization lines, like [FeII], are produced when the bow-shock sweeps the denser molecular disk. This explains why low-ionization lines are exclusively detected as redshifted (blueshifted) emission on the northern (southern) side of the disk. The amplitude of velocity shifts are larger for the high-ionization lines, reaching ~3000 km s^-1 in the northeastern cone, while they remain ≤500 km s^-1 for the low-ionization lines. Furthermore, the sign of the velocity shifts for the strongest emission components of the high-ionization lines are noticeably reversed with respect to low-ionization lines.

The molecular outflow detected by ALMA stands out as a redshifted (blueshifted) kinematic component on the northern (southern) side of the disk, as discussed in Sects. 6.1 and 6.2, consistent with the low-ionization lines of Tecza et al. (2001). In this scenario, depicted in Fig. 17, the molecular outflow is launched when the ionization cone of the NLR sweeps the disk in the CND and farther out north in the bow-shock arc region where Wilson & Ulvestad (1987) found that the radio-lobe nebulosity becomes limb-brightened and highly polarized, the signature of a bow-shock in the disk. The amplitudes of the velocity shifts in the molecular outflow are significantly smaller than for the high-ionization lines. There is also evidence that the molecular outflow becomes decelerated at r ~ 400 pc: the terminal velocities of the outflow in the bow-shock arc region are a factor 2 smaller than in the CND, as shown in Fig. 16.

The detection of emission from dense gas tracers at the extreme velocities of the molecular outflow in NGC 1068 (≥150 km s^-1) raises the question of how molecular clouds can survive during the launching of the outflow, since shocks of velocities greater than ~ 50 km s^-1 are known to be fast enough to destroy molecules (Hollenbach & McKee 1989; Neufeld & Dalgarno 1989a,b). However, the chemistry in both dissociative J-type and non-dissociative C-type shocks, which involve dust-mantle disruption and a high-temperature environment, favors an efficient fast reformation of molecules in the gas. Molecules can reform in the post-shocked gas on timescales ≤hundreds of yrs. There is ample observational evidence that the emission of dense gas tracers (≥10^{5 − 6} cm^-3) can coexist with fast shocks (~50–100 km s^-1) in galactic bipolar outflows of Young Stellar Objects (YSOs) (e.g., Arce et al. 2007 and references therein). These observations show that the abundance of some molecular species can undergo spectacular enhancements in the outflow gas, due to the onset of shock chemistry. Besides classical shock tracers such as SiO, whose abundance can be enhanced by factors of about 10⁴, molecular species which are less specific to shock environments, such as HCN or CS, can also undergo significant order-of-magnitude enhancements (e.g., Tafalla et al. 2010; Tafalla 2013). The detection of strong emission from some of these tracers, such as HCN or CS, at the velocities identified as abnormal in NGC 1068 suggests that molecular clouds survive in the outflow and that shock chemistry is likely at work in this component.

7. Molecular line ratios and environment

7.1. The SB ring versus the CND

Figure 18 shows the overlay of the HST Paα emission map on the CO(3–2)/CO(1–0) (hereafter R32 / 10) and CO(6–5)/ CO(3–2) (hereafter R65 / 32) line ratio maps. Line ratios were derived in T_mb units. To derive the 3–2/1–0 brightness temperature ratio, shown in panels a and b, we degraded the 3–2 map to the spatial resolution of the 1–0 observations of Schinnerer et al. (2000) (). The CO(6–5)/CO(3–2) brightness temperature ratio map, shown in panel c, was derived at the spatial resolution of the 3–2 ALMA observations (). Line ratios were obtained assuming a common 3σ clipping on the integrated intensities to assure image reliability.

The R32 / 10 ratio changes significantly depending on the particular environment of the disk. The average ratio is ~1.2 ± 0.02 in the SB ring yet R32 / 10 shows a large dispersion of values: R32 / 10 ~ 0.7 − 1 in the less actively star-forming regions of the SB ring, which are characterized by low or undetected Paα emission, while the strongest Paα emitting regions show ~a factor 3–4 higher ratios (R32 / 10 ~ 2 − 3).

Overall, the R32 / 10 ratio in the CND, ~2.7 ± 0.1, is higher than in the SB ring, in agreement with the global values measured with the PdBI and the SMA by Krips et al. (2011) and Tsai et al. (2012). In summary, the R32 / 10 ratios measured in NGC 1068, most particularly those of the CND, are at the high end of the typical values observed in the centers of nearby normal and starburst galaxies: ~0.5–1.4 (Devereux et al. 1994), ~0.2–0.7 (Mauersberger et al. 1999), ~0.2–1.9 (Mao et al. 2010), an indication that the excitation of molecular gas is extreme in the CND (Krips et al. 2011; Paper II).

7.2. Line ratio changes in the CND: the footprint of AGN feedback?

At the spatial resolution used in Figs. 18a and b, the R32 / 10 ratio shows a range of values in the CND: R32 / 10 goes from ~1 at the southern end of the ring to ~5 at the northern end, and it is ~3–4 at the AGN and at the E and W CO knots. The R32 / 10 ratio shows a north-south gradient which runs in parallel with the orientation of the bipolar AGN/jet nebulosities. Overall, the R65 / 32 ratio in the CND is ~0.8, with a significant range of values: R65 / 32 goes from ~0.7 at the southern end of the CND ring to a maximum value ~2.0 at the AGN. At the higher spatial resolution of Fig. 18c, the excitation of CO in the CND probed by the R65 / 32 ratio unveils the footprint of AGN feedback: the R65 / 32 ratio is enhanced by the degree of illumination of the molecular gas by the photons of the ionized gas outflow traced by Paα emission.

Other molecular line ratios show dramatic changes of up to an order of magnitude inside the CND on the spatial scales probed by ALMA (~35 pc). Figure 19 shows the overlay of the HST Paα emission map on a set of four molecular line ratio maps (in T_mb units) obtained at the resolution of ALMA observations in Band 7: (a) HCN(4–3)/ CO(3–2) (R_{HCN / CO}); (b) HCO⁺(4–3)/CO(3–2) (R_{HCO⁺/ CO}); and (c) HCN(4–3)/HCO⁺(4–3) (R_{HCN / HCO⁺}).

A visual inspection of Fig. 19 indicates that, as for CO, the R_{HCN / CO} and R_{HCO⁺/ CO} ratios also show the footprint of AGN feedback: these line ratios go from ~0.02–0.03 at the outer radii of the CND to ~0.8 at the AGN, showing an enhancement which runs in parallel with the irradiation of molecular gas by the photons of the bipolar ionized gas nebulosity.

The R_{HCN / HCO⁺} ratio is globally very high in the CND with an average value of ~2.5. These high values have been widely used as a diagnostic ratio to identify an AGN-like environment (Kohno et al. 2001; Usero et al. 2004; García-Burillo et al. 2006; Graciá-Carpio et al. 2006; 2008; Imanishi et al. 2007; 2013; Krips et al. 2008; 2011). High HCN/HCO⁺ ratios measured in galaxy nuclei have also been interpreted as the signature of significant mechanical heating in shock/mechanically dominated regions (MDRs) (Kazandjian et al. 2012). Nevertheless, the ALMA maps show that the lowest value of this ratio is found precisely at the AGN locus, where R_{HCN / HCO⁺} ~ 1.3, an indication that a simplistic interpretation of these ratios can be misleading (see Paper II). The comparatively lower value of R_{HCN / HCO⁺} at the AGN might reflect departures from chemical equilibrium: the HCN/HCO⁺ mass ratio is strongly time-dependent with variations that can reach orders of magnitude for times ≥10⁴ years (Meijerink et al. 2013).

As discussed in García-Burillo et al. (2010), the morphology of hard-X ray emission in the 6–8 keV band obtained by Chandra (Young et al. 2001; Ogle et al. 2003) can be used to search for observational evidence of AGN feedback on the excitation/chemistry of molecular gas in the CND. Emission in the 6–8 keV band observed in NGC 1068 is dominated by reflection of X rays by cold neutral (presumably mostly molecular) gas (Iwasawa et al. 1997). Figure 20 shows the R_{HCN / CO} and R_{HCO⁺/ CO} ratios in the CND as a function of the irradiation of molecular gas by hard X-rays normalized by gas column density; the latter is derived from the ratio of X-ray flux in the 6–8 keV band to the CO(3–2) intensities (X_hard/CO). All quantities have been measured at a common spatial resolution of 0.5″. With an admittedly large scatter, Fig. 20 provides quantitative evidence of a correlation between the R_{HCN / CO} and R_{HCO⁺/ CO} ratios and X_hard/CO. This is similar to the correlation found for the SiO/CO and CN/CO ratios as function of hard X-ray illumination of molecular gas in the CND, discussed by García-Burillo et al. (2010).

The existence of an AGN-driven outflow proves that the feedback of activity is shaping the gas kinematics in the inner r ~ 400 pc of NGC 1068. The change in molecular line ratios in the CND described above can be taken as independent evidence that radiative and/or mechanical feedback from the AGN is at work. The feedback of activity on the excitation and chemistry of molecular gas can be radiative, in PDR/XDR environments, or mechanical, in MDRs. In Paper II we use the line ratio maps derived in this work, complemented by those derived from interferometric observations of NGC 1068 done in other molecular tracers, to model the excitation and chemistry of molecular gas in the different environments of the CND and the SB ring and find a link between the different sources of energy present in the disk.

8. Summary and conclusions

We have used ALMA to map the emission of a set of dense molecular gas tracers (CO(3–2), CO(6–5), HCN(4–3), HCO⁺(4–3), and CS(7–6)) and their underlying continuum emission in the central r ~ 2 kpc of NGC 1068 with spatial resolutions ~0.3″ − 0.5″ (20–35 pc). These observations have greatly improved the sensitivity and spatial resolution of any previous interferometric study of this Seyfert, allowing us to make a significant progress in our understanding of the fueling and the feedback of activity in NGC 1068.

We summarize the main results of our study as follows:

• The CO(3–2) line and underlying continuum image the distribution of the dense molecular gas and dust emission from three main regions: the CND, the bar, and the SB ring. We also report the detection of CO(3–2) emission at different locations over the interarm region.

• The CND, which is fully resolved in the CO(6–5) map, is a highly-structured closed elliptical ring of 350 pc-size. Similar to the dust continuum ring, the CO ring is noticeably off-centered relative to the AGN.

• Outside the CND, the CO(3–2) emission is detected along two elongated offset lanes that run mostly parallel to the bar major axis. We also detect a gas component with anomalous velocities on the northeastern side of the disk at distances ~400 pc from the AGN. This CO feature has a counterpart in dust continuum emission and is interpreted as the signature of a bow-shock in the molecular disk.

• As for dust emission, most of the CO(3–2) flux in the disk is detected in the SB ring: a two-arm spiral structure that starts from the ends of the stellar bar and unfolds in the disk over ~180° in azimuth forming a pseudo-ring.

• In stark contrast with the CO(3–2) map, most of the emission in HCO⁺(4–3), HCN(4–3), and CS(7–6) stems from the CND. This reflects the different line ratios measured in the CND and in the SB ring. The line ratio maps at the spatial resolution of ALMA show dramatic changes of up to an order of magnitude inside the CND. These changes are tightly correlated with the UV/X-ray illumination by the AGN, which varies significantly on these spatial scales, hinting at the footprint of AGN feedback.

• We detect dust continuum and molecular line emission at the AGN. The emission peak is nevertheless shifted to the northeast in a structure that connects the AGN with the CND. We used the ALMA fluxes in Bands 9 and 7 together with NIR/MIR data to constrain the properties of the putative torus using CLUMPY models. The fitted outer torus radius is ${\mathit{R}}_{\mathrm{o}}\mathrm{=}\mathrm{2}{\mathrm{0}}_{-10}^{\mathrm{+}\mathrm{6}}\hspace{0.17em}$pc.; the torus mass is M_torus = 2.1( ± 1.2) × 10⁵ M_⊙. The ALMA fluxes lie above the model predictions, an indication that a non-negligible fraction of cold dust emission not associated with the torus was included in the 20 pc – resolution data of ALMA.

• The Fourier decomposition of the gas velocity field derived from the CO(3–2) data indicates that rotation is perturbed by an inward radial flow in the SB ring and the bar region, as expected if we are inside corotation of the perturbations in these regions. However, the kinematics of the CND and the bow-shock arc are shaped by an outward radial pattern superposed to rotation. There is a tight spatial correlation between the AGN ionized nebulosity, the radio jet plasma and the molecular outflow signature identified in the CO velocity field from r ~ 50 pc to r ~ 400 pc, which suggests that the outflow is AGN driven. The molecular outflow is launched when the ionization cone of the NLR sweeps the disk in the CND and farther out north in the bow-shock arc region.

• The estimated outflow rate in the CND, $\mathrm{d}\mathit{M}\mathit{/}\mathrm{d}\mathit{t}\mathrm{~}\mathrm{6}{\mathrm{3}}_{-37}^{\mathrm{+}\mathrm{21}}{\mathit{M}}_{\mathrm{\odot }}$ yr^-1, an order of magnitude higher than the SFR measured at these radii, confirms that the outflow is AGN driven. The molecular mass load rate implies a very short gas depletion timescale of ≤1 Myr in the CND.

The AGN-driven molecular outflow in NGC 1068 could quench star formation in the inner r ~ 400 pc of the galaxy on short timescales and at the same time regulate gas accretion in the CND. However, the molecular gas reservoir in the CND is likely replenished on longer timescales by efficient gas inflow from the outer disk. The signature of gas inflow has been identified in the velocity field at larger radii and has been attributed to the combined action of the bar and the spiral arms. Self-regulation of star formation and gas accretion would be thus possible by a combination of density wave-driven inflows and AGN-driven outflows. A similar scenario of self-regulation has been invoked to be at work in gas-rich high-redshift disk galaxies, analyzed in the numerical simulations of Gabor & Bournaud (2014).

While the gas kinematics are dominated by the molecular outflow in the CND, on theoretical grounds it is nevertheless expected that the gas should be fueling the AGN at smaller radii. The NIR data of Müller-Sánchez et al. (2009) showed elliptical streamers that bridge the CND and the central engine. They were interpreted as a tentative evidence of ongoing AGN fueling. However, with the spatial resolution of the observations presented in this paper we cannot spatially resolve the kinematics of the dense molecular gas in the inner r ~ 50 pc. Obtaining higher-resolution data is important to studying the inflow/outflow signatures in the region that connects the CND with the AGN. A significant improvement in spatial resolution would also allow us to resolve with quasi-thermal molecular lines the region where the jet is known to be interacting with the gas, image the connection to the inner r ~ 0.7 pc H₂O megamaser disk (Gallimore et al. 2001) and spatially resolve the emission from the putative molecular torus.

---

Acknowledgments

We acknowledge the staff of ALMA in Chile and the ARC-people at IRAM-Grenoble in France for their invaluable help during the data reduction process. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2011.0.00083.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. We used observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive. S.G.B. and I.M. acknowledge support from Spanish grants AYA2010-15169 and from the Junta de Andalucia through TIC-114 and the Excellence Project P08-TIC-03531. S.G.B., A.L., and A.F. acknowledge support from MICIN within program CONSOLIDER INGENIO 2010, under grant “Molecular Astrophysics: The Herschel and ALMA Era–ASTROMOL” (ref CSD2009-00038). S.G.B., A.U., L.C., and P.P. acknowledge support from Spanish grant AYA2012-32295. FC acknowledges the European Research Council for the Advanced Grant Program Num. 267399-Momentum. A.A.H. acknowledges support from the Universidad de Cantabria through the Augusto G. Linares programme and from the Spanish Plan Nacional grants AYA2009-05705-E and AYA2012-31447. C.R.A. is supported by a Marie Curie Intra European Fellowship within the 7th European Community Framework Programme (PIEF-GA-2012-327934). C.R.A. also ackowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) through project PN AYA2010-21887-C04.04 (Estallidos).

References

Appendix A: CLUMPY torus models

Over the past few years, a number of studies have demonstrated clumpy torus models (Nenkova et al. 2008a,b; Schartmann et al. 2008; Hoenig & Kishimoto 2010) represent well the nuclear infrared emission of local Seyfert galaxies, provided there is no contamination from extended dust components not related to the torus (e.g., foreground absorbing dust, optically thin emitting dust in ionization cones, etc.). More specifically, these models reproduce the non-stellar NIR and MIR SEDs and MIR spectroscopy of local Seyfert galaxies (Ramos Almeida et al. 2009; 2011a; Hoenig et al. 2010; Alonso-Herrero et al. 2011; Lira et al. 2013) and PG quasars (Mor et al. 2009). By modeling the nuclear infrared data of AGN we can constrain the torus geometry, dust properties, and distribution, and the AGN bolometric luminosity. Moreover, adding far-infrared nuclear photometric points or even upper limits has proven to be useful to constrain the torus sizes, especially if the dusty clumps have a relatively uniform radial distribution along the torus (Ramos Almeida et al. 2011b; Asensio Ramos & Ramos Almeida 2013; Ramos Almeida et al. 2014).

Alonso-Herrero et al. (2011) fitted the nuclear NIR and MIR SED and MIR spectroscopy of NGC 1068 using the so-called CLUMPY torus models of Nenkova et al. (2008a,b). The parameters in the CLUMPY models are the torus size Y defined as the ratio between the outer radius R_o and the inner radius⁵R_sub, the torus angular size σ, the viewing angle i, the number of clouds along the equatorial direction N₀, the optical depth of the clouds τ_V, and the index q of the radial distribution of the clouds ∝r^{− q}. Using a Bayesian approach to fit the data with the BayesClumpy tool (Asensio Ramos & Ramos Almeida 2009) the torus model parameters of NGC 1068 were well constrained. In their fitting they restricted the torus size to small values based on the 12 μm interferometric sizes of a few parsecs inferred for nearby AGN (Burtscher et al. 2013).

We can use the ALMA Band 9 and Band 7 continuum thermal fluxes at 435 μm and 860 μm, respectively, to investigate whether we can set further constraints on the torus properties of NGC 1068. To do so we also added the NGC 1068 SED and MIR spectroscopy presented in Alonso-Herrero et al. (2011). As in Alonso-Herrero et al. (2011) we used the detection of a maser in the nuclear region of NGC 1068, which implies a close to edge-on view of the AGN accretion disk, to set the following prior for the viewing angle i = 60 − 90°. Because we are including the ALMA far-infrared photometry we allowed the full range for the torus size Y = 5 − 100. We used the new version of the BayesClumpy tool (Asensio Ramos & Ramos Almeida 2009) which now interpolates linearly between the CLUMPY models.

Based on the fit, we can also estimate the gas mass in the torus. In the models of Nenkova et al. (2008b) the total mass in torus clouds can be written as: ${\mathit{M}}_{\mathrm{torus}}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{m}}_{\mathrm{H}}\mathrm{\left(}\sin \mathit{\sigma }\mathrm{\right)}{\mathit{N}}_{\mathrm{torus}}^{\mathrm{eq}}{\mathit{R}}_{\mathrm{sub}}^{\mathrm{2}}\mathit{Y}{\mathit{I}}_{\mathit{q}}\mathrm{\left(}\mathit{Y}\mathrm{\right)}$(A.1)where the function I_q(Y) = 1 for q = 2 and ${\mathit{N}}_{\mathrm{torus}}^{\mathrm{eq}}$ is the column density along the equatorial direction. The inner radius of the torus R_sub was computed from the AGN bolometric luminosity inferred from the fit and we assumed a gas-to-A_V ratio of N_H2/A_V = 1 × 10²¹ mol cm^-2 mag^-1, taken from Bohlin et al. (1978). The typical scatter reported by Bohlin et al. (1978) for this factor in H₂ dominated gas clouds is 30%.

The results of the fit in NGC 1068 and their implications are discussed in Sect. 4.1.2.

All Figures