© GRAVITY Collaboration 2021

1. Introduction

Recent advancements in infrared interferometric observations have allowed significant progress in understanding the gas and dust structure and dynamics around active galactic nuclei (AGNs). With the new capabilities of GRAVITY (GRAVITY Collaboration 2017), the second-generation instrument at the Very Large Telescope Interferometer (VLTI), the broad line region (BLR) has been resolved and modelled to provide a measurement of the supermassive black hole mass (GRAVITY Collaboration 2018, 2020a), the hot dust has been imaged revealing a thin ring at the dust sublimation radius (GRAVITY Collaboration 2020c), and hot dust sizes have been measured for an increasing number of AGNs (GRAVITY Collaboration 2020b). NGC 3783 provides a unique laboratory for studying not only these aspects but also the coronal line region (CLR) in a single object, and can be used to build a comprehensive picture of the nuclear and circumnuclear region around an AGN.

NGC 3783 hosts one of the most luminous local AGNs, with a bolometric AGN luminosity off log L_AGN ∼ 44.5 erg s⁻¹ for a distance of 38.5 Mpc (Tully & Fisher 1988; Davies et al. 2015)¹, and has been studied extensively, most notably for its ionised outflows (especially X-ray warm absorbers) and variability. As one of the most luminous AGNs, NGC 3783 has been intensely monitored via reverberation mapping (RM) and is one of the few AGNs for which simultaneous UV (Reichert et al. 1994) and optical (Stirpe et al. 1994) reverberation results have been obtained. NGC 3783 is also one of the AGNs that clearly demonstrates the virial relationship between emission line lag and width (Onken & Peterson 2002). Through RM, the black hole mass of NGC 3783 has been estimated to be ∼3 × 10⁷ M_⊙ assuming a virial factor of f_σ = 4−5 (Peterson et al. 2004; Bentz & Katz 2015), which implies that the AGN is radiating at ∼0.1 L_Edd. For a velocity dispersion of 130 km s⁻¹ as measured from both the Ca II triplet and the CO 2–0 bandhead (Caglar et al. 2020), this puts the object very close to the M_BH − σ* relation (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Ferrarese et al. 2001; Nelson et al. 2004; Onken et al. 2004; Gültekin et al. 2009; McConnell & Ma 2013). The black hole mass also matches that derived from the Fe Kα line at 6.4 keV (Brenneman et al. 2011; Capellupo et al. 2017) when modelled together with the spin of the black hole.

The X-ray-warm absorber was initially modelled by Netzer et al. (2003), and was found to require three different ionisation components at two different velocities, and at distances within limits of 0.2−25 pc, with a total column of N_H ∼ 4 × 10²² cm⁻². With more detailed spectra, this was later expanded to three velocities in the range −460 to −1600 km s⁻¹ by Mao et al. (2019), comparable to those identified previously in UV spectra of −550, −720, and −1370 km s⁻¹ (Kraemer et al. 2001). An X-ray obscuration event lasting about a month, which led to a significantly reduced flux at energies ≲5 keV and a modified ionisation structure, was reported by Mehdipour et al. (2017). Their modelling suggested this component had a column of N_H ∼ 10²³ cm⁻² and a covering factor of ∼0.5, and simultaneous UV observations (Kriss et al. 2019) indicated it was moving out at a velocity of ∼2000 km s⁻¹. In addition, from the implied density of 3 × 10⁹ cm⁻³ and similarity of the radial location at only ten light-days to the size of the BLR from RM (Onken & Peterson 2002), these authors argued that it is an obscuring wind originating in the outer part of the BLR. Subsequent analysis of archival data by Kaastra et al. (2018) suggested that such obscuring events, with columns exceeding 5 × 10²¹ cm⁻², may be rather common for NGC 3783.

Ionised outflows have also been observed on larger scales of tens to hundreds of parsecs (Rodríguez-Ardila et al. 2006). Guided by the almost unresolved (< 200 pc diameter) appearance of the optical [OIII] image, Fischer et al. (2013) modelled the kinematics in a slit spectrum as being consistent with an ionisation cone that is oriented 15° from face-on with inner and outer opening angles of 45 and 55°, respectively, and a maximum outflow velocity of 130 km s⁻¹. In contrast, based on the biconical appearance of the near-infrared (near-IR) [SiVI] distribution at 16 pc resolution, Müller-Sánchez et al. (2011) modelled the line kinematics as an ionisation cone that is 60° from face-on with inner and outer opening angles of 27 and 34° and a maximum outflow velocity of 400 km s⁻¹. Neither of these are completely satisfactory: the former is more consistent with what is expected for a Seyfert 1 while the latter matches the expected outflow velocity better and explains the [SiVI] morphology. However, face-on orientations appear to be ruled out by mid-IR interferometric data which show strong elongation attributed to polar dust at a position angle of −50 to −60° (Hönig et al. 2013; Burtscher et al. 2013; López-Gonzaga et al. 2016). Fitting the near- to mid-IR spectral energy distribution (SED) using a disc+wind model, Hönig & Kishimoto (2017) reached a conclusion that is plausibly consistent with both, favouring an inclination of 30° and an opening angle of 38°. This would imply that in NGC 3783, our line of sight (LOS) is close to the edge of the ionisation cone, which could perhaps explain the slight extinction to the BLR (A_V ∼ 0.1 mag, Schnorr-Müller et al. 2016a), the classification as a Sy 1.2−1.5, and the frequent X-ray obscuration events.

These results indicate a polar axis that is oriented several tens of degrees west of north, consistent with that deduced from optical polarisation measurements. Applying an equatorial disc model, Smith et al. (2002, 2004) concluded from polarisation data that the polar axis is at −45°. However, there are several issues that make this conclusion uncertain: NGC 3783 does not show the expected position angle rotation across the Hα line profile, the depolarisation in the line core may partially be due to the narrow lines, and the polarisation peak in the line wing is not unique to this specific model. Lira et al. (2020) also noted that NGC 3783 has unusual polarisation characteristics, with the position angle showing a more M-like profile possibly indicative of a radially outflowing scatterer. Further, both the mid-IR and polarisation measurements could be biased towards the edges of the outflow cones if this is the location of most of the mass in the cones and this bias has been hinted at recently by observations of Circinus (e.g., Stalevski et al. 2017, 2019). Unfortunately, radio maps do not help in determining the orientation of the polar axis because they are unresolved at scales of 0.2−0.6″ at 8.5 GHz as well as 10−30 mas at 1.6 GHz (Schmitt et al. 2001; Orienti & Prieto 2010). Ironically, in both cases the beam is elongated in roughly the same direction as the polar axis, making it harder to assess whether there is a small-scale radio jet.

To compound the issue surrounding the inner geometry of NGC 3783, the orientation of the host galaxy disc is such that the kinematic major axis, measured from both the stellar and H₂ 1–0 S(1) velocity fields on scales of 30−300 pc, are consistent with −30 to −40° (Davies 2007; Hicks et al. 2009; Lin et al. 2018) which matches the kinematic axis on kiloparsec scales measured from Hα in a recent VLT/MUSE observation (den Brok et al. 2020). These are very similar to the −45° kinematic major axis on tens of kiloparsecs measured from HI (García-Barreto et al. 1999), suggesting that over most of the galaxy disc there is little warping. This means that the host galaxy kinematic axis and the outflow direction, at least on scales of a few tens of parsecs and more, are both oriented north to northwest, and cannot provide a guide for disentangling the innermost geometry of this galaxy.

Here, we present new data from the VLTI spectrometer GRAVITY (GRAVITY Collaboration 2017) which combines the light from all four Unit Telescopes (UTs) of the VLT, and can resolve structures on scales smaller than a few milliarcseconds and provide spectro-astrometric relative precision of tens of microarseconds. This resolution is essential to resolve the BLR and the hot dust distribution immediately around it. Analysis of the H- and K-band variability by Lira et al. (2011) revealed a 70 day lag with respect to the optical continuum, indicative of a radius of 0.06 pc for the hot dust distribution. The size measured from K-band interferometry is very comparable at 0.1 pc (Weigelt et al. 2012; GRAVITY Collaboration 2020b). The new data we report here present a new look at this innermost region of NGC 3783.

This work adopts the following parameters for a ΛCDM cosmology: Ω_m = 0.308, Ω_Λ = 0.692, and H₀ = 67.8 km s⁻¹ Mpc⁻¹ (Planck Collaboration XIII 2016). Using this cosmology and our adopted distance of 38.5 Mpc, 1 pc subtends 5.36 mas on the sky and 1 μas corresponds to 0.22 light-days.

2. Observations and data reduction

2.1. GRAVITY

We observed NGC 3783 with GRAVITY (GRAVITY Collaboration 2017) over 6 nights spanning 3 years through a series of Open Time programmes and our Large Programme² aimed at measuring the size of the BLR and the mass of the supermassive black hole (SMBH). We used single-field on-axis mode with combined polarisation and medium (R ∼ 500) spectral resolution in the science channel for all observations. In single-field on-axis mode, both the fringe tracker (FT) and science channel (SC) fibres are centred on the same object with each receiving 50% of the light. A single exposure lasted 5 min and consists of (1) coherent 3.3 ms integrations with the FT fibre and (2) a series of ten exposures with detector on-chip integration times (DITs) of 30 s with the SC fibre. Exposures were done in sequence with only intermittent sky and calibrator star observations. Table 1 provides a log of the GRAVITY observations used for our analyses in this paper including the range of seeing and coherence time throughout the observations.

We reduced all data using the latest version of the GRAVITY pipeline (Lapeyrere et al. 2014). For the FT continuum visibility data we applied the default settings. For the SC data, we followed GRAVITY Collaboration (2020a, hereafter GC20a) and GRAVITY Collaboration (2020b, hereafter GC20b) and chose to retain all SC DITs regardless of FT signal-to-noise ratio (S/N) (snr-min-ft = 0) and estimated visibility loss (vfactor-min-sc = 0) which we found results in a substantial improvement in the phase noise.

After the pipeline reduction of the FT data, we further applied the selection method first presented in GC20b where only FT DITs with a group delay < 3 μm are selected before averaging together and then calibrating using the calibrator star observations. This improves the averaged visibility squared data and removes its dependence on the Strehl ratio.

