Open Access
Issue
A&A
Volume 634, February 2020
Article Number A1
Number of page(s) 12
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201936255
Published online 28 January 2020

© GRAVITY Collaboration et al. 2020

Licence Creative Commons
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Open Access funding provided by Max Planck Society.

1. Introduction

NGC 1068 is often regarded as an archetypical Seyfert 2 galaxy (e.g. Bland-Hawthorn et al. 1997). In particular, observations of NGC 1068 played an important role in the early development of unified models of active galactic nuclei (AGN). In this model, opaque material in the equatorial plane obscures the supermassive black hole, accretion disc, and broad emission lines, so that the AGN can only be observed directly from polar directions (Miller & Antonucci 1983; Antonucci & Miller 1985). Following the inception of the torus concept (Antonucci 1982, 1984) and the first efforts to explain its physical properties (Krolik & Begelman 1988), there have been numerous observations of the central structures in this and other AGN from across the electromagnetic spectrum and covering many different tracers. Ideas about the geometry of the central obscurer have evolved and modified significantly over the years. Indeed, in the literature the term “torus” is used to mean a variety of different things. While it was always understood that this structure might be clumpy (Krolik & Begelman 1988), early calculations adopted a smooth geometry (Pier & Krolik 1992). The subsequent incorporation of a clumpy structure led to a suite of models that provide very good fits to the near-to-mid infrared spectral energy distribution (SED) of many AGN (Nenkova et al. 2008; Hönig et al. 2006; Stalevski et al. 2012), and which are also consistent with dust reverberation measurements (Koshida et al. 2014). We refer to Netzer (2015) for a detailed review. Models and observations have achieved a level of detail where it is possible to apply a wide range of tests to the properties and assumptions used in creating them. Additionally, the data we present here allow a new test to be applied. It concerns the observed morphology of the hot dust close around the AGN. To provide a context for this, we highlight a few key observations pertinent to understanding the central few parsecs around the black hole in NGC 1068 as follows.

The detection of broad permitted emission lines in polarised light (Antonucci & Miller 1985) represents a cornerstone of the AGN unification scheme (Antonucci 1993), which postulates that the observed properties of an AGN mainly depend on inclination. Cameron et al. (1993) pointed out that other structures must also be present on small scales in NGC 1068 since the mid-infrared (MIR) continuum is extended on scales of ∼1″ rather than dominated by a central source. Using higher resolution data, Bock et al. (2000) showed there was a central < 0.2″ core, which they associated with the torus, and a tongue extending ∼1″ along a north-south axis that mimics the inner part of the radio jet (Gallimore et al. 1996a), which is thought to be launched in a northerly direction and bends to the north-east when it encounters a dense molecular cloud.

Early K-band speckle and interferometry work found hot dust emission coming from multiple components on scales of 400 mas to a few mas (Weigelt et al. 2004; Wittkowski et al. 2004). Observations with MIDI at 10 μm indicated the presence of multiple dust components with a central 600–800 K elongated dust core. The first component is a 800 K, disc-like structure about 1.4 × 0.5 pc in size at a position angle (PA) of approximately −45° (Jaffe et al. 2004; Raban et al. 2009; López-Gonzaga et al. 2014). A second, cooler ∼300 K component, 3 × 4 pc in size, was identified with the “body” of the torus (Raban et al. 2009).

The presence of a complex structure consisting of multiple components is supported by X-ray spectra, indicating a compact high column density NH = 1025 cm−2 X-ray reflector with lower column NH = 1023 cm−2 clouds on more extended scales (Bauer et al. 2015). This does not necessarily imply anything about a dusty structure, since the X-ray absorption could be through dust-free clouds in the optically hidden broad line region (BLR). The lack of Brα detection led Lutz et al. (2000) to put a lower limit on the extinction corresponding to AV ∼ 100 mag, or an equivalent column of nH ∼ 1023 cm−2 for a standard dust-to-gas ratio.

A maser disc is seen extending up to ∼15 mas (1 pc) from the AGN, with velocities reaching ±330 km s−1 (Gallimore et al. 1996b; Greenhill & Gwinn 1997; Gallimore et al. 2004). In addition, Greenhill & Gwinn (1997) noted that the ∼ − 45° position angle of the disc rotation axis differs from that of the jet by 30–40°. As such, there is clearly some warping on sub-pc scales in NGC 1068, and any description of the geometry needs to be associated with a particular structure or radial scale. The maser disc can be described by a thin edge-on (i ∼ 80°) disc. Studies based on SED modelling also favour torus inclinations between 70° (Hönig et al. 2007; Lopez-Rodriguez et al. 2018) and 90° (Hönig et al. 2008). Attempts to model the three dimensional orientation of the system found that the northern narrow line region (NLR) cone is directed towards the observer and the southern part away from the observer (Packham et al. 1997; Kishimoto 1999).

Near-infrared (NIR) integral field adaptive optics observations revealed streamers in the H2 1-0 S(1) line, tracing hot molecular gas, which were interpreted as infalling material on a few tens of parsecs scales (Müller Sánchez et al. 2009). They also provided the first glimpse of a molecular structure on a ∼10 pc scale.

Gratadour et al. (2015) obtained polarimetric imaging with SPHERE. This showed a bicone structure perpendicular to the maser disc, probably related to outflows, as well as an elongated patch aligned with the central molecular structure but extending further out (R ≈ 25 pc), which they interpreted as the less dense outskirt of the putative torus.

ALMA observations of the 432 μm continuum and CO(6-5) molecular line resolved a 7–10 pc thick-disc like structure (Gallimore et al. 2016; García-Burillo et al. 2016). As expected for a Seyfert 2 galaxy, observational data argue for a system close to edge-on. Based on the linewidth of the spatially resolved HCN and HCO+ lines in more recent observations of this structure, Imanishi et al. (2018) argue for turbulent motions with σ = 70–80 km s−1, of similar order as the rotational velocity. Impellizzeri et al. (2019) find similar line widths close to the nucleus. García-Burillo et al. (in prep.) detected the molecular disc in CO(2-1) and (3-2) at a PA ∼115 ° with an extension of 30 pc.

Taken together, these observations, while still being consistent with and requiring the existence of a nuclear obscuring structure, are incompatible with geometrically thick clumpy torus models that attempt to account for all of the nuclear near-to-mid infrared continuum. Vollmer et al. (2018) propose a resolution to this discrepancy by constructing a model which has a thin disc on scales ≲1 pc and a thick disc on scales of 1–10 pc, the structure of which is determined by inflow from larger scales in the host galaxy. In their model, a magnetocentrifugal wind is launched at the boundary between the discs, and accounts for the elongated polar structures seen in many MIR interferometric measurements (López-Gonzaga et al. 2016). This multi-component model was able to account for a wide variety of detailed observations when applied to NGC 1068. A similar approach has been described by Hönig (2019) who posits a dusty disc with small to moderate scale height, through which gas flows inwards and is then unbound in a dusty wind launched at the inner puffed-up rim by radiation pressure from the AGN.

2. Prerequisites and observations

2.1. Bolometric luminosity and sublimation distance

For consistency with most of the NGC 1068 literature, we adopt the “standard” distance of 14.4 Mpc (Bland-Hawthorn et al. 1997) but note that d ≈ 16.5 Mpc would arise from Virgocentric infall velocity 1117 km s−1 (Mould et al. 2000, NED1) and the Hubble constant (H0 = 67.8 ± 0.9 km s−1Mpc−1) of Planck Collaboration XIII (2016).

The estimated black hole mass of up to MBH ≈ 1.7 × 107M (see Table 1) corresponds, for a solar metallicity gas, to an Eddington luminosity of LEdd ≈ 2.5 × 1045 erg s−1.

Table 1.

Basic parameters of NGC 1068 taken from the literature.

The literature reports bolometric luminosities between Lbol ≈ 0.5 × 1045 erg s−1 (Lopez-Rodriguez et al. 2018) based on SED modelling and Lbol ≈ 1.6 × 1045 erg s−1 based on MIR data (Raban et al. 2009). Prieto et al. (2010) on the other hand found a significantly lower luminosity of ∼9 × 1043 erg s−1.

Our estimates of Lbol are based on various X-ray and optical observations of the source. When we compared typical luminosity probes, we found that the X-ray based Lbol estimates are significantly lower than the estimates from ionised lines. Based on detailed modelling, Bauer et al. (2015) infer an X-ray luminosity of L2 − 10 keV ≈ 2 × 1043 erg s−1, which includes a number of X-ray components and a detailed SED fitting of the source. Using broadband X-ray data Ricci et al. (2017) find an intrinsic L14 − 195 keV ∼ 5 × 1042 erg s−1. Moreover, recent observations of Marinucci et al. (2016) find X-ray variability, due to column density variations. The authors infer an intrinsic L2 − 10 keV ≈ 7 × 1043 erg s−1. The X-ray estimates combined with the Marconi et al. (2004) bolometric correction factor, give Lbol = 0.4–1.4 × 1045 erg s−1. An estimate based on the nuclear luminosity of the 12 μm continuum can be made using the relations of Asmus et al. (2015), together with additional relations from Winter et al. (2012). These yield Lbol ∼ 1.8 × 1045 erg s−1. For the optical data we use estimates based on two combinations of three narrow emission lines: Hβλ4861, [OIII] λ5007, [OI] λ6300. The expressions are taken from Netzer (2009) and make use of reddening corrected line fluxes and galactic type extinction. The line fluxes are taken from Storchi-Bergmann et al. (1995). The resulting bolometric luminosities are 2.3 × 1045 erg s−1 based on the [OIII] and [OI] lines, and 4.7 × 1045 erg s−1 based on the Hβ and [OIII] lines. An early estimate, using scattered light from an off-nuclear cloud that arguably sees the hidden nucleus, also favoured luminosities rather in excess of 1045 erg s−1 (Miller et al. 1991). The large range of luminosity estimates Lbol = 0.4–4.7 × 1045 erg s−1 reflects the extreme extinction of the source. In any case, NGC 1068 is radiating close to the Eddington limit.

