Molecular line emission in NGC1068 imaged with ALMA. I An AGN-driven outflow in the dense molecular gas
S. Garcia-Burillo, F. Combes, A. Usero, S. Aalto, M. Krips, S. Viti, A. Alonso-Herrero, L. K. Hunt, E. Schinnerer, A. J. Baker, F. Boone V. Casasola, L. Colina, F. Costagliola, A. Eckart, A. Fuente, C. Henkel, A. Labiano, S. Martin, I. Marquez, S. Muller, P. Planesas, C. Ramos Almeida, M. Spaans, L. J. Tacconi, P. P. van der Werf
AAstronomy & Astrophysics manuscript no. ngc1068-paper-I-inpress-pdflatex c (cid:13)
ESO 2014June 16, 2014
Molecular line emission in NGC 1068 imaged with ALMA (cid:63)
I. An AGN-driven outflow in the dense molecular gas
S. García-Burillo , F. Combes , A. Usero , S. Aalto , M. Krips , S. Viti , A. Alonso-Herrero , L. K. Hunt ,E. Schinnerer , A. J. Baker , F. Boone , V. Casasola , L. Colina , F. Costagliola , A. Eckart , A. Fuente ,C. Henkel , , A. Labiano , , S. Martín , I. Márquez , S. Muller , P. Planesas , C. Ramos Almeida , ,M. Spaans , L. J. Tacconi , and P. P. van der Werf (A ffi liations can be found after the references) Received —; accepted —-
ABSTRACT
Aims.
We investigate the fueling and the feedback of star formation and nuclear activity in NGC 1068, a nearby ( D =
14 Mpc) Seyfert2 barred galaxy, by analyzing the distribution and kinematics of the molecular gas in the disk. We aim to understand if and how gasaccretion can self-regulate.
Methods.
We have used the Atacama Large Millimeter Array (ALMA) to map the emission of a set of dense molecular gas ( n (H ) (cid:39) − cm − ) 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 ∼ ∼ . (cid:48)(cid:48) − . (cid:48)(cid:48) ( ∼ −
35 pc for the assumed distance of D =
14 Mpc).
Results.
The sensitivity and spatial resolution of ALMA give an unprecedented detailed view of the distribution and kinematics of thedense molecular gas ( n (H ) ≥ − cm − ) in NGC 1068. Molecular line and dust continuum emissions are detected from a r ∼
200 pco ff -centered circumnuclear disk (CND), from the 2.6 kpc-diameter bar region, and from the r ∼ . + , HCN, and CS stems from the CND. Molecular line ratios show dramatic order-of-magnitude changes insidethe CND that are correlated with the UV / X-ray illumination by the AGN, betraying ongoing feedback. We used the dust continuumfluxes measured by ALMA together with NIR / MIR data to constrain the properties of the putative torus using CLUMPY modelsand found a torus radius of 20 + − pc. The Fourier decomposition of the gas velocity field indicates that rotation is perturbed by aninward radial flow in the SB ring and the bar region. However, the gas kinematics from r ∼
50 pc out to r ∼
400 pc reveal a massive( M mol ∼ . + . − . × M (cid:12) ) outflow in all molecular tracers. The tight correlation between the ionized gas outflow, the radio jet, andthe occurrence of outward motions in the disk suggests that the outflow is AGN driven. Conclusions.
The molecular outflow is likely launched when the ionization cone of the narrow line region sweeps the nuclear disk.The outflow rate estimated in the CND, d M / d t ∼ + − M (cid:12) yr − , is an order of magnitude higher than the star formation rate at theseradii, confirming that the outflow is AGN driven. The power of the AGN is able to account for the estimated momentum and kineticluminosity of the outflow. The CND mass load rate of the CND outflow implies a very short gas depletion timescale of ≤ ffi cient gas inflow from the outer disk. Key words.
Galaxies: individual: NGC 1068 – Galaxies: ISM – Galaxies: kinematics and dynamics – Galaxies: nuclei – Galaxies:Seyfert – Radio lines: galaxies
1. Introduction
The study of the content, distribution, and kinematics of inter-stellar gas is a key to understanding the origin and maintenanceof nuclear activity in galaxies. The processes involved in the fu-eling of active galactic nuclei (AGNs) encompass a wide rangeof scales, both spatial and temporal (Combes 2003, 2006; Jo-gee 2006). Current mm-interferometers are instrumental in pro-viding a sharp view of the distribution and kinematics of molec-ular gas, the dominant gas phase in galaxy nuclei, through ex-tensive CO line mapping. Combined with high-resolution nearinfrared images, interferometric CO maps are used to derive theangular momentum transfer budget in the circumnuclear disks ofAGNs (e.g., García-Burillo et al. 2005; Haan et al. 2009; Meidtet al. 2013). Maps of CO in nearby AGNs unveil a wide rangeof large-scale and embedded m = , (cid:63) Based on observations carried out with ALMA in Cycle 0. clei of Galaxies (NUGA) project, a CO interferometric surveyof a sample of 25 nearby low-luminosity AGNs (García-Burilloet al. 2003), indicate that molecular gas is frequently stalled inrings, which are the signposts of gravity torque barriers, and thatonly ∼ / m = ff erent spatial scales ofgalaxy disks (Hopkins et al. 2010a, 2011, 2012).Furthermore, the use of molecular tracers specific to thedense gas phase can probe the feedback of activity on the chem-istry and energy balance / redistribution of the interstellar mediumof galaxies. Observations suggest that the excitation and chem-istry of the main molecular species in AGNs are di ff erent withrespect to those found in purely star-forming galaxies (Tacconi Article number, page 1 of 24 a r X i v : . [ a s t r o - ph . GA ] J un & A proofs: manuscript no. ngc1068-paper-I-inpress-pdflatex et al. 1994; Kohno et al. 2001; Usero et al. 2004; Graciá-Carpioet al. 2008; Krips et al. 2008, 2011; García-Burillo et al. 2010;Imanishi et al. 2007, 2013; Aladro et al. 2013). Although the-oretical models have tried to explain these di ff erences, the un-derlying reasons for the apparent AGN specificity are still de-bated (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 consideredas a footprint of the mechanical feedback of activity, are beingdiscovered in a growing number of nearby active galaxies in-cluding ultra luminous infrared galaxies (ULIRGs), radio galax-ies, and Seyferts (Feruglio et al. 2010; Sturm et al. 2011; Alat-alo et al. 2011; Chung et al. 2011; Aalto et al. 2012; Dasyra &Combes 2012; Combes et al. 2013; Morganti et al. 2013; Ciconeet al. 2012, 2014). Radiative and mechanical feedback is ofteninvoked 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 bespatially resolved, are thus instrumental if we are to understandif and how gas accretion can self-regulate in galaxies.
NGC 1068 is a prototypical nearby ( D =
14 Mpc; Bland-Hawthorn et al. 1997) Seyfert 2 galaxy. It has a large-scale ovaland a nuclear bar with a pseudo-bulge which is overly mas-sive with respect to its central black hole (e.g., Kormendy &Ho 2013). NGC 1068 has been the subject of numerous cam-paigns using molecular line observations to study the fueling andthe feedback of activity. Schinnerer et al. (2000) used the Plateaude Bure Interferometer (PdBI) to map the emission of molecu-lar gas in the central r ∼ J = − J = − ∼ − . r ∼
200 pc circumnu-clear disk (CND) that surrounds the AGN. The CO( J = − r ≤
50 pc, the kinematicsof the molecular gas revealed by the 2.12 µ m H ff erent pic-ture. These data, sensitive to hot ( T K (cid:39) K) and moderatelydense ( n (H ) (cid:39) cm − ) molecular gas, show elliptical stream-ers that bridge the CND and the central engine. Müller-Sánchezet al. (2009) suggest that these structures correspond to gas feed-ing the AGN. Although on di ff erent spatial scales, inflowingand 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 ) (cid:39) − cm − ; Sternberg et al. 1994; Usero et al. 2004;Krips et al. 2008, 2011; Pérez-Beaupuits et al. 2007, 2009). Theuse of molecular tracers which are more representative of the to-tal H content is thus essential to derive the mass and energeticsassociated with the di ff erent inflowing / outflowing gas compo-nents at the CND.Interferometric images of NGC 1068, obtained in tracersspecific to the dense molecular gas, have also revealed the ex-istence of a strong chemical di ff erentiation in the disk of thisSeyfert. Essential diagnostic line ratios are di ff erent in the SBring compared to the CND. Tacconi et al. (1994) derived a highHCN / 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 abun-dance of HCN relative to CO is globally enhanced in the CND:HCN / CO ∼ − (Sternberg et al. 1994; Usero et al. 2004; Kripset al. 2011; Kamenetzky et al. 2011). The detection of molecu-lar ions like HOC + and H O + has been interpreted as the signa-ture of X-ray processing (Usero et al. 2004; Aalto et al. 2011).García-Burillo et al. (2010) analyzed the likely drivers of chem-ical di ff erentiation inside the CND with high-resolution obser-vations of CN and SiO. The abundances of SiO and CN areenhanced at the extreme velocities of gas associated with non-circular motions / shocks close to the AGN ( r <
70 pc). On theother hand, the correlation of CN / CO and SiO / CO ratios withhard X-ray irradiation suggests that the CND is a giant X-ray-dominated region (XDR). Although these results imply a strongradiative and mechanical feedback in the CND, the mechanismthat drives the excitation and chemistry of molecular gas is yetto be elucidated.
We use here the Atacama large millimeter array (ALMA) tomap the emission of a set of molecular gas tracers (the J = − J = − J = − + , and the J = − r ∼ ∼ . (cid:48)(cid:48) − . (cid:48)(cid:48) (20–35 pc). These line transitions span a range of critical densities( n crit [CO(3 − ∼ a few 10 cm − , n crit [CO(6 − ∼ a few10 cm − , n crit [HCO + (4 − ∼ cm − , n crit [CS(7 − ∼ afew 10 cm − , and n crit [HCN(4 − ∼ cm − ). The ac-tual average densities of molecular gas probed by these linesin NGC 1068 are seen to be ∼ − cm − due to subthermalexcitation of the higher J-lines and radiative trapping (Krips etal. 2011; Viti et al. 2014, hereafter, paper II). The line and contin-uum maps of the dense molecular gas and dust emission greatlyimprove the sensitivity and spatial resolution of any previous in-terferometric study of NGC 1068 in the (sub)millimeter range.We analyze the distribution, kinematics, and excitation of themolecular gas and study the processes associated with the fuel-ing and the feedback of activity in this Seyfert. In paper II we usethe line ratio maps derived in this work to model the excitationand chemistry of molecular gas in the di ff erent environments ofthe CND and SB ring. Star formation laws will be analyzed in afuture paper (García-Burillo et al. in prep.; paper III).We describe in Sect. 2 the ALMA observations and the an-cillary data used in this work. Section 3 presents the contin-uum maps obtained at 349 GHz and 689 GHz. Section 4 dis-cusses dust masses derived from continuum observations andillustrates the use of the CLUMPY torus models (Nenkova etal. 2008a, 2008b) to constrain the parameters of the torus. The Article number, page 2 of 24. García-Burillo et al.: Molecular line emission in NGC 1068 imaged with ALMA distribution of the molecular gas derived from the CO, HCN,HCO + , and CS line maps is discussed in Sect. 5. We describethe kinematics of the molecular gas and derive the main proper-ties of the outflow component in Sect. 6. A first description ofline ratio maps is presented in Sect. 7. The main conclusions ofthis work are summarized in Sect. 8.
2. Data
We observed the CO( J = −
2) emission and the CO( J = − CASA while the calibrated uv-tables were subsequentlyexported to GILDAS where the mapping and cleaning wereperformed as detailed below. Hereafter we adopt a distance toNGC 1068 of D ∼
14 Mpc (Bland-Hawthorn et al. 1997); thisimplies a spatial scale of ∼
70 pc / (cid:48)(cid:48) . In order to cover the CND and the SB ring of NGC 1068 weused an eleven-field mosaic in Band 7 with a field-of-view of17 (cid:48)(cid:48) per mosaic pointing. In total four tracks were observed be-tween June and August 2012 resulting in initial on-source times(i.e., before flagging) of ∼ −
39 minutes per track or a to-tal of 138 minutes. Between 18 and 27 antennas were availableduring the observations with projected baselines ranging from17 m to 400 m. Weather conditions were good with mediansystem temperatures of T sys = −
220 K and peak valuesnot exceeding 300 K. Four spectral windows with a bandwidthof 1.875 GHz each and a spectral resolution of 488 kHz wereplaced in Band 7, two in the lower side band (LSB) and twoin the upper sideband (USB); the centers of the two sidebandsare separated by 12 GHz. The four windows were centered onthe following sky frequencies: 341.955 GHz and 343.830 GHzin the LSB and 353.830 GHz and 354.705 GHz in the USB.This setup allowed us to simultaneously observe CO( J = − + ( J = −
3) (356.734 GHz atrest) in the LSB bands, and HCN( J = −
3) (354.505 GHz atrest) and CS( J = − +
017 was used to calibrate the amplitudes and phases intime. To set the absolute flux scale, Uranus and Callisto were ob-served using the Butler-JPL-Horizons 2012 catalogue to modeltheir visibilities as both were resolved at the given resolution ofthe observations. We estimate that the absolute flux accuracy isabout 10–15%.The angular resolution obtained using natural weighting was ∼ (cid:48)(cid:48) . × (cid:48)(cid:48) . ∼ ◦ in all the line and con-tinuum data cubes. The conversion factor between Jy beam − and K is 37 K Jy − beam. The line data cubes were binned to acommon frequency resolution of 14.65 MHz (equivalent to ∼ − in Band 7). The point source sensitivities in the linedata cubes were derived selecting areas free from emission in allchannels. They range from 2 mJy beam − (in the CS line) up to2.8 mJy beam − (in the CO line) in channels of 12.4-12.8 km s − width. Images of the continuum emission were obtained by av-eraging in each of the four sub-bands those channels free of http // casa.nrao.edu / http: // / IRAMFR / GILDAS line emission. The resulting maps were averaged to produce animage of the continuum emission of the galaxy at 349 GHz.The corresponding point source sensitivity for the continuum is0.14 mJy beam − .As our observations do not contain short-spacing correction,we expect that a significant amount of flux will be filtered outon scales ≥ (cid:48)(cid:48) for continuum emission in Band 7. A com-parison with the single-dish fluxes measured by Papadopoulos& Seaquist (1999) with the James Clerk Maxwell Telescope(JCMT) in apertures of ∼ (cid:48)(cid:48) indicates that we are filteringup to (cid:39) / ff ectmostly the low level emission that may extend on large-scalesin the interarm and spiral arm region. However, fluxes estimatedin small ( ∼ − (cid:48)(cid:48) ) apertures centered on the brightest spots ofemission of the CND and the SB ring are not expected to besignificantly underestimated. The percentage of missing flux inthe continuum maps in Band 7 represents the worst case sce-nario: the same estimate points to significantly lower missingfluxes in the molecular line ALMA maps. In order to estimatethe missing flux due to the lack of short-spacing correction in theCO(3–2) ALMA map we used as references the fluxes measuredin two studies of the galaxy done with single-dish telescopes:the single-pointed observations done at 22 (cid:48)(cid:48) -angular resolutionwith the Heinrich Hertz Telescope (HHT) by Mao et al. (2010),and the fully-sampled map of the central 1 (cid:48) region of NGC 1068done at 14 (cid:48)(cid:48) -angular resolution with the JCMT by Israel (2009a).Compared on these apertures (14–22 (cid:48)(cid:48) ), the ALMA map recov-ers 80–90% of the CO(3–2) flux at the CND and over the SBring. This percentage is lowered to 60–70% at the edges of thearea mapped by the ALMA mosaic and over the interarm region.As expected, the clumpy distribution of gas aided by the veloc-ity structure of the emission helps us recover the bulk of the fluxover the map on these spatial scales. Fluxes estimated on smaller( ∼ − (cid:48)(cid:48) ) apertures centered on the brightest emission spots ofthe CND and the SB ring are thus not expected to be significantlya ff ected.We set the phase tracking center of the central field of themosaic to α = h m . s , δ = − ◦ (cid:48) . (cid:48)(cid:48) ,which is the galaxy’s center according to SIMBAD (takenfrom the Two Micron All Sky Survey–2MASS survey; Strut-skie et al. 2006). By default, ( ∆ α , ∆ δ )-o ff sets are relative tothis position. Rest frequencies for all the lines were correctedfor the recession velocity initially assumed to be v o (LSR) = − = v o (HEL) = − . Systemic velocity is nev-ertheless re-determined in this work (Sect. 6.1) as v sys (LSR) = − = v sys (HEL) = − . Relative velocitiesthroughout the paper refer to v sys . Given the small field-of-view ( ∼ (cid:48)(cid:48) ) and the di ffi culty of Band 9observations, we used only one field to cover the circumnucleardisk in NGC 1068. In total four di ff erent tracks were observedduring July and August 2012 of which one had to be discardedbecause of poor quality. Weather conditions during the remain-ing three tracks were excellent with median system tempera-tures of T sys =
600 K for a range of 400–1200 K. This re-sults in an initial (i.e., prior to flagging) on-source time of 7–29 minutes per track or a total of ∼
52 minutes. Between 21 and27 antennas were used during each track with projected base-lines ranging between 17 m and 400 m. Three spectral windowswith a spectral bandwidth of 1.875 GHz and a spectral resolu-tion of 488 kHz were placed in the USB at sky frequencies of
Article number, page 3 of 24 & A proofs: manuscript no. ngc1068-paper-I-inpress-pdflatex
Fig. 1. a) ( Left panel ) The continuum emission map of NGC 1068 obtained with ALMA at 349 GHz. The map is shown in color scale (in Jy beam − units as indicated) with contour levels 3 σ , 5 σ , 10 σ , 15 σ , 20 σ , 30 σ to 120 σ in steps of 15 σ , where 1 σ = .
14 mJy beam − . ( ∆ α , ∆ δ )-o ff sets arerelative to the location of the phase tracking center adopted in this work: α = h m . s , δ = − ◦ (cid:48) . (cid:48)(cid:48) . We highlight the locationof the regions and components of the emission described in Sect. 3: the CND, the bar (identified by a representative isophote of the NIR K-bandimage of 2MASS), the bow-shock arc, and the SB ring. We also plot the edge of the eleven-field mosaic (gray dashed contour). The filled ellipseat the bottom right corner represents the beam size at 349 GHz (0 (cid:48)(cid:48) . × (cid:48)(cid:48) . PA = ◦ ). b) ( Upper right panel ) Same as a) but zooming in on theCND region. The position of the AGN ([ ∆ α , ∆ δ ] ≈ [–0.9 (cid:48)(cid:48) , –0.1 (cid:48)(cid:48) ] = [ α = h m . s , δ = − ◦ (cid:48) . (cid:48)(cid:48) ]) is highlighted by the starmarker. c) ( Lower right panel ) Continuum emission in the CND region at 689 GHz. Color scale is given in Jy beam − units. Contour levels are3 σ , 5 σ , 7 σ , and 10 σ to 25 σ in steps of 5 σ , where 1 σ = .
95 mJy beam − . The filled ellipse at the bottom right corner represents the beam size at689 GHz (0 (cid:48)(cid:48) . × (cid:48)(cid:48) . PA = ◦ ). ∼
16 GHz in Band 9. This setup allowed us to cover theCO( J = −
5) line entirely and having enough line-free channelsto determine the continuum emission. Bandpass calibration wasperformed on either Uranus or Ganymede while phase and am-plitude calibration in time was done on J0339-017. The absoluteflux scale was set by using models from the Butler-JPL-Horizons2012 catalogue for Uranus, Ganymede, and / or Callisto. For mosttracks, we had at least two solar bodies available so that we couldcrosscheck the flux calibration. However, we assume a conser-vative accuracy of ∼
30% for the absolute flux scale. We per-formed a self-calibration on the final maps using the strongestpeak emission, located in the eastern knot of the CND, whichpurposely coincides within 0.1 (cid:48)(cid:48) with the adopted phase trackingcenter. Self-calibration helped to improve the fidelity and moreimportantly the dynamic range of the final images.The synthesized beam obtained using natural weighting is ∼ (cid:48)(cid:48) . × (cid:48)(cid:48) . ∼ ◦ . The conversionfactor between Jy beam − and K is 27 K Jy − beam. The linedata cube was binned to a frequency resolution of 29.30 MHz(equivalent to ∼ . − ). The point source sensitivity inthe line data cube, derived by selecting areas free of emissionis 23 mJy beam − in channels of 12.8 km s − width. Imagesof the continuum emission were obtained by averaging in eachof the two nearby sub-bands at 686.899 GHz and 690.884 GHzthose channels free of line emission. The resulting maps werecombined to produce an image of the continuum emission of the galaxy at 689 GHz. The corresponding point source sensitivityfor the continuum is 1.95 mJy beam − .As our observations do not include zero-spacings, we expectto filter a non-negligible amount of flux on scales ≥ (cid:48)(cid:48) for con-tinuum emission at 689 GHz. The total integrated flux measuredinside the 9 (cid:48)(cid:48) field-of-view amounts to ∼
815 mJy. A compari-son with the fluxes measured with the JCMT by Papadopoulos& Seaquist (1999), and with Herschel by Hailey-Dunsheath etal. (2012), on apertures of ∼ (cid:48)(cid:48) indicates that we are filtering atmost (cid:39) −
45% of the total flux on these scales. However, fluxesestimated on small ( ∼ − (cid:48)(cid:48) ) apertures centered on the brightestspots of emission of the CND are not expected to be significantlya ff ected. As it happens in Band 7 observations, a similar esti-mate would also suggest that significantly lower missing fluxesare likely to be filtered out in the CO(6–5) maps, due to the highclumpiness and velocity structure of the emission.The phase tracking center of the central field and the refer-ence velocity are the same as the ones assumed in Band 7 (seeSect. 2.1.1 and Fig. 1). 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 theHLA with a fully revamped calibration pipeline which includestemperature-dependent dark frames, improved temperature mea-surements from bias levels, and other routines to reduce theimpact of image artifacts and improve the calibrated data. Thepixel size of the HLA images is 0 (cid:48)(cid:48) .
10 square. The redshift of
Article number, page 4 of 24. García-Burillo et al.: Molecular line emission in NGC 1068 imaged with ALMA
NGC 1068 is su ffi ciently small that Pa α falls within the F187Nfilter. 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 thetwo images at the nucleus, where the F190N emission peaks.Because the Pa α line is not centered at the filter central wave-length, we corrected the flux and converted the data numbers inthe image to flux densities using the IRAF task synphot .We also use the hard X-ray image obtained in the 6–8 keVband by Chandra. The X-ray map was obtained by applying anadaptive smoothing on the Chandra archive data available in thiswave-band for NGC 1068, following the same procedure de-scribed by Usero et al. (2004). This image is identical to thehard X-ray map published by Ogle et al. (2003) (see also Younget al. 2001).
3. Continuum emission
Figure 1a shows the continuum emission map obtained in Band 7(at 349 GHz). The emission detected inside the inner r ∼ (cid:48)(cid:48) (1.7 kpc) region imaged by ALMA, with a total integratedflux ∼
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 of4 (cid:48)(cid:48) × (cid:48)(cid:48) . (cid:48)(cid:48) . × (cid:48)(cid:48) . ∼
370 pc ×
200 pc),derived assuming the disk geometry fit of Sect. 6.1 ( PA = ± ◦ and i = ± ◦ ), reinforces the idea that the ring is of ellipticalshape. The ring, which shows a rich substructure, with two knotsof emission located east and west of the nucleus, is remarkablyo ff -centered relative to the location of the AGN. The latter isidentified as the strongest emission peak in the continuum mapat [ ∆ α , ∆ δ ] ≈ [–0.9 (cid:48)(cid:48) , –0.1 (cid:48)(cid:48) ] = [ α = h m . s , δ = − ◦ (cid:48) . (cid:48)(cid:48) ]. This maximum coincides within the errors withthe position of the AGN core ( S
1) determined from VLBI radio-continuum maps of the galaxy (Gallimore et al. 1996, 2004). Thestriking o ff -centered ring morphology of the CND is echoed bythe molecular line emission maps discussed in Sect. 5. This par-ticular geometry questions the simplest scenario that interpretsthe gas ring as the signature of the ILR region located at r ≤ (cid:48)(cid:48) (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 ∼ ± ◦ up to a corotation radius R COR ∼ ± (cid:48)(cid:48) ( ∼ . ∼ − (cid:48)(cid:48) from the AGN (depro-jected radii r ∼ −
650 pc). The source, hereafter referred toas the bow-shock arc , has a V-shaped morphology on its far side.This feature coincides in position with the northern edge of theAGN nebulosity identified in ionized gas emission, and with thenortheast radio lobe. Wilson & Ulvestad (1987) interpreted thelimb-brightened shape and the high polarization of the northeastradio lobe as the signature of a bow-shock in the galaxy disk.Leaving aside some isolated clumps, no significant emission isdetected anywhere else inside the bar region.
3. The SB ring:
Most of the emission in the disk ( ≥
72% ofthe flux integrated within the ALMA mosaic) is detected froma two-arm spiral structure that starts from the ends of the stellarbar and unfolds in the disk over ∼ ◦ in azimuth forming a Fig. 2.
The 689 GHz–to–349 GHz continuum flux ratio in the CND atthe spatial resolution of observations in Band 7 (0 (cid:48)(cid:48) . × (cid:48)(cid:48) . PA = ◦ ,shown by the filled ellipse in the lower right corner). Levels go from 5to 35 in steps of 5. pseudo-ring at r ∼ (cid:48)(cid:48) ( ∼ . Figure 1c shows the continuum emission map obtained in Band 9(at 689 GHz) in the inner ∼ (cid:48)(cid:48) field-of-view imaged by ALMA,which covers completely the CND region . Similar to the pic-ture drawn from observations in Band 7, the CND emission inBand 9 stems from a clumpy ring that is o ff -centered relativeto the AGN. The prominent east and west knots identified inBand 7 are spatially resolved into several clumps, due to thehigher spatial resolution at 689 GHz. Furthermore, while emis-sion is also detected at the AGN locus, like in Band 7 observa-tions, the strongest peak coincides with the position of the eastknot. There is also a secondary peak of emission in Band 9 whichis shifted to the northeast forming a spatially resolved ∼
35 pc–size elongated structure that connects the AGN with the CNDalong PA ∼ ◦ . This feature is reminiscent of the ∼ north-southelongated source present in the mid-infrared (MIR) maps of thecentral r ∼ − (cid:48)(cid:48) of the galaxy. This component has been at-tributed to emission of hot dust from narrow-line-region (NLR)clouds (Bock et al. 1998, 2000; Alloin et al. 2000; Tomono etal. 2001, 2006; Galliano et al. 2005; López-Gonzaga et al. 2014).Gallimore et al. (2001) mapped the subarcsecond radio jet andfound evidence of a jet-cloud interaction around 0 (cid:48)(cid:48) .
30 (20 pc)north and PA ∼ ◦ relative to the AGN locus, based on the de-tection of strong H O and OH maser emission and the bendingof the radio jet at this location. This is very close to the positionof a dust cloud identified in our maps, as shown in Fig. 3.
4. Dust masses
We combined observations of continuum emission in Bands 7and 9 of ALMA with observations done in two PACS bands The bow-shock arc region, not shown in Fig. 1c, lies beyond the halfpower of the primary beam of ALMA at 689 GHz.Article number, page 5 of 24 & A proofs: manuscript no. ngc1068-paper-I-inpress-pdflatex
Fig. 3. a) ( Upper panel ) The nuclear SED of the dust continuum emis-sion in NGC 1068 derived using NIR and MIR continuum and spec-troscopy data (blue squares) from Alonso-Herrero et al. (2011), andthe ALMA data points in Bands 7 and 9 (red stars). The SED was de-rived in apertures centered at the AGN that range from 0 (cid:48)(cid:48) . ∼
14 pc) forNIR data to 0 (cid:48)(cid:48) . ∼
35 pc) for MIR and Band-7 ALMA data. The bestCLUMPY model fit to the observations (curve) and the 1 σ uncertaintyrange of the fit (gray shaded region) are superposed onto the data points. b) ( Lower panel ) A close-up view of the dust continuum emission inBand 9. Levels and markers are as in Fig. 1. The (0 (cid:48)(cid:48) . × (cid:48)(cid:48) . at 70 µ m and 160 µ m by Hailey-Dunsheath et al. (2012), andwith high-resolution observations done at near- and mid-infrared(NIR and MIR) wavelengths by Alonso-Herrero et al. (2011) toconstruct the spectral energy distributions (SED) and derive themass of the dust in two di ff erent 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). r ≤ pc Figure 2 shows the map of the 689 GHz–to–349 GHz continuumflux ratio (hereafter, B / B
7) in the CND derived at the (lower)spatial resolution of observations in Band 7, after convolutionof the 689 GHz continuum image with a Gaussian beam withFWHM ∼ . (cid:48)(cid:48) . The B / B (cid:39) r ≤ (cid:48)(cid:48) .
5) up to a maximum of (cid:39) −
25 in aring-like region of radius r (cid:39) . (cid:48)(cid:48) . The highest ratios can be ex-plained if we assume a steep dependence of the dust emissivity κ ν , which scales as ∼ ν β , on frequency ( β ≥ . T dust ≥
100 K in the ring. On the otherhand, the low value of B / B ff erentmechanisms of non-thermal emission at the vicinity of the AGN:electron-scattered synchrotron emission, pure synchrotron emis-sion, and thermal free-free emission. The inclusion of the twoALMA band fluxes (integrated in 1 (cid:48)(cid:48) -apertures) into the AGNSED confirms that pure synchrotron emission is a poor repre-sentation of the SED in the submillimeter range, as also arguedby Krips et al. (2011): this scenario over-predicts by 50% thetotal flux measured in Band 7, i.e., well beyond the associateduncertainties ( < ≤ − We used the ALMA Band 9 and Band 7 continuum thermalfluxes at 435 µ m and 860 µ m, respectively, to investigate whetherwe can set further constraints on the properties of the putativetorus of NGC 1068 using CLUMPY torus models (Nenkova etal. 2008a, 2008b). To do so we took into account the NGC 1068SED and MIR spectroscopy presented in Alonso-Herrero etal. (2011). We added the ALMA Band 9 measurement at the fullresolution of 0 (cid:48)(cid:48) .
3, which is similar to that of the MIR data. TheALMA Band 7 continuum measurement has a slightly worse res-olution (0 (cid:48)(cid:48) .
5) and a higher contribution from non-thermal emis-sion as discussed in Sect. 4.1.1, and therefore we used the es-timated thermal component as an upper limit. For the ALMABand 9 flux uncertainties we included both the error in the photo-metric calibration ( ∼ ∼ ∼
10 larger value, which corresponds to a torus outer ra-dius of 20 + − pc, derived using the AGN bolometric luminosityof ∼ . + . − . × erg s − estimated from scaling the best fitmodel to the data. The torus size is not well constrained becausethe large value of the fitted index of the radial distribution of theclouds, which implies that most of the clouds are located close tothe inner radius of the torus. For comparison the modeled 12 µ minterferometric half radius of the resolved and unresolved com-ponents of NGC 1068 is 1.6 pc (Burtscher et al. 2013). Thesedi ff erences might be explained because the NIR and MIR emis-sion is probing warm dust that is on average closer to the AGN,whereas in the submillimeter range we are more sensitive to colddust distributed over larger distances from the AGN.The best CLUMPY model fit and the 1 σ uncertainty to thenuclear emission of NGC 1068 is presented in Fig. 3 a . It is clearfrom this figure that the measured continuum ALMA fluxes areabove the CLUMPY torus fit. This could be explained if cold Article number, page 6 of 24. García-Burillo et al.: Molecular line emission in NGC 1068 imaged with ALMA
Fig. 4. a) ( Left panel ) The CO(3–2) integrated intensity map obtained with ALMA using an eleven-field mosaic in the disk of NGC 1068. The mapis shown in color scale with contour levels 5 σ , 10 σ , 20 σ , 30 σ , 45 σ , 70 σ , 100 σ to 500 σ in steps of 50 σ , and 600 σ to 800 σ in steps of 100 σ , where1 σ = .
22 Jy km s − beam − . The filled ellipse at the bottom right corner represents the CO(3-2) beam size (0 (cid:48)(cid:48) . × (cid:48)(cid:48) . PA = ◦ ). b) ( Upper rightpanel ) Same as a) but zooming in on the circumnuclear disk (CND) region. c) ( Lower right panel ) Same as b) but for the CO(6–5) line, obtainedwith a single field mosaic. Contour levels are: 5 σ , 10 σ , 20 σ , 30 σ , 40 σ , 70 σ , 90 σ , 120 σ to 240 σ in steps of 40 σ , where 1 σ = − beam − .The filled ellipse at the bottom right corner represents the CO(6-5) beam size (0 (cid:48)(cid:48) . × (cid:48)(cid:48) . PA = ◦ ). dust not associated with the torus was included in the ALMABand 7 and 9 flux measurements. Figure 3 b shows a close-upof the nuclear region of NGC 1068 at 435 µ m at full resolu-tion. Clearly there is much cold dust in this region including theionization cone, as discussed in Sect. 3.2, but more importantlythere is no point source associated with the position of the AGN.This seems to confirm our suspicion that even at the 0 (cid:48)(cid:48) . × (cid:48)(cid:48) . . + . − . (cid:48)(cid:48) . The lower limit is consistent withnot having resolved the torus in Band 9 with the current angularresolution. Incidentally, the NIR data points are also above theCLUMPY torus model fit. The extra flux at these wavelengthsmight come from the dust at the base of the ionization cone, ascan be the case with other Seyfert galaxies (see, e.g., Hoenig etal. 2013).We can finally estimate the gas mass in the torus ofNGC 1068 based on the fit, according to Eq.A.1 of Appendix A.We obtained: M torus = . ± . × M (cid:12) , where we consid-ered that the main uncertainties come from the relatively un-constrained torus size and from the scatter around the adopted N H / A V scaling ratio taken from Bohlin et al. (1978) (see Ap-pendix A). This mass estimate is comfortingly similar to the es-timated molecular gas mass detected inside the central r =
20 pcaperture derived from the CO(3–2) emission, as discussed inSect. 5.4. r ≤ pc region: the CND and the bow-shock arc ALMA observations in Bands 7 and 9 were combined withPACS observations obtained at 70 µ m and 160 µ m in compa- rable apertures by Hailey-Dunsheath et al. (2012) to constrainthe overall SED of the dust emission in the central r =
400 pc ofthe galaxy. The estimate of M gas based on this fit and the impliedconversion factor for CO ( X CO ), discussed in Sect. 5.4, are usedin Sect. 6 to derive the mass load of the outflow identified in thisregion.The SED was fit using a modified black-body (MBB) modelto derive the dust temperature ( T dust ), dust mass ( M dust ) andemissivity index ( β ) in this region. In this approach, the mea-sured fluxes, S ν , can be expressed as S ν = M dust × κ ν × B ν ( T dust ) / D , where the emissivity of dust, κ ν , scales as ∼ κ
352 GHz × ( ν [GHz] / β , with κ
352 GHz = .
09 m kg − (a valuerounded up from κ
352 GHz = . kg − used by Klaas etal. 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 correctedfor the contamination by non-thermal emission in the central 1 (cid:48)(cid:48) aperture at the AGN, as estimated in Sect. 4.1.1. The best MBBfit is found for M dust = (8 ± × M (cid:12) , T dust = ± β = . ± .
2. The errorbars on the parameters of the fit reflectthe uncertainties due to the estimated range of missing flux inthe ALMA bands, derived in Sects. 3.1 and 3.2.The value of M dust can be used to predict the associated neu-tral gas mass budget in this region. Applying the linear dust / gasscaling ratio of Draine et al. (2007) (see also Sandstrom etal. 2013) to the gas-phase oxygen abundances measured in thecentral 2 kpc of NGC 1068 ( ∼ + log (O / H) ∼ .
8; Pilyuginet al. 2004, 2007), which yields a gas–to–dust mass ratio of ∼ + − , we estimate that the total neutral gas mass in the cen-tral 10 (cid:48)(cid:48) aperture is M gas = (5 ± × M (cid:12) . This number isa good proxy for the molecular gas mass in this region becausethe HI distribution in the disk of NGC 1068, studied by Brinks etal. (1997), shows a bright ring between 30 (cid:48)(cid:48) and 80 (cid:48)(cid:48) with a truecentral hole. Article number, page 7 of 24 & A proofs: manuscript no. ngc1068-paper-I-inpress-pdflatex
Fig. 5.
Zoom in of the inner disk of NGC 1068 to show the overlay oftwo CO(3–2) velocity channel maps obtained for v − v sys = −
50 (blackcontours: 5, 10, 20 to 100% in steps of 20% of the peak value = − (white contours: 5, 10, 20 to 100% in steps of 20% ofthe peak value = . v sys (HEL) = ± − , on the NIRK-band image (in color scale and arbitrary units) of NGC 1068 obtainedby the Two Micron All Sky Survey (2MASS). We highlight the locationof gas emission along the bar’s leading edges as well as the existence ofan anomalous component associated with the outflow. The orientationof the stellar bar’s major axis along PA = ◦ is marked by a dashedline.
5. Molecular line emission
In this section we analyze our newly acquired ALMA maps andcompare the morphology of the dense-gas tracers.
Figure 4 shows the CO(3–2) and CO(6–5) velocity-integrated in-tensity maps of NGC 1068 obtained with ALMA. The CO mapsreveal the distribution of the dense molecular gas in NGC 1068with unprecedented high-dynamic range capabilities: ≥
600 and ≥
200 for the CO(3–2) and CO(6–5) maps, respectively. As forthe dust emission, we identify three main regions in the disk ofNGC 1068:
1. The CND:
The brightest CO(3–2) line emission comesfrom the CND, which appears as a closed asymmetric ellipticalring of 6 (cid:48)(cid:48) × (cid:48)(cid:48) –size as shown in Fig. 4 b . The substructure of theCO(3–2) ring reveals two strong emission peaks located ∼ (cid:48)(cid:48) east and ∼ . (cid:48)(cid:48) west of the AGN. Similar to the dust contin-uum ring, the CO(3–2) ring is noticeably o ff -centered relative tothe location of the AGN: the two emission knots are bridged bylower-level emission that completes the ring north and south ofthe AGN and leaves a gas-emptied region to the southwest. How-ever, unlike for Band 7 continuum emission, CO(3–2) emissiondoes not peak at the AGN. The morphology of the ALMA mapof the CND is to a large extent similar to that of the SMA map ofKrips et al. (2011). Nevertheless, the order of magnitude higherdynamic range of the ALMA image of the CND helps reveal thatthe ring closes south of the AGN, i.e., similarly to the molecularring imaged with the Very Large Telescope (VLT) in the 2.12 µ mH line by Müller-Sánchez et al. (2009). The brightest CO(3-2)emission features are also detected and spatially resolved in theCO(6–5) map of Fig. 4 c .
2. The bar:
The CND appears connected to lower-level emis-sion which extends farther out in the disk into the stellar bar re-gion. On these scales the CO(3–2) emission is detected alongtwo lanes, o ff set by 2 − (cid:48)(cid:48) , which run mostly parallel to the bar’smajor axis ( PA = ± ◦ ; Scoville et al. 1988; Bland-Hawthornet al. 1997; Schinnerer et al. 2000; Emsellem et al. 2006). Theseare reminiscent of the gas leading edge morphology that is ex-pected to prevail between the corotation of the bar and the ILRregion. This is illustrated in Fig. 5, which shows the emissionof CO(3–2) at two velocity channels symmetrically located rel-ative to v sys : v − v sys = −
50 km s − and +
50 km s − , with v sys (HEL) = ± − as determined in Sect. 6.1.1. Gasemission along the bar’s leading edges appears at the velocityrange expected for gas falling to the nucleus: at redshifted veloc-ities on the southwestern (near) side and blueshifted velocities onthe northeastern (far) side. We nevertheless detect a gas compo-nent with anomalously redshifted velocities on the northeasternside of the disk at distances ∼ − (cid:48)(cid:48) from the AGN. This CO fea-ture, closely associated with the bow-shock arc identified in dustcontinuum emission, is interpreted in Sect. 6 as the signature ofa molecular outflow.
3. The SB ring:
Similar to the dust emission, most of theCO(3–2) flux in the disk mapped by ALMA ( (cid:39)
63% of the grandtotal) is detected in the SB ring. The two spiral arms, which un-fold over ∼ ◦ in azimuth, are spatially resolved. The (trans-verse) widths of the arms go typically from 100 pc to 350 pc.Figure 6 shows that the SB ring concentrates most of the ongo-ing massive star formation in the disk identified by strong Pa α emission. The emission of CO(3–2) over the SB ring is highlyclumpy and it appears to be organized as coming from molecu-lar cloud associations of ≥
50 pc-size which are also identifiedin the dust continuum map, as shown in Fig. 7.We also report the detection of CO(3–2) emission at di ff erentlocations throughout the interarm region. These interarm com-plexes, which remained undetected in dust emission due to thelower dynamic range of the Band 7 continuum map, are orga-nized into a network of filaments that extends out to the edge ofour mapped region. + , and CS maps Figure 8 shows the integrated intensity maps of NGC 1068 ob-tained 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 emis-sion in these likely denser gas tracers that are characterized by ∼ a factor 100 comparatively higher critical densities, stems fromthe CND. Notwithstanding, a few ( ∼
12) isolated clumps are de-tected in the HCO + (4–3) and HCN(4–3) lines at significant lev-els towards the SB ring (Figure 8 a shows the HCO + emission inthe SB ring). The di ff erent distributions of the CO(3–2), HCN(4–3), and HCO + (4–3) line emission in the disk is reflected in thedi ff erent line ratios measured in the CND and in the SB ring (seeSect. 7 and paper II). Other molecular tracers also show remark-ably di ff erent distributions in the SB ring and the CND (e.g.,Takano et al. 2014).The overall morphology of the CND in HCO + , and HCNis similar to that of the ALMA CO maps, i.e., an asymmetricclosed molecular ring, which is noticeably o ff -centered relativeto the AGN. However, the lower S / N ratio of the ALMA CS(7–6)map restricts the detected emission in this line mostly to the twoprominent knots of the CND. The superior capabilities of ALMAdramatically improve the picture drawn for the distribution of thedense molecular gas traced by the HCN and HCO + lines in the Article number, page 8 of 24. García-Burillo et al.: Molecular line emission in NGC 1068 imaged with ALMA
Fig. 6. a) ( Left panel ) Overlay of the CO(3–2) ALMA intensity contours (levels as in Fig. 4 a ) on the Pa α emission HST map (color scale asshown in counts s − pixel − ). b) ( Upper right panel ) Same as a) but zooming in on the CND region. c) ( Lower right panel ) Overlay of the CO(6–5)ALMA intensity contours (levels as in Fig. 4 c ) on the HST Pa α emission map (color scale as shown). The filled ellipses at the bottom right cornersrepresent the CO beam sizes. Fig. 7. a) ( Left panel ) Overlay of the CO(3–2) intensity contours (levels as in Fig. 4 a ) on the dust continuum emission at 349 GHz (color scalein Jy beam − units as indicated). b) ( Upper right panel ) Same as a) but zooming in on the CND region. c) ( Lower right panel ) Overlay of theCO(6–5) intensity contours (levels as in Fig. 4 c ) on the dust continuum emission at 689 GHz (color scale in Jy beam − units as indicated). Thefilled ellipses at the bottom right corners represent the CO beam sizes. CND compared to the previous SMA data of Krips et al. (2011):while most of the emission detected in the HCO + (4–3) SMAmap is restricted to the western knot, the ALMA map shows theemission to be widespread over the whole CND. The integratedflux of this line over the CND is a factor ∼ . ±
10 Jy km s − ) compared to the SMA map(60 ±
10 Jy km s − ). With the exception of the CS(7–6) line, emission is detected inall the molecular gas tracers probed by ALMA at the positionof 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 ex-pected, molecular emission is within the errors centered around v = v sys (HEL) = ± − (determined in this work: Article number, page 9 of 24 & A proofs: manuscript no. ngc1068-paper-I-inpress-pdflatex
Fig. 8. a) ( Left panel ) The HCO + (4–3) integrated intensity map obtained for NGC 1068. Contour levels are: 5 σ , 10 σ to 70 σ in steps of 10 σ , and85 σ . The filled ellipse shows the HCO + beam size, similar to that of CO(3–2). Panels b) ( upper right panel ) and c) ( middle right panel ) show azoomed-in view of the HCO + and HCN images, respectively, with 3 σ contours added to the list of displayed levels. d) ( Lower right panel ) Sameas b) and c) but for the CS(7–6) line. Contour levels are: 3 σ , and 5 σ to 40 σ in steps of 5 σ . The assumed value of 1 σ common for all lines is ∼ .
20 Jy km s − beam − . Fig. 9.
We show from left to right the CO(3–2), CO(6–5), HCN(4–3), and HCO + (4–3) emission line profiles towards the position of the AGN.The corresponding apertures are 0 . (cid:48)(cid:48) × . (cid:48)(cid:48) (40 pc) and 0 . (cid:48)(cid:48) × . (cid:48)(cid:48) (20 pc) for observations in Band 7 (CO(3–2), HCN(4–3) and HCO + (4–3))and Band 9 (CO(6–5)), respectively. The red lines identify the 1 σ level for each transition in order to illustrate the reliability of detections. Theundetected CS(7–6) line is included in the last panel (green histograms). Velocities refer to v sys (HEL) = ± − . see Sect. 6.1.1) and extends in contiguous channels across 200-300 km s − over significant levels ( > σ ) in all tracers. Gaus-sian fits to the line profiles extracted from the same apertures( ∼
40 pc) indicate that the higher density tracers (HCN(4–3) andHCO + (4–3)) show significantly wider lines (FWHM ∼ ±
10 km s − ) compared to CO(3–2) (FWHM ∼ ± − ).This result indicates that the excitation of molecular gas, esti-mated 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 alllikelihood correspond to gas lying closer to the central engine.This trend runs in parallel with the observed tendency to findhigher velocity-integrated ratios at the smallest radii inside theCND (see discussion in Sect. 7). As shown in Fig. 9 the half width of the CO(6–5) line derivedusing the velocity channels that show contiguous emission abovea 1 σ -level is ∼
125 km s − . Adopting an inclination angle i = − ◦ (Bland-Hawthorn et al. 1997; Sect. 6.1), the impliedspherical mass enclosed at r ∼ (cid:48)(cid:48) .
15 (10 pc) is M dyn ∼ − × M (cid:12) . This is a factor ∼ − − . × M (cid:12) ) estimated from the H O maser kinematicsin the inner r ∼ . r ∼ (cid:48)(cid:48) .
15 (10 pc) is contained in the central stellar cluster.
Article number, page 10 of 24. García-Burillo et al.: Molecular line emission in NGC 1068 imaged with ALMA conversion factor Using the CO(1–0) map of Schinnerer et al. (2000), we findthat the CO(1–0) flux ( S CO(1 − ) integrated in the central10 (cid:48)(cid:48) (700 pc)-aperture of the galaxy is ∼
100 Jy km s − . A com-parison with the CO(1–0) flux derived using a similar aper-ture on the BIMA map of Helfer et al. (1995), which, unlikethe Schinnerer et al. map, includes zero spacings, suggests thatthe PdBI map recovers (cid:39)
80% of the total flux in this re-gion. The implied conversion factor X CO = N (H ) / I CO(1 − thatis required to match the value of M gas derived from the dustSED fitting discussed in Sect. 4.2 is ∼ / (4 + − ) × X MWCO , where X MWCO = × cm − (K km s − ) − is the average CO conver-sion factor assumed to hold in the molecular clouds of the MilkyWay disk (Strong et al. 1988; Strong & Mattox 1996; Dame etal. 2001). A similar conclusion, pointing to an X CO ∼ a factorof 4–8 lower than X MWCO in the CND of NGC 1068, was reachedby Usero et al. (2004) who used LVG models to fit the CO lineratios observed in this region. This result is also in line with theobservational work of Israel (2009a, 2009b) and Sandstrom etal. (2013), who found that X CO can be up to a factor 4–10 lowerthan X MWCO in the central ∼ ∼ a factor of 10 lower than X MWCO intheir models of dense PDRs ( n (H ) ≥ cm − ) of starburst nu-clei. Lower X CO factors have also been commonly assumed forULIRGs; however, it is still debated whether this is a commonfeature of extreme starburst systems (e.g., see discussion in Pa-padopoulos et al. 2012).An estimate of the molecular gas mass detected inside thecentral r =
20 pc aperture ( M gas [AGN]) can be made basedon the detected CO(3–2) emission, provided that we assumea conversion factor for this line. Taking the average value de-rived above for X CO , which is representative for the 1–0 linein the CND, the conversion factor for the 3–2 line has to bescaled down by an additional factor 3, based on the 3–2 / ∼ mb units) derived in the neighborhood ofthe central engine of NGC 1068 (see Sect. 7). This implies that M gas [AGN] = . + . − . × M (cid:12) , i.e., the typical mass of a GMCand a value compatible within the errors with the gas mass de-rived from the dust emission in Sect. 4.1. This is equivalent to anaverage H column density N H = × mol cm − .Similarly, if we use the CO(6–5) integrated intensities in theALMA aperture for Band 9, i.e., r =
10 pc, and the factor of ∼ − X CO for this line estimated fromthe 6–5 / ∼ − mb units) measured at the AGN, wederive M gas [AGN] = . + . − . × M (cid:12) . Considering the smalleraperture used in the 6–5 line, this value is, not surprisingly,slightly below the estimate derived from 3–2. This is equivalentto an average H column density N H = . × mol cm − .
6. Gas kinematics: the molecular outflow
Figures 10 and 11 show the mean-velocity field of molecular gasin 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 arederived by integrating the emission above 5 σ -levels throughoutthe disk in all tracers. The clipping is lowered to 3 σ –levels whenwe zoom in on the CND region for HCN(4–3) and HCO + (4–3)in Fig. 11.The CO(3–2) isovelocities of Fig. 10 a , which sample thekinematics of molecular gas throughout the central ∼ (cid:48)(cid:48) (2.8 kpc)–region, show the expected pattern for a rotating diskcharacterized by an overall east-west orientation of its kinematic major axis ( PA = ◦ ± ◦ ; Bland-Hawthorn et al. 1997; Schin-nerer et al. 2000 ). At close sight, the orientation of the majoraxis is seen to change from the CND region ( PA ∼ ◦ ) to thebar region ( PA ∼ ◦ ) and farther out to the spiral arm region( PA ∼ ◦ ). The PA of the CND is close to the orientation ofthe kinematic major axis derived for the much smaller r ∼ . O maser emission published by Green-hill et al. (1996) ( PA = ◦ ). The PA values for the CO disk,derived from a qualitative fit of isovelocities, are confirmed bythe kinematic modeling of Sect. 6.1.1. If we assume that gas or-bits are roughly coplanar, the observed trend can be interpretedas the footprint of non-circular motions on the gas flow. As dis-cussed in Sect. 6.1.1, these distortions of the gas kinematics arecaused by di ff erent mechanisms which are at work on di ff erentspatial scales: the spiral arm structure in the outer disk, the stellarbar on intermediate scales and the nuclear outflow in the CND.All molecular tracers show a similar tilt of the CND kinematicmajor axis relative to the large-scale disk as shown in Fig. 11.The observed line widths (FWHM), displayed in Fig. 10 d ,go from 15–20 km s − in the well detected ( > σ ) molecularcomplexes of the interarm region, to an average value of 35–40 km s − in the SB ring. The corresponding velocity dispersionestimates would be σ v ∼ − − in the interarm regionand σ v ∼ −
17 km s − in the SB ring. The measured widthsat the CND, where FWHM values range from 50 to 200 km s − ,likely overestimate σ v for molecular gas in this region as theyreflect the superposition of di ff erent velocity components in-side the beam associated with rotation but also with an under-lying strong outflow pattern (see Sect. 6.1). Figures 10 e and 10 f show that the region in the CND where FWHM values are above ∼
100 km s − has a ring-like morphology. The widespread dis-tribution of high line width values in this region suggests thatthe outflow, which is found in Sect. 6.1 to be a mode superposedonto rotation, spatially extends over most of the CND.We analyze below the evidence supporting the existence of amassive molecular outflow in NGC 1068 based on the modelingof 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 alsofound in the molecular gas traced at higher densities by CO(6–5), HCN(4–3), and HCO + (4–3). The possible powering sourcesof the molecular outflow are discussed in Sect. 6.3. We explore inSect. 6.4 if other alternative scenarios can explain the observedkinematics. Section 6.5 compares the properties derived for themolecular outflow detected by ALMA with those derived fromprevious studies of NGC 1068. The general description of the two-dimensional line-of-sight ve-locity field of a galaxy disk can be expressed as v los ( x , y ) = v sys + v θ ( x , y ) cos ψ sin i + v R ( x , y ) sin ψ sin i , (1)where ( v R , v θ ) is the velocity in polar coordinates, ψ is the phaseangle measured in the galaxy plane from the receding side of theline of nodes, and i is the inclination angle restricted to the range0 < i < π/
2. With this convention v θ > v R > v R <
0) indicates outflow (inflow) for counterclockwise rota-tion and inflow (outflow) for clockwise rotation, respectively.We note that in the case of NGC 1068, where the receding side The position angle of the kinematic major axis is measured east fromnorth for the receding side. Article number, page 11 of 24 & A proofs: manuscript no. ngc1068-paper-I-inpress-pdflatex
Fig. 10.
Upper panels : a) Overlay of CO(3–2) isovelocity contours that span the range (–180 km s − , 180km s − ) in steps of 30 km s − on afalse-color velocity map (linear color scale as shown). Velocities refer to v sys (HEL) = − . b) Same as a) but zooming in on the CNDregion with a 20 km s − -velocity spacing. c) Same as b) but derived from the CO(6–5) data. Lower panels : d) Overlay of the CO(3–2) line widths(FWHM) shown in contours (10, 30, 50 to 200 km s − in steps of 30 km s − ) on a false-color width map (logarithmic color scale as shown). e) Same as d) but zooming in on the CND region. f) Same as e) but derived from the CO(6–5) data. of the major axis is to the west, rotation should be counterclock-wise 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 dividedinto a number of elliptical ring profiles defined by ( PA , i ) for agiven radius r , and v los ( r , ψ ) is decomposed as a Fourier serieswith harmonic coe ffi cients c j ( r ) and s j ( r ), where v los ( r , ψ ) = c + n (cid:88) j = [ c j ( r ) cos( j ψ ) + s j ( r ) sin( j ψ )] . (2)In the most general case, c reflects the contribution from cir-cular rotation while all remaining terms represent contributions from noncircular motions (Schoenmakers et al. 1997; Schoen-makers 1999). It is known that expanding the series of Eq. 2out to n = 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 c sin i and all remain-ing terms can be neglected. In the case of a pure (axisymmetric)radial flow ( v R ), c = v sys and s = v R sin i , with the rest of coef-ficients being 0.We derived the Fourier terms that best describe the v los extracted from the CO(3–2) mean-velocity field presented inSect. 5. We adopted an iterative three-step process using the soft-ware package kinemetry developed by Krajnovic et al. (2006). Article number, page 12 of 24. García-Burillo et al.: Molecular line emission in NGC 1068 imaged with ALMA
Fig. 11.
Mean-velocity fields derived from the HCO + (4–3) ( a) left panel ) data set. Panels b) ( upper right panel ) and c) ( lower right panel ) show aclose-up of the CND isovelocities derived from HCO + (4–3) and HCN(4–3), respectively. Contour levels, color scales and velocity reference in allpanels as in Fig. 10. See Sect. 6 for details on the di ff erent clipping thresholds used. Fig. 12. a) ( Left panel ) The magnitude of the c term of the Fourier decomposition of the velocity field of NGC 1068 described in Sect. 6.1 as afunction of radius (black curve). The c term represents the best fit of the (projected) axisymmetric circular component of the velocity field ( v circ ).We also plot the radial variation of the overall magnitude of the (projected) non-circular motions ( v nonc ) derived from the Fourier decompositiontill third order (red curve). b) ( Middle panel ). The variation of the v nonc /v circ ratio as function of radius. c) ( Right panel ) The variation of the s / c ratio as function of radius. The s term represents the best fit of the (projected) axisymmetric radial component of the velocity field. We used in the fit 28 ellipses with semi major axes covering thedisk from r = (cid:48)(cid:48) . r = (cid:48)(cid:48) , with a spacing ∆ r = . (cid:48)(cid:48) inthe inner 4 (cid:48)(cid:48) and ∆ r = (cid:48)(cid:48) farther out. At each step, kinemetry finds the best fitting ellipses by minimizing the contribution ofnoncircular motions ( v nonc ), evaluated as v nonc ( r ) = (cid:113) s ( r ) + s ( r ) + c ( r ) + s ( r ) + c ( r ) . (3)In a first step we left the position angle PA ( r ) and the inclina-tion i ( r ) as free parameters in the fit and assumed that the dynam-ical center coincides with the AGN. We then derived the averagevalues of PA ( r ) and i ( r ): (cid:104) PA (cid:105) = ± ◦ and (cid:104) i (cid:105) = ± ◦ (excluding outliers). In a second step, we fixed (cid:104) i (cid:105) = ± ◦ andre-determined the PA ( r ) profile. This profile shows a changing PA value from the CND region ( PA ∼ ◦ ± ◦ at r < (cid:48)(cid:48) ) tothe bar ( PA ∼ ◦ ± ◦ at 5 (cid:48)(cid:48) < r < (cid:48)(cid:48) ) and the spiral armregion ( PA ∼ ◦ ± ◦ at r > (cid:48)(cid:48) ). From this second itera-tion, we derived (cid:104) PA (cid:105) = ± ◦ . These values of (cid:104) PA (cid:105) and (cid:104) i (cid:105) are compatible within the errors with previous determinations ofthe disk orientation available in the literature ( PA = ± ◦ and i = ± ◦ , e.g., see compilation by Bland-Hawthorn et al. 1997).Finally, we derived the Fourier decomposition of the NGC 1068velocity field fixing PA = ± ◦ and i = ± ◦ at all radii.Implicit in this assumption is that the gas kinematics at all radiican be satisfactorily modeled by coplanar orbits (alternative sce-narios are described in Sect. 6.4). As an output parameter of thisfinal iteration we obtained the best fit for the systemic velocity, v sys (HEL) = ± − . This coincides within the errorswith the value of v sys inferred from the kinematics of the water Article number, page 13 of 24 & A proofs: manuscript no. ngc1068-paper-I-inpress-pdflatex
Fig. 13.
Comparison of s and s terms normalized by the circular ro-tation term c derived from the adopted best fit of the CO(3-2) veloc-ity field. The continuous line represents the least-squares fit to the datapoints, and the dashed line represents the expected warp line locationpredicting a relation between the s and s terms for an error in the po-sition angle adopted throughout the disk ( PA = ◦ ) (see discussionin Wong et al. 2004). The sign of s is taken as a signature of inflow oroutflow as shown, for the assumed geometry of the galaxy. Black sym-bols correspond to points at radii > (cid:48)(cid:48) , while red symbols correspondto radii ≤ (cid:48)(cid:48) . Fig. 14.
Overlay of the integrated intensity map of CO(3–2) (incontours: 15 σ , 30 σ , 60 σ , 100 σ , 200 σ , and 300 σ ; with 1 σ = .
22 Jy km s − beam − ) on the residual mean-velocity field ( (cid:104) V res (cid:105) , infalse-color scale spanning the range: –90 km s − to 90 km s − ) ob-tained after subtraction of the best fit rotation component, as describedin Sect. 6.1.1. Velocities refer to v sys (HEL) = − . maser emission detected in the central r ∼ . ab compare the magnitude of the (projected) cir-cular component of the CO(3–2) velocity field ( v circ = c ) ob-tained in the fit with that of the non-circular motions ( v nonc ) de-fined in Eq. 3. While the value of v nonc /v circ stays within the range ∼ . − . (cid:48)(cid:48) < r < (cid:48)(cid:48) ), reflecting the expected order of magnitude of the contribu- tion from the bar and the spiral structure to v nonc in this region,the same ratio reaches extremely high values ( ∼ . − .
30) in theinner r ≤ (cid:48)(cid:48) of the galaxy. As further illustrated in Fig. 12 c , themain contribution to v nonc comes from the s term. In particular,the sign ( >
0) and normalized strength of s ( s / c > .
3) indi-cates that strong outflow motions prevail at r ≤ (cid:48)(cid:48) . The s termchanges sign at r ∼ (cid:48)(cid:48) and stays moderately strong and negative( − . < s / c < − .
1) farther out in the disk; this behavior fairlyreflects the expected influence of the bar and the spiral structureon the gas flow provided that we are inside corotation of therelevant perturbation (bar or spiral) at these radii (Schinnerer etal. 2000; Wong et al. 2004; Emsellem et al. 2006).Figure 13, where we compare the s and s terms of theFourier decomposition of v los , also supports that bar-and-spiralinduced streaming motions (inside corotation) are the simplestexplanation for the s profile in the outer disk: the dominanceof the s term over the s term indicates that we remain insidecorotation in this region. This supports the conclusion that themolecular gas pseudo-ring at r ∼ (cid:48)(cid:48) (1.3 kpc) is not part ofthe inner bar pattern, but that it constitutes an independent wavefeature characterized by a lower pattern speed (e.g., see Rand& Wallin 2004). By contrast, this mostly excludes the stream-ing motion scenario to explain the strong outward radial mo-tions identified in the inner disk, which we rather interpret as amolecular outflow: values of s (cid:29) outside corotation for the assumed geometry of NGC 1068.Figure 14 shows the residual mean-velocity field ( (cid:104) V res (cid:105) ) ob-tained after subtraction of the rotation component derived in theanalysis above. The residuals clearly show a pair of approaching-receding regions in the outer disk. This morphology follows theexpected 2D-pattern produced by the combined action of the barand the spiral on the gas flow provided that we are inside corota-tion of the perturbations on these scales (Canzian 1993; Sempereet al. 1995, Colombo et al. 2014). This is in agreement with thepredominance of an inward radial flow that characterizes the fitto non-circular motions in the outer disk.Closer to the nucleus the comparison between the resid-ual velocity field, shaped by outward radial motions, and thegas / dust distribution in the inner disk suggests that a significantfraction of the gas in this region is associated with the outflow,as shown in Figs. 15 ab . Most remarkably, Figs. 15 cd show a no-ticeable spatial correlation between the AGN ionized nebulosity,traced by Pa α emission, the radio jet plasma, traced by the radiocontinuum emission at 22 GHz, and the molecular outflow sig-nature identified in the CO velocity field. This close association,which applies to a large range of radii and to a wide angle inthe disk, suggests that a sizable fraction of the dense moleculargas traced by CO(3–2) is being entrained due to AGN feedbackin the CND, and farther out in the disk, in the bow-shock arcregion (located at r ∼
400 pc).
To derive the mass load of the outflow we need to estimate thecharacteristic mass ( M mol ) as well as the projected values of theradial size ( R out ) and velocity ( V out ) of the outflow, and also as-sume a certain geometry (i.e., a certain angle α between the out-flow 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 (d M / d t ) can be estimatedfrom the expressiond M d t = × V out × M mol / R out × tan( α ) . (4) Article number, page 14 of 24. García-Burillo et al.: Molecular line emission in NGC 1068 imaged with ALMA
As discussed by Maiolino et al. (2012), if instead of a ho-mogeneous multi-conical outflow, we assume a single shell-like geometry of thinness dR out (where generally dR out (cid:28) R out ), theabove estimate should be replaced byd M d t = V out × M mol / dR out × tan( α ) . (5)In the following we use Eq. 4, as it provides a more conser-vative lower limit to the outflow rate estimated globally over theCND.The mass ( M mol ) has been calculated from the CO(3–2) datacube, after subtraction of the projected rotation curve (derivedin Sect. 6.1), by integrating the emission of the line outside avelocity range (cid:104) V res (cid:105) = [ − , +
50] km s − , which encompassesmost of the expected virial motions around rotation. We deter-mine that (cid:39)
50% of the total CO(3–2) flux in the CND stemsfrom the outflow component. A similar percentage is derived forthe bow-shock arc region, where we also find the signature ofthe outflow at larger radii. This translates into a total molecularmass M mol ∼ . + . − . × M (cid:12) for the CND, including the massof helium, if we assume a CO conversion factor ∼ / (4 + − ) of theMW value (see discussion in Sect. 5.4). We note that this esti-mate is rather conservative as the X CO conversion factor for thefraction of molecular gas that participates in the outflow couldbe higher if this component consists of an ensemble of opticallythick dense clumps embedded in a di ff use medium.The average projected radial extent of the outflow in theCND is R out ∼ . (cid:48)(cid:48) (100 pc) and the projected radial velocities V out are close to ∼
100 km s − according to the residual velocityfield shown in Figs. 14, 15 and 16.The implied outflow rate given by Eq. 4 is d M / d t ∼ + − × tan( α ) M (cid:12) yr − , where α reflects the unknown angle betweenthe outflow and the line-of-sight. If we assume that the outflowis coplanar with the main disk, tan( α ) = / tan( i ), where i = ◦ ,then d M / d t ∼ + − M (cid:12) yr − . This is significantly above themass load rate estimated on similar spatial scales for the ion-ized gas outflow: ∼ (cid:12) yr − (Müller-Sánchez et al. 2011).The molecular mass load rate implies a very short gas depletiontimescale of ≤ ff erent approach, Krips et al. (2011) also concluded that asignificant fraction ( ≥ µ mH line are suggestive of an outflow with typical projected veloc-ities V out ∼
100 km s − , i.e., in quantitative agreement with thevalues derived in this work. While the NIR line traces a minorfraction of the total H content and a distinctly di ff erent molec-ular gas phase to that probed by the ALMA data, the similarpicture drawn from these di ff erent data sets supports the outflowscenario in the CND.The same estimate applied to the outflow detected in thiswork in the bow-shock arc region, with V out ∼
75 km s − , M mol ∼ + − . × M (cid:12) and R out ∼ (cid:48)(cid:48) ( ∼
400 pc), gives d M / d t ∼ + − . M (cid:12) yr − for the uniformly filled model. The bow-shockarc component was not detected in the SMA map of Krips etal. (2011), nor is there a signature of the outflow on these scalesin the VLT map of Davies et al. (2008). n (H ) ≥ − cm − ) In the following we will show that the analysis of the kinemat-ics derived from CO(6–5), HCN(4–3), and HCO + (4–3) confirmsthat the molecular outflow detected in the CO(3–2) line has ahigher density ( n (H ) ≥ − cm − , according to paper II) coun-terpart in NGC 1068.As discussed in Sect. 5, the mean-velocity field of the CNDshows in all tracers a similar pattern: the kinematic major axisof the CND is tilted by ≥ ◦ relative to that of the large-scaledisk traced by CO(3–2), as shown in Fig. 11. This distinct kine-matic feature has been modeled as an outflow in CO(3–2) (seeSect. 6.1), and as such can be similarly interpreted as a signa-ture of outflow in the dense gas tracers. Figure 16 shows theposition-velocity (p-v) plots taken along the minor axis, deter-mined in Sect. 6.1 ( PA = ◦ ), for various dense moleculartracers. 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 / outflowmotions. An inspection of this figure shows that the outflow sig-natures identified in CO(3–2) out to radii r ∼
400 pc are echoedin CO(6–5), HCN(4–3), and HCO + (4–3) on the scales of theCND ( r ∼ −
200 pc). A sizable fraction of the emission ofthese molecular tracers lies outside the expected range of virialmotions attributable to rotation and dispersion: on average, emis-sion is ≥
50 km / s-redshifted on the northern side of the CND,while it is ≥
100 km / s-blueshifted on the southern side. Thisreflects the sign and the right order of magnitude of the mean-velocity field deviations seen in the CO(3–2) map of Fig. 15 atthe CND.The limited velocity coverage of the CO(3–2) data at the bluevelocity end ( (cid:104) V res (cid:105) < −
170 km s − ) implies that we underesti-mate the most extreme outflow velocities on the southern sideof the CND in this line. We note that outflow velocities up to (cid:104) V res (cid:105) ∼ +
280 km s − are detected in CO(3–2) on the northernside. The reality of these highly redshifted velocities in the out-flow is confirmed by the HCN(4–3) data, which show emissionfrom (cid:104) V res (cid:105) = –260 to +
260 km s − in the CND, as shown inFig. 16. This implies that the outflow mass load rate estimated inSect. 6.1.2 from CO(3–2) should be taken as a lower limit. Evidence of a young stellar population in the CND has been re-cently found by Storchi-Bergmann et al. (2012), who located theyoung star formation (SF) episode inside the expanding molec-ular ring at r ∼
100 pc. While the presence of young starsis well established, the total energy possibly injected by SFin the CND is much lower compared to the AGN. Davies etal. (2007) (see also discussion in Hailey-Dunsheath et al. 2012)estimated that only a fraction of the total FIR continuum lumi-nosity ( L FIR ) from the CND can be attributed to ongoing / recentSF. 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 . − × L (cid:12) ; this implies an integrated star formationrate of S FR ∼ . − . (cid:12) yr − . This estimate is similar tothe nuclear S FR measured from the 11.3 µ m PAH luminosity inthe inner r ∼
12 pc of NGC 1068 by Esquej et al. (2014), whoestimated that
S FR nuclear ∼ . (cid:12) yr − . The SFR estimated byEsquej et al. (2014) for the circumnuclear region out to a radius r ∼ (cid:48)(cid:48) (140 pc) is about 1 M (cid:12) yr − . The total SFR in the CND istherefore about an order of magnitude lower than the estimated Article number, page 15 of 24 & A proofs: manuscript no. ngc1068-paper-I-inpress-pdflatex
Fig. 15.
Overlay of the residual mean-velocity field( (cid:104) V res (cid:105) ) of Fig. 14, in color scale, with the contours rep-resenting: the integrated intensity of CO(3–2) ( a) up-per left panel ; contours as in Fig. 14), the 349 GHzcontinuum emission ( b) lower left panel ; contours asin Fig. 1), the HST Pa α emission ( c) upper rightpanel : with contours 0.2%, 0.5%, 1% to 5% in stepsof 1%, 10% to 90% in steps of 20% of the peakvalue = − pixel − ), and the 22 GHz VLAmap of Gallimore et al. (1996) ( d) lower right panel :with contours 1%, 2%, 5%,10%, 15%, 20%, 30%, 50%,70%, and 90% of the peak value =
36 mJy beam − ). mass load rate of the molecular outflow. Taken at face value thisdiscrepancy suggests that SF is not able to drive the molecularoutflow in the CND of NGC 1068 (see Murray et al. 2005 andVeilleux et al. 2005 for a general discussion).The kinetic luminosity of the CND outflow can be derivedfrom the expression L kin = / × d M d t × (cid:32) V out cos( α ) (cid:33) . (6)Assuming that the gas in the outflow is coplanar, L kin ∼ + . − × erg s − . AGN feedback models require that a sig-nificant fraction of the radiated luminosity ( ∼ L bol ; di Matteoet al. 2005) should be coupled to the ISM to produce an outflow.This fraction is lowered to ∼ .
5% in the two-phase feedbackmodel of Hopkins & Elvis (2010b). The bolometric luminosityof 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 ≥ − erg s − (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 theAGN can power the outflow in the CND of NGC 1068.The momentum flux of the CND outflow can be computedfrom the expressiond P out d t = d M d t × V out cos( α ) . (7)In the case of coplanar gas, Eq. 7 yields d P out / d t ∼ + − . × g cm s − . Compared to the momentum provided by the AGNphotons, derived as L bol / c ∼ (0 . − × g cm s − , d P out / d t isa moderate factor 1–27 larger. The required boost is only ∼ L bol = . × erg s − from Sect. 4.1.2. The molec-ular outflow could thus be in the momentum-conserving regime (where cooling is fast ), and radiation pressure could be the driv-ing mechanism. This is also within the range of values of the mo-mentum boost factors predicted by AGN feedback models underthe assumption that molecular outflows are energy-conserving(adiabatic Sedov phase): (d P out / d t ) / ( L bol / c ) ∼ −
50 (Faucher-Giguère & Quataert 2012).Alternatively, the radio jet of NGC 1068 could also inject therequired outflow power. The jet power ( W jet ) can be estimatedfrom the monochromatic luminosity at 1.4 GHz, according toBîrzan et al. (2008). In the case of NGC 1068, the spatially inte-grated flux of the jet at 1.4 GHz, including the main components(known in the literature as NE , C , and S ), is ∼
840 mJy fromthe radio continuum map of Gallimore et al. (1996). This im-plies W jet = . × erg s − . Since W jet ∼ (30 − × L kin , weconclude that the jet can drive the molecular outflow in the CNDeven assuming a low coupling e ffi ciency. The velocity field of the outer disk from r ∼
400 pc out to ∼ . ff -centered ring, which is apparently rotating but witha di ff erent kinematic axis, as shown in Fig. 10. As an alterna-tive to the outflow scenario described in Sects. 6.1 and 6.2, theabrupt shift of ≥ ◦ degrees in the kinematic PA of the twodisks can be interpreted as due to the CND being a non-coplanardisk which is far from being dynamically relaxed. Two types ofnon-coplanar instabilities can be invoked to account for the de- Article number, page 16 of 24. García-Burillo et al.: Molecular line emission in NGC 1068 imaged with ALMA
Fig. 16.
Position-velocity (p-v) plots takenalong the kinematic minor axis ( PA = ◦ )of NGC 1068 help identify outflow signaturesin several molecular tracers: a) ( Upper leftpanel ) CO(3–2) (color scale in Jy beam − andcontours: 3 σ , 5 σ , 10 σ , 20 σ , 30 σ , 45 σ , 70 σ ,and 95 σ ; 1 σ = . − ). b) ( Upperright panel ) CO(6–5) (color scale in Jy beam − and contours: 2.5 σ , 4 σ , 6 σ , 8 σ ,10 σ , 14 σ , and18 σ ; 1 σ =
23 mJy beam − ). c) ( Lower leftpanel ) HCN(4–3) (color scale in Jy beam − and contours: 2.5 σ , 4 σ , 6 σ , 8 σ , 12 σ , 16 σ ,20 σ , 26 σ , and 30 σ ; 1 σ = . − ). d) ( Lower right panel ) HCO + (4–3) (color scalein Jy beam − and contours: 2.5 σ , 4 σ , 6 σ , 8 σ ,10 σ , and 12 σ ; 1 σ = . − ). Thevelocity scale ( (cid:104) V res (cid:105) ) is identical to that ofFig. 14. The spatial scale ( ∆ y) along the minoraxis refers to the AGN locus; positive o ff setson the northern side. The black dashed linesat (cid:104) V res (cid:105) = ±
50 km s − identify the expectedrange of virial motions along the minor axis.The red solid lines in panels a) and c) indicatethe edges of the bands beyond which data havebeen flagged ( (cid:104) V res (cid:105) < −
170 km s − in CO(3–2)and (cid:104) V res (cid:105) >
260 km s − in HCN(4–3)). Sim-ilarly we also identify the 9 (cid:48)(cid:48) field-of-view inthe CO(6–5) p-v plot of panel b ). coupled kinematics observed in the CND: a nuclear warp or anon-coplanar lopsided instability. Schinnerer et al. (2000) explored a nuclear warp scenario in anattempt to model the kinematics of molecular gas in the nu-clear region of NGC 1068, based on CO(2–1) observations ob-tained with the PdBI. They concluded that gas motions could beequally fit with either a warp or a nuclear bar in the CND. Oneof the main predictions of their warp model was that the CO diskshould become edge-on at a radius of ∼
70 pc. However, we donot see the signature of an edge-on disk on these spatial scales inthe higher-resolution ALMA maps presented in this work. Fur-thermore the internal kinematics of the CND do not show thetypical
S-shaped distortion attributable to a nuclear warp insta-bility.The diagram shown in Fig. 13, originally introduced byWong et al. (2004), can also be used as a diagnostic tool to iden-tify nuclear warp signatures in the velocity field of galaxies in theparticular case where the hypothesized tilted orbits share a com-mon center. An inspection of Fig. 13 indicates that the continu-ous line, which represents the least-squares fit to the NGC 1068data points, shows no correspondence to the expected locationof the warp line , which would relate the s and s terms in thecase of a warp having the AGN as a common center of the tiltednon-coplanar orbits. This disagreement, more evident in the in-ner disk, leads us to conclude that the simplest nuclear warp sce- nario described above is not a satisfactory explanation for thevelocity residuals observed in this region. If we abandon the restriction of having the AGN as common or-biting center for both the CND and the outer disk, the scenarioof non-coplanarity remains nevertheless viable. In this scenario,a non-coplanar
CND would be orbiting around a secondary nu-cleus at ( ∆ α , ∆ δ ) ∼ ( − . (cid:48)(cid:48) , − (cid:48)(cid:48) ), i.e., 1 (cid:48)(cid:48) =
70 pc o ff set to thesouthwest relative to the AGN.Two types of mechanisms have been described as potentialtriggers of lopsided non-coplanar gas instabilities in galactic nu-clei:
1. External trigger:
The postulated secondary nucleus couldbe the footprint of a recent minor merger with a nucleated satel-lite. In this scenario dynamical friction would make the satel-lite quickly sink toward the nucleus of the host dragging the ac-creted gas disk to the central region (Taniguchi & Wada 1996;Taniguchi 1999, 2013). Minor mergers could thus explain a ran-dom orientation of circumnuclear gas disks in the NLR of Seyfertgalaxies, as the orbital plane of the resulting nuclear disk wouldbe determined by the orbital parameters of the accreted satellite(Taniguchi 2013). As the CND is noticeably o ff -centered, thiswould imply that NGC 1068 would be at an early stage in themerging process where the secondary nucleus has not yet sunktoward the nucleus. However, the absence of any signature of asecondary nucleus in the high-resolution HST / NICMOS images
Article number, page 17 of 24 & A proofs: manuscript no. ngc1068-paper-I-inpress-pdflatex ! ! !" ! " %& ’( " ) * + " ’( % , " (- . , ’ / % & ’ ( " ) * + " ’ ( % , " ( - . , ’ / % & ’ ( " ) * + " ’ ( % , " ( - . ! " " $ % & ’ ( ) ! " " $ % & ’ ( ) ! " %& ’( " ) * + " ’( % , " (- . %&’()*’+,-./01 ’ . ! " $ " % % % α % % % % " % Fig. 17.
A revised version of the kinematic model of the NLR, first proposed by Tecza et al. (2001) and Cecil et al. (2002), which now accounts forthe molecular outflow (denoted in figure as
CO outflow ) detected by ALMA in the CND (seen in projection in N and S, for northern and southernknots) and farther north in the molecular disk. The figure shows a crosscut of the NLR as viewed from inside the galaxy disk along the projecteddirection of the radio jet ( PA ∼ ◦ ; shown by the green line). We highlight the extent of the ionized gas outflow (light purple shade) and theassumed geometry defined by angles α and i . of the CND, which should be at the position of the apparent guiding center of the non-coplanar disk, makes the minor mergerscenario implausible in NGC 1068.
2. Internal trigger:
Alternatively, a lopsided instabilitytriggered on the gas could have an internal origin (Jog &Combes 2009). The CND could be the result of a recent episodeof gas infall produced by the feedback of a star formation burstin the disk. This feedback action would be able to eject gas per-pendicular to the plane, at heights of 100-200 pc, as frequentlyseen in numerical simulations (e.g., Emsellem et al. 2014). Thegas fountain would fall back while settling initially in an inclineddisk, the perpendicular orientation being preferred since di ff er-ential precession is then canceled. This would explain the o ff setbetween the AGN and the guiding center of the instability. Whilesettling back to the disk, the nuclear gas would trigger a slowand long-lasting lopsided perturbation, which could remain dur-ing more than 100 dynamical times at this radius, i.e., 400 Myrfor the CND (e.g., Jog & Combes 2009). In this scenario it isdi ffi cult to predict what evolutionary stage during the settling ofthe disk best represents the case of NGC 1068. However, at theend of the process, the m = ffi ciently interact with molecular gas in a coplanar CND. Overall, the outflow hypothesis is the simplest explana-tion for the distorted kinematics of molecular gas at the CND, butalso for the bow-shock arc region, the signature of the interac-tion at larger radii. In this scenario the remarkable o ff -centeringof the molecular ring would be triggered before the jet or the ion-ized gas wind hits the disk and would launch the outflow ratherthan being the result of it. Furthermore, the lopsided morphol-ogy in a coplanar CND, which is characterized by a strong near side / far side asymmetry in NGC 1068, could also be enhanceddue to opacity e ff ects in molecular lines, as discussed by Booneet al. (2011). NGC 1068 harbors a wide-angle (FWHM ∼ − ◦ ) bicon-ical outflow of ionized gas (Macchetto et al. 1994; Arribaset al. 1996; Crenshaw et al. 2000; Tecza et al. 2001; Cecilet al. 2002; Mueller-Sánchez et al. 2011). The orientation ofthe outflow is not perpendicular to the galaxy disk: i out f lo w ∼ − ◦ , while i disk ∼ ◦ . This particular geometry favors in-teraction with the molecular disk out to a galactocentric radius of r ∼
400 pc. An interaction between the jet and the ISM was al-ready identified close to the AGN ( r ∼ −
40 pc) by Gallimoreet al. (1996, 2001) (see also Bicknell et al. 1998).Figure 17 shows a simple kinematic model of the NLR, firstproposed by Tecza et al. (2001) (see also: Cecil et al. 2002; Daset al. 2006, 2007), here revised to account for the molecularoutflow detected by ALMA. The figure shows a scaled cross-cut of the NLR along the projected direction of the radio jet( PA ∼ ◦ ). The observer’s line-of-sight, which makes an angle i disk ∼ ◦ relative to the direction orthogonal to the moleculardisk, 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 rightof the figure. The CND, represented by two molecular knots lo-cated asymmetrically north and south of the AGN (N and S infigure), is embedded in the molecular disk.According to Tecza et al. (2001) (see also Cecil et al. 2002and Mueller-Sánchez et al. 2011), high-ionization lines, like[Si VI ] or [O III ], are produced where the radio lobes encounterdi ff use material located outside the galaxy plane. The interac-tion of the radio lobes with this medium generate two bow-shock Article number, page 18 of 24. García-Burillo et al.: Molecular line emission in NGC 1068 imaged with ALMA
Fig. 18. a) ( Left panel ) Overlay of the HST Pa α emission map (in contours as in Fig. 15 c ) on the CO(3–2) / CO(1–0) brightness temperature ratiomap (in color scale and T mb units) derived at the spatial resolution of the 1–0 observations of Schinnerer et al. (2000) (2 (cid:48)(cid:48) × (cid:48)(cid:48) . PA = ◦ ;ellipse in lower right corner). b) ( Upper right panel ) Same as a) but showing a zoom in on the CND region. c) ( Lower right panel ) Same as a) but showing the CO(6–5) / CO(3–2) brightness temperature ratio map (in color scale and T mb units) derived at the spatial resolution of the 3-2observations (0 (cid:48)(cid:48) . × (cid:48)(cid:48) . PA = ◦ ; ellipse in lower left corner). fronts that are responsible for the simultaneous detection of red-shifted and blueshifted emission in both ionization cones, asshown in Fig. 17. The low-ionization lines, like [Fe II ], are pro-duced when the bow-shock sweeps the denser molecular disk.This explains why low-ionization lines are exclusively detectedas redshifted (blueshifted) emission on the northern (southern)side of the disk. The amplitude of velocity shifts are larger for thehigh-ionization lines, reaching ∼ − in the northeast-ern cone, while they remain ≤
500 km s − for the low-ionizationlines. Furthermore, the sign of the velocity shifts for the strongestemission components of the high-ionization lines are noticeablyreversed with respect to low-ionization lines.The molecular outflow detected by ALMA stands out as aredshifted (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 islaunched when the ionization cone of the NLR sweeps the diskin the CND and farther out north in the bow-shock arc regionwhere Wilson & Ulvestad (1987) found that the radio-lobe neb-ulosity becomes limb-brightened and highly polarized, the sig-nature of a bow-shock in the disk. The amplitudes of the velocityshifts in the molecular outflow are significantly smaller than forthe high-ionization lines. There is also evidence that the molec-ular outflow becomes decelerated at r ∼
400 pc: the terminalvelocities of the outflow in the bow-shock arc region are a factor2 smaller than in the CND, as shown in Fig. 16.The detection of emission from dense gas tracers at theextreme velocities of the molecular outflow in NGC 1068 ( ≥
150 km s − ) raises the question of how molecular clouds cansurvive during the launching of the outflow, since shocks of ve-locities greater than ∼
50 km s − are known to be fast enoughto destroy molecules (Hollenbach & McKee 1989; Neufeld &Dalgarno 1989a, 1989b). However, the chemistry in both dis-sociative J-type and non-dissociative C-type shocks, which in-volve dust-mantle disruption and a high-temperature environ- ment, favors an e ffi cient fast reformation of molecules in thegas. Molecules can reform in the post-shocked gas on timescales ≤ hundreds of yrs. There is ample observational evidence thatthe emission of dense gas tracers ( ≥ − cm − ) can coexist withfast shocks ( ∼ −
100 km s − ) in galactic bipolar outflows ofYoung Stellar Objects (YSOs) (e.g., Arce et al. 2007 and refer-ences therein). These observations show that the abundance ofsome molecular species can undergo spectacular enhancementsin the outflow gas, due to the onset of shock chemistry. Be-sides classical shock tracers such as SiO, whose abundance canbe enhanced by factors of about 10 , molecular species whichare 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 strongemission from some of these tracers, such as HCN or CS, atthe velocities identified as abnormal in NGC 1068 suggests thatmolecular clouds survive in the outflow and that shock chemistryis likely at work in this component.
7. Molecular line ratios and environment
Figure 18 shows the overlay of the HST Pa α emission map onthe CO(3–2) / CO(1–0) (hereafter R /
10) and CO(6–5) / CO(3–2) (hereafter R /
32) line ratio maps. Line ratios were derivedin T mb units. To derive the 3–2 / (cid:48)(cid:48) × (cid:48)(cid:48) . / CO(3–2) brightness temperature ratiomap, shown in panel c, was derived at the spatial resolution of the3–2 ALMA observations (0 (cid:48)(cid:48) . × (cid:48)(cid:48) . σ clipping on the integrated intensities toassure image reliability.The R /
10 ratio changes significantly depending on the par-ticular environment of the disk. The average ratio is ∼ . ± . Article number, page 19 of 24 & A proofs: manuscript no. ngc1068-paper-I-inpress-pdflatex
Fig. 19. a) ( Upper panel ) Overlay of the HST Pa α emission map (incontours as in Fig. 15 c ) on the HCN(4–3) / CO(3–2) brightness temper-ature ratio map (in color scale and T mb units). b) ( Middle panel ) Sameas a) but showing the HCO + (4–3) / CO(3–2) brightness temperature ratiomap. c) ( Lower panel ) Same as a) but showing the HCN(4–3) / HCO + (4–3) brightness temperature ratio map. All the ratio maps have been de-rived at the spatial resolution of the CO(3-2) observations (0 (cid:48)(cid:48) . × (cid:48)(cid:48) . PA = ◦ ; ellipses in lower left corners). in the SB ring yet R /
10 shows a large dispersion of values: R / ∼ . − α emission, while the strongest Pa α emitting regions show ∼ a fac-tor 3–4 higher ratios ( R / ∼ − R /
10 ratio in the CND, ∼ . ± .
1, is higherthan in the SB ring, in agreement with the global values mea-sured with the PdBI and the SMA by Krips et al. (2011) andTsai et al. (2012). In summary, the R /
10 ratios measured inNGC 1068, most particularly those of the CND, are at the highend of the typical values observed in the centers of nearby nor-mal and starburst galaxies: ∼ . − . ∼ . − . ∼ . − . At the spatial resolution used in Figs. 18 ab , the R /
10 ratioshows a range of values in the CND: R /
10 goes from ∼ ∼ ∼ − R / / jet nebulosities. Overall, the R /
32 ratio in the CND is ∼ .
8, with a significant range ofvalues: R65 /
32 goes from ∼ . ∼ . c , the excitation of CO in the CND probedby the R /
32 ratio unveils the footprint of AGN feedback: the R /
32 ratio is enhanced by the degree of illumination of themolecular gas by the photons of the ionized gas outflow tracedby Pa α emission.Other molecular line ratios show dramatic changes of upto an order of magnitude inside the CND on the spatial scalesprobed by ALMA ( ∼
35 pc). Figure 19 shows the overlay ofthe HST Pa α emission map on a set of four molecular line ra-tio maps (in T mb units) obtained at the resolution of ALMAobservations 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 AGNfeedback: these line ratios go from ∼ . − .
03 at the outerradii of the CND to ∼ R HCN / HCO + ratio is globally very high in the CND withan average value of ∼ / HCO + ratios measuredin galaxy nuclei have also been interpreted as the signatureof significant mechanical heating in shock / mechanically domi-nated regions (MDRs) (Kazandjian et al. 2012). Nevertheless,the ALMA maps show that the lowest value of this ratio is foundprecisely at the AGN locus, where R HCN / HCO + ∼ .
3, an indica-tion that a simplistic interpretation of these ratios can be mislead-ing (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 vari-ations that can reach orders of magnitude for times ≥ years(Meijerink et al. 2013).As discussed in García-Burillo et al. (2010), the morphol-ogy of hard-X ray emission in the 6–8 keV band obtained byChandra (Young et al. 2001; Ogle et al. 2003) can be used tosearch for observational evidence of AGN feedback on the exci-tation / chemistry of molecular gas in the CND. Emission in the6–8 keV band observed in NGC 1068 is dominated by reflec-tion 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 ofmolecular gas by hard X-rays normalized by gas column den-sity; 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 quanti-ties have been measured at a common spatial resolution of 0.5 (cid:48)(cid:48) .With an admittedly large scatter, Fig. 20 provides quantitativeevidence of a correlation between the R HCN / CO and R HCO + / CO ra-tios and X hard / CO. This is similar to the correlation found for the
Article number, page 20 of 24. García-Burillo et al.: Molecular line emission in NGC 1068 imaged with ALMA log ( X hard CO ) l og ( R HCNC O ) ● ●●●●●●● ●●●●●● ● ● ●●●●●●●● ● ●●●●●●●●●● ● ● ●●●●●●●● ●●●●●●●● ● ●●●●● ● ●●●●●●●●●●●● ● ● ● ●●●●●●●●●●● ● ● ● ● ● ●●●●●●●●● ● ● ● ● ● ● ●●●●●● ● ●●●●● ● ●●● ● ●●●●●●● ● ●●●●●●●● ●●●● ●●●●●●● ● ●● ●●●●●●●●●●●●● ●● ● ● ●●●●●●●●●● ●●●● ● ● ●●●●●●●●●●●●●●●●● ● ●●●●●●●●●●●●●● ●●● ● ●●●●●●●●●●●●● ●●● ● ● ●●●●●●●●●●●●●●●● ● ● ● ●●●●●●●●●●●●●●● ● ●●●●●●● ● ●●●●●●● ●●●●●●● ● ●●●●●●●●●●●●●●● ● ●●● ●●●●●●●●●●●●●● ●●●●●●●●●●●● ●●●●●●●●●●●● ●●●●●●●●●● ●●●●●●●● ● ●●●●● ●● ● −3 −2 −1 0 − − y BS = + log ( X hard CO ) l og ( R HC O + C O ) ● ●●●●●● ● ●●●●●● ● ●● ●●●●●●● ● ● ●●●●●●●●● ● ●●●●●●●●●●●●●●● ● ● ●●●●●●●●● ● ● ● ● ●●●●●●●●●● ● ● ● ● ● ●●● ● ●●●●● ● ●●● ●●●● ● ●●●●●●●● ●●●●● ● ●● ●●●●●●●●●● ●● ● ● ●●●●●●●●●● ●● ●●●●●●●●●●●● ●●●●●●●●●●●●●● ● ●●●●●●●●●●● ● ● ●●●●●●●●●●●● ● ● ● ●●●●●●●●●●●●●●●●● ● ●●●●●●● ● ● ●●●●●●● ●●●●●●● ● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●● −3 −2 −1 0 − − y BS = −0.2 + Fig. 20. a) ( Left panel ) The HCN(4–3) / CO(3–2) velocity-integrated intensity ratio ( R HCN / CO , in logarithmic scale) versus the ratio of X-ray flux inthe 6–8 keV band (in counts s − ) to the CO(3-2) intensities (in K km s − ) ( X hard / CO, in logarithmic scale) measured at a common spatial resolutionof 0.5 (cid:48)(cid:48) over the CND (open circles). Isodensity contours illustrate the distribution of values in this parameter space for 10%, 25%, and 50% of thedata points. The straight line represents the bisector linear fit to the data in logarithmic scale [log( X hard / CO), log( R HCN / CO )]. The yellow-shaded areashows the range allowed by the ordinary least squares fits of y as a function of x and viceversa. The parameters for the bisector fit are indicated. b) ( Right panel ) Same as a) but particularized for the HCO + (4–3) / CO(3–2) ratio ( R HCO + / CO ) versus X hard / CO. The linear correlation coe ffi cients ofthe regressions are r ∼ . SiO / CO and CN / CO ratios as function of hard X-ray illumina-tion of molecular gas in the CND, discussed by García-Burilloet al. (2010).The existence of an AGN-driven outflow proves that thefeedback of activity is shaping the gas kinematics in the inner r ∼
400 pc of NGC 1068. The change in molecular line ratios inthe CND described above can be taken as independent evidencethat radiative and / or mechanical feedback from the AGN is atwork. The feedback of activity on the excitation and chemistryof molecular gas can be radiative, in PDR / XDR environments,or mechanical, in MDRs. In paper II we use the line ratio mapsderived in this work, complemented by those derived from inter-ferometric observations of NGC 1068 done in other moleculartracers, to model the excitation and chemistry of molecular gasin the di ff erent environments of the CND and the SB ring andfind a link between the di ff erent sources of energy present in thedisk.
8. Summary and conclusions
We have used ALMA to map the emission of a set of densemolecular gas tracers (CO(3–2), CO(6–5), HCN(4–3), HCO + (4–3), and CS(7–6)) and their underlying continuum emission inthe central r ∼ ∼ . (cid:48)(cid:48) − . (cid:48)(cid:48) (20–35 pc). These observations have greatly im-proved the sensitivity and spatial resolution of any previous in-terferometric study of this Seyfert, allowing us to make a signif-icant progress in our understanding of the fueling and the feed-back 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 dis-tribution of the dense molecular gas and dust emission fromthree main regions: the CND, the bar, and the SB ring. Wealso report the detection of CO(3–2) emission at di ff erent lo-cations over the interarm region. – The CND, which is fully resolved in the CO(6–5) map, is ahighly-structured closed elliptical ring of 350 pc-size. Sim-ilar to the dust continuum ring, the CO ring is noticeablyo ff -centered relative to the AGN. – Outside the CND, the CO(3–2) emission is detected alongtwo elongated o ff set lanes that run mostly parallel to the barmajor axis. We also detect a gas component with anomalousvelocities on the northeastern side of the disk at distances ∼
400 pc from the AGN. This CO feature has a counterpart indust continuum emission and is interpreted as the signatureof a bow-shock in the molecular disk. – As for dust emission, most of the CO(3–2) flux in the disk isdetected in the SB ring: a two-arm spiral structure that startsfrom the ends of the stellar bar and unfolds in the disk over ∼ ◦ in azimuth forming a pseudo-ring. – In stark contrast with the CO(3–2) map, most of the emis-sion in HCO + (4–3), HCN(4–3), and CS(7–6) stems fromthe CND. This reflects the di ff erent line ratios measured inthe CND and in the SB ring. The line ratio maps at thespatial resolution of ALMA show dramatic changes of upto an order of magnitude inside the CND. These changesare tightly correlated with the UV / X-ray illumination by theAGN, which varies significantly on these spatial scales, hint-ing at the footprint of AGN feedback. – We detect dust continuum and molecular line emission at theAGN. The emission peak is nevertheless shifted to the north-east in a structure that connects the AGN with the CND.We used the ALMA fluxes in Bands 9 and 7 together withNIR / MIR data to constrain the properties of the putativetorus using CLUMPY models. The fitted outer torus radius is R o = + − pc.; the torus mass is M torus = . ± . × M (cid:12) .The ALMA fluxes lie above the model predictions, an indi-cation that a non-negligible fraction of cold dust emissionnot associated with the torus was included in the 20 pc–resolution data of ALMA. – The Fourier decomposition of the gas velocity field derivedfrom the CO(3–2) data indicates that rotation is perturbed by
Article number, page 21 of 24 & A proofs: manuscript no. ngc1068-paper-I-inpress-pdflatex an inward radial flow in the SB ring and the bar region, asexpected if we are inside corotation of the perturbations inthese regions. However, the kinematics of the CND and thebow-shock arc are shaped by an outward radial pattern su-perposed to rotation. There is a tight spatial correlation be-tween the AGN ionized nebulosity, the radio jet plasma andthe molecular outflow signature identified in the CO velocityfield from r ∼
50 pc to r ∼
400 pc, which suggests that theoutflow is AGN driven. The molecular outflow is launchedwhen the ionization cone of the NLR sweeps the disk in theCND and farther out north in the bow-shock arc region. – The estimated outflow rate in the CND, d M / d t ∼ + − M (cid:12) yr − , an order of magnitude higher than the SFRmeasured at these radii, confirms that the outflow is AGNdriven. The molecular mass load rate implies a very shortgas depletion timescale of ≤ r ∼
400 pc of the galaxyon short timescales and at the same time regulate gas accretionin the CND. However, the molecular gas reservoir in the CNDis likely replenished on longer timescales by e ffi cient gas in-flow from the outer disk. The signature of gas inflow has beenidentified in the velocity field at larger radii and has been at-tributed to the combined action of the bar and the spiral arms.Self-regulation of star formation and gas accretion would be thuspossible by a combination of density wave-driven inflows andAGN-driven outflows. A similar scenario of self-regulation hasbeen invoked to be at work in gas-rich high-redshift disk galax-ies, analyzed in the numerical simulations of Gabor & Bour-naud (2014).While the gas kinematics are dominated by the molecularoutflow in the CND, on theoretical grounds it is nevertheless ex-pected that the gas should be fueling the AGN at smaller radii.The NIR data of Müller-Sánchez et al. (2009) showed ellipti-cal streamers that bridge the CND and the central engine. Theywere interpreted as a tentative evidence of ongoing AGN fuel-ing. However, with the spatial resolution of the observations pre-sented in this paper we cannot spatially resolve the kinematics ofthe dense molecular gas in the inner r ∼
50 pc. Obtaining higher-resolution data is important to studying the inflow / outflow signa-tures in the region that connects the CND with the AGN. A sig-nificant improvement in spatial resolution would also allow usto resolve with quasi-thermal molecular lines the region wherethe jet is known to be interacting with the gas, image the con-nection to the inner r ∼ . O megamaser disk (Gallimoreet al. 2001) and spatially resolve the emission from the putativemolecular torus.
Acknowledgements.
We acknowledge the sta ff of ALMA in Chile and theARC-people at IRAM-Grenoble in France for their invaluable help during thedata reduction process. This paper makes use of the following ALMA data:ADS / JAO.ALMA / NRAO, andNAOJ. We used observations made with the NASA / ESA Hubble Space Tele-scope, and obtained from the Hubble Legacy Archive.SGB and IM acknowl-edge support from Spanish grants AYA2010-15169 and from the Junta de An-dalucia through TIC-114 and the Excellence Project P08-TIC-03531. SGB, AL,and AF acknowledge support from MICIN within program CONSOLIDER IN-GENIO 2010, under grant ’Molecular Astrophysics: The Herschel and ALMAEra–ASTROMOL’ (ref CSD2009-00038). SGB, AU, LC, and PP acknowledgesupport from Spanish grant AYA2012-32295. FC acknowledges the EuropeanResearch Council for the Advanced Grant Program Num. 267399-Momentum.AAH acknowledges support from the Universidad de Cantabria through theAugusto G. Linares programme and from the Spanish Plan Nacional grantsAYA2009-05705-E and AYA2012-31447. CRA is supported by a Marie Curie Intra European Fellowship within the 7th European Community FrameworkProgramme (PIEF-GA-2012-327934). CRA also ackowledges financial supportfrom the Spanish Ministry of Science and Innovation (MICINN) through projectPN AYA2010-21887-C04.04 (Estallidos).
References
Aalto, S., Costagliola, F., van der Tak, F., & Meijerink, R. 2011, A&A, 527, A69Aalto, S., Garcia-Burillo, S., Muller, S., et al. 2012, A&A, 537, A44Aladro, R., Viti, S., Bayet, E., et al. 2013, A&A, 549, A39Alatalo, K., Blitz, L., Young, L. M., et al. 2011, ApJ, 735, 88Alloin, D., Pantin, E., Lagage, P. O., & Granato, G. L. 2000, A&A, 363, 926Alonso-Herrero, A., Ramos Almeida, C., Mason, R., et al. 2011, ApJ, 736, 82Arce, H. G., Shepherd, D., Gueth, F., et al. 2007, Protostars and Planets V, 245Arribas, S., Mediavilla, E., & Garcia-Lorenzo, B. 1996, ApJ, 463, 509Asensio Ramos, A., & Ramos Almeida, C. 2009, ApJ, 696, 2075Asensio Ramos, A., & Ramos Almeida, C. 2013, MNRAS, 428, 195Baker, A. J. 2000, Ph.D. ThesisBell, T. A., Viti, S., & Williams, D. A. 2007, MNRAS, 378, 983Bicknell, G. V., Dopita, M. A., Tsvetanov, Z. I., & Sutherland, R. S. 1998, ApJ,495, 680Binney, J., & Tremaine, S. 1987, Princeton, NJ, Princeton University Press, 1987,747 p.Bîrzan, L., McNamara, B. R., Nulsen, P. E. J., Carilli, C. L., & Wise, M. W. 2008,ApJ, 686, 859Bland-Hawthorn, J., Gallimore, J. F., Tacconi, L. J., Brinks, E., Baum, S. A.,Antonucci, R. R. J., & Cecil, G. N. 1997, Ap&SS, 248, 9Bock, J. J., Marsh, K. A., Ressler, M. E., & Werner, M. W. 1998, ApJ, 504, L5Bock, J. J., Neugebauer, G., Matthews, K., et al. 2000, AJ, 120, 2904Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132Boone, F., García-Burillo, S., Combes, F., et al. 2011, A&A, 525, A18Brinks, E., Skillman, E. D., Terlevich, R. J., & Terlevich, E. 1997, Ap&SS, 248,23Burtscher, L., Meisenheimer, K., Tristram, K. R. W., et al. 2013, A&A, 558,A149Canzian, B. 1993, ApJ, 414, 487Cecil, G., Dopita, M. A., Groves, B., et al. 2002, ApJ, 568, 627Chung, A., Yun, M. S., Naraynan, G., Heyer, M., & Erickson, N. R. 2011, ApJ,732, L15Cicone, C., Feruglio, C., Maiolino, R., et al. 2012, A&A, 543, A99Cicone, C., Maiolino, R., Sturm, E., et al. 2014, A&A, 562, A21Colombo, D., Meidt, S. E., Schinnerer, E., et al. 2014, ApJ, 784, 4Combes, F. 2003, Active Galactic Nuclei: From Central Engine to Host Galaxy,290, 411Combes, F. 2006, Astrophysics Update 2, 159Combes, F., García-Burillo, S., Casasola, V., et al. 2013, A&A, 558, A124Contopoulos, G. 1971, ApJ, 163, 181Crenshaw, D. M., & Kraemer, S. B. 2000, ApJ, 532, L101Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792Das, V., Crenshaw, D. M., Kraemer, S. B., & Deo, R. P. 2006, AJ, 132, 620Das, V., Crenshaw, D. M., & Kraemer, S. B. 2007, ApJ, 656, 699Dasyra, K. M., & Combes, F. 2012, A&A, 541, L7Davies, R. I., Müller Sánchez, F., Genzel, R., et al. 2007, ApJ, 671, 1388Davies, R., Genzel, R., Tacconi, L., Sánchez, F. M., & Sternberg, A. 2008, Map-ping the Galaxy and Nearby Galaxies, 144Devereux, N., Taniguchi, Y., Sanders, D. B., Nakai, N., & Young, J. S. 1994, AJ,107, 2006Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604Di Matteo, T., Colberg, J., Springel, V., Hernquist, L., & Sijacki, D. 2008, ApJ,676, 33Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866Emsellem, E., Fathi, K., Wozniak, H., Ferruit, P., Mundell, C. G., & Schinnerer,E. 2006, MNRAS, 365, 367Emsellem, E., Renaud, F., Bournaud, F., Elmegreen, B. & Combes, F. 2014,A&A, in prep.Esquej, P., Alonso-Herrero, A., González-Martín, O., et al. 2014, ApJ, 780, 86Faucher-Giguère, C.-A., & Quataert, E. 2012, MNRAS, 425, 605Feruglio, C., Maiolino, R., Piconcelli, E., et al. 2010, A&A, 518, L155Gabor, J. M., & Bournaud, F. 2014, arXiv:1402.4482Galliano, E., & Alloin, D. 2002, A&A, 393, 43Galliano, E., Pantin, E., Alloin, D., & Lagage, P. O. 2005, MNRAS, 363, L1Gallimore, J. F., Baum, S. A., & O’Dea, C. P. 2004, ApJ, 613, 794Gallimore, J. F., Baum, S. A., O’Dea, C. P., & Pedlar, A. 1996, ApJ, 458, 136Gallimore, J. F., Henkel, C., Baum, S. A., et al. 2001, ApJ, 556, 694García-Burillo, S., & Combes, F. 2012, Journal of Physics Conference Series,372, 012050García-Burillo, S., Combes, F., Hunt, L. K., et al. 2003, A&A, 407, 485
Article number, page 22 of 24. García-Burillo et al.: Molecular line emission in NGC 1068 imaged with ALMA
García-Burillo, S., Combes, F., Schinnerer, E., Boone, F., & Hunt, L. K. 2005,A&A, 441, 1011García-Burillo, S., Graciá-Carpio, J., Guélin, M. et al. 2006, ApJ, 645, L17García-Burillo, S., Usero, A., Fuente, A., et al. 2010, A&A, 519, A2Graciá-Carpio, J., García-Burillo, S., Planesas, P., & Colina, L. 2006, ApJ, 640,L135Graciá-Carpio, J., García-Burillo, S., Planesas, P., Fuente, A., & Usero, A. 2008,A&A, 479, 703Greenhill, L. J., & Gwinn, C. R. 1997, Ap&SS, 248, 261Greenhill, L. J., Gwinn, C. R., Antonucci, R., & Barvainis, R. 1996, ApJ, 472,L21Hailey-Dunsheath, S., Sturm, E., Fischer, J., et al. 2012, ApJ, 755, 57Haan, S., Schinnerer, E., Emsellem, E., et al. 2009, ApJ, 692, 1623Harada, N., Thompson, T. A., & Herbst, E. 2013, ApJ, 765, 108Helfer, T. T., & Blitz, L. 1995, ApJ, 450, 90Hollenbach, D., & McKee, C. F. 1989, ApJ, 342, 306Hönig, S. F., & Kishimoto, M. 2010, A&A, 523, A27Hönig, S. F., Prieto, M. A., & Beckert, T. 2008, A&A, 485, 33Hönig, S. F., Kishimoto, M., Gandhi, P., et al. 2010, A&A, 515, A23Hönig, S. F., Kishimoto, M., Tristram, K. R. W., et al. 2013, ApJ, 771, 87Hopkins, P. F., & Elvis, M. 2010b, MNRAS, 401, 7Hopkins, P. F., & Quataert, E. 2010a, MNRAS, 407, 1529Hopkins, P. F., & Quataert, E. 2011, MNRAS, 415, 1027Hopkins, P. F., Hayward, C. C., Narayanan, D., & Hernquist, L. 2012, MNRAS,420, 320Imanishi, M., Nakanishi, K., Tamura, Y., Oi, N., & Kohno, K. 2007, AJ, 134,2366Imanishi, M., & Nakanishi, K. 2013, AJ, 146, 91Israel, F. P. 2009a, A&A, 493, 525Israel, F. P. 2009b, A&A, 506, 689Iwasawa, K., Fabian, A. C., & Matt, G. 1997, MNRAS, 289, 443Jog, C. J., & Combes, F. 2009, Phys. Rep., 471, 75Jogee, S. 2006, Physics of Active Galactic Nuclei at all Scales, 693, 143Kamenetzky, J., Glenn, J., Maloney, P. R., et al. 2011, ApJ, 731, 83Kazandjian, M. V., Meijerink, R., Pelupessy, I., Israel, F. P., & Spaans, M. 2012,A&A, 542, A65Klaas, U., Haas, M., Müller, S. A. H., et al. 2001, A&A, 379, 823Kohno, K., Matsushita, S., Vila-Vilaró, B., Okumura, S. K., Shibatsuka, T., Ok-iura, M., Ishizuki, S., & Kawabe, R. 2001, The Central Kiloparsec of Star-bursts and AGN: The La Palma Connection, 249, 672Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511Krajnovi´c, D., Cappellari, M., de Zeeuw, P. T., & Copin, Y. 2006, MNRAS, 366,787Krips, M., Neri, R., García-Burillo, S., Martín, S., Combes, F., Graciá-Carpio, J.,& Eckart, A. 2008, ApJ, 677, 262Krips, M., Martín, S., Eckart, A., et al. 2011, ApJ, 736, 37Laurent, O., Mirabel, I. F., Charmandaris, V., et al. 2000, ISO Beyond the Peaks:The 2nd ISO Workshop on Analytical Spectroscopy, 456, 249Lepp, S., & Dalgarno, A. 1996, A&A, 306, L21Lira, P., Videla, L., Wu, Y., et al. 2013, ApJ, 764, 159López-Gonzaga, N., Ja ff e, W., Burtscher, L., Tristram, K. R. W., & Meisen-heimer, K. 2014, arXiv:1401.3248Macchetto, F., Capetti, A., Sparks, W. B., Axon, D. J., & Boksenberg, A. 1994,ApJ, 435, L15Maiolino, R., Gallerani, S., Neri, R., et al. 2012, MNRAS, 425, L66Maloney, P. R., Hollenbach, D. J., & Tielens, A. G. G. M. 1996, ApJ, 466, 561Mao, R.-Q., Schulz, A., Henkel, C., et al. 2010, ApJ, 724, 1336Matt, G., Fabian, A. C., Guainazzi, M., et al. 2000, MNRAS, 318, 173Mauersberger, R., Henkel, C., Walsh, W., & Schulz, A. 1999, A&A, 341, 256Meidt, S. E., Schinnerer, E., García-Burillo, S., et al. 2013, ApJ, 779, 45Meijerink, R., & Spaans, M. 2005, A&A, 436, 397Meijerink, R., Spaans, M., & Israel, F. P. 2007, A&A, 461, 793Meijerink, R., Spaans, M., Kamp, I., et al. 2013, Journal of Physical ChemistryA, 117, 9593Mor, R., Netzer, H., & Elitzur, M. 2009, ApJ, 705, 298Morganti, R., Frieswijk, W., Oonk, R. J. B., Oosterloo, T., & Tadhunter, C. 2013,A&A, 552, L4Müller-Sánchez, F., Davies, R. I., Genzel, R., Tacconi, L. J., Eisenhauer, F.,Hicks, E. K. S., Friedrich, S., & Sternberg, A. 2009, ApJ, 691, 749Müller-Sánchez, F., Prieto, M. A., Hicks, E. K. S., et al. 2011, ApJ, 739, 69Murray, N., Quataert, E., & Thompson, T. A. 2005, ApJ, 618, 569Nenkova, M., Sirocky, M. M., Ivezi´c, Ž., & Elitzur, M. 2008a, ApJ, 685, 147Nenkova, M., Sirocky, M. M., Nikutta, R., Ivezi´c, Ž., & Elitzur, M. 2008b, ApJ,685, 160Neufeld, D. A., & Dalgarno, A. 1989a, ApJ, 340, 869Neufeld, D. A., & Dalgarno, A. 1989b, ApJ, 344, 251Ogle, P. M., Brookings, T., Canizares, C. R., Lee, J. C., & Marshall, H. L. 2003,A&A, 402, 849Papadopoulos, P. P., & Seaquist, E. R. 1999, ApJ, 514, L95 Papadopoulos, P. P., van der Werf, P., Xilouris, E., Isaak, K. G., & Gao, Y. 2012,ApJ, 751, 10Pérez-Beaupuits, J. P., Aalto, S., & Gerebro, H. 2007, A&A, 476, 177Pérez-Beaupuits, J. P., Spaans, M., van der Tak, F. F. S., et al. 2009, A&A, 503,459Pilyugin, L. S., Vílchez, J. M., & Contini, T. 2004, A&A, 425, 849Pilyugin, L. S., Thuan, T. X., & Vílchez, J. M. 2007, MNRAS, 376, 353Planesas, P., Scoville, N., & Myers, S. T. 1991, ApJ, 369, 364Prieto, M. A., Reunanen, J., Tristram, K. R. W., et al. 2010, MNRAS, 402, 724Raban, D., Ja ff e, W., Röttgering, H., Meisenheimer, K., & Tristram, K. R. W.2009, MNRAS, 394, 1325Ramos Almeida, C., Levenson, N. A., Rodríguez Espinosa, J. M., et al. 2009,ApJ, 702, 1127Ramos Almeida, C., Levenson, N. A., Alonso-Herrero, A., et al. 2011a, ApJ,731, 92Ramos Almeida, C., Sánchez-Portal, M., Pérez García, A. M., et al. 2011b, MN-RAS, 417, L46Ramos Almeida, C., Alonso-Herrero, A., Levenson, N. A., et al. 2014, MNRAS,444Rand, R. J., & Wallin, J. F. 2004, ApJ, 614, 142Sandstrom, K. M., Leroy, A. K., Walter, F., et al. 2013, ApJ, 777, 5Scoville, N. Z., Matthews, K., Carico, D. P., & Sanders, D. B. 1988, ApJ, 327,L61Sempere, M. J., Garcia-Burillo, S., Combes, F., & Knapen, J. H. 1995, A&A,296, 45Schartmann, M., Meisenheimer, K., Camenzind, M., et al. 2008, A&A, 482, 67Schinnerer, E., Eckart, A., Tacconi, L. J., Genzel, R., & Downes, D. 2000, ApJ,533, 850Schoenmakers, R. H. M. 1999, Ph.D. ThesisSchoenmakers, R. H. M., Franx, M., & de Zeeuw, P. T. 1997, MNRAS, 292, 349Sternberg, A., Genzel, R., & Tacconi, L. 1994, ApJ, 436, L131Storchi-Bergmann, T., Ri ff el, R. A., Ri ff el, R., et al. 2012, ApJ, 755, 87Strong, A. W., & Mattox, J. R. 1996, A&A, 308, L21Strong, A. W., Bloemen, J. B. G. M., Dame, T. M. et al. 1988, A&A, 207, 1Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163Sturm, E., González-Alfonso, E., Veilleux, S., et al. 2011, ApJ, 733, L16Tacconi, L. J., Genzel, R., Blietz, M., Cameron, M., Harris, A. I., & Madden, S.1994, ApJ, 426, L77Tafalla, M. 2013, Astronomical Society of the Pacific Conference Series, 476,177Tafalla, M., Santiago-García, J., Hacar, A., & Bachiller, R. 2010, A&A, 522, A91Takano, S., Nakajima, T., Kohno, K., et al. 2014, arXiv:1406.0782Taniguchi, Y. 1999, ApJ, 524, 65Taniguchi, Y. 2013, Astronomical Society of the Pacific Conference Series, 477,265Taniguchi, Y., & Wada, K. 1996, ApJ, 469, 581Tecza, M., Thatte, N., & Maiolino, R. 2001, Galaxies and their Constituents atthe Highest Angular Resolutions, 205, 216Thompson, R. I., Chary, R.-R., Corbin, M. R., & Epps, H. 2001, ApJ, 558, L97Tomono, D., Doi, Y., Usuda, T., & Nishimura, T. 2001, ApJ, 557, 637Tomono, D., Terada, H., & Kobayashi, N. 2006, ApJ, 646, 774Toomre, A. 1981, Structure and Evolution of Normal Galaxies, 111Trachternach, C., de Blok, W. J. G., Walter, F., Brinks, E., & Kennicutt, R. C., Jr.2008, AJ, 136, 2720Tsai, M., Hwang, C.-Y., Matsushita, S., Baker, A. J., & Espada, D. 2012, ApJ,746, 129Usero, A., García-Burillo, S., Fuente, A., Martín-Pintado, J., & Rodríguez-Fernández, N. J. 2004, A&A, 419, 897Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005, ARA&A, 43, 769Viti, S., García-Burillo, S, Combes, F. et al. 2014, A&A, submitted (paper II)Wilson, A. S., & Ulvestad, J. S. 1987, ApJ, 319, 105Wong, T., Blitz, L., & Bosma, A. 2004, ApJ, 605, 183Yamada, M., Wada, K., & Tomisaka, K. 2007, ApJ, 671, 73Young, A. J., Wilson, A. S., & Shopbell, P. L. 2001, ApJ, 556, 6 Article number, page 23 of 24 & A proofs: manuscript no. ngc1068-paper-I-inpress-pdflatex Observatorio Astronómico Nacional (OAN)-Observatoriode Madrid, Alfonso XII, 3, 28014-Madrid, Spain e-mail: [email protected] Observatoire de Paris, LERMA, CNRS, 61 Av. de l’Observatoire,75014-Paris, France Department of Earth and Space Sciences, Chalmers University ofTechnology, Onsala Observatory, 439 94-Onsala, Sweden Institut de Radio Astronomie Millimétrique (IRAM), 300 rue dela Piscine, Domaine Universitaire de Grenoble, 38406-St.Martind’Hères, France Department of Physics and Astronomy, UCL, Gower Place, LondonWC1E 6BT, UK Instituto de Física de Cantabria, CSIC-UC, E-39005 Santander,Spain. Augusto G. Linares Senior Research Fellow. INAF-Osservatorio Astrofisico di Arcetri, Largo Enrico Fermi 5,50125-Firenze, Italy Max-Planck-Institut für Astronomie, Königstuhl, 17, 69117-Heidelberg, Germany Department of Physics and Astronomy, Rutgers, The State Univer-sity of New Jersey, Piscataway, NJ 08854, USA Université de Toulouse, UPS-OMP, IRAP, 31028, Toulouse, France INAF - Istituto di Radioastronomia, via Gobetti 101, 40129,Bologna, Italy Centro de Astrobiología (CSIC-INTA), Ctra de Torrejón a Ajalvir,km 4, 28850 Torrejón de Ardoz, Madrid, Spain Instituto de Astrofísica de Andalucía (CSIC), Apdo 3004, 18080-Granada, Spain I. Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77,50937, Köln, Germany Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69,53121, Bonn, Germany Astronomy Department, King Abdulazizi University, P. O. Box80203, Jeddah 21589, Saudi Arabia Institute for Astronomy, Department of Physics, ETH Zurich, CH-8093 Zurich, Switzerland Instituto de Astrofísica de Canarias, Calle Vía Láctea, s / n, E-38205La Laguna, Tenerife, Spain Departamento de Astrofísica, Universidad de La Laguna, E-38205,La Laguna, Tenerife, Spain Kapteyn Astronomical Institute, University of Groningen, PO Box800, NL-9700 AV Groningen Max-Planck-Institut für extraterrestrische Physik, Postfach 1312,85741-Garching, Germany Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Lei-den, Netherlands
Appendix A: CLUMPY torus models
Over the past few years, a number of studies have demonstratedclumpy torus models (Nenkova et al. 2008a, 2008b; Schart-mann et al. 2008; Hoenig & Kishimoto 2010) represent wellthe nuclear infrared emission of local Seyfert galaxies, providedthere is no contamination from extended dust components notrelated to the torus (e.g., foreground absorbing dust, opticallythin emitting dust in ionization cones, etc). More specifically,these models reproduce the non-stellar NIR and MIR SEDs andMIR spectroscopy of local Seyfert galaxies (Ramos Almeida etal. 2009, 2011a; Hoenig et al. 2010; Alonso-Herrero et al. 2011;Lira et al. 2013) and PG quasars (Mor et al. 2009). By model-ing the nuclear infrared data of AGN we can constrain the torusgeometry, dust properties, and distribution, and the AGN bolo-metric luminosity. Moreover, adding far-infrared nuclear photo-metric points or even upper limits has proven to be useful toconstrain the torus sizes, especially if the dusty clumps havea 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 MIRSED and MIR spectroscopy of NGC 1068 using the so-calledCLUMPY torus models of Nenkova et al. (2008a, 2008b). Theparameters in the CLUMPY models are the torus size Y definedas the ratio between the outer radius R o and the inner radius R sub , the torus angular size σ , the viewing angle i , the numberof clouds along the equatorial direction N , the optical depth ofthe clouds τ V , and the index q of the radial distribution of theclouds ∝ 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 basedon the 12 µ m interferometric sizes of a few parsecs inferred fornearby AGN (Burtscher et al. 2013).We can use the ALMA Band 9 and Band 7 continuum ther-mal fluxes at 435 µ m and 860 µ m, respectively, to investigatewhether we can set further constraints on the torus propertiesof NGC 1068. To do so we also added the NGC 1068 SED andMIR spectroscopy presented in Alonso-Herrero et al. (2011). Asin Alonso-Herrero et al. (2011) we used the detection of a maserin the nuclear region of NGC 1068, which implies a close toedge-on view of the AGN accretion disk, to set the followingprior for the viewing angle i = − ◦ . Because we are includ-ing the ALMA far-infrared photometry we allowed the full rangefor the torus size Y = − M torus = π m H (sin σ ) N eqtorus R Y I q ( Y ) (A.1)where the function I q ( Y ) = q = N eqtorus is the col-umn density along the equatorial direction. The inner radius ofthe torus R sub was computed from the AGN bolometric lumi-nosity inferred from the fit and we assumed a gas-to- A V ratioof N H / A V = × mol cm − mag − , taken from Bohlin etal (1978). The typical scatter reported by Bohlin et al. (1978) forthis factor in H dominated gas clouds is 30%.The results of the fit in NGC 1068 and their implications arediscussed in Sect. 4.1.2.5