For the SC data, we applied the additional steps outlined in GC20a to improve the SC differential phase spectra. This involves first fitting and subtracting an instrumental phase model that removes residual phase features introduced by the use of a third-order polynomial in the pipeline; this latter polynomial is insufficient to flatten the phase spectra below the 1° level. As a final step, we fit and subtract a local first-order polynomial around the Brγ line to further refine the flattening. Science channel differential amplitude spectra are produced by fitting and dividing out a second-order polynomial from the pipeline reduced visibility amplitude spectra. We estimate the uncertainty for the differential phase and amplitude spectra by calculating the root-mean-square (rms) in the line-free regions.

To still further improve the S/N of the SC data, we chose to bin and stack the data based on their location in the uv plane. We split the data into three angular uv bins such that each bin covers one-third of the full uv track for each baseline. We stack all SC data within a uv bin using the inverse of the squared rms uncertainty as the weight. Figure 1 shows the average uv coordinates for each bin with extensions in the radial direction due to the spectral coverage of GRAVITY. Appendix A shows the full uv-binned spectra, but we note that the differential phase spectra have had the ‘continuum phase’ signal subtracted (see Sect. 3.2).

Fig. 1. uv coverage of the binned GRAVITY observations. The radial span of each stripe is due to the wavelength range of the data. uv bins 1, 2, and 3 are ordered in clockwise direction for each baseline.

2.2. SINFONI

We observed NGC 3783 on April 20, 2019, using SINFONI with adaptive optics in service mode at the VLT (Eisenhauer et al. 2003; Bonnet et al. 2004). The observations were obtained at the 25 mas pix⁻¹ scale and the K grating (R ≃ 4000). A total of six exposures of 100 s were taken in a dithering object-sky-object sequence. Data were reduced with our standard pipeline and wavelength calibration scheme (SPRED, Schreiber et al. 2004; Abuter et al. 2006). Sky frames were subtracted from the images to correct for instrumental and atmospheric background. We then reconstructed 12.5 mas pix⁻¹ images with a spectrum at each pixel, applying flat-fielding, bad-pixel, distortion, and cosmic-ray hit corrections. The wavelength calibration was done using emission line lamps and is further tuned on the atmospheric OH lines in the raw frames. The individual data cubes were then combined, and a telluric correction was applied using molecfit (Smette et al. 2015; Kausch et al. 2015). Finally, to improve the S/N in the outer regions of the field of view (FOV) to derive emission line maps, in particular of the [Ca VIII] line, we smoothed the data cube with a six-pixel full-width at half-maximum (FWHM) Gaussian in the spatial directions, equal to the 75 mas FWHM point spread function (PSF) measured from a fit to the continuum³, and a two-pixel FWHM Gaussian in the spectral direction, matching the instrumental dispersion.

Figure 2 shows example spectra (left panel) from the smoothed SINFONI cube including (1) the full integrated spectrum, (2) a spectrum from a nuclear spaxel at the position of the green cross in the continuum image (right panel), and (3) a spectrum from an off-nuclear region indicated by the orange circle in the continuum image. We further show the expected locations of the [Si VI], H₂ (1–0) S(1), Brγ, and [Ca VIII] emission lines.

Fig. 2. Left: example spectra extracted from the SINFONI cube including the full integrated spectrum (blue), a spectrum from a nuclear spaxel (green), and an integrated spectrum from an off-nuclear region (orange). All spectra have been normalised to their median flux and slightly offset to improve visualisation. Right: K-band continuum image created by integrating the SINFONI cube between 2.25 and 2.31 μm. The green cross indicates the spaxel used for the nuclear spectrum in the right panel. The orange circle indicates the aperture used for the off-nuclear spectrum.

2.3. Normalised Brγ and [Ca VIII] profiles

Besides the interferometric observables (e.g., differential phase and amplitude, closure phases, etc.), an important component of our analysis is the emission line flux profile normalised to the continuum which is also measured by GRAVITY. Following GC20a, to construct the normalised line profiles, we only used observations where an early-type star (Feb. 16, 2019 and Mar. 08, 2020) was observed as a calibrator to correct for telluric features. Using an early-type star avoids the complicated stellar absorption features that are prominent in late-type stars above 2.3 μm where we expect the faint [Ca VIII] line. However, early-type stars have strong Brγ absorption which occurs in the blue wing of the broad Brγ profile of NGC 3783. Fortunately, these calibrators were also observed on the same night as IRAS 09149−6206, whose Brγ emission does not overlap with the Brγ absorption of the calibrator. Thus, we fit a line profile to the absorption of the calibrator in the IRAS 09149−6206 spectrum from each night and used these fits to remove the stellar Brγ feature from the NGC 3783 spectra. The spectra from each of the two nights were then averaged together, weighted by their uncertainties to produce the final GRAVITY flux spectrum. Figure 3 shows the final GRAVITY mean spectrum for both Brγ and [Ca VIII].

Fig. 3. Normalised Brγ and [Ca VIII] line profiles observed by GRAVITY (blue) and SINFONI (orange).

With higher spectral resolution and S/N, the SINFONI data also provide a high-quality K-band spectrum that can be compared with the GRAVITY spectrum. We extracted a nuclear spectrum using a circular aperture centred on the peak of the continuum and a diameter of 9 pixels that matches the PSF FWHM of the smoothed cube. The region around the Brγ and [Ca VIII] lines were then normalised to a local fit of the continuum. Figure 3 shows the normalised Brγ and [Ca VIII] profiles along with the GRAVITY profiles.

For both Brγ and [Ca VIII], the SINFONI spectra show slightly higher peak normalised fluxes. This can be explained by the SINFONI spectrum covering a much larger physical area compared to GRAVITY. Whereas the FOV of GRAVITY is ∼60 mas, the SINFONI spectrum is produced by integrating over a circle with a diameter of 112.5 mas. Narrow Brγ emission, which is largely contributing to the peak of the Brγ profile, and [Ca VIII] emission is expected to occur over a large range of size scales. The additional area covered by SINFONI then increases the relative flux of the narrow Brγ and [Ca VIII] lines because the continuum is much more compact.

Our choice of line profile depends on the specific analysis we want to perform. For our BLR modelling (Sect. 3.2), we choose to use the SINFONI profile because of its much higher spectral resolution. Small deviations from a smooth line profile can be important and signal BLR substructure. For our coronal line region (CLR) analysis, we use the GRAVITY profiles because they represent the correct line-to-continuum ratio that matches the interferometric data.

For the BLR analysis, we further have to account for the narrow Brγ component. To be as accurate as possible, we use the [O III]λ5007 Å line profile as a template for the narrow line. To produce the template, we fit the [O III] line from previous X-shooter spectra (Schnorr-Müller et al. 2016a) with four Gaussian components that were needed to accurately model the line. The velocities, widths, and relative amplitudes of these Gaussian components were then fixed to create the narrow line template. The template was broadened to account for the spectral resolution difference between X-shooter and SINFONI. This broadened template along with two Gaussian components to describe the broad component were fit to the SINFONI Brγ spectrum with only an overall scale factor and velocity shift for the narrow line template allowed to vary. Figure 4 shows the results of our line fitting and shows that this model reproduces the Brγ profile very well. The narrow component was then subtracted from the spectrum to produce the profile to be used for our BLR analysis.

Fig. 4. Normalised Brγ profile from SINFONI (blue) with the best-fit line profile model (orange). The model is the combination of a narrow line template determined by the [O III]λ5007 Å line from X-shooter spectra (green) and a broad component consisting of two Gaussians (red). The residuals (purple) show that our fit very closely matches the observed Brγ profile.

3. Structure of the BLR and the SMBH mass

The baseline-averaged differential visibility amplitude spectra plotted in Fig. 5 definitively show we have detected the BLR. The positive bumps (ΔV > 1) seen in the longest baselines (UT4−UT2, UT4−UT1, and UT3−UT1) over the Brγ line indicate a smaller BLR size compared to the continuum, as expected. While the differential phase spectra are noisier, we also observe positive bumps (Δϕ > 0) in the UT4−UT3 and UT4−UT2 baselines and an ‘S-shape’ signal in the UT4−UT1 and UT3−UT1 baselines as shown in Fig. 6. As seen for IRAS 09149−6206 (GC20a), a differential phase signal following the line profile (i.e. ‘continuum phase’) is produced from a difference in interferometric phase between the BLR and hot dust continuum. An S-shape signal on the other hand can be produced by ordered rotation such as that detected in 3C 273 (GRAVITY Collaboration 2018).

Fig. 5. Baseline averaged differential visibility amplitude spectra spanning the Brγ line (coloured points and grey error bars). The grey profile in each panel plots the same normalised Brγ spectrum from GRAVITY. The black curves represent the best-fit model from our CLR analysis in Sect. 5.

However, in contrast to IRAS 09149−6206, we are able to directly constrain the ‘continuum phase’ contribution through our image reconstruction of the hot dust continuum (see Sect. 4.1). To do this, we first simulate differential phase GRAVITY data for each observation based on our best image reconstruction by running it through our in-house built GRAVITY simulator. This produces the complex phase signal associated with the hot dust continuum structure. We then calculate the differential phase signal by multiplying the complex phase by −2πf_λ/(1 + f_λ) where f_λ is the line flux at wavelength λ relative to a continuum level of unity. These ‘continuum phase’ spectra are finally subtracted from the original differential phase spectra. The subtracted differential phase spectra are shown in Fig. A.1. In the following sections we use these spectra to measure photocentres and fit a dynamic BLR model.

3.1. Photocentre fitting

A first test of our qualitative assessment of the BLR structure is calculating the photocentres of the spectral channels where Brγ line emission dominates. Following our analysis of IRAS 09149−6206 (GC20a), we fit the continuum-phase-subtracted, uv-binned differential phase spectra with the following equation,

(1)

where u is the uv coordinate of the baseline, and x_BLR, λ is the BLR coordinate of the photocentre at wavelength λ with respect to the photocentre of the continuum.

