© GRAVITY Collaboration et al. 2020

1. Introduction

The spectral energy distributions (SEDs) of AGN show an excess in the near-infrared (NIR) due to thermal emission from hot dust grains with a color temperature T_color ≃ 1200−1900 K. The dust luminosity and temperature is explained as reprocessed emission from the central accretion disk (e.g., Rieke 1978). In the standard unification paradigm of AGN, dust distributed in a flattened parsec-scale torus with a large covering factor obscures the view of the broad emission lines and central accretion disk for large viewing angles (Antonucci 1993; Urry & Padovani 1995). The NIR radiation might then arise at the inner edge near the sublimation radius, where irradiation from the center is strongest and the dust temperature highest. The physical origin of such a structure remains unclear. Support by gas pressure would require a sound speed far in excess of the Keplerian speed at that distance from the central black hole. Alternatives include radiation or magnetic pressure support, usually in the form of an outflow (e.g., Krolik & Begelman 1988).

Mid-infrared interferometry presents a challenge to the torus paradigm. Many objects show an unresolved core, consistent with an origin in the inner part of the torus. At the same time, a significant fraction of the total flux may originate in the polar region on parsec and larger scales, attributed to dusty outflows from the center (e.g., López-Gonzaga et al. 2016; Hönig & Kishimoto 2017). It remains unclear whether this material produces significant obscuration of the central source. The first resolved image of hot dust found from GRAVITY observations of NGC 1068 shows a size consistent with that expected for the sublimation radius, but in a geometrically thin ring geometry (GRAVITY Collaboration 2020). An additional obscuring structure is required to explain the absence of broad emission lines seen in polarized (scattered) light (Antonucci & Miller 1985).

Continuum reverberation experiments find correlated variability between the optical and NIR emission with a lag that is consistent with reprocessing. The inferred emission radius scales with luminosity as R  ∼  L^1/2 (Suganuma et al. 2006; Koshida et al. 2014), as expected if hot dust radiation peaks near the sublimation radius where the central engine radiation is weak enough for dust to survive (e.g., Barvainis 1987; Netzer 2015). This location depends on dust composition and grain size. The observed normalization of the relation is smaller than predicted for standard ISM dust composition and grain sizes, which may imply the presence of large graphite dust grains near AGN (Kishimoto et al. 2007).

The inferred sub-parsec (≲ milliarcsecond, mas) scales for the brightest type 1 AGN are too compact to be spatially resolved with single telescopes. NIR interferometry, mostly with the Keck Interferometer (Swain et al. 2003; Kishimoto et al. 2009, 2011a) and the Very Large Telescope (VLT) Interferometer (VLTI) instrument AMBER (Weigelt et al. 2012) have measured compact sizes associated with partially resolved sources. The size measurements are consistent with an origin at the sublimation radius, as found from reverberation.

The second-generation VLTI instrument GRAVITY has vastly improved sensitivity and coverage as a result of combining light from all four UT telescopes, resulting in a six-baseline array (GRAVITY Collaboration 2017). We use data from our ongoing AGN observing program to measure sizes for a sample of eight of the brightest type 1 AGN, nearly doubling the sample for which NIR interferometry is available. We describe the data acquisition and selection procedure (Sect. 2) and fitting methods used to measure sizes from both continuum and differential visibilities (Sect. 3). All AGN are partially resolved, with consistent results from both methods (Sect. 4). We find similar angular sizes for objects of similar flux but spanning four orders of magnitude in luminosity, meaning that the physical size of the hot dust emission increases with luminosity. The two luminous quasars in our sample are more compact than the Seyfert 1 galaxies that have been observed. All measured sizes are broadly consistent with the radius-luminosity relation determined using previous NIR interferometry measurements. We discuss the implications of our results in terms of dust emissivity, composition, and the relation of the hot dust emission to that of the broad emission line region (Sect. 5).

2. Observations, data reduction, and data selection

2.1. Observations

The main science goal of our GRAVITY AGN observing program is to spatially resolve the broad emission line region. This has recently been achieved with observations of the quasar 3C 273 (GRAVITY Collaboration 2018). Targets are selected as the brightest Type 1 AGN on the sky visible from the VLTI and above the GRAVITY sensitivity limit (K ≲ 11). For broad-line region (BLR) science we require deep integrations and repeated observations over many nights, with less emphasis on observing calibrators. The observations are designed to probe continuum size and structure from the same data sets. A radius-luminosity relation R ∼ L^1/2 predicts constant angular size on sky for a constant optical flux (Elvis & Karovska 2002). For our targets spanning V ∼ 12−15 we therefore expect a relatively narrow range of continuum angular size ≃0.3−1 mas.

In the past two years we have successfully observed eight type 1 AGN over 15 epochs¹. Details of the targets and observations are given in Table 1, and for more information on the GRAVITY instrument and operations, see GRAVITY Collaboration (2017). In short, for each observation, we first closed the loop of the MACAO visible adaptive optics (AO) systems on the visible AGN source with each of the four UT telescopes. The AGN was then acquired on the GRAVITY acquisition camera, and we placed both the GRAVITY fringe tracking (FT) and science channel (SC) fibers on the AGN (observing on-axis), and split the light between the two detectors. After fringes were acquired, we collected exposures with coherent integrations of 3.3 ms on the FT and 30 s on the medium spectral resolution SC detector (R ∼ 500). We recorded a sequence of exposures of ten SC DITs each, interrupted by occasional sky exposures or calibrator observations. Even the brightest AGN are relatively faint in V, resulting in a poor AO correction even in exceptional conditions. The fringe tracking signal-to-noise ratio is generally low (≲3) at 300 Hz but remains stable due to the excellent performance of the GRAVITY fringe tracker (Lacour et al. 2019). This allows us to simultaneously study continuum (FT) and spectral differential (SC) source properties.

2.2. Data reduction