Theoretical and observational work argue for a strong correlation between the NIR size and the AGN luminosity (Suganuma et al. 2006; Kishimoto et al. 2011; Koshida et al. 2014; GRAVITY Collaboration 2019). In the standard picture, the NIR emission originates from the dust sublimation region, which is determined by the dust sublimation temperature Tsub. The temperature at which dust sublimates depends on the grain size and composition and ranges from Tsub ≈ 1500 K for silicate (Si) grains to Tsub ≈ 2000 K for graphite (C) grains (Baskin & Laor 2018). By assuming a grain size distribution typical for the interstellar medium (ISM), the mean sublimation distance can be calculated for graphite and silicate dust as [pc] with aC = 0.5 and aSi = 1.3 (e.g. Barvainis 1987; Netzer 2015; for grain size dependence see Baskin & Laor 2018). However, there is evidence both from SED fitting (Mor & Netzer 2012) and from dust reverberation mapping (Kishimoto et al. 2007; Koshida et al. 2014) that, close to the AGN, one may find only larger graphite grains. Since, at the same distance from the AGN, large grains are cooler than small grains, one might naturally expect to find a distribution of large graphite grains at the smallest radii. In addition changes in the brightness of the AGN may lead to changes in the dust sublimation radius (Kishimoto et al. 2013) or instead, depending on the dust distribution, its temperature (Schnülle et al. 2015). As such, the sublimation radius is really rather a sublimation region according to grain size and species (e.g. Hönig & Kishimoto 2017; Baskin & Laor 2018). Empirical studies have found R2.2 μm ≈ 0.4 pc × (LAGN/1046 erg s−1)1/2 (assuming the V-band to Rsub relation from Koshida et al. 2014, and using the conversion LAGN ∼ 8L5500 Å from Netzer 2015), which corresponds to a good approximation to the sublimation radius of graphite. For NGC 1068 the observed relation from Koshida et al. (2014) predicts a dust sublimation radius in the range R2.2 μm ≈ 0.08–0.27 pc, for the luminosity range quoted in Table 1.

2.2. Measured column densities

Numerous studies have tried to estimate the column density obscuring the line-of-sight to the central region of NGC 1068. Early CO/HCN measurements of the central few arcsec suggested column densities of NH ∼ 0.2 − 1 × 1023 cm−2 (Planesas et al. 1991; Tacconi et al. 1994; Sternberg et al. 1994). Because the beam is relatively large, these measurements may be affected by the complex structures in the nuclear region, for example the bright CO emission about 1″ west of the nucleus. More recent higher resolution observations focus on the central molecular structure. Imanishi et al. (2018) found gas masses of 2 × 106M within 14 pc based on ALMA observations of HCN 3-2 and HCO tracers, which imply a column density of NH ∼ 2 × 1023 cm−2. ALMA 432 μm continuum observations by García-Burillo et al. (2016) inferred gas masses of Mgas = 105M within a structure of 7 pc diameter, suggesting a column density of 8 × 1022 cm−2.

In the infrared, ISO upper limits on the Brackett α broad-line emission led Lutz et al. (2000) to place a lower limit on the extinction ABrα, 4.05 ≥ 2.4, which corresponds to AK ≥ 6. For an ISM-like gas-to-dust ratio, the optical extinction can be converted into a column density. Based on X-ray calibrations, the relationship between the hydrogen column density NH and the optical extinction AV is NH (cm−2) = a × 1021AV (mag), where a has a value ranging from about 1.79–2.21 (Predehl & Schmitt 1995; Güver & Özel 2009). Therefore, the observed extinction implies a hydrogen column density of NH ≳ (1.8 ± 0.6)×1023 cm−2.

The gas masses found by Müller Sánchez et al. (2009) imply a column density of NH ∼ 3 × 1024 cm−2 assuming a source size of 0.2″ and a hot-to-cold gas mass ratio of . Lopez-Rodriguez et al. (2018) estimated NH = 5 × 1023 cm−2 based on 2–400 μm clumpy torus SED fitting over a 0.4″ aperture.

Generally larger column densities are found in the X-rays, which could partially originate from dust-free gas. Various studies found consistently a column density of NH ∼ 1 × 1025 cm−2 (with a 50% coverage) based on 3–100 keV data (Matt et al. 1997; Bauer et al. 2015) as well as based on broad-band 0.3–150 keV X-ray model fitting (Ricci et al. 2017).

Overall the different measurements argue for a mean column density of NH ∼ 1023–1024 cm−2 within the central few parsecs to tens of parsecs. The higher column density towards the X-ray emitting corona in the very centre may largely originate in dust-free gas (Risaliti et al. 2002; Netzer 2015; Davies et al. 2015).

2.3. Extinction correction

Throughout the paper we have assumed the same extinction law as can be found in the highly absorbed Galactic Centre region (Fritz et al. 2011). For wavelength < 2.6 μm we assumed an extinction slope of A(λ)∝λ−2.05. At larger wavelengths the extinction shows spectral features, which we interpolated based on the extinction published by Fritz et al. (2011). Some useful conversions are: AV/ABrγ ≈ 20, ABrγ/ABrα = 2.5, ABrγ/A8.76 μm = 1.22 and ABrγ/A12.4 μm = 1.86.

2.4. Interferometric observations

The observations were carried out with the GRAVITY instrument, which interferometrically combines either four 8 m Unit Telescopes (UTs), or four 1.8 m Auxiliary Telescopes (ATs), of the European Southern Observatory (ESO) Very Large Telescope (VLT). Observations of NGC 1068 included both of these, and provided 3 milli-arcsec (mas) resolution K-band (2.2 μm) continuum imaging. For a detailed description of the instrument and the data analysis we refer to GRAVITY Collaboration (2017). In brief, the light of the four telescopes is extracted into mono-mode fibres (Pfuhl et al. 2012, 2014) for two positions on the sky and then interfered in the beam combiner for all six baselines of the interferometer. The fibre diameters of 56 mas and 250 mas are matched to the 2.2 μm diffraction limited beams of an 8 m UT and 1.8 m AT respectively. In a “normal” single-field observation, the light of the science target is split 50/50 and fed into two channels, the fringe-tracker and the science spectrometer. The fringe-tracker is used to correct the fast atmospheric optical path variations with typical integration times of 1–10 milli-seconds (Lacour et al. 2019). This allows coherent integration times with the science spectrometer up to ∼100 s. The interferometric observations of NGC 1068 are challenging due to the complex nature of the object. The central 1″ shows extended continuum as well as line emission. The bright core has been partially resolved with speckle data obtained at the 6-m telescope of the Special Astrophysical Observatory in Russia (Weigelt et al. 2004) as well as with adaptive optics assisted 8-m telescope data from the VLT (e.g. Rouan et al. 2019). Overall there is a significant contribution of coherent flux on all spatial scales, which can be probed with the VLT interferometer (VLTI). The visibilities range from 44% on the shortest AT baselines to ∼0.1% on the longest UT baselines. In order to reach sufficient signal-to-noise, we decided to feed all the light to the fringe-tracker channel. This provided us with a factor two more signal on the fringe-tracker at the expense of spectral resolving power (The fringe-tracker has five spectral channels across the K-band, corresponding to R ∼ 20).

The observations2 were taken on 20 November 2018 (1:00–5:20 UTC) with the UT array and on 23 December 2018 (1:30–3:30 UTC) with the AT compact array. The UT observations on November 20th had excellent conditions with seeing measured by the differential image motion monitor (DIMM) between 0.47–0.65″ and 7–13 ms (optical) coherence time. The conditions during the AT observations were almost as good with the seeing ranging from 0.43 to 0.75″ and a coherence time between 4 and 8 ms. In total we obtained 45 min with the UTs and 30 min with the ATs of integration time on source. The science observations were bracketed by observations of close calibrator stars HIP 54 (mK = 7.95), HI P16739 (mK = 8.93) and HIP 17272 (mK = 8.34). We used the standard pipeline (Lapeyrere et al. 2014; GRAVITY Collaboration 2017) to reduce and calibrate the data.

2.5. Adaptive optics based photometry