We fit spectral channels with f_λ > 1.04 resulting in 13 model-independent photocentres that are shown in Fig. 7a. Interestingly, even though the continuum phase has already been subtracted, we still see a slight residual systematic offset of the photocentres from the origin. This is likely caused by the relatively large pixel scale (100 μas) in our image reconstruction. Any slight shift of the continuum from the centre on scales less than 1 pixel would not be captured in the reconstructed continuum phase. Indeed, the residual offset is much less than 100 μas. Despite the small systematic shift, the photocentres still show a general NW–SE velocity gradient which would cause the ‘S-shape’ signal in other baselines.

We further test the robustness of the velocity gradient by fitting for a single photocentre for all five blueshifted channels and all eight redshifted channels relative to the line centre of 2.1866 μm, our ‘2-pole’ model. The blue and redshifted ‘poles’ are shown in Fig. 7b as blue and red dots respectively. As in the individual photocentre fitting, we still find the general velocity gradient in the same direction. We use an F-test to estimate the significance of the gradient by comparing to a null hypothesis where all of the spectral channels are located at a single photocentre (i.e. ‘null’ model). The best-fit photocentre of the null model is shown as a black dot in Fig. 7b. Comparing the 2-pole and null models, we find the detection of the velocity gradient is significant at > 8σ. Using the separation of the blue and red poles, we can also place a rough estimate on the size of the BLR with R_BLR ∼ 30 μas, consistent with the approximately ten light-day time-lags measured from RM (Stirpe et al. 1994; Onken & Peterson 2002; Peterson et al. 2004; Zu et al. 2011) However, to truly measure the size of the BLR and its other properties we need to fit the data with a physical model for the BLR.

3.2. Modelling the broad line region

Given the high significance of the velocity gradient and to constrain the properties of the BLR, we choose to fit our data with a physical model of the BLR. We use the model first developed in Pancoast et al. (2014) and updated for GRAVITY data in GRAVITY Collaboration (2018) and GC20a. As GC20a contains a detailed description of the BLR model, we only provide here a brief description necessary for our analysis. The underlying assumption of the model is that the BLR is composed of a large number of non-interacting clouds under the gravitational influence of a central SMBH with mass, M_BH. The radial distribution of the clouds is described by a shifted gamma distribution with a hard lower limit of the Schwarzschild radius and free parameters; β, controlling the shape of the distribution; F, the fractional inner radius; and R_BLR the mean BLR radius. Clouds are then distributed randomly both around the rotation axis and above the midplane up to the angular thickness of the disc, θ₀.

Three different parameters (κ, γ, ξ) introduce asymmetry in the BLR emission: κ controls the fraction of line emission emitted from individual clouds towards or away from the central source and hence the fractional emission in the observer’s direction from the near and far sides of the BLR; γ controls whether or not the outer surface of the BLR preferentially emits line emission; and ξ controls the level of midplane transparency. Together, these three parameters control the relative weight of each cloud in determining the total emission.

Each cloud is placed on a bound elliptical orbit to model the kinematics. The tangential (v_ϕ) and radial (v_r) velocities are randomly chosen based on a distribution centred around v_ϕ = v_circ and v_r = 0 with v_circ equal to the Keplerian circular velocity given the radial position of the cloud and M_BH of the SMBH. The distribution follows an ellipse described by Eq. (6) of GC20a and has a Gaussian shape with standard deviations, σ_Θ, circ along the ellipse, and σ_ρ, circ perpendicular to the ellipse. A fraction of the clouds (1 − f_ellip) can be placed on highly elongated orbits dominated by radial motion. Whether these clouds are inflowing or outflowing is controlled by f_flow where f_flow < 0.5 indicates inflow and f_flow > 0.5 indicates outflow.

The cloud distribution is then rotated on the sky by an inclination angle, i, and position angle, PA, and translated by an offset, (x₀, y₀). Line-of-sight velocities are calculated including both the full relativistic Doppler effect and gravitational redshift. Finally, clouds are binned into spectral channels according to their LOS velocity. The flux for each spectral channel is the sum of the weights for all clouds within the bin and the photocentre is the weighted average position on the sky of all the clouds within the bin. The model flux profile is normalised such that the maximum is f_peak and differential phases are then calculated according to Eq. (1). Both the model flux profile and differential phase spectra are compared to the observed flux profile and differential phase spectra and the model parameter posterior distributions are sampled using nested sampling with the dynesty Python package (Speagle 2020). Priors on each parameter are the same as given in GC20a.

Table 2 lists the best-fit parameters and uncertainties from our BLR modelling and Fig. A.2 shows the posterior distributions. Best-fit parameter values are determined from the maximum a posteriori values of the joint posterior distribution. Uncertainties on these values represent the 95% credible interval estimated from the marginalised posterior distributions.

Our best-fit BLR model well reproduces both the differential phase and Brγ line profile with a reduced chi-square, . Figure 6 shows the combined best-fit BLR model added to the continuum phase signal for the baseline-averaged differential phase spectra (see Fig. A.1 for the comparison to the uv-binned data). A representation of the on-sky cloud distribution is plotted in Figs. 7c and 8a shows the corresponding model photocentres which agree well with the orientation and size of the observed photocentres. We note that the PA of the model photocentres does not match the best-fit PA. This is due to the moderate fraction of inflowing clouds which twists the PA on the sky away from the PA of the rotating disc.

Fig. 6. Baseline averaged differential phase spectra spanning the Brγ line (coloured points and grey error bars). The grey profile in each panel plots the same normalised Brγ spectrum from SINFONI with the narrow component removed. The black curves represent the best-fit BLR model summed with the continuum phase calculated from our image reconstruction. The model line profile appears noisy because of the continuum phase contribution being calculated based on the observed line profile.

Fig. 7. Panel a: best-fit photocentres for the 13 brightest spectral channels covering the Brγ line to the differential phase data. The colour indicates wavelength and the ellipses outline the uncertainty on the photocentres. The red cross plots the (0,0) position. Panel b: best fit ‘2-pole’ photocentre model where all blueshifted (blue point) and redshifted (red point) channels are assumed to have the same photocentre with the addition of the wavelength-independent offset. The black point indicates the best-fit position of the ‘null’ model where all channels are assumed to have the same photocentre. Panel c: photocentres produced from our best-fitting BLR model (see Sect. 3.2) which we calculated by first producing mock differential phase spectra before fitting them in the same way as panel a. Colours are the same as in panel a.

Fig. 8. Panel a: cloud distribution of the best-fit BLR model. Each circle represents one cloud colour-coded by the LOS velocity. The green ellipse at the origin illustrates the uncertainty of the offset of the BLR centre. Panel b: observed averaged differential phase from baselines UT4−UT2, UT4−UT1, and UT3−UT1 after removing the residual ‘continuum phase’ signal (blue points) compared to the averaged differential phase from the best-fit BLR model (lower red line). The black points and overlaid red line show the observed and best model line flux profile respectively. Panel c: averaged differential phase data (points) and the best-fit models (lines) of the three baselines that show the strongest signal of the BLR component (dashed lines). The phase in panel b is calculated by averaging the phases of these three baselines after subtracting the best-fit residual continuum phases (dotted lines).

With a best-fit inclination of i = 23° and θ₀ = 24°, we find a relatively face-on and moderately thick disc describes the BLR of NGC 3783 well, as expected for a type 1 AGN and similar to both 3C 273 (GRAVITY Collaboration 2018) and IRAS 09149−6206 GC20a. The centre of the BLR is offset from the centre of the hot dust by (−0.5, −19) μas similar to the offset in the photocentres. In Fig. 8b we show the intrinsic BLR differential phase signal compared to the data by averaging the continuum-phase-subtracted differential phase spectra from the longest three baselines shown in Fig. 8c. Both the model and data show the characteristic ‘S-shape’ expected of a BLR with kinematics dominated by Keplerian rotation. Figure 8c also shows the decomposition of the model differential phases into the BLR and continuum components. Because we had already removed most of the continuum phase signal, the primary component is the BLR component.

We find very little anisotropy in the emission of the BLR clouds. κ is consistent with zero, indicating that neither the near nor far side of the clouds is preferentially emitting. ξ is ∼0.7, suggesting little mid-plane obscuration. Finally, γ is close to 1.0, indicating the clouds on the outer surface of the BLR are not radiating more than the inner clouds.

The motion of the clouds is split relatively evenly between bound elliptical orbits and radial motion with f_ellip ∼ 0.46. The clouds with radial motion are inflowing with P_inflow = 0.65⁴. Our finding of a significant fraction of inflowing clouds agrees with the simple velocity-resolved time-lag measurements from Bentz et al. (2021) which showed slightly shorter lags in the red wing of the line profile.

An interesting feature of our best-fit model is the radial distribution of the clouds. With β = 1.7, which is moderately high, the distribution is more sharply peaked near R_min with a tail towards larger radii compared to 3C 273 and IRAS 09149−6206 which each had a β ∼ 1.2. The best-fit distribution produces a mean radius of R_BLR = 71 μas, which corresponds to 16 light-days and a minimum radius of R_min = 32 μas corresponding to 7 light-days. In our model, there is no maximum radius at which a cloud can exist.

We can ask whether these various radii make physical sense. Reverberation mapping time-lags for high ionisation lines such as He II, C IV, and Si IV have been measured between 1 and 4 days (Reichert et al. 1994; Onken & Peterson 2002; Bentz et al. 2021) indicating gas at smaller radii than the R_min measured here. At large radii, both observations and models suggest a maximum radius of the BLR at the dust sublimation radius (Laor & Draine 1993; Netzer & Laor 1993; Korista et al. 1997; Baskin et al. 2014; Schnorr-Müller et al. 2016b; Suganuma et al. 2006) which for NGC 3783 should occur around 80 light-days. Our best-fit model only contains < 1% of clouds at radii larger than this, and so indeed for NGC 3783 it does seem the model BLR is confined within the dust sublimation radius.

As a test of the sensitivity of the model to the inner and outer radii, we further decided to fit the data with a model with fixed R_min = 4 light-days and R_max = 80 light-days. We find that by fixing these radii we can still achieve a fit that is almost as good as our fiducial one. The main change to the BLR model is that β decreases to 1.2, κ decreases to −0.3, and f_ellip increases to 0.83. All other parameters remain unchanged, including M_BH. This highlights both the extreme flexibility of the BLR model and in particular the choice of the shifted gamma distribution to describe the radial distribution of the clouds as well as the insensitivity of the current data set to R_min. We will explore this issue in more detail in upcoming publications that will also include a joint RM and GRAVITY analysis of NGC 3783, which will help to constrain these parameters. For the rest of this paper, we remain with the BLR parameters as presented in Table 2.