The data were reduced with the standard GRAVITY pipeline (Lapeyrere et al. 2014; GRAVITY Collaboration 2017) in two separate modes. For the continuum FT visibility data we used the default pipeline settings. The low signal-to-noise ratio and reduced SC visibility amplitude (loss of coherence) often cause the pipeline to flag SC DITs or entire exposures. In analyzing these data, we have found a substantial improvement in residual phase noise (rms scatter) by retaining all data independent of fringe tracker signal-to-noise ratio or V factor and then averaging all DITs by the inverse variance of their differential phases (GRAVITY Collaboration 2018). The differential phase was constructed by removing a mean and slope across the spectrum. We used a method that computes the mean and slope separately for each spectral channel, excluding that channel itself from the measurement (Millour et al. 2008), as implemented in the GRAVITY pipeline. This method improves our differential phase precision by ≃10−20%. The differential phase was constructed in a similar fashion by normalizing the visibility amplitude of each baseline and exposure to its median value, flattening slopes across the wavelength range by filtering out low-frequency modes, and averaging over time with weights measured empirically in each exposure.

2.3. Acquisition camera images and photometry

We estimated the AGN H-band magnitude using the GRAVITY acquisition camera. We fit a 2D Gaussian model to the image from each telescope in each exposure, and estimated the source flux as that integrated over the Gaussian model. We then flux-calibrated using the same acquisition camera measurements of calibrator stars of known magnitude. The results are listed in Table 1. Based on observing the same source on successive nights or using different calibrators within the same night, we estimate the uncertainty in these measurements to be 0.2 mag.

2.4. Data selection

The fringe tracker (FT) keeps the record of fringe measurements in every DIT of 3.3 ms. The group delay and phase delay were calculated in real time to track the variation of the optical path differences of the baselines. The fringe tracker stabilizes the fringe to allow long integrations of the science channel. When the tracking source was bright, most (≳90%) of the FT group delays were ≲3 μm. However, the group delay of a considerable fraction (e.g., ≳50%) of the AGN fringe-tracking data was usually much longer, resulting in considerable visibility loss. As an example, we show in Fig. 1 that when we average all of the FT data, the visibility is much lower than unity. The correlation of the visibility squared with the Strehl ratio implies that the visibility loss is likely affected by the atmosphere. When we only select the DITs with group delay < 3 μm (rejecting 50−80% of DITs), the averaged visibility squared data are visibly improved. The correlation of visibility squared with the Strehl ratio also vanishes. Although other effects might still keep the visibility of the shortest baseline (UT3–UT2) below unity, the data selection based on the group delay apparently helps to improve the data.

3. Visibility model fitting

The typical size scales of hot dust in the brightest type 1 AGN are ≲1 mas, while the maximum VLTI/UT baseline resolution is λ/B ≃ 3 mas. As a result, all of our targets are only partially resolved. In this limit, all source models predict a similar trend of visibility amplitude (or V²) with uv distance (spatial frequency). The observable property is then the characteristic size, related to the normalized second moment of the image (e.g., Lachaume 2003; Johnson et al. 2018). We adopted Gaussian source models throughout to measure the full width at half-maximum (FWHM) angular size of the hot dust emission region using both FT and SC data. We report angular measurements as Gaussian FWHM to avoid assumptions about the hot dust emission, for instance, that it comes from a thin ring near the sublimation radius, as expected for an obscuring torus.

3.1. Continuum visibility fitting

Sample FT continuum V² data for two sources are shown in Fig. 2, color-coded by baseline. In both cases the source is partially resolved, that is, the visibilities fall with increasing baseline length. For NGC 1365, the visibility at zero baseline reaches unity, as it should for a model of a single compact source. In IRAS 09149-6206, it reaches only V²(0) ≃ 0.8. The zero-baseline visibility varies between sources but also between exposures and nights for the same target. We attribute this to a likely coherence loss, although some fraction could also result from extended nuclear K-band continuum emission. To compensate for this effect, we fit for the zero-baseline visibility along with the source size in each exposure for each object over all nights. The function fit is then

Fig. 1. Sample V2 data of PDS 456 on 2018-08-27 before (top) and after (bottom) selecting data with a group delay of < 3 μm. Without the data selection, the averaged visibility squared data are ≲0.65. There is also a trend that the visibility squared increases as the Strehl ratio increases for each baseline. In contrast, the visibility squared reaches ≃0.7−0.8 after the data selection, with no clear trend with Strehl ratio.

Fig. 2. Sample calibrated V2 GRAVITY FT data as a function of baseline length for single-epoch observations of two sources. The data pointed are color-coded by baseline. In both cases the sources are partially resolved, e.g., the visibilities fall with baseline length. For NGC 1365 on 2018 January 8, the V2 approaches 1 at short baselines. For IRAS 09149-6206 on 2018 November 20, the short baselines remain significantly below V2 = 1. The black lines are the best -fitting 1D Gaussian models for each epoch.

$$\begin{array}{c}\hfill V^2=V_0^{2}exp(-2{\pi }^2{r}_{\mathrm{uv}}^2{\sigma }^2),\end{array}$$(1)

for zero-baseline visibility V₀, size parameter σ in radians, and r_uv the baseline length in units of the observed wavelength λ. We report sizes after converting σ into FWHM measured in mas. Sample fits are shown as solid lines in Fig. 2, plotted against all data from one night for each source.

Best-fitting zero baseline $V_0^2$ and FWHM were found individually for each exposure using the scipy function curve_fit. We determined average sizes and uncertainties from the median and rms scatter of all measured sizes over all exposures. To account for correlated systematic errors in calibration, which might be related to AO performance and seeing conditions, we chose an uncertainty in the mean reduced by $\sqrt{N_{\mathrm{epochs}}}$ (nights) rather than by the total number of exposures.