The combination of four telescopes allows us to simultaneously obtain the coherent as well as the incoherent (absolute) flux of the object. At the same time, the acquisition camera (Anugu et al. 2018) provides an estimate of the Strehl ratio (in H-band). By comparing the observed flux of NGC 1068 with observations of known calibrator stars HIP 17272 (mK = 8.34), HIP 16739 (mK = 8.93) and HIP 54 (mK = 7.95) we derived an effective magnitude of NGC 1068. As previously described, NGC 1068 is resolved even on scales of a single telescope. This means that the fibre-injected flux depends on the size of the fibre mode on sky. With the UTs, we find a total flux of 310 ± 80 mJy within the fibre beam (Dfibre ≈ 56 mas). This corresponds to a magnitude of mK = 8.33 ± 0.25. Overall, our flux measurement is in good agreement with the previously reported flux of the central component 350 mJy by Weigelt et al. (2004) and with Rouan et al. (2019), who found a central mK = 8.65.

Weigelt et al. (2004) used speckle data to derive the flux distribution of the central few 100 mas. They modelled the data with two elongated Gaussian components, a central component with a full width at half maximum (FWHM) of 18/39 mas contributing F1 = 350 mJy at PA ≈ −20°, and a larger structure with a FWHM of 400 mas that had a quite uncertain flux estimate of F2 = 30–230 mJy. The first component has a FWHM significantly smaller than the fibre mode (56/250 mas). Therefore the component is fully injected into our GRAVITY fibres, which agrees with our flux measurement. The 400 mas component however is only partially injected. Less than 30% (9–69 mJy) of the large component is injected into the fibres using ATs and only ∼2% (< 5 mJy) is injected into the fibres using the UTs.

3. Interferometric data and image reconstruction

The interferometric data reveal significant structures on all spatial scales. In Fig. 1 we show the calibrated visibilities from the AT and UT observations. For comparison we show the early speckle data from Weigelt et al. (2004) obtained with a single telescope. The observed visibilities cannot be approximated with a simple symmetric model like a Gaussian or a uniform disc. Also, as shown in the figure, a ring model fails to reproduce the data. This is however not surprising, given the large asymmetries in the source intensity distribution indicated by the highly non-zero closure phases (see Fig. 2). The complexity of the visibility and closure phase data prevent us from using simple model fitting to interpret the observations. Encouraged by the excellent quality of the data and the large uv coverage, we therefore turned to standard image reconstruction codes to analyse the data.

thumbnail Fig. 1.

Squared visibilities measured for NGC 1068. The light coloured filled data points were obtained with the AT array in compact configuration (spatial frequencies 5 <  f <  17 Mλ.). The high spatial frequency data were taken with the UT array. The small inset shows the corresponding UV coordinates. The open blue and red symbols (spatial frequencies f <  4 Mλ) represent speckle data from Weigelt et al. (2004), which are shown only for comparison; they were not used for the image reconstruction. The data point and error bar at a spatial frequency of 21 Mλ, lying above our new data, shows the only previous NIR VLTI detection of NGC 1068 (Wittkowski et al. 2004). For comparison, the grey continuous line shows the visibility of a thin ring-shaped emission region with a radius of 3.3 mas and a width of 0.5 mas scaled to a maximum visibility of 0.3. While this reflects the general shape of the visibility distribution, it is a poor fit and indicates that more complex structures are present.

Open with DEXTER

thumbnail Fig. 2.

Closure phase data of NGC 1068. See Fig. 1 for the explanations of the symbols. The highly non-zero closure phases indicate strong source asymmetries on all spatial scales.

Open with DEXTER

3.1. An image of the dust sublimation region

We used the publicly available MiRA tool (Thiébaut 2008) to reconstruct an image based on the visibility and closure phase data. A “grey image” was reconstructed assuming that the intensity distribution on sky I(α, β) is independent of the spectral channel across the observed K-band. We normalised the visibility by a factor 0.5 to match the visibilities at the smallest spatial frequencies. This means that we reconstruct 50% of the photometric flux in our image (≈155 mJy). The remaining flux originates at larger scales and/or is smoothly distributed. Figure 3 shows the reconstructed image using MiRA with a l2l1 smoothness prior (threshold 5 × 10−4) on a 200 × 200 pixel grid with a spatial scale of 0.4 mas/pixel. As an initial guess, we used a single Gaussian image with 5 mas FWHM. The final reconstructed image was then convolved with a Gaussian beam of FWHM 3.1 × 1.1 mas and PA = 48.6°, which corresponds to 0.8 times the (UT array) interferometric beam. In order to asses the systematic uncertainty of the image, we repeated the reconstruction using different penalty functions (for more information, refer to Thiébaut 2008; Thiébaut & Young 2017) such as a smoothness prior, (i.e. a penalty on sharp transitions in the image), a compact prior (i.e. penalty on large radii) and a l2l1 smoothness prior (i.e. a prior able to maintain sharp transitions in the image). We stepwise optimised the gain and hyper parameters, where applicable. The resulting images varied but the central structure hardly changed in shape and amplitude. In a second step, we determined the reliability of individual detections in the image by repeating the reconstruction with temporal and bandwidth restricted subsets of the data (with and without the bluest 1.995 μm channel which is partially contaminated by metrology laser light). Spurious sources, which were not present throughout the repeated imaging attempts likely correspond to the reconstruction of noise in the data. This noise floor has a root-mean-square (rms) of ∼7 mJy beam−1. In addition, the image reconstruction methods used here, where calibrated phases are not available, assume positive fluxes. Therefore, a reconstructed image never contains negative values, unlike single-dish dark- or sky-subtracted images. Numerous spurious noise sources can be seen in Fig 3. The central ring-like structure is highly significant at 5–7σ. The surrounding features to the north-east and north-west are significant at 3–4σ. Figure 4 represents a zoom on the significant central structure.

thumbnail Fig. 3.

Reconstructed image of the central 5.1 × 5.1 pc of NGC 1068. The elliptical beam, which was used to convolve the image, is shown in the lower left. Details of the reconstruction, which only allows positive sources, are given in Sect. 3.1. The central structures are robust to changing the reconstruction parameters. The fainter features towards the edge of the field, with fluxes close to the noise floor of ∼7 mJy beam−1 rms (i.e. dark cyan in this colour scale), are uncertain.

Open with DEXTER

thumbnail Fig. 4.