In the following section, we address and compare our mean radius with the time-lags measured through RM and discuss the SMBH mass and the position of NGC 3783 on the R–L relation.

3.3. Time-lags, the R–L relation, and SMBH mass

Our BLR mean radius is a factor of 1.6 larger than the radius measured through RM. However, in order to compare the model size from interferometry with RM, we need to take into account that both techniques are sensitive to different parts of the cloud distribution in the BLR. Reverberation mapping determines a response-weighted mean radius while the interferometry modelling measures a deprojected brightness-weighted size of the distribution, and so differences between the two are not unexpected.

To test this, we generated a mock emission-line light curve based on the real continuum light curve from Bentz et al. (2021). From our best-fit BLR model, we calculate a transfer function based on the distribution of time-lags associated with each BLR cloud. This was then used to convolve the continuum light curve and produce an emission-line light curve. Figure 9a shows the input continuum light curve and Fig. 9b shows the emission-line light curve corresponding to our best-fit BLR model. Following reverberation studies, we finally calculated the time-lag between the two light curves using the cross correlation function (CCF). The results of the CCF are shown in Fig. 9c.

Fig. 9. Panel a: observed continuum light curve from Bentz et al. (2021) (black points) with the finer sampled interpolated continuum light curve over-plotted (grey line). Panel b: observed broad emission-line light curve (black points) with the emission-line light curve produced by reverberating the input continuum light curve off the clouds in our best-fit BLR model (red). Panel c: cross-correlation function between the continuum and broad emission-line light curves as a function of time-lag. The dashed line indicates the peak cross-correlation time-lag of 10.2 days.

We find a peak time-lag of days, very consistent with the time-lags measured in previous studies. This suggests that, at least for NGC 3783, there is a modest difference between the physical mean radius of the BLR and the observed time-lag. In Fig. 10, we place NGC 3783 on the relation measured by Bentz et al. (2013) using a 5100 Å continuum luminosity, log λL₅₁₀₀ = 42.93 erg s⁻¹, from Bentz et al. (2021), which is adjusted for our luminosity distance of 38.5 Mpc and both R_BLR determined from our modelling (filled star) and the measured time-lag from the literature (open star). While the 5100 Å continuum luminosity from Bentz et al. (2021) was measured at about the same time as our last set of GRAVITY observations in March 2020, our full set of observations span several years. Long-term monitoring of NGC 3783 shows that the AGN can vary by a factor of two over several years (Lira et al. 2011). Therefore, in Fig. 10 we include an uncertainty of 0.3 dex in the 5100 Å continuum luminosity. GRAVITY-based results for IRAS 09149−6206 and 3C 273 are also plotted along with two RM samples from Du & Wang (2019) and Grier et al. (2017). We note here that both IRAS 09149−6206 and 3C 273 do not show the same difference between the interferometrically and RM measured size. This is due to the lower β values for their cloud distributions and therefore we only plot one point for them in Fig. 10.

Fig. 10. Relationship between BLR radius and 5100 Å AGN luminosity with the best-fit line from Bentz et al. (2013) (dashed line) and the RM samples from Du & Wang (2019) (black points) and Grier et al. (2017) (grey points). Previous GRAVITY BLR radius measurements of IRAS 09149−6206 and 3C 273 are shown along with the BLR mean radius for NGC 3783 as filled red stars. The time-lag measurement of ∼10 days for NGC 3783 from the literature is shown as an open star. The horizontal error bar indicates a factor of two uncertainty on the luminosity of NGC 3783 based on variability.

As seen before, the literature time-lag for NGC 3783 lies perfectly on the R−L relation (e.g., Peterson et al. 2004; Bentz et al. 2013). However, the physical mean radius is above it by 0.2 dex. We note that the discrepancy between NGC 3783 and the R−L relation is different from the increased scatter that both Du & Wang (2019) and Grier et al. (2017) have reported. Du & Wang (2019) specifically showed that shorter time-lags are observed in AGNs with high Eddington ratios while Grier et al. (2017) point to selection effects or a change in the BLR structure at higher luminosities rather than accretion rate. NGC 3783, on the other hand, is at lower luminosity and lower Eddington ratio (∼0.1 based on the modelled black hole mass and log L_bol = 44.52 from GC20b) and we see a consistent time-lag with the R−L relation and rather find an offset mean radius. This highlights the fact that the R−L relation is fundamentally a time-lag luminosity relation and the conversion to a radius from a time-lag relies on simple assumptions about the geometry and structure of the BLR. White & Peterson (1994) predicted that if the true BLR cloud distribution is significantly extended, RM time-lags could be biased towards the inner radius.

The discrepancy we observe could also be related to the same trend we observed with the dust sizes. GC20b showed that less luminous AGNs have larger interferometric dust sizes compared to RM-measured dust sizes. It is possible that both the BLR and hot dust geometry are changing in similar ways with increasing AGN luminosity. In particular, this could be related to the recent result that the obscuration fraction of AGNs significantly decreases with increasing Eddington ratio due to the AGNs driving gas and dust out of the central few parsecs (Ricci et al. 2017). If the gas and dust are able to accumulate in the inner regions at low luminosity this would lead to the discrepancies we observe, while at higher luminosity, gas and dust are driven out primarily from the inner regions and reduce the difference between RM and physical radii. This is certainly speculative and based primarily on two objects. We need more interferometric BLR measurements to test this explanation.

Our BLR modelling also constrains the SMBH mass with log M_BH = 7.68. Interestingly, this is quite consistent with the masses based on the time-lag and using a virial factor, f ≈ 4−5, which is needed to match the M − σ relation. Thus, the inferred virial factor for NGC 3783 is consistent with the average value used in RM studies.

4. Structure of the hot dust continuum

Beyond the BLR, we can also study, using the FT data, the surrounding hot dust structure which produces the underlying near-IR continuum. Figure 11 shows a general trend of decreasing V² with increasing baseline length which indicates a partially resolved primary source of the hot dust continuum. GC20b used this data to measure a Gaussian FWHM size of 0.82 mas. However, GC20b also noted that NGC 3783 was the only AGN within their sample to show strong signatures of asymmetry evidenced by non-zero closure phases. The right panel of Fig. 11 plots the closure phases which show consistently negative values between −2° and −10°. Symmetric structures such as a circular Gaussian cannot produce non-zero closure phases and therefore higher-order structures must be present. In the following sections, we investigate the nature of the hot dust continuum using two methods: (1) model-independent image reconstruction and (2) uv plane model fitting.

Fig. 11. GRAVITY FT squared visibilities (left) and closure phases (right) as a function of spatial frequency for all nights of observation. The significant non-zero closure phases indicate resolved asymmetric hot dust structure in NGC 3783.

4.1. Continuum image reconstruction

The high quality and large uv coverage of the FT data allows us to use image reconstruction codes to spatially map the hot dust. As we did for our imaging of NGC 1068 (GRAVITY Collaboration 2020c), we used MiRA⁵ (Multi-aperture Image Reconstruction Algorithm Thiébaut 2008) to reconstruct the K-band image of NGC 3783 based on the V² and closure phase data as shown in Fig. 11. The V² uncertainties shown have been multiplied by ten to match the S/N of the closure phases and avoid overweighting V². This also takes into account calibration uncertainty that is not included in the pipeline-generated data.

In general, image reconstruction from interferometric data is an ill-posed problem due to sparse coverage of the uv plane, especially for optical/near-IR interferometry. Therefore, to reduce the number of possible solutions (i.e. images) that can equally reproduce the data, image reconstruction codes use priors to constrain the brightness distribution. Within MiRA we imposed both a positivity prior (i.e. flux must be ≥0) and the hyperbolic regularisation which is an edge-preserving smoothness prior. The hyperbolic regularisation favours solutions where the flux is smooth inside the structure but contains sharp edges. Two hyperparameters can be tuned with this regularisation, μ and τ. μ is simply the weight given to the regularisation in determining the total likelihood, L = f_data + μf_prior, where f_data is a function comparing the model to the data and f_prior is a function comparing the model to the constraints given by the prior.

τ is a hyperparameter specifically associated with the hyperbolic regularisation and is an edge threshold that controls how sharp the edges are expected to be. Overly small values of τ will produce a cartoon-like image with very smooth regions that are then sharply cut off. Large values of τ instead lead to a compactness regularisation which favours a single centrally concentrated source.

Thiébaut & Young (2017) outline best practices for choosing optimum values of μ and τ. Following these, we ran MiRA over a grid of values for μ and τ. For each value of τ, we chose the value of μ that corresponds to the elbow of the ‘L-curve’ which is a plot of f_data against f_prior. Choosing μ at the elbow is a compromise between over- and under-regularisation. We then visually inspected all of the images associated with each μ-optimised τ value and chose the image that avoided the cartoon-like effects but also many spurious compact sources. We found optimum values of μ and τ of 5 × 10⁵ and 10⁻³ respectively.

Figure 12 shows our final image reconstruction at a pixel scale of 0.1 mas, a FOV of 25.6 mas, and an initial image of a Dirac delta function. As expected, the hot dust image is dominated by a central, marginally resolved source, however immediately noticeable is the presence of a fainter smaller source to the SW. Simple Gaussian fitting directly on the image finds a size of 1.8 × 1.2 mas and PA ∼ −44° (east of north) for the bright central source and 1.5 × 0.8 mas and PA ∼ −32° for the fainter SW source (hereafter referred to as the offset cloud). The offset cloud is 3.3 mas (0.6 pc) away from the centre at a PA of ∼ − 96° and contains 5% of the flux.

Fig. 12. MiRA image reconstruction of the K-band continuum for NGC 3783. We used the hyperbolic regularisation with a FOV of 25.6 mas and a pixel size of 0.1 mas. The image is normalised to the maximum pixel flux. The GRAVITY interferometric beam is shown as a white ellipse in the bottom left corner for comparison. North is up and east is to the left.