We also explored possible biases in size measurements from continuum FT visibility data resulting from the calibration step using a faint star. Specifically, we tried to measure sizes from unresolved calibrator stars using data from the night of 2018 November 20, where three calibrator stars of K = 8−9 with seven total exposures were taken. We clearly see structure in the reduced V² data that varies over time for the same calibrator and between calibrators with V² ≃ 0.9−1. The variations are largest between objects and over long time separations. Treating each star as a science target, we measured their sizes as described above. The median size inferred can be as large as ≃0.5 mas, but given the scatter between frames, it is consistent with zero in each case. This exercise is conservative in the sense that during science observing, we chose calibrators that were are as close as possible to the science target in time and location on sky, while here the different calibrators were widely separated. As shown below, we also find consistent results for hot dust sizes measured from differential amplitudes, where no calibrator star was used. The calibration step is therefore unlikely to contribute a dominant error to our inferred source sizes.

3.2. Analysis of the differential visibility amplitude

We also independently estimated hot dust continuum visibility amplitudes from the differential amplitude measured in the continuum and broad emission lines (GRAVITY Collaboration 2018, low-redshift Br γ in Seyfert galaxies or Pa α for quasars at z >  0.1).

We constructed the differential visibility by normalizing the visibility amplitude of each baseline. At a spectral channel with line flux f relative to a continuum level of 1, this is

$$\begin{array}{c}\hfill v=\frac{1+f{V}_{\mathrm{l}}/V_{\mathrm{c}}}{1+f},\end{array}$$(2)

where V_c and V_l are the continuum and line visibilities. The differential amplitude v at the line depends on the ratio V_l/V_c. In the partially resolved limit, this in turn measures the quadrature size difference between the line and continuum (Waisberg et al. 2017). For a larger (smaller) spectral line emission region, |v| decreases (increases) at the line. We used these data previously to show that the hot dust emission size is larger than that of the Pa α BLR in 3C 273 (GRAVITY Collaboration 2018). The continuum visibility is

$$\begin{array}{c}\hfill \frac{V_{\mathrm{c}}}{V_{\mathrm{l}}}=\frac{f}{v(1+f)-1}·\end{array}$$(3)

We measured continuum visibilities V_c from the data (v, f) by assuming V_l = 1, for instance, that the broad emission lines are unresolved. The values f were taken from the photometric flux spectra, flattened to remove the instrument profile and averaged over the four UT telescopes and all exposures. They were then normalized to the continuum level. For 3C 273, we found a BLR size from modeling differential phase data of R_BLR = 46 ± 10 μas (GRAVITY Collaboration 2018). This size corresponds to a visibility V_l ≃ 0.996 at a long VLTI baseline of 120 m, justifying our approximation. We note that in general, the BLR is found from reverberation to be a factor of ≳2 smaller than the hot dust continuum (e.g., Koshida et al. 2014.) In this case, our approximation will tend to overestimate continuum visibilities, leading to a ≃10−20% underestimation of the continuum size.

These measurements require deep integrations, but do not suffer from systematic errors related to calibration or coherence loss. We measured V_c from this method in four sources (3C 273, PDS 456, NGC 3783, and IRAS 09149−6206). We expect V_c to be independent over the line spectral channel. From Eq. (2), we then predict an amplitude signature |v| that follows the shape of the spectral line, with a peak at the peak of the line where the contrast between line and continuum images is highest. An example is shown in Fig. 3 for PDS 456, where the visibility amplitude averaged over long baselines shows a significant increase at the spectral line, following the expected behavior.

Sample comparisons of our independent FT continuum and SC differential measurements of the continuum visibility amplitude V_c are shown in Fig. 4. They are generally consistent. Following our analysis of the FT data, we also performed single 1D Gaussian source fits to the SC differential data for each source to infer a characteristic FWHM size.

Fig. 3. Differential visibility amplitude vs. wavelength for PDS 456 from 2018 August 26 averaged over long baselines (blue circles and error bars) compared with the photometric flux. The orange curve shows a model of an unresolved BLR and a Gaussian continuum. The rise of the differential amplitude following the shape of the line is consistent with a partially resolved continuum, and its amplitude allows an independent measurement of continuum hot dust size.

Fig. 4. Sample squared visibilities measured for FT (small circles colored by baseline) compared to those derived from SC differential amplitudes (black circles and their averages as large squares) for 3C 273 (left) and IRAS 09149−6206 (right). The FT data have been divided by the best-fitting zero-baseline visibility, so that V2(0) = 1. With this scaling, the two independent measurements are generally consistent.

4. Results

4.1. Elongation and asymmetry

We do not find evidence for 2D (elongated) or asymmetric structure in either the SC or FT data, except in NGC 3783, which will be analyzed in more detail elsewhere. Figure 5 shows distributions of closure phases combined for 3C 273, PDS 456, IRAS 09149−6206, and Mrk 509 (our sources with the most and highest precision data). The closure phase is formed by summing the visibility phase over baseline triangles, and is immune to telescope-based phase errors. The distributions show a typical closure phase rms of ±1° with a median consistent with zero, indicating a symmetric structure down to < 0.1 mas scales² The closure phases measured for NGC 3783 are always < 0° (shown separately), and are indicative of asymmetry. For consistency with the other sources, we report size measurements for NGC 3783 from 1D Gaussian model fits, even though the model is inconsistent with the observed nonzero closure phases.

Fig. 5. Distributions of measured closure phases for NGC 3783 (orange) and other AGN targets with sufficient data (3C 273, PDS 456, Mrk 509, and IRAS 09149−6206).

4.2. Hot dust is more extended than the BLR

In all sources where we detect a differential amplitude signature, the amplitude increases at the line, showing robustly that the broad emission line region is more compact than the hot dust continuum. In the remaining targets we were limited by sensitivity, that is, the size difference is not constraining. This finding is consistent with previous reverberation and spectral measurements (e.g., Netzer 2015, and references therein). In 3C 273 we also have an interferometric BLR radius measured from kinematic modeling of the detected velocity gradient in differential phase data (GRAVITY Collaboration 2018), which is a factor ≃3 smaller than that of the hot dust.