Reconstructed image of the inner 2.1 × 2.1 pc of NGC 1068 (i.e. the central region of Fig. 3). The best fitting circular ring is shown in dashed white (r = 0.24 pc (3.5 mas). The dust sublimation sizes for the range of bolometric luminosities (∼0.08–0.27 pc) are shown in orange dotted ellipses (Lbol = 0.4–4.7 × 1045 erg s−1). The filled black circle indicates the position of the AGN, corresponding to the kinematic centre of the masers. The radio continuum (grey contours) and the maser emission (filled coloured circles) were extracted from Gallimore et al. (2004) and aligned based on the re-analysis of Gallimore & Impellizzeri (2019). The radio and IR images are aligned based on the continuum peaks. The green dashed lines indicate the bipolar outflow cone observed in IR polarisation on scales of few hundred mas (Gratadour et al. 2015). The orientation of the polarisation cone is very similar to the large scale ionisation cone found by Das et al. (2006). The grey dashed ellipse indicates the observed 10 μm emission size (taken from Raban et al. 2009). The elliptical beam, used for the reconstruction, is shown in the lower left corner.

Open with DEXTER

The image shown in Fig. 4 features an extended structure with a central hole, which can be described as a 3/4 ring with a radius of r ≈ 0.24 ± 0.03 pc (3.5 ± 0.3 mas), position angle PA = −50  ±  4° and an inclination of i = 70  ±  5°. These parameters are summarised in Table 2. The extent of the NIR structure resembles in size and orientation the observed 10 μm emission, yet at much higher resolution. The size of the inner bright ring-like structure matches remarkably well the expected dust sublimation distance of graphite grains in the radiation field of a central engine. The south-western ridge is about a factor two brighter than the north-eastern side of the ring-structure. The north-east orientation of the jet and the ionised outflow, which are approaching the observer (Gallimore et al. 2004), suggests that the south-western ridge corresponds to the near-side of the emission ring (but see Sect. 4.4 for an alternative interpretation). The three dimensional orientation is reinforced if the radio continuum emission (Gallimore et al. 2004) is aligned close to the NIR emission. We find a striking resemblance between the maser emission locations and the bright south western rim of the K-band image. Since maser emission is beamed towards the observer, the emitting clouds must lie between the source (AGN) and the observer, and hence the maser emission should originate from the near-side. This alignment, discussed in Sect. 3.3, supports the association of the ring-like structure with the dust sublimation region around the AGN, and it means that the near-side of the dusty ring is not only observable but also brighter than the far-side. Those findings are very difficult to reconcile with geometrically thick clumpy torus models, which can only reproduce a NIR ring-like structure in systems that are relatively face-on, and struggle to make the near side of the ring brighter than the far side (the near and far sides are given by other observations, notably Kishimoto 1999 and Gallimore et al. 2016, as indicated in Sect. 1).

Table 2.

Parameters of dust sublimation ring.

3.2. Disc height

The observed width of the ring structure is consistent with it being unresolved. We can place a firm upper limit on the width of the ring of < 0.5 mas or 0.035 pc. If this is interpreted in terms of a vertical thickness, then the scale height of the ring must be h/r <  0.14.

3.3. Cospatial dust and maser disc

Maser emission requires a high velocity coherence of the maser emitting gas as well as a population inversion of the masering molecules (Lo 2005). Energy and momentum conservation imply that the induced photon has the same frequency and direction as the stimulating photon (Elitzur 1990, 2002). There is a broad consensus that both the maser emission and the NIR emission are powered by the central AGN (e.g. Neufeld & Melnick 1991). The X-rays are able to heat the gas above the 400 K threshold where the water abundance is enhanced. To achieve a state in which atomic and molecular phases can co-exist, it is important that the medium is optically thin to the non-masing far-infrared photons, otherwise the trapping of these photons by the maser molecules would lead to a thermalisation of the population. Ultimately, the maser emission rate is limited by the loss or destruction rate of this accompanying far-infrared radiation (Rank et al. 1971; Neufeld et al. 1994). This constrains the geometry of the maser source because the surface area for cooling should be optimised, favouring an elongated or flat structure (Lo 2005).

The fact that mega-masers are almost exclusively found in Seyfert 2 galaxies (Ramolla et al. 2011) is consistent with the picture that disc maser emission is produced and beamed close to the plane of the obscuring material. In NGC 1068, nuclear H2O 22 GHz maser emission is detected, spanning a velocity range of 600 km s−1 nearly symmetric about the systemic velocity of the galaxy (Greenhill & Gwinn 1997; Gallimore et al. 2001). The maser spots extend out to ∼0.8 pc, aligned along a position angle of roughly −45°. While symmetric in velocity, the “red masers” at positive velocities are a factor ∼4–5 brighter than the “blue masers” (Gallimore et al. 2001). In particular H2O, 22 GHz emission requires hydrogen densities of ∼107–1011 cm−3 and dust temperatures above ∼400 K, but not too far in excess of 1000–1500 K (Neufeld 2000). The pump rate of different water emission lines (22, 183, 321 and 325 GHz) depends sensitively on the local density and temperature (Neufeld 2000). The detection of these higher transitions in the future (e.g. with ALMA) can therefore provide valuable information on the conditions close to the sublimation region.

Unfortunately the absolute astrometric accuracy of the infrared emission (discussed in Sect. 5.1) is not sufficient to align the radio and the infrared emission on the scales of interest for this paper, although the radio continuum and the maser emission can be aligned to ∼1.3 mas accuracy (Gallimore et al. 2004). Given the lack of precise positioning information, one option for the relative alignment would be, given their striking similarities, to match the bright side of the ring-like structure in the K-band continuum with the high density of masers spots. However, this does not take into account the implied location of the central black hole. If one does so, a more natural choice is to align the centre of the ring with the kinematic centre of the masers (Gallimore & Impellizzeri 2019). The overall morphological match between the maser spots, the radio-continuum, and the south-western ridge (side facing the observer) of the K-band continuum is remarkable (see Fig. 4). This suggests that the maser and the K-band emission originate from the same region, the dust sublimation ring. The fact that maser emission is only detected on the side facing the observer is consistent with the directional beaming of the maser emission (masing clouds are located between the central source and observer). We only find infrared counterparts of the systemic masers and red masers. The weaker blue masers show no infrared counterpart.

4. Infrared spectral energy distribution

Our K-band image of NGC 1068 shows the central emitting region at an unprecedented spatial resolution. This allows us to estimate the brightness temperature based on the observed NIR and MIR fluxes. The GRAVITY fringe-tracker samples the K-band with five spectral elements between 1.95 and 2.35 μm. The observed flux of NGC 1068 shows a very red colour across the K-band. Based on the observed slope we can constrain the effective temperature and the extinction. Figure 5 shows the observed K-band flux. We also included the published MIR (8–13 μm) SED for the central component (taken from Raban et al. 2009; López-Gonzaga et al. 2014).

thumbnail Fig. 5.

K-band flux density observed with GRAVITY is shown in filled circles. The MIR flux density (incl. uncertainty) taken from Raban et al. (2009) is shown in grey (central 20 mas component only) and in blue (total MIR flux observed with MIDI within few 10 pc). In black, we show three torus models from Stalevski et al. (2012) with different inclinations of i = 60° (dash-dotted), i = 70° (solid) and i = 80° (dashed). The two single temperature black body models are overplotted. Model 2 in orange (Teff = 720 K, r = 0.24 pc and dr = 0.045 pc.). Model 3 in green (Teff = 1500 K, AK = 5.5, r = 0.24 pc and dr = 0.026 pc). The bolometric luminosities of the different models are indicated in units of L45 = Lbol/1045 erg s−1.

Open with DEXTER

In the following subsections, we discuss several models that might explain the observed K-band emission. We first provide a short description of the key points of each model, where it can succeed, and where it fails. We focus on whether each of the models can explain the ring-like structure seen in the interferometric data, whether they match the slope of the NIR continuum, and whether they are consistent with other observational constraints (e.g. extinction, luminosity and covering factor). We conclude on which model we favour in Sect. 5, where we provide a heuristic view of the central regions.

Model 1 is a geometrically thick clumpy torus composed of dusty clouds. This class of models is inconsistent with several aspects of the GRAVITY data. Such models, when viewed as close to edge-on as the orientation in NGC 1068 requires, cannot reproduce a ring-like structure. In addition the NIR continuum slope is rather shallower than that observed.

Model 2 enables us to explain both the observed K-band and MIR emission with a single, rather cool dust structure of Teff ≈ 720 K. The required NIR luminosity is rather low, and hence so is the implied covering factor, which is naturally the case for a co-planar thin disc. The major difficulty is that this also implies too little foreground extinction.

Model 3 assumes hot dust of Teff ≈ 1500 K absorbed by significant foreground extinction. This is mainly motivated by the observed emission size, which suggests hot dust close to the sublimation region. This model requires that the bulk of the MIR continuum originates in separate structures (consistent with the model proposed by Vollmer et al. 2018). Its positive aspects are that it can match the observed steep red slope of the NIR continuum, and that the required foreground extinction is consistent with that expected from other observations. The main difficulty is that the high NIR luminosity requires a large covering factor. But this can be resolved if the accretion disc is tilted with respect to the maser disc as hinted at in point (4) of the introduction.

Model 4 is an alternative to Model 3, in which all of the bright NIR continuum observed by GRAVITY originates from the far side, rather than near-side, of the inner disc. Rather than recognising the K-band ring as a physical structure as in Model 3, this model assigns the NIR morphology to chance alignment due to clouds along the line of sight. In particular, it requires unrelated patchy obscuration to create an apparent ring-like structure with a position angle, aspect ratio, and size that match the geometry given by the maser disc and dust sublimation radius.

4.1. Model 1: geometrically thick clumpy torus

We have already seen in Sect. 3 that the ring-like structure, especially with the near side brighter than the far side, is hard to reconcile with a geometrically thick clumpy torus model that is close to edge-on. Similarly, these models require a careful matching of the vertical density structure to remain consistent with the presence of the thin maser disc at a radius ≲1 pc.

Here we consider the photometric aspect. We compare our infrared SED to the torus models of Stalevski et al. (2012). These models use a two-phase medium of high density clumps and low-density gas filling the space in between. The model assumes standard silicate and graphite ISM-type grains and does not include a pure graphite grain component at small distances. The models have been calculated for a primary source of 1011L (Lbol = 0.4 × 1045 erg s−1) with inner and outer radii of 0.5 and 15 pc, respectively. We scaled the models to the distance of NGC 1068 and to the observed K-band flux. In Fig. 5, models with three different inclinations are compared to the observed SED. We find the best matching models to be at inclinations between 60° and 70°. Models with higher inclinations (e.g. 80°) tend to over-predict the MIR emission. Models with lower inclinations, on the other hand, under-predict the MIR emission. These scaled models would imply a bolometric luminosity between Lbol = 0.4–1.0 × 1045 erg s−1 for i = 60°–70°, in good agreement with earlier estimates of the bolometric luminosity (see Sect. 2.1). However, they do not match the red slope of the NIR continuum (with the caveat that the model is missing a pure graphite-grain component).

This model has difficulties reproducing the ring-like dust structure, accounting for the maser disc, and a NIR continuum slope that is too shallow. We therefore rule out a geometrically thick clumpy torus model for NGC 1068.

4.2. Model 2: cool dust disc

MIR interferometry (Jaffe et al. 2004; Raban et al. 2009; López-Gonzaga et al. 2014) was used to model the central component with 600–800 K dust. In the cool dust model, we try to explain the observed K-band and MIR emission with a single source and temperature. The observed SED can be well approximated by a black body with an effective temperature of Teff = 720 ± 30 K seen through a screen of extinction of AK ≈ 0.9 ± 0.3 (see Fig. 5). It is somewhat surprising that it is possible to explain the 2–13 μm SED with a single absorbed black body. In particular, the strong silicate absorption dip at 10 μm depends sensitively on the amount of dust and the dust composition. The very good agreement of the model and the observed silicate dip would indicate that the dust composition in NGC 1068 is very similar to the dust in our Galactic Centre.

The minimum emission surface would be at least ∼0.068 pc2, which would correspond to a ring with r = 0.24 pc and a thickness of δr = 0.045 pc (0.65 mas). Such a large thickness would be partially resolved and is therefore barely consistent with the data. A single black body ring would produce 1 percent of the estimated bolometric luminosity (L = σAT4 = 1 × 1043 erg s−1). This covering factor is consistent with a thin dusty disc having h/r <  0.14, that is co-planar with the accretion disc even under the assumption of anisotropic emission.

The cool dust model would allow us to explain the near-to-mid infrared SED with a surprisingly simple model. However it predicts too little extinction, which is inconsistent with the Brα constraint from ISO (ABrα >  2.4, Lutz et al. 2000) as well as the millimetre constraints on the column density (Tacconi et al. 1994; Sternberg et al. 1994; García-Burillo et al. 2017; Imanishi et al. 2018). The low extinction value results in a remarkably low bolometric luminosity of the dust. Moreover, it is hard to explain a lack of a dust emission component with a much higher temperature given the distance of the central ring from the very luminous nucleus. We therefore disfavour the cool dust and low extinction model.

4.3. Model 3: hot dust disc

The size of the observed structure in Fig. 3 suggests a temperature close to the dust sublimation temperature of Teff = 1500 K–1800 K. In order to reproduce the K-band slope, we therefore assumed a black body of T = 1500 K, which is highly absorbed (by foreground dust) with a mean K-band extinction of AK = 5.5 ± 1.7. Figure 5 shows the observed fluxes and the proposed model. The corresponding optical extinction is AV ≈ 90 ± 30. This is fully consistent with the extinction to the BLR derived from limits on the Brα broad line emission by Lutz et al. (2000), as well as other estimates of the column density towards the nucleus noted in Sect. 2.2.

In order to explain the observed K-band flux with the hot dust model, the emission surface needs to be at least ∼0.042 pc2. If the emission is not a black body but a grey body then the emission region can be correspondingly larger. The minimum required surface is equivalent to a ring with r = 0.24 pc and a width of δr = 0.026 pc (0.38 mas). This means that the observed ring-like structure is sufficient to explain the observed K-band flux. The width of the ring would be unresolved, which is also consistent with the observations.

The hot dust model fails to produce the observed MIR emission. This is, in principle, not a problem, but it means that the bulk of the 8–13 μm emission must originate from significantly cooler dust at larger radii, which is directly exposed to the central continuum and could be related to the dusty outflow (Cameron et al. 1993; Bock et al. 2000; Vollmer et al. 2018).

The model also requires about 10 percent of the estimated AGN bolometric luminosity to be reprocessed by dust and re-emitted at NIR wavelengths (L = σAT4 = 1.1 × 1044 erg s−1). This would be inconsistent with a system in which the accretion disc and thin dusty maser disc are co-planar due to the small covering factor. However, there are clear indications that these structures are misaligned: the relative position angle of the radio jet and maser disc (see Fig. 8 of Gallimore et al. 2004 for an illustration), as well as the inclination of the maser disc with respect to the orientation of the outflowing clouds (Kishimoto 1999). Adopting anisotropic emission from the accretion disc as described by Netzer (2015)L ∝ cos(θ)[1 + b cos(θ)] with b ≃ 2, we find that about 10% of the AGN luminosity is intersected by a surrounding disc with h/r <  0.14 (see Sect. 3.2) when the accretion disc is tilted by 40°. This tilt is consistent with position angle and inclination discrepancies above. A weaker anisotropy, for example L ∝ cos(θ), requires a somewhat smaller tilt. This means that the “covering factor” issue, which is often taken to imply a geometrically thick structure can, for the NIR continuum, be entirely accounted for by a warped inner disc at scales ≲0.1 pc.

In the scenario of a thin disc screened by foreground extinction, the visibility of the “central engine” (accretion disc and BLR) inside the sublimation hole is an obvious question. The K-band contribution of the accretion disc can be estimated based on the bolometric luminosity Lbol ≈ 1 × 1045 erg s−1, the distance 14.4 Mpc, a λ5100Å bolometric correction factor ∼9 (Netzer 2015) and a SED slope of the accretion disc of Fν ∝ να with α = −0.44 (Vanden Berk et al. 2001). Also considering the AK ∼ 5.5 foreground extinction, we infer a K-band accretion disc flux of ∼9 mJy. The empirical slope of Vanden Berk et al. (2001) is significantly redder than the canonical Shakura & Sunyaev (1973) slope of α = 0.33. The K-band accretion disc emission for the canonical slope would be much weaker, and so the inferred contribution is probably a conservative upper limit. Other evidence for yet bluer slopes as found by Kishimoto et al. (2007, 2008) would further enhance that argument. The estimated accretion disc contribution of ∼9 mJy corresponds to 6% or less of the flux in the reconstructed image (Fig. 4). This means that a point-like accretion disc could easily be hidden in the image including the centre of the dust ring.

Both aspects of (1) accretion disc and (2) BLR visibility are consistent with the hot dust scenario. However, in both cases expected signatures are not far from the observational results/limits.

4.4. Model 4: optically thick hot dusty disc

Model 3 above places the AGN at the centre of the K-band hot dust ring that has the expected position angle and aspect ratio for the inner geometry of NGC 1068 and a size comparable to the dust sublimation radius, and identifies this position with the maser kinematic centre. It implies the NIR continuum comes from both sides of the central disc, with the brightest on the near side while the far side is fainter. This implicitly requires that the re-radiating clouds are optically thin to the hot dust emission. Applying a standard gas-to-dust ratio as described in Sect. 2.2 and the extinction law from Sect. 2.3 mean that τ2.2 μm = 1 occurs at a column density of 4 × 1022 cm−2. This is the maximum column density allowed for a typical cloud in the dust disc for those models. In order for the far side of the disc to be fainter, there must be some additional obscuration to that line of sight. The modest factor 4 difference in intensity is plausibly equivalent to 1–2 such clouds. Confirming this alignment requires precise astrometry, which is addressed in Sect. 5.1.

An alternative alignment would be implied if the column density through individual clouds were higher, NH = 1023 cm−2 or more, as inferred from physical considerations of cloud stability in the vicinity of an AGN (Krolik & Begelman 1988; Vollmer et al. 2004; Hönig & Beckert 2007), so that they are optically thick at NIR wavelengths. In this scenario, the NIR radiation is anisotropic and primarily emerges from the hot, AGN-facing side of each cloud. As a consequence, a ring consisting of such clouds would be brighter on the far side than the near side. We would be exposed to the directly heated surface of clouds on the far side, while on the near side we would mostly observe the clouds’ cool back sides. The hot dust disc is still co-incident with the maser dust, but these phenomena trace different phases of the dusty gas within the disc, so that the alignment shown in Fig. 4 is lost. Models for NGC 1068 constructed in this way could be similar to the disc-wind scenario proposed by Hönig & Kishimoto (2017) for NGC 3783, which, in contrast to conventional geometrically thick torus models, indicate that the bulk of the MIR continuum arises in a wind. Such a model would posit that the NIR continuum comes from the τ2.2 μm = 1 surface of the hot disc that is still geometrically thin, although offset from the equatorial plane according to h/r, which could be in the range 0.1-0.2; and that at longer wavelengths pertinent to the masers, the disc is optically thin. In this model, the observed NIR continuum would still extend inwards to the dust sublimation radius, but the brightest spot observed in the K-band by GRAVITY would coincide with the far side of the disc. The ring morphology with radius matching the expected dust sublimation radius, the similarity between the maser disc and the K-band emission, as well as the alignment of the maser and radio continuum with the K-band would all be chance associations. The right panel of Fig. 6 shows the alignment for the optically thick hot dust disc model. The alignment of the GRAVITY data to the disc mid-plane has been implied by radiative transfer simulations using the CAT3D-WIND model (Hönig & Kishimoto 2017) with an assumed inclination of i = 80deg, a scale height of the disc of h/r = 0.1, and the optical depth of an individual cloud of τV = 50.

thumbnail Fig. 6.

Comparison of the alignment of the hot dusty disc models. The models themselves are not shown, only the implied interpretation of the observed continuum features. Left: Model 3 is the favoured interpretation, in which the ring-like structure in the K-band continuum arises from (partially) optically thin clouds (same as in Fig.4). In this case, it traces a disc that has a similar position angle and aspect ratio to the inner geometry of NGC 1068, is co-planar with the maser emission and traces the dust sublimation ring. Right: Model 4 comprises optically thick clouds, and the alignment of the mid-plane (indicated by the large dashed ring) to the NIR continuum observed with GRAVITY is given by radiative transfer simulations using the CAT3D-WIND model (Hönig & Kishimoto 2017). In Model 4, the near side of the mid-plane as traced by the masers, but is completely obscured in the K-band. The K-band emission originates from the τ = 1 surface of the hot disc, that is still geometrically thin. But all the structure in the image is assigned to chance effects of patchy obscuration along the line of sight.

Open with DEXTER

For this model, the resulting structures originate with the random distribution of clouds, and hence are due to variations in obscuration along the line of sight. As such, each instance of this model that is created will look different. This makes the model seem arbitrary because one needs to create many such instances in order to pick one that approximates the characteristics of the observed continuum distribution.

5. Discussion

5.1. Astrometry

The required accuracy to make a decisive statement on the astrometry of the K-band continuum with respect to the masers or radio continuum is on the order ≤2–3 mas. Without reference features, which can be identified simultaneously in the K-band and in the radio, absolute astrometry is required to relate the radio and infrared measurements. The radio position is sufficiently well known with an uncertainty σ <  1.5 mas (Gallimore et al. 2004). Unfortunately the astrometric uncertainty of the infrared data is a factor 10 worse (σ ∼ 20 mas). Future interferometric observations should be able to exploit nearby stars from the Gaia catalogue as astrometric references. Then, by constructing a local baseline solution using these reference stars and by logging the optical path difference between the stars and NGC 1068, a sufficiently accurate astrometric solution might be achieved.

Without such an astrometric reference yet available, we favour Model 3 which provides a natural explanation for several aspects of the observed ring-like structure:(i) the match of its position angle and aspect ratio to those of the maser disc; (ii) the match of its radius to the dust sublimation region around the AGN; and (iii) the alignment between its near-side and the masers, when the centre of the ring is aligned to the black hole location derived from the maser kinematics (see Sect. 3.3). The hot dust disc, with clouds that are optically thin in the NIR, is consistent with observations and provides a rationale for the similarity of the masers, radio continuum and the hot dust K-band emission. Based on that model we propose a simple disc-like geometry, which is obscured by foreground (within a few parsec) molecular gas.

5.2. A heuristic model for the central region of NGC 1068

With a wealth of observations across the electromagnetic spectrum at hand, we try to construct a simple heuristic picture of the central few parsecs of NGC 1068 (see Fig. 7). The K-band emission traces the hot dust sublimation radius of a thin disc, seen at an inclination of 70°. The inner masers are cospatial with the K-band emission. The masers extend out to radii of ∼1 pc (Gallimore et al. 2001), similar to the size of the warm dust structure (∼600 K) seen at MIR wavelengths (Jaffe et al. 2004; Raban et al. 2009; López-Gonzaga et al. 2016). This suggests that the warm dust originates from the outer region of the same thin disc. Beyond ∼1 pc distance no masers have been observed, probably because the conditions disfavour masing (e.g. too low dust temperature, < 400 K). In contrast to the thin disc at scales of < 1 pc, ALMA observations by Imanishi et al. (2018) reveal a significant amount of molecular gas, which is rotating with vc ≈ 40–50 km s−1 (notably in the opposite sense as the maser disc) at 5–10 pc distance. The comparably large velocity dispersion of σ ≈ 30–40 km s−1 and the correspondingly high σ/vc ≃ 1 suggests a thick disc with a large scale height of hz/R ≃ 1. The disc contains enough gas mass to reach column densities of NH ≈ 1023..24 cm−2 that screen the central region by AV ≈ 90 (AK ≈ 5.5 ± 1.7). The central few parsecs are fed by two gas streamers (Müller Sánchez et al. 2009), which come very close to the nucleus. The northern and southern streamers reach pericentre distances of ≈5 pc and ≈1 pc repectively. Their opposite sense of rotation could contribute to the turbulence observed in the molecular disc.

thumbnail Fig. 7.

Sketch of the observed central structures. The K-band emission traces the inner rim of a thin disc of hot gas and dust, at or close to the dust sublimation radius r ≈ 0.24 pc. The inner water masers are cospatial with the hot K-band dust. The masers stretch out to r ≈ 1 pc (Gallimore et al. 2001). MIR observations show warm dust on similar scales as the outer masers (e.g. Raban et al. 2009), likely originating from the disc periphery. ALMA observations of HCN and HCO+ show a turbulent structure, which rotates in the opposite direction as the maser disc (Imanishi et al. 2018). The vc/σ of the molecular gas structure argues for a thick disc, which contains enough gas mass to reach column densities of NH ≈ 1023 cm−2 that screen the central region by AV ≈ 90 (AK ≈ 5.5) from the observer.

Open with DEXTER

6. Conclusion

We present new 0.2 pc resolution observations from the GRAVITY interferometer on the VLT, which spatially resolve the hot dust continuum in the central parsec of NGC 1068. The following conclusions assume a distance of 14.4 Mpc. We find that:

  • The dust structure is dominated by a thin, ring-like structure with a radius of r = 0.24 ± 0.03 pc, a scale height h/r <  0.14, an inclination i = 70° ±5° and a position angle of PA = −50° ±4°.

  • The spatial scale of the ring agrees with the predicted sublimation distance of graphite dust illuminated by a Lbol ≈ 0.4–4.7 × 1045 erg s−1 AGN (radiating at or in excess of its Eddington limit). We therefore associate this with the dust sublimation region.

  • The NIR continuum structures show a striking resemblance to the maser disc. We suggest that both types of emission originate from the same region at ∼0.4 pc distance from the central source. The maser disc is co-spatial with the south-western rim of the dust sublimation region, the side closer to the observer (as expected due to the maser beaming).

  • Based on the observed ring-like structure and steep NIR continuum slope, we disfavour geometrically thick clumpy torus models (which also do not account for the maser disc).

  • The dust structure and photometry are consistent with a simple model of hot dust at T ∼ 1500 K that is behind AK ∼ 5.5 (AV ≈ 90) mag of foreground extinction. This amount of screen extinction is expected from an upper limit to the Brα broad line and could be provided by the dense and turbulent gas distribution observed on scales of 1–10 pc. The model implies that the NIR and MIR emitting dust are in spatially separate locations (the maser disc and the outflow respectively).


1

The NASA/IPAC Extragalactic Database (NED) is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

2

ESO Telescopes at the La Silla Paranal Observatory programme IDs 0102.B-0667, 0102.C-0205 and 0102.C-0211.

Acknowledgments

We thank the referees for their careful reading of the manuscript and their suggestions that have helped to improve it. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This work has made extensive use of the software QFitsView, which has been developed by Thomas Ott, Max Planck Institute for Extraterrestrial Physics (http://www.mpe.mpg.de/~ott/QFitsView/). F. E. and O. P. acknowledge support from ERC synergy grant No. 610058 (Black-HoleCam). J. D. was supported by a Sofja Kovalevskaja award from the Alexander von Humboldt foundation and in part by NSF grant AST 1909711. A. A. and P. G. acknowledge funding from Fundação para a Ciência e a Tecnologia through grants UID/FIS/00099/2013 and SFRH/BSAB/142940/2018. P.O.P acknowledges support from the Programme National Hautes Energies (PNHE) of the Centre National de la Recherche Scientifique (CNRS).

References

  1. Antonucci, R. 1993, ARA&A, 31, 473 [NASA ADS] [CrossRef] [Google Scholar]
  2. Antonucci, R. R. J. 1982, Nature, 299, 605 [NASA ADS] [CrossRef] [Google Scholar]
  3. Antonucci, R. R. J. 1984, ApJ, 278, 499 [NASA ADS] [CrossRef] [Google Scholar]
  4. Antonucci, R. R. J., & Miller, J. S. 1985, ApJ, 297, 621 [NASA ADS] [CrossRef] [Google Scholar]
  5. Anugu, N., Amorim, A., Gordo, P., et al. 2018, MNRAS, 476, 459 [NASA ADS] [CrossRef] [Google Scholar]
  6. Asmus, D., Gandhi, P., Hönig, S. F., Smette, A., & Duschl, W. J. 2015, MNRAS, 454, 766 [NASA ADS] [CrossRef] [Google Scholar]
  7. Barvainis, R. 1987, ApJ, 320, 537 [NASA ADS] [CrossRef] [Google Scholar]
  8. Baskin, A., & Laor, A. 2018, MNRAS, 474, 1970 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bauer, F. E., Arévalo, P., Walton, D. J., et al. 2015, ApJ, 812, 116 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bland-Hawthorn, J., Gallimore, J. F., Tacconi, L. J., et al. 1997, Ap&SS, 248, 9 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bock, J. J., Neugebauer, G., Matthews, K., et al. 2000, AJ, 120, 2904 [NASA ADS] [CrossRef] [Google Scholar]
  12. Cameron, M., Storey, J. W. V., Rotaciuc, V., et al. 1993, ApJ, 419, 136 [NASA ADS] [CrossRef] [Google Scholar]
  13. Das, V., Crenshaw, D. M., Kraemer, S. B., & Deo, R. P. 2006, AJ, 132, 620 [NASA ADS] [CrossRef] [Google Scholar]
  14. Davies, R. I., Burtscher, L., Rosario, D., et al. 2015, ApJ, 806, 127 [NASA ADS] [CrossRef] [Google Scholar]
  15. Elitzur, M. 1990, ApJ, 363, 628 [NASA ADS] [CrossRef] [Google Scholar]
  16. Elitzur, M. 2002, in Cosmic Masers: From Proto-Stars to Black Holes, eds. V. Migenes, & M. J. Reid, IAU Symp., 206, 452 [NASA ADS] [Google Scholar]
  17. Fritz, T. K., Gillessen, S., Dodds-Eden, K., et al. 2011, ApJ, 737, 73 [NASA ADS] [CrossRef] [Google Scholar]
  18. Gallimore, J. F., & Impellizzeri, C. M. V. 2019, ApJ, submitted [Google Scholar]
  19. Gallimore, J. F., Baum, S. A., & O’Dea, C. P. 1996a, ApJ, 464, 198 [NASA ADS] [CrossRef] [Google Scholar]
  20. Gallimore, J. F., Baum, S. A., O’Dea, C. P., Brinks, E., & Pedlar, A. 1996b, ApJ, 462, 740 [NASA ADS] [CrossRef] [Google Scholar]
  21. Gallimore, J. F., Henkel, C., Baum, S. A., et al. 2001, ApJ, 556, 694 [NASA ADS] [CrossRef] [Google Scholar]
  22. Gallimore, J. F., Baum, S. A., & O’Dea, C. P. 2004, ApJ, 613, 794 [NASA ADS] [CrossRef] [Google Scholar]
  23. Gallimore, J. F., Elitzur, M., Maiolino, R., et al. 2016, ApJ, 829, L7 [NASA ADS] [CrossRef] [Google Scholar]
  24. García-Burillo, S., Combes, F., Ramos Almeida, C., et al. 2016, ApJ, 823, L12 [NASA ADS] [CrossRef] [Google Scholar]
  25. García-Burillo, S., Viti, S., Combes, F., et al. 2017, A&A, 608, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Gratadour, D., Rouan, D., Grosset, L., Boccaletti, A., & Clénet, Y. 2015, A&A, 581, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. GRAVITY Collaboration (Abuter, R., et al.) 2017, A&A, 602, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. GRAVITY Collaboration (Dexter, J., et al.) 2019, A&A, submitted, ArXiv e-prints [arXiv:1910.00593] [Google Scholar]
  29. Greenhill, L. J., & Gwinn, C. R. 1997, AP&SS, 248, 261 [NASA ADS] [CrossRef] [Google Scholar]
  30. Greenhill, L. J., Gwinn, C. R., Antonucci, R., & Barvainis, R. 1996, ApJ, 472, L21 [NASA ADS] [CrossRef] [Google Scholar]
  31. Güver, T., & Özel, F. 2009, MNRAS, 400, 2050 [NASA ADS] [CrossRef] [Google Scholar]
  32. Hönig, S. F. 2019, ApJ, 884, 171 [NASA ADS] [CrossRef] [Google Scholar]
  33. Hönig, S. F., & Beckert, T. 2007, MNRAS, 380, 1172 [NASA ADS] [CrossRef] [Google Scholar]
  34. Hönig, S. F., & Kishimoto, M. 2017, ApJ, 838, L20 [NASA ADS] [CrossRef] [Google Scholar]
  35. Hönig, S. F., Beckert, T., Ohnaka, K., & Weigelt, G. 2006, A&A, 452, 459 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Hönig, S. F., Beckert, T., Ohnaka, K., & Weigelt, G. 2007, in The Central Engine of Active Galactic Nuclei, ed. L. C. Ho, ASP Conf. Ser., 373, 487 [NASA ADS] [Google Scholar]
  37. Hönig, S. F., Prieto, M. A., & Beckert, T. 2008, A&A, 485, 33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Huré, J.-M. 2002, A&A, 395, L21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Imanishi, M., Nakanishi, K., Izumi, T., & Wada, K. 2018, ApJ, 853, L25 [NASA ADS] [CrossRef] [Google Scholar]
  40. Impellizzeri, C. M. V., Gallimore, J. F., Baum, S. A., et al. 2019, ApJ, 884, L28 [NASA ADS] [CrossRef] [Google Scholar]
  41. Jaffe, W., Meisenheimer, K., Röttgering, H. J. A., et al. 2004, Nature, 429, 47 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  42. Kishimoto, M. 1999, ApJ, 518, 676 [NASA ADS] [CrossRef] [Google Scholar]
  43. Kishimoto, M., Hönig, S. F., Beckert, T., & Weigelt, G. 2007, A&A, 476, 713 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Kishimoto, M., Antonucci, R., Blaes, O., et al. 2008, Nature, 454, 492 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  45. Kishimoto, M., Hönig, S. F., Antonucci, R., et al. 2011, A&A, 527, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Kishimoto, M., Hönig, S. F., Antonucci, R., et al. 2013, ApJ, 775, L36 [NASA ADS] [CrossRef] [Google Scholar]
  47. Koshida, S., Minezaki, T., Yoshii, Y., et al. 2014, ApJ, 788, 159 [NASA ADS] [CrossRef] [Google Scholar]
  48. Krolik, J. H., & Begelman, M. C. 1988, ApJ, 329, 702 [NASA ADS] [CrossRef] [Google Scholar]
  49. Lacour, S., Dembet, R., Abuter, R., et al. 2019, A&A, 624, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Lapeyrere, V., Kervella, P., & Lacour, S. 2014, in Optical and Infrared Interferometry IV, Proc. SPIE, 9146, 91462D [Google Scholar]
  51. Lo, K. 2005, Annu. Rev. Astron. Astrophys., 43, 625 [NASA ADS] [CrossRef] [Google Scholar]
  52. Lodato, G., & Bertin, G. 2003, A&A, 398, 517 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. López-Gonzaga, N., Jaffe, W., Burtscher, L., Tristram, K. R. W., & Meisenheimer, K. 2014, A&A, 565, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. López-Gonzaga, N., Burtscher, L., Tristram, K. R. W., Meisenheimer, K., & Schartmann, M. 2016, A&A, 591, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Lopez-Rodriguez, E., Fuller, L., Alonso-Herrero, A., et al. 2018, ApJ, 859, 99 [NASA ADS] [CrossRef] [Google Scholar]
  56. Lutz, D., Genzel, R., Sturm, E., et al. 2000, ApJ, 530, 733 [NASA ADS] [CrossRef] [Google Scholar]
  57. Marconi, A., Risaliti, G., Gilli, R., et al. 2004, MNRAS, 351, 169 [NASA ADS] [CrossRef] [Google Scholar]
  58. Marinucci, A., Bianchi, S., Matt, G., et al. 2016, MNRAS, 456, L94 [NASA ADS] [CrossRef] [Google Scholar]
  59. Matt, G., Guainazzi, M., Frontera, F., et al. 1997, A&A, 325, L13 [NASA ADS] [Google Scholar]
  60. Miller, J. S., & Antonucci, R. R. J. 1983, ApJ, 271, L7 [NASA ADS] [CrossRef] [Google Scholar]
  61. Miller, J. S., Goodrich, R. W., & Mathews, W. G. 1991, ApJ, 378, 47 [NASA ADS] [CrossRef] [Google Scholar]
  62. Mor, R., & Netzer, H. 2012, MNRAS, 420, 526 [NASA ADS] [CrossRef] [Google Scholar]
  63. Mould, J. R., Huchra, J. P., Freedman, W. L., et al. 2000, ApJ, 529, 786 [NASA ADS] [CrossRef] [Google Scholar]
  64. Müller Sánchez, F., Davies, R. I., Genzel, R., et al. 2009, ApJ, 691, 749 [NASA ADS] [CrossRef] [Google Scholar]
  65. Nenkova, M., Sirocky, M. M., Nikutta, R., Ivezić, Ž., & Elitzur, M. 2008, ApJ, 685, 160 [NASA ADS] [CrossRef] [Google Scholar]
  66. Netzer, H. 2009, MNRAS, 399, 1907 [NASA ADS] [CrossRef] [Google Scholar]
  67. Netzer, H. 2015, ARA&A, 53, 365 [NASA ADS] [CrossRef] [Google Scholar]
  68. Neufeld, D. A. 2000, ApJ, 542, L99 [NASA ADS] [CrossRef] [Google Scholar]
  69. Neufeld, D. A., & Melnick, G. J. 1991, ApJ, 368, 215 [NASA ADS] [CrossRef] [Google Scholar]
  70. Neufeld, D. A., Maloney, P. R., & Conger, S. 1994, ApJ, 436, L127 [NASA ADS] [CrossRef] [Google Scholar]
  71. Packham, C., Young, S., Hough, J. H., Axon, D. J., & Bailey, J. A. 1997, MNRAS, 288, 375 [NASA ADS] [CrossRef] [Google Scholar]
  72. Pfuhl, O., Haug, M., Eisenhauer, F., et al. 2012, in Optical and Infrared Interferometry III, Proc. SPIE, 8445, 84451U [CrossRef] [Google Scholar]
  73. Pfuhl, O., Haug, M., Eisenhauer, F., et al. 2014, in Optical and Infrared Interferometry IV, Proc. SPIE, 9146, 914623 [CrossRef] [Google Scholar]
  74. Pier, E. A., & Krolik, J. H. 1992, ApJ, 401, 99 [NASA ADS] [CrossRef] [Google Scholar]
  75. Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Planesas, P., Scoville, N., & Myers, S. T. 1991, ApJ, 369, 364 [NASA ADS] [CrossRef] [Google Scholar]
  77. Predehl, P., & Schmitt, J. H. M. M. 1995, A&A, 293, 889 [NASA ADS] [Google Scholar]
  78. Prieto, M. A., Reunanen, J., Tristram, K. R. W., et al. 2010, MNRAS, 402, 724 [NASA ADS] [CrossRef] [Google Scholar]
  79. Raban, D., Jaffe, W., Röttgering, H., Meisenheimer, K., & Tristram, K. R. W. 2009, MNRAS, 394, 1325 [NASA ADS] [CrossRef] [Google Scholar]
  80. Ramolla, M., Haas, M., Bennert, V. N., & Chini, R. 2011, A&A, 530, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Rank, D. M., Townes, C. H., & Welch, W. J. 1971, Science, 174, 1083 [NASA ADS] [CrossRef] [Google Scholar]
  82. Ricci, C., Trakhtenbrot, B., Koss, M. J., et al. 2017, ApJS, 233, 17 [NASA ADS] [CrossRef] [Google Scholar]
  83. Risaliti, G., Elvis, M., & Nicastro, F. 2002, ApJ, 571, 234 [NASA ADS] [CrossRef] [Google Scholar]
  84. Rouan, D., Grosset, L., & Gratadour, D. 2019, A&A, 625, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Schnülle, K., Pott, J. U., Rix, H. W., et al. 2015, A&A, 578, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  87. Stalevski, M., Fritz, J., Baes, M., Nakos, T., & Popović, L. Č. 2012, MNRAS, 420, 2756 [NASA ADS] [CrossRef] [Google Scholar]
  88. Sternberg, A., Genzel, R., & Tacconi, L. 1994, ApJ, 436, L131 [NASA ADS] [CrossRef] [Google Scholar]
  89. Storchi-Bergmann, T., Kinney, A. L., & Challis, P. 1995, ApJS, 98, 103 [NASA ADS] [CrossRef] [Google Scholar]
  90. Suganuma, M., Yoshii, Y., Kobayashi, Y., et al. 2006, ApJ, 639, 46 [NASA ADS] [CrossRef] [Google Scholar]
  91. Tacconi, L. J., Genzel, R., Blietz, M., et al. 1994, ApJ, 426, L77 [NASA ADS] [CrossRef] [Google Scholar]
  92. Thiébaut, E. 2008, in Optical and Infrared Interferometry, Proc. SPIE, 7013, 70131I [NASA ADS] [CrossRef] [Google Scholar]
  93. Thiébaut, É., & Young, J. 2017, J. Opt. Soc. Am. A, 34, 904 [NASA ADS] [CrossRef] [Google Scholar]
  94. Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001, AJ, 122, 549 [NASA ADS] [CrossRef] [Google Scholar]
  95. Vollmer, B., Beckert, T., & Duschl, W. J. 2004, A&A, 413, 949 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Vollmer, B., Schartmann, M., Burtscher, L., et al. 2018, A&A, 615, A164 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Weigelt, G., Wittkowski, M., Balega, Y. Y., et al. 2004, A&A, 425, 77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Winter, L. M., Veilleux, S., McKernan, B., & Kallman, T. R. 2012, ApJ, 745, 107 [NASA ADS] [CrossRef] [Google Scholar]
  99. Wittkowski, M., Kervella, P., Arsenault, R., et al. 2004, A&A, 418, L39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Tables

Table 1.

Basic parameters of NGC 1068 taken from the literature.

Table 2.

Parameters of dust sublimation ring.

All Figures

thumbnail Fig. 1.

Squared visibilities measured for NGC 1068. The light coloured filled data points were obtained with the AT array in compact configuration (spatial frequencies 5 <  f <  17 Mλ.). The high spatial frequency data were taken with the UT array. The small inset shows the corresponding UV coordinates. The open blue and red symbols (spatial frequencies f <  4 Mλ) represent speckle data from Weigelt et al. (2004), which are shown only for comparison; they were not used for the image reconstruction. The data point and error bar at a spatial frequency of 21 Mλ, lying above our new data, shows the only previous NIR VLTI detection of NGC 1068 (Wittkowski et al. 2004). For comparison, the grey continuous line shows the visibility of a thin ring-shaped emission region with a radius of 3.3 mas and a width of 0.5 mas scaled to a maximum visibility of 0.3. While this reflects the general shape of the visibility distribution, it is a poor fit and indicates that more complex structures are present.

Open with DEXTER
In the text
thumbnail Fig. 2.

Closure phase data of NGC 1068. See Fig. 1 for the explanations of the symbols. The highly non-zero closure phases indicate strong source asymmetries on all spatial scales.

Open with DEXTER
In the text
thumbnail Fig. 3.

Reconstructed image of the central 5.1 × 5.1 pc of NGC 1068. The elliptical beam, which was used to convolve the image, is shown in the lower left. Details of the reconstruction, which only allows positive sources, are given in Sect. 3.1. The central structures are robust to changing the reconstruction parameters. The fainter features towards the edge of the field, with fluxes close to the noise floor of ∼7 mJy beam−1 rms (i.e. dark cyan in this colour scale), are uncertain.

Open with DEXTER
In the text
thumbnail Fig. 4.

Reconstructed image of the inner 2.1 × 2.1 pc of NGC 1068 (i.e. the central region of Fig. 3). The best fitting circular ring is shown in dashed white (r = 0.24 pc (3.5 mas). The dust sublimation sizes for the range of bolometric luminosities (∼0.08–0.27 pc) are shown in orange dotted ellipses (Lbol = 0.4–4.7 × 1045 erg s−1). The filled black circle indicates the position of the AGN, corresponding to the kinematic centre of the masers. The radio continuum (grey contours) and the maser emission (filled coloured circles) were extracted from Gallimore et al. (2004) and aligned based on the re-analysis of Gallimore & Impellizzeri (2019). The radio and IR images are aligned based on the continuum peaks. The green dashed lines indicate the bipolar outflow cone observed in IR polarisation on scales of few hundred mas (Gratadour et al. 2015). The orientation of the polarisation cone is very similar to the large scale ionisation cone found by Das et al. (2006). The grey dashed ellipse indicates the observed 10 μm emission size (taken from Raban et al. 2009). The elliptical beam, used for the reconstruction, is shown in the lower left corner.

Open with DEXTER
In the text
thumbnail Fig. 5.

K-band flux density observed with GRAVITY is shown in filled circles. The MIR flux density (incl. uncertainty) taken from Raban et al. (2009) is shown in grey (central 20 mas component only) and in blue (total MIR flux observed with MIDI within few 10 pc). In black, we show three torus models from Stalevski et al. (2012) with different inclinations of i = 60° (dash-dotted), i = 70° (solid) and i = 80° (dashed). The two single temperature black body models are overplotted. Model 2 in orange (Teff = 720 K, r = 0.24 pc and dr = 0.045 pc.). Model 3 in green (Teff = 1500 K, AK = 5.5, r = 0.24 pc and dr = 0.026 pc). The bolometric luminosities of the different models are indicated in units of L45 = Lbol/1045 erg s−1.

Open with DEXTER
In the text
thumbnail Fig. 6.

Comparison of the alignment of the hot dusty disc models. The models themselves are not shown, only the implied interpretation of the observed continuum features. Left: Model 3 is the favoured interpretation, in which the ring-like structure in the K-band continuum arises from (partially) optically thin clouds (same as in Fig.4). In this case, it traces a disc that has a similar position angle and aspect ratio to the inner geometry of NGC 1068, is co-planar with the maser emission and traces the dust sublimation ring. Right: Model 4 comprises optically thick clouds, and the alignment of the mid-plane (indicated by the large dashed ring) to the NIR continuum observed with GRAVITY is given by radiative transfer simulations using the CAT3D-WIND model (Hönig & Kishimoto 2017). In Model 4, the near side of the mid-plane as traced by the masers, but is completely obscured in the K-band. The K-band emission originates from the τ = 1 surface of the hot disc, that is still geometrically thin. But all the structure in the image is assigned to chance effects of patchy obscuration along the line of sight.

Open with DEXTER
In the text
thumbnail Fig. 7.

Sketch of the observed central structures. The K-band emission traces the inner rim of a thin disc of hot gas and dust, at or close to the dust sublimation radius r ≈ 0.24 pc. The inner water masers are cospatial with the hot K-band dust. The masers stretch out to r ≈ 1 pc (Gallimore et al. 2001). MIR observations show warm dust on similar scales as the outer masers (e.g. Raban et al. 2009), likely originating from the disc periphery. ALMA observations of HCN and HCO+ show a turbulent structure, which rotates in the opposite direction as the maser disc (Imanishi et al. 2018). The vc/σ of the molecular gas structure argues for a thick disc, which contains enough gas mass to reach column densities of NH ≈ 1023 cm−2 that screen the central region by AV ≈ 90 (AK ≈ 5.5) from the observer.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.