The PA and extent of the individual sources follow the PA of the GRAVITY beam, and therefore it is possible these properties are an artefact of the reconstruction and rather reflect the limited uv coverage of our data. Importantly, the distance and direction of the offset source places it well outside the beam of the central source and strengthens the reliability of its detection. In Appendix B, we show that the offset source is robust against the choice of regularisation, choice of image reconstruction algorithm, and random removal of data, while the other much fainter sources are likely spurious.

The bottom panels of Fig. 11 show the V² and closure phases of our image reconstruction. The closure phases match the observed ones quite well and accurately reproduce all of the key features including the rise to 0° towards smaller spatial frequencies and the vertical streaks prominent in the UT3−UT2−UT1 triangle. V², however, is only moderately well matched. In particular, the gradient of V² seems steeper for the image reconstruction which leads to a larger primary source. This leads to the factor of approximately two difference between the image reconstructed size and the size found in GC20b (0.82 mas) from fitting a Gaussian model to only V² data. Therefore, as a further test of the detection and specific properties of each source, in the following section we apply the same model fitting.

4.2. Visibility model fitting

We fit the individual V² and closure phase data of each night with a model composed of a central 2D Gaussian, an offset point source, and an unresolved background. As for the image reconstruction, we inflated the V² errors by a factor of ten to match the S/N of the closure phases and partially account for calibration uncertainty.

Table 3 lists the best-fit parameters for each night and Fig. 14 shows an example fit for January 7, 2018. In all nights, we find good qualitative agreement with the features of the data, and in particular, the offset point source provides a good match to the observed non-zero closure phases as was found in the image reconstruction. Between nights, we also find very good consistency in the best-fit results with an average central source FWHM of 0.72 ± 0.1 mas, and an offset point source at (−2.9 ± 0.2, −1.0 ± 0.2) mas relative to the central Gaussian, emitting an average fractional flux of 0.043 ± 0.01. This places the offset point source 3.1 mas (0.57 pc) away from the central Gaussian at a PA of −109°. All values except for the size of the central Gaussian are in excellent agreement with the features seen in our reconstructed image.

Fig. 13. Same as Fig. 12 but with lines indicating important angles and regions. Blue lines outline the ionisation cone as produced by the BLR polar axis (red line) and thickness. The yellow line indicates the position angle from polarimetery (Lira et al. 2020) and the green line indicates the position angle from mid-IR interferometry (Hönig et al. 2013).

Fig. 14. Sample 2D Gaussian plus point-source model fit to the observed V2 and closure phase data for NGC 3783 from Jan 7, 2018. The best-fitting Gaussian FWHM geometric mean size of ≃0.6 mas is consistent with that found from 1D Gaussian fitting (GC20b). The Gaussian PA and point source flux fraction and offset are consistent with those found from image reconstruction.

While our model fitting agrees well with GC20b, we suspect the choice of regularisation is causing the much larger size in the reconstructed image. Indeed, both increasing τ and switching to a compactness regularisation reduced the size of the central Gaussian to values similar to those in Table 3. This illustrates how sensitive image reconstruction can be to the choice of regularisation, especially as we are at sizes well below the diffraction limit, Δθ ≪ λ/B. For this reason, we choose to use the average Gaussian FWHM found through our model fitting as the size of the main hot dust continuum source for this paper. Both methods, image reconstruction and visibility model fitting, robustly detect the offset cloud and we discuss possible origins in the following section.

4.3. Origin of the offset hot dust source

We explore two possibilities for the origin of the offset dust cloud: (1) an orbiting secondary SMBH and (2) a dust cloud heated by the central AGN. In both cases, we use a projected radial distance from the primary AGN of 0.6 pc and a K-band flux density of 3.3 mJy. This flux density was calculated starting with a total K-band fibre magnitude of 10 mag (GC20b) and using the measured ∼5% fractional flux of the offset cloud.

In scenario 1, we could be observing a secondary SMBH orbiting the primary AGN. The hot dust emission would then suggest that the secondary SMBH is also a faint AGN that is heating the dust around it and would also contain its own BLR. The secondary would then have a near-IR luminosity of 10⁴² erg s⁻¹ and a bolometric luminosity of 10^42.7 erg s⁻¹ using the near-IR–X-ray relation from Burtscher et al. (2015) and converting to a bolometric luminosity using the relation in Winter et al. (2012). Assuming the secondary AGN is at the Eddington luminosity places a lower limit on the black hole mass of 4 × 10⁴ M_⊙, and assuming it is accreting at 10% Eddington, similar to the primary, gives M_BH = 4 × 10⁵ M_⊙. As a first-order calculation of the expected orbital velocity of the secondary, we can simply use with M = 10⁸ M_⊙ found from our BLR modelling because the secondary must be much less massive than the primary. At a radius of 0.6 pc, we would then expect v ≈ 900 km s⁻¹ which would induce a shift in the secondary BLR and produce a double peak in the combined spectrum, albeit at only 5% the strength of the primary broad emission lines and assuming there is little obscuration towards the secondary BLR. This velocity is also the maximum velocity we would expect given a fully inclined orbit and the secondary currently at its maximum projected distance from the primary. We attempted to fit the normalised Brγ profiles with the addition of a second broad emission component but were not able to achieve a reasonable fit. Another potential signal of a binary SMBH would be a periodic light curve, however given the distance of the secondary and the mass of primary, the expected period of the orbit is 4400 years. Thus, based on the current data on NGC 3783, we cannot fully rule out a SMBH binary as the origin of the offset cloud.

A simpler explanation is scenario 2, where the offset cloud is a massive cloud of gas either inflowing towards the AGN from the circumnuclear disc or potentially outflowing away from it. The dust emission would then be heated, not internally, but externally by the central AGN. Indeed, as shown in Fig. 13, the cloud does lie near the edge of the ionisation cone that is produced given the geometry of our best-fit BLR, which further matches the geometry inferred from polarimetry (Lira et al. 2020) and mid-IR interferometry (Hönig et al. 2013). This would give the cloud a direct view towards the AGN and its radiation.

To test this scenario, we can use the observed luminosity of the AGN to constrain the dust temperature and measure a dust mass under the assumption that it is emitting as a modified blackbody. If a non-physical amount of dust would be needed to produce the near-IR emission observed then we can safely rule out this scenario. We calculate the dust temperature using the simple scaling relation that the temperature decreases as a power law with a maximum dust temperature, T_sub, at the dust sublimation radius, r_sub, (e.g., Hönig & Kishimoto 2010),

(2)

Both α and r_sub depend on the specific grain distribution. We use the ‘ISM large grains’ model from Hönig & Kishimoto (2010) which consists of 47% graphites and 53% silicates and a Mathis et al. (1977) size distribution between 0.1 μm and 1 μm. For this dust model, r_sub = 0.5 pc for T_sub = 1500 K and an AGN bolometric luminosity of 10⁴⁶ erg s⁻¹ and α = −2.1. We adjust for the luminosity of the AGN in NGC 3783 (log L_bol = 44.5) by scaling with and find r_sub = 0.09 pc, consistent with the radius measured for the bright component (0.07 pc) in our visibility model fitting. At the location of the offset cloud, we then calculate a dust temperature of 604 K. We note that other grain models, such as those based on pure graphite grains in the innermost part of the obscuring structure, give similar r_sub.

With only one SED data point, we use a simple modified blackbody model to calculate the dust mass:

(3)

where M_d is the dust mass, D_L is the luminosity distance, c is the speed of light, h is the Planck constant, and k is the Boltzmann constant. We use κ₀ = 167 m² kg⁻¹, β = 1.25, and ν₀ = 137 THz (2.19 μm) from Draine (2003). With T_d = 604 K from Eq. (2) and S_K = 3.3 mJy, we estimate M_d = 0.2 M_⊙. Assuming a standard dust-to-gas ratio of 100, the gas mass is 20 M_⊙. This cloud mass does not seem completely unreasonable, especially on sub-parsec scales around an AGN. Assuming a radius of ∼0.2 pc for the cloud directly measured from the image, the density of the cloud would be ∼10⁴ cm⁻³ which is at the high end of the observed densities in the narrow line region (NLR) of local AGNs (Davies et al. 2020). Therefore, while we cannot specifically dismiss the binary SMBH origin of the offset cloud, we prefer the primary AGN heated dust scenario for simplicity.

MATISSE, the new mid-IR instrument at the VLTI could potentially help distinguish between our two proposed explanations. While we showed that we expect a dust temperature T_d ∼ 600 K if heated by a single central AGN, if instead it is heated by a secondary it should have a more complex SED, likely with an inner hot dust component at ∼1400 K and cooler dust at further radii (e.g., Kishimoto et al. 2011). Using the simple modified blackbody, we estimate flux densities of 26, 45, and 28 mJy for the L, M, and N bands accessible by MATISSE. Based on the exposure time calculator, closure phase uncertainties of 1° should be achievable with a few hours of observations in the L and M bands.

5. Nuclear size of the coronal line region

Our final analysis from this rich data set involves, for the first time, measuring the size of the coronal line region (CLR) on nuclear scales. In addition to the BLR and NLR, a subset of AGNs also show lines from highly ionised atoms, the so-called coronal lines, for their first observation in the solar corona. These lines all have ionisation potentials ≥100 eV and are collisionally excited forbidden transitions with relatively high critical densities (10⁷ − 10¹⁰ cm⁻³). Because of the high critical densities and observations that many coronal lines have FWHM intermediate between BLR and NLR lines (e.g., Appenzeller & Wagner 1991; Veilleux 1991; Rodríguez-Ardila et al. 2006, 2011; Lamperti et al. 2017) it is thought that the CLR lies between the BLR and NLR. A distinct increase in line FWHM with ionisation potential for some AGNs is further evidence of its intermediate location (e.g., Wilson 1979; Penston et al. 1984; Thompson 1995; Rodríguez-Ardila et al. 2002, 2011) under the assumption of photoionisation as the primary ionisation mechanism (e.g., Penston et al. 1984; Ferguson et al. 1997; Mazzalay et al. 2010; Rodríguez-Ardila et al. 2011). This has also led to the hypothesis that the CLR originates at the inner region of the central obscuring structure (Pier & Voit 1995; Ferguson et al. 1997; Murayama & Taniguchi 1998). Nevertheless, arguments against this hypothesis include the fact that coronal lines are observed equally in type 1 and type 2 AGNs (e.g., Rodríguez-Ardila et al. 2011) and that significant coronal line emission is seen in even the most heavily Compton-thick AGNs like Circinus, Centaurus A, and NGC 1068 (e.g., Moorwood et al. 1996; Reunanen et al. 2003; Rodríguez-Ardila et al. 2006; Mazzalay et al. 2010; Müller-Sánchez et al. 2011).