4.3. Hot dust size measurements

The FWHM 1D Gaussian size measurements for each source are shown in Fig. 6 and are listed in Table 2 for the FT and where available, also the SC measurements. All sources are partially resolved with sizes ≲1 mas, and as compact as ≃0.3 mas (3C 273). These small sizes relative to the VLTI interferometric beam are still detected due to the high sensitivity and improved uv-coverage of GRAVITY. The sizes measured by the SC and FT methods are generally consistent. The scatter in our size measurements is dominated by systematic calibration errors, which cause scatter in the sizes measured in individual exposures that are much larger than expected by the signal-to-noise ratio of the individual V² measurements. This is reflected in the scatter of their zero-baseline visibilities V₀.

Fig. 6. Average source angular sizes (Gaussian FWHM) and their uncertainties measured from both continuum (FT) and spectral differential (SC) data. The latter are measured in deep integrations and are only possible for four objects so far.

Previous NIR interferometry measured sizes using a thin ring (delta function) intensity distribution rather than a Gaussian, as we have done here. Keck interferometer observations found a hot dust radius of the quasar 3C 273 of 0.25 ± 0.1 mas assuming a thin ring model (Kishimoto et al. 2011a), which is equivalent to a Gaussian HWHM of 0.21 ± 0.09 mas and consistent with our Gaussian HWHM result of 0.15 ± 0.03 mas. VLTI/AMBER observations found a thin ring radius of 0.74 ± 0.23 mas for NGC 3783 (Gaussian HWHM of 0.56 ± 0.17 mas), consistent with our result of 0.42 ± 0.06 mas. We note that the Keck observations used only one baseline (compared to our six baselines), while the VLTI/AMBER data used three baselines but at lower signal-to-noise ratio. In analyzing the AMBER data, Weigelt et al. (2012) further assumed a fixed zero-baseline visibility of unity.

4.4. BLR size estimates

We have inferred continuum visibilities and sizes from SC differential data by assuming a point source BLR (V_l = 1). Alternatively, we can take the measured continuum size from the FT data and use the SC differential data to infer the BLR size. For 3C 273, this gives R_BLR = 75 ± 65 μas, consistent with both GRAVITY modeling of the SC differential phases (GRAVITY Collaboration 2018) and various estimates from reverberation mapping (Kaspi et al. 2000; Peterson et al. 2004; Zhang et al. 2019). For IRAS 09149−6206, the result is an upper limit of R_BLR ≲ 100 μas. For PDS 456 and NGC 3783, taken at face value, we would obtain large BLR sizes R_BLR = 200 ± 70 μas and R_BLR = 320 ± 100 μas. The uncertainties are large because the SC differential amplitude is weakly sensitive to the BLR size for such small angular sizes relative to the baseline resolution. We also note that any systematic offsets between our FT and SC results would likely change the results by a large amount compared to the statistical uncertainty. Still, combining FT and SC data with sufficient sensitivity in principle provides a measurement of BLR size that does not rely on adopting a kinematic model of the broad emission line profile.

5. Discussion

We have measured the K-band continuum sizes of eight bright type 1 AGN using NIR interferometry with the VLTI instrument GRAVITY. We used continuum visibilities as measured by the GRAVITY fringe tracker V², and from the spectral differential visibility amplitudes. Each source was partially resolved, with FWHM sizes ≃0.3−0.8 mas smaller than the VLTI beam ≃3 mas. The narrow range of angular sizes is expected given the similar optical fluxes of our targets (V ≃ 12−15). We find no evidence for significant elongation or asymmetry in seven of eight targets observed (see Appendix B). We confirm that the BLR is more compact than the hot dust emission from direct measurements of the spectral differential visibility amplitude. Here we discuss the implications of our results for the size, structure, and physical properties of hot dust near AGN.

5.1. Hot dust surface emissivity

We measured the hot dust surface emissivity ϵ_ν ≡ I_ν/B_ν (product of grain emissivity and surface filling factor) for objects in our sample, assuming a true dust temperature T = 1500 K (Table 2). This metric is equivalent to the surface brightness reported by Kishimoto et al. (2011b). Because our range of angular sizes is mostly consistent with previous work, the dust surface emissivity values are as well, ϵ_ν ≃ 0.1−0.2. The compact angular size of 3C 273 leads to a higher emissivity of ϵ_ν = 0.7 ± 0.3.

Because we did not resolve the emission region, our inferred values correspond to averages over the effective VLTI beam and are lower limits to the peak hot dust emissivity. The correction could be large if a significant amount of the K-band emission originates in a narrow ring (Kishimoto et al. 2009; GRAVITY Collaboration 2020). Such high surface emissivities close to blackbody emission can be expected if the dust emission is dominated by large grains ≳0.2 μm as compared to a standard ISM grain size distribution with mean grain size < 0.1 μm. This is consistent with differential grain sublimation in the AGN vicinity (Hönig & Kishimoto 2017).

5.2. Hot dust radius-luminosity relation

Near-infrared reverberation (Glass 1992; Suganuma et al. 2006; Koshida et al. 2014) and interferometry (Kishimoto et al. 2009, 2011a; Weigelt et al. 2012) observations have used hot dust continuum sizes to measure a radius-luminosity relation, R ∼ L^1/2, consistent with the scaling for hot dust originating near the sublimation radius. At the same time, the normalization of this relation indicates an emission region that is more compact than expected, which might be explained by the presence of large graphite dust grains (Kishimoto et al. 2007). A relationship R ∼ L^1/2 predicts an equal angular size at constant observed flux (e.g., Elvis & Karovska 2002). Here all objects are K ≃ 10 and V ≃ 12−15. Their relatively constant angular sizes ≃0.3−0.8 mas over four orders of magnitude in bolometric luminosity show that the physical radius indeed increases with luminosity.

Our observations nearly double the size of the NIR interferometry sample, and are measured with significantly improved uv-coverage and sensitivity. In the partially resolved limit, the source size corresponds to the normalized second moment of the image. For comparison with previous work, we measured physical radii by converting our Gaussian HWHM measurements into thin ring angular radii (in the partially resolved limit, divide by a factor of $\sqrt{ln2}\simeq 0.8$). We then converted into a physical radius using the angular diameter distances to our targets. Significant point-source contributions to the NIR source (from the accretion disk and/or jet) could cause us to underestimate the hot dust size. Then the apparent source size for a point source plus extended model can be written

$$\begin{array}{cc}& {\sigma }_{\mathrm{app}}^2=f_{\mathrm{pt}}{\sigma }_{\mathrm{pt}}^2+(1-f_{\mathrm{pt}}){\sigma }_{\mathrm{true}}^2,\hfill \end{array}$$(4)

$$\begin{array}{cc}& \frac{{\sigma }_{\mathrm{app}}}{{\sigma }_{\mathrm{true}}}\simeq \sqrt{1-f_{\mathrm{pt}}},\hfill \end{array}$$(5)

where σ are normalized second moments for the point-source component (σ_pt = 0), the total apparent source, and the true hot dust component. The fractional point source flux, f_pt, varies per object. We adopted a typical constant value f_pt = 0.2 in converting into physical radius (Kishimoto et al. 2007), which results in a size increase of ≃10%. The bolometric luminosities of our objects and those from the literature were estimated uniformly as described in Appendix A. The physical thin ring radii and bolometric luminosity values we used are listed in Table 2.

Figure 7 shows the resulting physical hot dust emission radius versus bolometric luminosity from previous and our new measurements. Fitting a power law $R\sim L_{\mathrm{bol}}^{1/2}$ to reverberation measurements from the literature (Suganuma et al. 2006; Kishimoto et al. 2007; Koshida et al. 2014), we obtain the dashed line with a normalization R₄₄ = 0.038 ± 0.02 pc, where R₄₄ is the radius cΔt at a bolometric luminosity of 10⁴⁴ erg s⁻¹ and Δt is the measured lag between the K and V bands (e.g., Koshida et al. 2014). Leaving the slope free, we find a best fit with a flatter luminosity dependence, $R\sim L_{\mathrm{bol}}^{0.40\pm 0.04}$. The individual data points are plotted as blue crosses. Combining the Seyfert galaxies in our sample (L_bol ≲ 10⁴⁶ erg s⁻¹) with previous KI/AMBER results, we find a consistent best-fitting radius-luminosity relation from NIR interferometry (slope of 0.38) as from reverberation (slope of 0.42). The two luminous quasars (L_bol >  10⁴⁶ erg s⁻¹, PDS 456 and especially 3C 273) show smaller hot dust angular sizes than the other sources. These sizes are robust, obtained with both SC and FT data, and precise (uncertainties ≲20%). The measurement for 3C 273 agrees with that reported by Kishimoto et al. (2009). We generally find that interferometry sizes are larger than those from reverberation at comparable bolometric luminosity, in agreement with previous work (e.g., Kishimoto et al. 2011a; Koshida et al. 2014).

Fig. 7. Radius−luminosity relation for literature reverberation and interferometry size measurements (blue and orange crosses) compared to our new GRAVITY measurements (filled circles). The dashed line is the best-fitting R ∼ L1/2 fit to reverberation measurements (Suganuma et al. 2006; Koshida et al. 2014). We use FT measurements and our own bolometric luminosity estimates in all cases. The radii reported here for all interferometry measurements is that of a thin ring with an assumed unresolved point-source flux fraction of 20% (see Sect. 5.2 for details).

Figure 8 shows NIR interferometry radius measurements normalized to the R − L relation from reverberation. The scatter between objects is larger than their errors, showing that there are differences in viewing geometry and/or physical structure (geometry or dust composition). At low luminosities, the interferometry sizes are systematically above the relation. Larger interferometric sizes might be expected as a result of an extended NIR emission region due to their weighting by total intensity, rather than by the response to a variable central source (e.g., Koratkar & Gaskell 1991). Kishimoto et al. (2011b) argued that the reverberation value corresponds to the sublimation radius, R_in, and that values R/R_in ≃ 1 imply a hot dust emissivity falling rapidly away from that location. Higher values instead correspond to a shallower emissivity. Figure 8 shows that R/R_in falls with increasing bolometric luminosity, if with significant scatter. One explanation is a variable radial dust emission profile, which is sharply peaked at the sublimation radius for more luminous objects and more extended at the lowest luminosities. Our data support this interpretation.

Fig. 8. Radius from GRAVITY and Keck/AMBER NIR interferometry, measured relative to the empirical radius-luminosity relation from reverberation mapping (Suganuma et al. 2006; Koshida et al. 2014). Reverberation sizes are shown for comparison, and the error bars correspond to 1σ in all cases. Both interferometry data sets show a larger median size R/Rin >  1 and decreasing relative size with increasing luminosity, particularly evident as R/Rin ≲ 1 for quasars.

5.3. Compact hot dust in quasars