Further, through long-slit and integral field spectroscopic observations, coronal line emission has also been found to be extended on ∼10−100 pc scales (e.g., Prieto et al. 2005; Mazzalay et al. 2010; Müller-Sánchez et al. 2011), but with the brightest emission still concentrated in the unresolved nuclear region. This more extended and more diffuse emission instead is likely to be powered by shocks that could be produced as radio jets propagate through the surrounding gas (e.g., Rodríguez-Ardila et al. 2002, 2006, 2017a; Reunanen et al. 2003; Müller-Sánchez et al. 2011). GRAVITY provides the first opportunity to measure the size of the nuclear CLR through the [Ca VIII] line.

The differential visibility amplitude spectra spanning the [Ca VIII] line (Fig. 15) show a dip near the peak of the emission line in all the baselines, in contrast to the Brγ spectra which show primarily a peak and indicate a smaller BLR size compared to the hot dust. A dip instead is caused by a larger emitting line region compared to the hot dust. We also observe a dip in differential visibility amplitude superimposed on the broad peak seen in the Brγ region, and the location of this dip further corresponds to the narrow peak in the flux profile, which suggests that the narrow Brγ line-emitting region is also larger than the hot dust continuum. The narrow Brγ line then could be tracing the base of the NLR or also originate in the CLR.

Fig. 15. Same as Fig. 5 but for the [Ca VIII] region.

The differential visibility amplitude across an emission line is described by the following equation,

(4)

where f is the flux line profile normalised to a continuum of 1, V_line is the visibility amplitude of the line, and V_c is the visibility amplitude of the continuum. However, Brγ is the combination of a broad and narrow line, each with a different line width and physical size relative to the continuum. The differential visibility amplitude for the Brγ line then follows

(5)

where f_t, f_b, and f_n are the total, broad, and narrow flux profiles normalised to a continuum of 1, and V_c, V_b, and V_n are the continuum, broad, and narrow visibility amplitudes.

To fit both Eqs. (4) and (5) we model all emitting regions as Gaussian sources. In the marginally resolved limit, the visibility of a Gaussian follows

(6)

where V₀ is the zero baseline visibility, r_uv is the baseline length in units of mas⁻¹, and FWHM is the size of the source in mas. While a single compact source should always have V₀ = 1, GC20b showed this is not the case for our AGN and can be caused either by coherence loss or unresolved background emission. We therefore choose to include this in our modelling. We further set V_b = 1 because the BLR is effectively unresolved. The measured size of 101 μas from spectroastrometry (see Sect. 3.2) results in V_b = 0.995 at the longest baseline observed.

Equations (4) and (5) also show that the differential visibility amplitude depends on the normalised flux profile. While for the [Ca VIII] line we were able to simply use the observed line profile, for Brγ we need to decompose the line into its narrow and broad components. To fold in the uncertainties related to this decomposition, we include in our model the shape of the narrow and broad components, parameterised as Gaussian lines with an amplitude, central wavelength, and FWHM. We further also fit for the FWHM Gaussian size of the narrow Brγ-emitting region. In total our model has 14 free parameters to fully describe the sizes of the CLR, hot dust, and narrow Brγ regions; the differential visibility amplitude spectra across the Brγ and [Ca VIII] lines; and the flux line profiles of Brγ and [Ca VIII].

Table 4 lists the fitted parameters and their best-fit values and uncertainties that were estimated using the median and 95% credible interval of their respective posterior distributions (see Fig. A.5). From this analysis, we measure a hot dust size of 0.73 mas corresponding to 0.13 pc and a CLR size of 2.2 mas corresponding to 0.4 pc. The hot dust size very closely matches the size measured from visibility model fitting of the FT data (see Sect. 4).

Our size firmly places the nuclear [Ca VIII]-emitting clouds beyond both the BLR and the near-IR-emitting hot dust. For NGC 3783, this argues against an origin of [Ca VIII] in the inner region of the central obscuring structure which has been suggested as the site of the CLR (Pier & Voit 1995; Ferguson et al. 1997; Murayama & Taniguchi 1998). However, this does not rule out an inner region origin for all coronal lines, especially as [Ca VIII] has one of the lowest ionisation potentials of the observed coronal lines (IP = 128 eV). Indeed, Mullaney et al. (2009) were able to model the flux and kinematics of multiple coronal lines with an AGN-driven outflow launched from the inner edge of a dusty obscuring structure. The high-velocity components of the highest ionisation lines all are produced near the dust sublimation radius while the lower ionisation lines are produced at larger radii. This could also explain the small velocity shift of the [Ca VIII] line and relatively low FWHM, because at larger radii the gravitational potential of the stellar bulge should reduce the velocity of the clouds.

The [Ca VIII] flux⁶ of 4 × 10⁻¹⁸ W m⁻² observed from a region of 0.92 mas in radius can be used to derive a crude estimate for the density of this inner part of the CLR. We adopt an estimated filling factor by [Ca VIII]-emitting gas of f = 0.1 in order to exclude regions either devoid of gas or with Ca in other ionisation stages. Using atomic data of Landi et al. (2004) and Saraph & Storey (1996) for a simple two-level analysis, and a Ca abundance of 2.e−6 by number (Landi et al. 2004), n_e ≈ 2 × 10⁵ cm⁻³ is needed to reproduce the observed [Ca VIII] flux. This is consistent with the notion that coronal lines arise in an inner and dense part of the NLR; see for example Davies et al. (2020) for a recent assessment of NLR densities. The required CLR density scales with f^−0.5, that is, the electron density in the [Ca VIII]- emitting region must be at least about 10⁵ cm⁻³ but could be clearly higher if arising in very low-filling-factor clouds or filaments.

6. Connecting the nuclear and circumnuclear regions

With these new GRAVITY observations of NGC 3783, we detected and measured the properties of three distinct components within a radius of 1 pc from the AGN:

1. A rotating BLR with a mean radius of 0.013 pc, an inclination of 23°, and PA of 295°.

2. A hot dust structure composed of a central bright source of 0.14 pc in size and a faint offset cloud 0.6 pc away at a PA of −109°.

3. A nuclear CLR with a size of 0.4 pc.

In this section, we aim to place these components in the context of the larger scale circumnuclear environment that has been well studied with previous observations. To help in this, we fit the [Si VI] and the ro-vibrational H₂ (1–0) S(1) lines detected in our SINFONI cube. We chose to also fit [Si VI] because it is a brighter coronal line than [Ca VIII] and has a similar ionisation potential (IP = 167 eV) and thus allows the ionised gas to be traced out to larger scales. The H₂ line traces hot molecular gas that is likely inflowing and feeding the AGN (e.g., Hicks et al. 2009, 2013; Davies et al. 2014).

Each emission line was fitted with a single Gaussian profile on top of a linear continuum to trace the bulk motion of the line-emitting gas. Spectral regions around each line were chosen to avoid other lines and regions strongly affected by telluric features (1.965−1.995 μm for [Si VI] and 2.128−2.157 for H₂). Only spaxels that had at least one spectral channel with S/N > 3 were fitted, where the noise was determined as the local line-free rms of the spectrum. For [Si VI], we masked the 1.975−1.978 μm region where the H₂ (1–0) S(3) line is expected. Velocities were allowed to be ±1000 km s⁻¹ and velocity dispersions were allowed to be 0−500 km s⁻¹. Finally, 100 Monte Carlo iterations of the fit were performed by adding Gaussian noise to the spectra to determine the uncertainties on the best-fit line parameters.

Figure 16 shows the results of our fits where velocities have been corrected for the systemic velocity of the host galaxy as given in NED (2917 km s⁻¹). In addition, Fig. 17 also shows example fits to two pixels. The location of the pixels are plotted as red crosses in Fig. 16. There is a clear difference in the structure and kinematics of the ionised and molecular gas. The flux distribution of [Si VI] seems to change PA from ∼ − 18° on small scales to ∼ + 10° on large scales. These PAs were measured by visually inspecting the contours shown in the [Si VI] flux map and therefore only represent estimates with uncertainties of 5°. Lines representing the PAs are also shown in the [Si VI] flux map as blue lines.

Fig. 16. Normalised flux (left), LOS velocity (middle), and velocity dispersion (right) maps of the [Si VI] line (top row) and H2 (1–0) S(1) line (bottom row). The maps were measured from the same SINFONI cube from which we derived the Brγ profile used in our BLR modelling. The red crosses indicate the locations of the spaxels used to show example line fits in Fig. 17. The blue lines in the [Si VI] flux map show a PA of −18° and +10° which were estimated from the contours of the flux distribution (white contours). The blue line in the H2 flux map shows a PA of +10° estimated from the contours of the flux distribution (white contours). The blue line in the H2 velocity map shows a PA of −20° estimated visually from the velocity field.

Fig. 17. Example fits to individual spectra for the [Si VI] (left panels) and H2 (1–0) S(1) line (right panels). The observed spectra are shown in blue while the best-fit Gaussian models are shown in orange. Top row: spectral regions from a nuclear pixel, bottom row: spectral regions from an off-nuclear pixel. The specific pixels used are shown as white crosses in Fig. 16. The grey shaded region in the [Si VI] panels show the masked region corresponding to the expected location of the H2 (1–0) S(3).

The [Si VI] kinematics show clear non-circular signatures with a strong redshifted component to the north and high-velocity dispersion. This matches the analysis of Müller-Sánchez et al. (2011) who interpreted and modelled the kinematics as an outflow with a small contribution from disc rotation. The hot molecular gas, in contrast, shows a flux distribution with a PA of ∼ + 10° on all scales which is shown as a blue line in the H₂ flux map of Fig. 16. The kinematics are more indicative of disc rotation with a kinematic major axis along a PA of ∼ − 20° (shown as a blue line in the H₂ velocity map of Fig. 16 and estimated visually from the gradient of the velocity field) matching previous SINFONI results (Davies et al. 2007; Hicks et al. 2009; Müller-Sánchez et al. 2011). The LOS velocities are also relatively low (±50 km s⁻¹) and the dispersions relatively high (for a rotating disc) suggesting a thick disc with low inclination.