The two luminous sources also have L/L_Edd  ∼  1 and could have particularly high point-source flux contributions. A bias by a factor of 2, needed to bring 3C 273 within the observed range of other objects, requires a point-source fraction f_pt = 3/4, that is, the observed K-band emission needs to be dominated by contaminating emission. Kishimoto et al. (2011a) estimated an accretion disk contribution for 3C 273 of f_pt ≃ 0.3, similar to our assumed f_pt  =  0.2. Nonthermal jet emission could also contribute in the NIR, as seen in rapidly variable K-band light curves (e.g., Robson et al. 1993; McHardy et al. 2007; Bewketu Belete et al. 2020), corresponding to flaring events. Our observations of 3C 273 have produced consistent acquisition camera H-band magnitudes and size measurements from July 2017 to May 2018. The acquisition camera H-band flux density of ≃40 mJy is consistent with the quiescent state of 3C 273 with no evidence of flare contributions (e.g., Robson et al. 1993; Soldi et al. 2008). Robson et al. (1993) estimated a jet contribution to the quiescent flux density of ≃30 mJy from extrapolating the power-law spectral slope seen at millimeter wavelengths. Soldi et al. (2008) found a significant accretion disk fraction but a much smaller jet contribution in SED fitting due to an imposed spectral cutoff to the synchrotron spectrum. Kishimoto et al. (2011a) estimated this contribution to be ≲8 mJy based on the low NIR polarization fraction and assuming an intrinsic ≃10% net polarization for the synchrotron component. These estimates produce a wide range of possible point-source flux fractions from jet emission (assuming a small jet emission size), f_jet ≃ 0.1−0.4. When we combine the highest inferred values for both the accretion disk and jet contributions, the total point-source contribution might well be large enough to explain our observed compact size. A total point-source contribution ≲0.5 seems more likely based on the clear thermal dust bump in the NIR SED (Soldi et al. 2008), and would increase our hot dust size by ≲20%. PDS 456 also shows a compact size and is radio quiet. We conclude that contaminating emission can contribute to the compact sizes we measure, particularly for 3C 273, but probably does not explain our finding of R/R_in ≲ 1 in two luminous quasars.

A further bias comes from errors in the bolometric luminosity estimates. Even random errors in L_bol will tend to produce an anticorrelation of R/R_in and L_bol because $R_{\mathrm{in}}\propto L_{\mathrm{bol}}^{1/2}$. By simulating many random data sets, we find that the degree of anticorrelation observed in GRAVITY and KI/AMBER data could be explained if the true luminosity errors are ≃0.6 dex. Our measured errors by comparing various methods appear to be somewhat smaller ≃0.3−0.5 dex (Appendix A). Still, this bias warrants caution in interpreting the results in terms of varying hot dust structure with source luminosity.

5.4. Geometric distance estimates

Previous measurements of the hot dust radius from reverberation have been made for NGC 3783 (Lira et al. 2011), Mrk 335, and Mrk 509 (Koshida et al. 2014). Geometric estimates of the angular diameter distance from combining time lags with interferometric sizes are given in Table 3. For these estimates we have followed Koshida et al. (2014) in assuming R = cΔt = D_Aθ_ring, where Δt is the reverberation lag and θ_ring is the angular radius of a thin ring model with a 20% unresolved point-source contribution. We find consistent angular diameter distances to the fiducial values given the redshift for Mrk 335, as well as for NGC 3783 when we use its differential SC angular size. The GRAVITY-inferred physical radius of Mrk 509 is a factor ≳2 larger than from reverberation because of its fiducial angular diameter distance. A similar discrepancy is found using the continuum FT size of NGC 3783, although this object clearly shows complex structure beyond a partially resolved thin ring. Understanding the discrepancy and the physical structure of the NIR emission region will be important for efforts to measure geometric distances by combining interferometric and reverberation sizes (Hönig et al. 2014).

6. Summary

We have reported NIR interferometry measurements of eight AGN made with the second-generation VLTI instrument GRAVITY. In all cases, we partially resolved the continuum hot dust emission region and were able to measure its size. In the four objects with sufficiently deep integrations, we used a new spectral differential method to show definitively that the hot dust continuum is much larger than the low-ionization Pa α or Br γ broad emission line region. In seven of eight objects, we found no clear evidence for elongation or asymmetry and in particular constrained the closure phases to ≲1°. The hot dust continuum sizes span 0.3−0.8 mas, mostly consistent with previous interferometry measurements with the Keck interferometer and the VLTI instrument AMBER. Our objects span four orders of magnitude in bolometric luminosity with similar V magnitude. Their roughly constant angular size implies an increasing physical radius with luminosity, consistent with previous reverberation and interferometry measurements. We find that at low luminosity, the hot dust sizes are systematically larger than those from reverberation, which might be related to emission from an extended region. At high luminosity, the two quasars in our sample are compact, with sizes consistent with sharply peaked emission at the sublimation radius. Accretion disk and/or jet contributions to the flux and biases from errors in estimated bolometric luminosities are probably not sufficient to fully explain this finding.

Acknowledgments

J. D., M. B., and A. J. R. were supported by a Sofja Kovalevskaja award from the Alexander von Humboldt foundation and in part by NSF grant AST 1909711. SFH acknowledges support by the EU Horizon 2020 framework programme via the ERC Starting Grant DUST-IN-THE-WIND (ERC-2015-StG- 677117). M. K. acknowledges support from JSPS (16H05731). F. E. and O. P. acknowledge support from ERC synergy grant No. 610058 (Black-HoleCam). 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

Appendix A: Bolometric luminosity

The bolometric luminosity of the AGNs are estimated mainly with X-ray measurements. Unless noted in particular, we take the 14–195 keV flux from the BAT AGN Spectroscopic Survey (Koss et al. 2017), taken directly from the 70-month Swift-BAT survey (Baumgartner et al. 2013). The hard X-ray emission is less strongly contaminated by the host galaxy and less affected by the absorption. We adopted the luminosity-dependent bolometric correction relation provided by Winter et al. (2012),

$$\begin{array}{c}\hfill log\left(\frac{L_{\mathrm{bol},14-195\mathrm{keV}}}{{\mathrm{erg}\mathrm{s}}^{-1}}\right)=1.1157log\left(\frac{L_{14-195\mathrm{keV}}}{{\mathrm{erg}\mathrm{s}}^{-1}}\right)-4.2280.\end{array}$$(A.1)

We also collected their optical (5100 Å) continuum from BAT catalog and 12 μm nuclear emission from Asmus et al. (2014) to confirm the consistency of different methods for deriving the bolometric luminosity. Following the 5100 Å bolometric correction of Trakhtenbrot et al. (2017), we used

$$\begin{array}{c}\hfill log\left(\frac{L_{\mathrm{bol},5100\mathrm{\AA }}}{{\mathrm{erg}\mathrm{s}}^{-1}}\right)=0.916log\left(\frac{\nu L_{\nu }(5100\mathrm{\AA })}{{\mathrm{erg}\mathrm{s}}^{-1}}\right)+4.596.\end{array}$$(A.2)

For the mid-IR (MIR) data, we first converted the 12 μm flux into 2–10 keV intrinsic flux according to Asmus et al. (2011),

$$\begin{array}{c}\hfill log\left(\frac{f_{2-10\mathrm{keV}}}{\mathrm{mJy}}\right)=0.89log\left(\frac{f_{12\mu \mathrm{m}}}{\mathrm{mJy}}\right)-12.81.\end{array}$$(A.3)

Then we converted the 2–10 keV luminosity into 14–195 keV luminosity according to Winter et al. (2009),

$$\begin{array}{c}\hfill log\left(\frac{L_{14-195\mathrm{keV}}}{\mathrm{erg}{\mathrm{s}}^{-1}}\right)=0.94log\left(\frac{L_{2-10\mathrm{keV}}}{{\mathrm{erg}\mathrm{s}}^{-1}}\right)+2.91,\end{array}$$(A.4)

and relied on the 14–195 keV bolometric correction (Winter et al. 2012) to obtain the final bolometric luminosity. The latter step was made to maintain consistency with our X-ray bolometric luminosity. The intrinsic scatter of the relation in Eq. (A.4) may introduce additional uncertainty. However, this is probably not as important considering the large uncertainty of the bolometric correction (see below). We compare the bolometric luminosities derived from the three methods in Fig. A.1. The scatters of the relations are ∼0.3–0.5 dex, which likely reflects the uncertainty of the bolometric correction methods at a single frequency. Considering the source variability and the difference of their intrinsic SEDs, the uncertainty of our bolometric correction is hardly below ∼0.5 dex (see Netzer 2019 and discussion therein). Therefore we conclude that the three methods produce mostly consistent results. A few objects are noted individually below, including the four objects (PDS 456, Mrk 231, IRAS 13349+2438, and PG 1202+281) that are not available in the BAT AGN sample.

Mrk 335. It is reported to be extremely variable in X-rays with an amplitude of a factor of ∼10 and a timescale of tens of days (Grupe et al. 2012). We therefore preferred to adopt the 5100 Å bolometric luminosity, which is consistent with the value reported by Woo & Urry (2002), who calculated the bolometric luminosity by integrating the observed flux over the SED.

3C 120. The bolometric luminosity derived from 5100 Å luminosity is about one order of magnitude lower than those with other methods. This might be due to variability. However, we note that the 5100 Å luminosity reported by Vestergaard & Peterson (2006) is ∼0.9 dex higher than that of the BAT AGN Spectroscopic Survey. The former provides a bolometric luminosity that is more consistent with our other methods. The X-ray derived bolometric luminosity is also consistent with that integrated over the SED from Woo & Urry (2002).

3C 273. The jet emission probably contributes significantly to the X-ray luminosity of 3C 273 (Dermer et al. 1997; Vasudevan & Fabian 2007). We therefore discarded the bolometric luminosities derived from X-ray data and adopted the value based on 5100 Å, which is found to be consistent with that based on 12 μm data.

PDS 456. This source is very luminous in the UV/optical. Its 5100 Å luminosity is ∼2 × 10⁴⁶ erg s⁻¹ (Simpson et al. 1999; Reeves et al. 2009), corresponding to L_{bol,5100 Å} = 1.0 × 10⁴⁷ erg s⁻¹. This is very well consistent with the bolometric luminosity derived from UV–IR SED (Simpson et al. 1999). PDS 456 is extremely variable in X-rays (Reeves et al. 2002), and the observed 2–10 keV luminosity is only 0.2% of the bolometric luminosity (Reeves et al. 2009). We therefore do not quote the bolometric luminosity estimated from X-rays.

Mrk 231. Recent NuSTAR measurements revealed Mrk 231 to be intrinsically X-ray weak (Teng et al. 2014). Therefore, X-ray measurements are not suitable to derive the bolometric luminosity for this target. We adopted the 12 μm measurement from Lopez-Rodriguez et al. (2017) and obtained a bolometric luminosity of ∼6.6 × 10⁴⁵ erg s⁻¹, which is consistent with the bolometric luminosity, 4−5 × 10⁴⁵ erg s⁻¹, derived from the decomposition of IR SED (Farrah et al. 2003; Teng et al. 2014).

IRAS 13349+2438. We derived the bolometric luminosity of IRAS 13349+2438, ∼1.8 × 10⁴⁶ erg s⁻¹, based on 12 μm data. This is within a factor of 2 consistent with the estimates based on quasar SEDs (Beichman et al. 1986; Lee et al. 2013).

IC 4329A. The bolometric luminosity derived from X-ray data is ∼0.9 dex higher than that from 5100 Å. This probably reflects the uncertainty of the different bolometric luminosity corrections. The 5100 Å luminosity of the BAT catalog is consistent with other measurements (e.g., Bentz et al. 2009). The X-ray bolometric luminosity is consistent with the values derived from the AGN SEDs (Woo & Urry 2002; Vasudevan et al. 2010). We therefore still adopted the X-ray derived bolometric luminosity in our analysis. It is worth noting that with the bolometric correction theoretically calculated by Netzer (2019), the bolometric luminosity based on 5100 Å luminosity is ∼0.4 dex higher than what we quote, and as a result is more consistent with that derived from X-rays.