Interestingly, neither the [Si VI] nor the hot H₂ axes match the kinematic axis of the BLR. The larger scale flux distributions are close but the kinematic axes of the [Si VI] and H₂ are very different from the BLR. The BLR is blueshifted to the northeast and redshifted to the southwest while [Si VI] is redshifted to the north and blueshifted to the south. H₂ is blueshifted to the north but much more to the northwest. The hot molecular gas disc has a measured inclination of ∼35° (Davies 2007; Hicks et al. 2009; Müller-Sánchez et al. 2011) which is twice the inclination of the BLR but matches the larger kiloparsec-scale inclination of the host galaxy. Therefore, there must be a warping of the gas disc as it flows from 50 pc down to subparsec scales. Such a warping has been observed in other AGNs, most notably NGC 1068 (Impellizzeri et al. 2019) which shows counter-rotation at larger scales compared to the subparsec maser disc.

Nevertheless, our BLR orientation is in relatively good agreement with the polar axis measured by Smith et al. (2002, 2004). We measure a BLR polar axis of −65° compared with a polar axis of −45° measured through polarisation. While this is a difference of 20°, the uncertainty on the BLR polar axis is large (+55°, −49°). This further matches the polar dust angle of −50 to −60° found through mid-IR interferometry (Hönig et al. 2013; Burtscher et al. 2013; López-Gonzaga et al. 2016). Our BLR inclination also well matches the inferred inclination from the disc+wind model of Hönig & Kishimoto (2017), providing more evidence in favour of the interpretation that the extended mid-IR emission is tracing a dusty outflow.

To connect to the larger scale outflow traced by [Si VI], there must be a gradual shifting of the orientation as the outflow has expanded and interacted with the host galaxy ISM. On parsec scales the outflow begins with a PA of ∼ − 60°, in line with the polar axis of the BLR and the extended mid-IR component. The outflow then seems to shift northward, first to −18° and finally +10° matching the host galaxy disc. We note this is all still consistent with the [Ca VIII] size measurement which is only probing the very nuclear regions. Coronal line emission is easily produced at small scales due to their high critical densities and thus its expected for the brightest emission to occur in the nucleus. At larger scales, where the density is lower, it is likely that shock excitation instead of photoionisation is producing the coronal lines and is expected to be fainter (Rodríguez-Ardila et al. 2017b; May et al. 2018). Indeed the [Si VI] map shows a dominant unresolved central core and fainter, more diffuse extended emission up to 100 pc. To quantitatively determine the excitation mechanism at every scale would require detailed line modelling which is out of the scope of this paper. However, Rodríguez-Ardila et al. (2006) was able to explain the nuclear Fe and Si coronal line emission of NGC 3783 with photoionisation alone, but required additional shock excitation at radii larger than 100 pc. These shocks could be produced by a radio jet interacting with the ISM of NGC 3783, however both VLA and VLBA observations show only an unresolved component on scales of less than ∼20 pc (Schmitt et al. 2001; Orienti & Prieto 2010). The shocks are more likely produced by the AGN-radiation-pressure-driven wind.

Of particular interest is that we do not observe a large velocity shift in the nuclear [Ca VIII] line while at 55 pc we observe a LOS velocity of 150 km s⁻¹. Assuming a constant inclination angle, this could indicate an acceleration of the ionised gas from subparsec scales to scales of tens of parsecs as seen in other AGNs (e.g., Müller-Sánchez et al. 2011). Another explanation could be that the outflow was launched with a distribution of velocities such that the material with the largest velocities is at the largest radii, which has been observed in the NLR of NGC 1068 (Miyauchi & Kishimoto 2020).

7. Summary

In this paper, we report our analysis of VLTI/GRAVITY observations of the nearby type 1 AGN NGC 3783. We investigated three distinct components of the nuclear region around the AGN: (1) the BLR, (2) the hot dust, and (3) the CLR. Our main conclusions are as follows:

• We detect and successfully model the BLR as a rotating, thick disc with a physical mean radius of 16 light-days. Due to the centrally peaked but heavy tailed distribution of clouds, this leads to a cross-correlation function measured peak time-lag of ∼10.2 light-days, fully consistent with RM results. We find that the key parameters of interest for this paper, such as M_BH, appear to be robustly determined, but others, in particular R_min, are not. We will return to this issue in a future paper exploring joint modelling of GRAVITY and RM data.

• We reconstruct an image of the hot dust, revealing the presence of an offset cloud of gas and dust 0.6 pc in projected distance away from the main central hot dust component. We measure a FWHM size of 0.14 pc for the main component, consistent with previous interferometric results and with the expected dust sublimation radius.

• We measure a gas mass for the offset cloud of 20 M_⊙ and interpret it primarily as an AGN-heated cloud outflowing within the ionisation cone.

• We measure a FWHM size for the nuclear [Ca VIII]- emitting CLR of 0.4 pc, firmly placing it at the very inner regions of the NLR. Combined with our VLT/SINFONI data, we show the CLR is composed of a bright compact nuclear component and a fainter extended component out to 100 pc with outflow kinematics.

• We combine our results with past mid-IR interferometric and our VLT/SINFONI integral field unit data to establish a comprehensive view of the nuclear and circumnuclear region of NGC 3783 which includes an extended dusty outflow originating along the polar axis of the BLR. The AGN sits within a thick molecular gas disc that is feeding the AGN. Either the outflow or the molecular gas disc could be the origin of the offset cloud seen in the image reconstruction.

In Fig. 18 we show all of the components together in a single picture in an effort to place them all in context with each other. We note the cartoon is not to scale, and so all inferred distances and sizes of components are not correct. We further infer a counterclockwise rotation direction for the gas given the winding direction of the larger scale spiral arms of NGC 3783 (see den Brok et al. 2020 for a recent image) which in most spiral galaxies trail the direction of rotation (e.g., Buta 2011). However, the relative orientations do match the description described here. Putting all of these observations together produces a comprehensive and coherent picture of the gas structure and dynamics from 0.01 to 100 pc around an AGN, which was only made achievable with the impressive capabilities of near-IR interferometry and VLTI/GRAVITY. With GRAVITY, we aim to perform a similar analysis for the brightest and nearest AGNs, but the upgrade to GRAVITY+ will make it possible to analyse fainter and higher redshift AGNs, allowing us to trace the evolution of gas around AGNs as a function of cosmic time.

Fig. 18. Cartoon of the nuclear and circumnuclear region as described in Sect. 6. Different coloured clouds within the ionisation cone correspond to coronal and narrow-emission-line-emitting clouds. Arrows for the BLR and molecular gas disc indicate the direction of rotation. The image is not to scale in order to be able to show all components together.

Acknowledgments

We thank the anonymous referee for their careful reading and suggestions that helped to improve this manuscript. J.D. was supported in part by NSF grant AST 1909711 and an Alfred P. Sloan Research Fellowship. A.A. and P.G. were supported by Fundação para a Ciência e a Tecnologia, with grants reference UIDB/00099/2020 and SFRH/BSAB/142940/2018. SH acknowledges support from the European Research Council via Starting Grant ERC-StG-677117 DUST-IN-THE-WIND. JSB acknowledges the full support from the UNAM PAPIIT project IA 101220. 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 research made use of ASTROPY (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018), NUMPY (Van Der Walt et al. 2011), SCIPY (Jones et al. 2001), and MATPLOTLIB (Hunter 2007).

References

Appendix A: uv-binned differential spectra and model fits

In this appendix, we show the uv-binned differential visibility phase and amplitude spectra along with the best-fit BLR and CLR models described in the main text. Also shown are the full joint and marginalised posterior distributions for our BLR and CLR fits.

Fig. A.1. Top row: normalised Brγ profile (black points) with the best-fit model line profile (red line). Bottom rows: uv-binned, continuum-phase-subtracted differential phase spectra for each baseline (coloured points) with the best-fit BLR model spectra (solid lines). Model spectra for the BLR component are shown as dashed lines while the offset component is shown as dotted lines.

Fig. A.2. Corner plot for our BLR model fit showing the joint and marginalised posterior distributions for each free parameter. The blue lines and orange crosses indicate the maximum a posteriori position used as the best-fit values given in Table 2.

Fig. A.3. Top row: normalised GRAVITY Brγ profile (black points) with best-fit model line profile (red line) from our coronal line region analysis. Bottom rows: uv-binned differential visibility amplitude spectra for each baseline (coloured points and lines) with the best-fit CLR model (solid lines).

Fig. A.4. Same as Fig. A.3 but for the [Ca VIII] line.

Fig. A.5. Corner plot for our CLR analysis showing the joint and marginalised posterior distributions for each free parameter. The blue lines and orange crosses indicate the maximum a posteriori position used as the best-fit values given in Table 4.

Appendix B: Image reconstruction robustness tests

In this appendix, we test the robustness of the features seen in our image reconstruction of the hot dust continuum for NGC 3783. We specifically tested the imaging as a function of (1) the choice of data, (2) the choice of regularisation, and (3) the choice of image reconstruction algorithm.

B.1. Bootstrapping

Our first test involves testing the robustness of our image against the choice of data used in the reconstruction. We generated 100 sets of GRAVITY data with each set having 30% of the original data randomly removed. We then ran MiRA on each of the 100 data sets to produce 100 images. From these 100 images, we calculated a minimum and median image. The minimum image is constructed using the minimum value of all 100 images for each pixel. The median is the median value of all 100 images for each pixel. These are shown in Fig. B.1.

Fig. B.1. Minimum (left) and median images from our bootstrapping analysis of our image reconstruction. Values in the minimum image are determined as the minimum value from all 100 images produced. Values in the median image are the median value of all 100 images.

While the median image largely shows all the same features as the original image reconstruction, the minimum image instead has removed nearly all of the fainter features and left primarily the central source and the offset cloud. This strongly suggests that the fainter features are likely due to noise in the data and/or are artefacts of the image reconstruction. The offset cloud however is clearly a strong feature of the data that persists even after removing 30% of the data.