PG 1202+281. We collected the 5100 Å data from Vestergaard & Peterson (2006). The derived the bolometric luminosity is well consistent with that derived from NIR-to-X-ray SED (Runnoe et al. 2012).

Mrk 744. Our X-ray derived bolometric luminosity is well consistent with that reported in Woo & Urry (2002) from interpolating the AGN SED.

Mrk 110. We collected the 5100 Å luminosity from Vestergaard & Peterson (2006).

IRAS 03450+0055. The 5100 Å luminosity comes from Afanasiev et al. (2019).

Fig. A.1. Comparison of AGN bolometric luminosities derived using the different methods listed in Table A.1. The median and standard deviation of the difference (y − x) are indicated in the upper left corner of each panel. The dashed lines indicate the one-to-one relation. a: 14–195 keV and 12 μm results show a possible offset, comparable to their scatter (≃0.3 dex). The 5100 Å bolometric luminosity tends to be better consistent with both 14–195 keV (b) and 12 μm (c), although the scatters of both relations are considerable. Considering the uncertainty of the bolometric correction is large (∼0.5 dex), the different methods are relatively well consistent with each other.

Appendix B: 2D visibility plots

Figure B.1 shows the calibrated V² data for all sources that we have used here.

Fig. B.1. 2D maps of V2 (colors) are shown for each source as a function of observed uv coordinates. The point size is proportional to its signal-to-noise ratio. All sources are partially resolved (on average, lighter colors closer to the origin), without clear evidence of elongation (systematic variation of V2 with position angle), except potentially for NGC 3783. The inner and outer dashed circles show baseline resolutions of 10 and 4 mas.

All Tables

All Figures

Fig. 1. Sample V2 data of PDS 456 on 2018-08-27 before (top) and after (bottom) selecting data with a group delay of < 3 μm. Without the data selection, the averaged visibility squared data are ≲0.65. There is also a trend that the visibility squared increases as the Strehl ratio increases for each baseline. In contrast, the visibility squared reaches ≃0.7−0.8 after the data selection, with no clear trend with Strehl ratio.

In the text

Fig. 2. Sample calibrated V2 GRAVITY FT data as a function of baseline length for single-epoch observations of two sources. The data pointed are color-coded by baseline. In both cases the sources are partially resolved, e.g., the visibilities fall with baseline length. For NGC 1365 on 2018 January 8, the V2 approaches 1 at short baselines. For IRAS 09149-6206 on 2018 November 20, the short baselines remain significantly below V2 = 1. The black lines are the best -fitting 1D Gaussian models for each epoch.

In the text

Fig. 3. Differential visibility amplitude vs. wavelength for PDS 456 from 2018 August 26 averaged over long baselines (blue circles and error bars) compared with the photometric flux. The orange curve shows a model of an unresolved BLR and a Gaussian continuum. The rise of the differential amplitude following the shape of the line is consistent with a partially resolved continuum, and its amplitude allows an independent measurement of continuum hot dust size.

In the text

Fig. 4. Sample squared visibilities measured for FT (small circles colored by baseline) compared to those derived from SC differential amplitudes (black circles and their averages as large squares) for 3C 273 (left) and IRAS 09149−6206 (right). The FT data have been divided by the best-fitting zero-baseline visibility, so that V2(0) = 1. With this scaling, the two independent measurements are generally consistent.

In the text

	Fig. 5. Distributions of measured closure phases for NGC 3783 (orange) and other AGN targets with sufficient data (3C 273, PDS 456, Mrk 509, and IRAS 09149−6206).
In the text

	Fig. 6. Average source angular sizes (Gaussian FWHM) and their uncertainties measured from both continuum (FT) and spectral differential (SC) data. The latter are measured in deep integrations and are only possible for four objects so far.
In the text

Fig. 7. Radius−luminosity relation for literature reverberation and interferometry size measurements (blue and orange crosses) compared to our new GRAVITY measurements (filled circles). The dashed line is the best-fitting R ∼ L1/2 fit to reverberation measurements (Suganuma et al. 2006; Koshida et al. 2014). We use FT measurements and our own bolometric luminosity estimates in all cases. The radii reported here for all interferometry measurements is that of a thin ring with an assumed unresolved point-source flux fraction of 20% (see Sect. 5.2 for details).

In the text

Fig. 8. Radius from GRAVITY and Keck/AMBER NIR interferometry, measured relative to the empirical radius-luminosity relation from reverberation mapping (Suganuma et al. 2006; Koshida et al. 2014). Reverberation sizes are shown for comparison, and the error bars correspond to 1σ in all cases. Both interferometry data sets show a larger median size R/Rin >  1 and decreasing relative size with increasing luminosity, particularly evident as R/Rin ≲ 1 for quasars.

In the text

Fig. A.1. Comparison of AGN bolometric luminosities derived using the different methods listed in Table A.1. The median and standard deviation of the difference (y − x) are indicated in the upper left corner of each panel. The dashed lines indicate the one-to-one relation. a: 14–195 keV and 12 μm results show a possible offset, comparable to their scatter (≃0.3 dex). The 5100 Å bolometric luminosity tends to be better consistent with both 14–195 keV (b) and 12 μm (c), although the scatters of both relations are considerable. Considering the uncertainty of the bolometric correction is large (∼0.5 dex), the different methods are relatively well consistent with each other.

In the text

Fig. B.1. 2D maps of V2 (colors) are shown for each source as a function of observed uv coordinates. The point size is proportional to its signal-to-noise ratio. All sources are partially resolved (on average, lighter colors closer to the origin), without clear evidence of elongation (systematic variation of V2 with position angle), except potentially for NGC 3783. The inner and outer dashed circles show baseline resolutions of 10 and 4 mas.

In the text