B.2. Regularisation

Our second test simply changes the specific regularisation used in the MiRA reconstruction. We instead applied the compactness regularisation which prioritises centrally located compact sources. This is in contrast to the hyperbolic regularisation which prioritises smooth, extended sources with sharp edges. Figure B.2 shows the reconstructed image under the compactness regularisation.

Fig. B.2. MiRA reconstructed image using the compactness regularisation.

As expected, the bright central source has decreased in size with a FWHM of 1.24 × 0.83 mas compared to 1.8 × 1.2 mas using the hyperbolic regularisation. However, the fainter offset cloud still remains at the same position relative the central source and contains the same 5% of the total flux. This shows that the offset cloud is robust against our choice of regularisation available in MiRA.

B.3. Image reconstruction algorithm

For our final test, we used another image reconstruction algorithm, SQUEEZE (Baron et al. 2010) to ensure our image features are robust against the choice of image reconstruction algorithm. In Fig. B.3, we show the final median image reconstruction produced by SQUEEZE. Just as with MiRA, we see two components, a bright central extended component and an offset fainter component. The location and brightness of the offset component is very consistent with MiRA and our model fitting. Therefore, we conclude that the properties of the offset cloud are robust against our choice of image reconstruction algorithm.

Fig. B.3. Image reconstruction of NGC 3783 using the SQUEEZE algorithm.

All Tables

All Figures

	Fig. 1. uv coverage of the binned GRAVITY observations. The radial span of each stripe is due to the wavelength range of the data. uv bins 1, 2, and 3 are ordered in clockwise direction for each baseline.
In the text

Fig. 2. Left: example spectra extracted from the SINFONI cube including the full integrated spectrum (blue), a spectrum from a nuclear spaxel (green), and an integrated spectrum from an off-nuclear region (orange). All spectra have been normalised to their median flux and slightly offset to improve visualisation. Right: K-band continuum image created by integrating the SINFONI cube between 2.25 and 2.31 μm. The green cross indicates the spaxel used for the nuclear spectrum in the right panel. The orange circle indicates the aperture used for the off-nuclear spectrum.

In the text

	Fig. 3. Normalised Brγ and [Ca VIII] line profiles observed by GRAVITY (blue) and SINFONI (orange).
In the text

	Fig. 4. Normalised Brγ profile from SINFONI (blue) with the best-fit line profile model (orange). The model is the combination of a narrow line template determined by the [O III]λ5007 Å line from X-shooter spectra (green) and a broad component consisting of two Gaussians (red). The residuals (purple) show that our fit very closely matches the observed Brγ profile.
In the text

	Fig. 5. Baseline averaged differential visibility amplitude spectra spanning the Brγ line (coloured points and grey error bars). The grey profile in each panel plots the same normalised Brγ spectrum from GRAVITY. The black curves represent the best-fit model from our CLR analysis in Sect. 5.
In the text

Fig. 6. Baseline averaged differential phase spectra spanning the Brγ line (coloured points and grey error bars). The grey profile in each panel plots the same normalised Brγ spectrum from SINFONI with the narrow component removed. The black curves represent the best-fit BLR model summed with the continuum phase calculated from our image reconstruction. The model line profile appears noisy because of the continuum phase contribution being calculated based on the observed line profile.

In the text

Fig. 7. Panel a: best-fit photocentres for the 13 brightest spectral channels covering the Brγ line to the differential phase data. The colour indicates wavelength and the ellipses outline the uncertainty on the photocentres. The red cross plots the (0,0) position. Panel b: best fit ‘2-pole’ photocentre model where all blueshifted (blue point) and redshifted (red point) channels are assumed to have the same photocentre with the addition of the wavelength-independent offset. The black point indicates the best-fit position of the ‘null’ model where all channels are assumed to have the same photocentre. Panel c: photocentres produced from our best-fitting BLR model (see Sect. 3.2) which we calculated by first producing mock differential phase spectra before fitting them in the same way as panel a. Colours are the same as in panel a.

In the text

Fig. 8. Panel a: cloud distribution of the best-fit BLR model. Each circle represents one cloud colour-coded by the LOS velocity. The green ellipse at the origin illustrates the uncertainty of the offset of the BLR centre. Panel b: observed averaged differential phase from baselines UT4−UT2, UT4−UT1, and UT3−UT1 after removing the residual ‘continuum phase’ signal (blue points) compared to the averaged differential phase from the best-fit BLR model (lower red line). The black points and overlaid red line show the observed and best model line flux profile respectively. Panel c: averaged differential phase data (points) and the best-fit models (lines) of the three baselines that show the strongest signal of the BLR component (dashed lines). The phase in panel b is calculated by averaging the phases of these three baselines after subtracting the best-fit residual continuum phases (dotted lines).

In the text

Fig. 9. Panel a: observed continuum light curve from Bentz et al. (2021) (black points) with the finer sampled interpolated continuum light curve over-plotted (grey line). Panel b: observed broad emission-line light curve (black points) with the emission-line light curve produced by reverberating the input continuum light curve off the clouds in our best-fit BLR model (red). Panel c: cross-correlation function between the continuum and broad emission-line light curves as a function of time-lag. The dashed line indicates the peak cross-correlation time-lag of 10.2 days.

In the text

Fig. 10. Relationship between BLR radius and 5100 Å AGN luminosity with the best-fit line from Bentz et al. (2013) (dashed line) and the RM samples from Du & Wang (2019) (black points) and Grier et al. (2017) (grey points). Previous GRAVITY BLR radius measurements of IRAS 09149−6206 and 3C 273 are shown along with the BLR mean radius for NGC 3783 as filled red stars. The time-lag measurement of ∼10 days for NGC 3783 from the literature is shown as an open star. The horizontal error bar indicates a factor of two uncertainty on the luminosity of NGC 3783 based on variability.

In the text

	Fig. 11. GRAVITY FT squared visibilities (left) and closure phases (right) as a function of spatial frequency for all nights of observation. The significant non-zero closure phases indicate resolved asymmetric hot dust structure in NGC 3783.
In the text

	Fig. 12. MiRA image reconstruction of the K-band continuum for NGC 3783. We used the hyperbolic regularisation with a FOV of 25.6 mas and a pixel size of 0.1 mas. The image is normalised to the maximum pixel flux. The GRAVITY interferometric beam is shown as a white ellipse in the bottom left corner for comparison. North is up and east is to the left.
In the text

	Fig. 13. Same as Fig. 12 but with lines indicating important angles and regions. Blue lines outline the ionisation cone as produced by the BLR polar axis (red line) and thickness. The yellow line indicates the position angle from polarimetery (Lira et al. 2020) and the green line indicates the position angle from mid-IR interferometry (Hönig et al. 2013).
In the text

	Fig. 14. Sample 2D Gaussian plus point-source model fit to the observed V2 and closure phase data for NGC 3783 from Jan 7, 2018. The best-fitting Gaussian FWHM geometric mean size of ≃0.6 mas is consistent with that found from 1D Gaussian fitting (GC20b). The Gaussian PA and point source flux fraction and offset are consistent with those found from image reconstruction.
In the text

	Fig. 15. Same as Fig. 5 but for the [Ca VIII] region.
In the text

Fig. 16. Normalised flux (left), LOS velocity (middle), and velocity dispersion (right) maps of the [Si VI] line (top row) and H2 (1–0) S(1) line (bottom row). The maps were measured from the same SINFONI cube from which we derived the Brγ profile used in our BLR modelling. The red crosses indicate the locations of the spaxels used to show example line fits in Fig. 17. The blue lines in the [Si VI] flux map show a PA of −18° and +10° which were estimated from the contours of the flux distribution (white contours). The blue line in the H2 flux map shows a PA of +10° estimated from the contours of the flux distribution (white contours). The blue line in the H2 velocity map shows a PA of −20° estimated visually from the velocity field.

In the text

Fig. 17. Example fits to individual spectra for the [Si VI] (left panels) and H2 (1–0) S(1) line (right panels). The observed spectra are shown in blue while the best-fit Gaussian models are shown in orange. Top row: spectral regions from a nuclear pixel, bottom row: spectral regions from an off-nuclear pixel. The specific pixels used are shown as white crosses in Fig. 16. The grey shaded region in the [Si VI] panels show the masked region corresponding to the expected location of the H2 (1–0) S(3).

In the text

	Fig. 18. Cartoon of the nuclear and circumnuclear region as described in Sect. 6. Different coloured clouds within the ionisation cone correspond to coronal and narrow-emission-line-emitting clouds. Arrows for the BLR and molecular gas disc indicate the direction of rotation. The image is not to scale in order to be able to show all components together.
In the text

	Fig. A.1. Top row: normalised Brγ profile (black points) with the best-fit model line profile (red line). Bottom rows: uv-binned, continuum-phase-subtracted differential phase spectra for each baseline (coloured points) with the best-fit BLR model spectra (solid lines). Model spectra for the BLR component are shown as dashed lines while the offset component is shown as dotted lines.
In the text

	Fig. A.2. Corner plot for our BLR model fit showing the joint and marginalised posterior distributions for each free parameter. The blue lines and orange crosses indicate the maximum a posteriori position used as the best-fit values given in Table 2.
In the text

	Fig. A.3. Top row: normalised GRAVITY Brγ profile (black points) with best-fit model line profile (red line) from our coronal line region analysis. Bottom rows: uv-binned differential visibility amplitude spectra for each baseline (coloured points and lines) with the best-fit CLR model (solid lines).
In the text

	Fig. A.4. Same as Fig. A.3 but for the [Ca VIII] line.
In the text

	Fig. A.5. Corner plot for our CLR analysis showing the joint and marginalised posterior distributions for each free parameter. The blue lines and orange crosses indicate the maximum a posteriori position used as the best-fit values given in Table 4.
In the text

	Fig. B.1. Minimum (left) and median images from our bootstrapping analysis of our image reconstruction. Values in the minimum image are determined as the minimum value from all 100 images produced. Values in the median image are the median value of all 100 images.
In the text

	Fig. B.2. MiRA reconstructed image using the compactness regularisation.
In the text

	Fig. B.3. Image reconstruction of NGC 3783 using the SQUEEZE algorithm.
In the text