The central parsec of NGC 3783: a rotating broad emission line region, asymmetric hot dust structure, and compact coronal line region

Using VLTI/GRAVITY and SINFONI data, we investigate the sub-pc gas and dust structure around the nearby type 1 AGN hosted by NGC 3783. The K-band coverage of GRAVITY uniquely allows a simultaneous analysis of the size and kinematics of the broad line region (BLR), the size and structure of the near-IR continuum emitting hot dust, and the size of the coronal line region (CLR). We find the BLR probed through broad Br$\gamma$ emission is well described by a rotating, thick disk with a radial distribution of clouds peaking in the inner region. In our BLR model the physical mean radius of 16 light days is nearly twice the 10 day time lag that would be measured, which matches very well the 10 day time lag that has been measured by reverberation mapping. We measure a hot dust FWHM size of 0.74 mas (0.14 pc) and further reconstruct an image of the hot dust which reveals a faint (5% of the total flux) offset cloud which we interpret as an accreting cloud heated by the central AGN. Finally, we directly measure the FWHM size of the nuclear CLR as traced by the [CaVIII] and narrow Br$\gamma$ line. We find a FWHM size of 2.2 mas (0.4 pc), fully in line with the expectation of the CLR located between the BLR and narrow line region. Combining all of these measurements together with larger scale near-IR integral field unit and mid-IR interferometry data, we are able to comprehensively map the structure and dynamics of gas and dust from 0.01--100 pc.


Introduction
Recent advancements in infrared interferometric observations has allowed significant progress in understanding the gas and dust structure and dynamics around active galactic nuclei (AGN). With the new capabilities of GRAVITY (Gravity Collaboration et al. 2017), the second generation instrument at the Very Large Telescope Interferometer (VLTI), the broad line region has been resolved and modelled to provide a measurement of the supermassive black hole mass (Gravity Collaboration et al. 2018, 2020a, the hot dust has been imaged revealing a thin ring at the dust sublimation radius (Gravity Collaboration et al. 2020c), and hot dust sizes have been measured for an increasing number of AGN (Gravity Collaboration et al. 2020b). NGC 3783 provides a unique laboratory to study not only these aspects but also the coronal line region (CLR) in a single object and build a comprehensive picture of the nuclear and circumnuclear region around an AGN.
GRAVITY is developed in a collaboration by the Max Planck Institute for extraterrestrial Physics, LESIA of Observatoire de Paris/Université PSL/CNRS/Sorbonne Université/Université de Paris and IPAG of Université Grenoble Alpes / CNRS, the Max Planck Institute for Astronomy, the University of Cologne, the CENTRA -Centro de Astrofisica e Gravitação, and the European Southern Observatory.
Corresponding author: T. Shimizu (e-mail: shimizu@mpe.mpg. de) NGC 3783 hosts one of the most luminous local AGN, with a bolometric AGN luminosity, log L AGN ∼ 44.5 erg s −1 for a distance of 38.5 Mpc (Tully & Fisher 1988;Davies et al. 2015) 1 , and is the subject of an extensive literature, most notably related to ionised outflows (especially X-ray warm absorbers) and variability. As one of the most luminous AGN, NGC 3783 has been intensely monitored via reverberation mapping and is one of the few AGN 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 AGN that clearly demonstrated the virial relationship between emission line lag and width (Onken & Peterson 2002). Through reverberation mapping, the black hole mass of NGC 3783 has been estimated to 1 We adopt a distance of 38.5 Mpc based on the Tully-Fisher (TF) relation and reported on the NASA Extragalactic Database (NED) which is significantly smaller than the luminosity or angular size distance inferred from the redshift (47-48 Mpc) and due to the peculiar velocity of the galaxy. This distance though could potentially be affected by AGN contamination and thus underestimated. However, Crook et al. (2007) found NGC 3783 to reside in a group of four galaxies at a distance of 37 Mpc and Kourkchi & Tully (2017) found NGC 3783 to reside in a group of nine galaxies at a distance of 42 Mpc. Therefore, we choose to use the TF based distance of 38.5 Mpc throughout our analysis and note the exact value of the distance within a few Mpc does not change the results.
The X-ray warm absorber was initially modelled by Netzer et al. (2003) as requiring three different ionization components at two different velocities, at distances within limits of 0.2-25 pc, and with a total column of N H ∼ 4 × 10 22 cm −2 . With more detailed spectra, this was later expanded to three velocities in the range −460 to −1600 km s −1 by Mao et al. (2019), comparable to those identified previously in UV spectra of −550, −720, and −1370 km s −1 (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 23 cm −2 , 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 −1 . In addition, from the implied density of 3 × 10 9 cm −3 and similarity of the radial location at only 10 light-days to the size of the broad line region (BLR) from reverberation mapping (Onken & Peterson 2002), these authors argued that it is an obscuring wind originating in the outer part of the BLR. A subsequent analysis of archival data by Kaastra et al. (2018) suggested that such obscuring events, with columns exceeding 5 × 10 21 cm −2 may be rather common for NGC 3783.
Ionized 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/outer opening angles of 45/55 • and a maximum outflow velocity of 130 km s −1 . In contrast, based on the biconical appearance of the near-infrared [SiVI] distribution at 16 pc resolution, Müller-Sánchez et al. (2011) modelled the line kinematics as an ionization cone that is 60 • from face-on with inner/outer opening angles 27/34 • and a maximum outflow velocity of 400 km s −1 . 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-infrared 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-infrared spectral energy distribution using a disk+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 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. 2016b), 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 disk model, Smith et al. (2002Smith et al. ( , 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 has been hinted at recently through observations of Circinus (e.g. Stalevski et al. 2017Stalevski et al. , 2019. Unfortunately, radio maps do not help in determining the orientation of the polar axis since 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 disk is such that the kinematic major axis, measured from both the stellar and H 2 1-0 S(1) velocity fields on scales 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 disk there is little warping. It 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 north-west, and cannot provide a guide for disentangling the innermost geometry of this galaxy.
In this paper we present new data from the VLTI spectrometer GRAVITY (Gravity Collaboration et al. 2017) which combines the light from all 4 UTs of the VLT, and can resolve structures on scales smaller than a few milliarcsec and provide spectro-astrometric relative precision of tens of microarsec. This resolution is essential to resolve the broad line region 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 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 et al. 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 0 = 67.8 km s −1 Mpc −1 (Planck Collaboration et al. 2016). Using this cosmology and our adopted distance of 38.5 Mpc, 1 pc subtends 5.36 mas on sky and 1 µas corresponds to 0.22 light days. black hole. 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.3ms integrations with the FT fibre and 2) a series of 10 exposures with detector on-chip integration times (DITs) of 30s 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 et al. (2020a, hereafter GC20a) andGravity Collaboration et al. (2020b, hereafter GC20b) and chose to retain all SC DITs regardless of FT SNR (snr-min-ft = 0) and estimated visibility loss (vfactor-min-sc = 0) which we found results in a substantial improvement in the phase noise.

Observations and Data Reduction
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 accounts for residual phase features introduced as a result of the 3rd order polynomial used in the pipeline being 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. SC differential amplitude spectra are produced by fitting and dividing out a 2nd order polynomial from the pipeline reduced visibility amplitude spectra. We estimate the uncertainty for the differential phase and amplitude spectra by calculating the RMS in the line-free regions.
To still further improve the SNR 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 : 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.
uncertainty as the weight. Fig. 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, however the differential phase spectra have had the "continuum phase" signal subtracted (see Sec. 3.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 −1 scale and the K grating (R 4000). A total of 6 exposures of 100s 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 reconstruct 12.5 mas pix −1 images with a spectrum at each pixel, applying flat-fielding, bad-pixel, distortion, and cosmic-ray hit corrections. The wavelength calibration is 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 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 6 pixel FWHM Gaussian in the spatial directions, equal to the 75 mas FWHM point spread function (PSF) measured from a fit to the continuum 3 , and a 2 pixel FWHM Gaussian in the spectral direction, matching the instrumental dispersion. Fig. 2 plots 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 2 (1-0) S(1), Brγ, and [Ca viii] emission lines.

The 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 NGC 3783's broad Brγ profile. Fortunately, these calibrators were also observed on the same night as IRAS 09149-6206, whose Brγ emission does not overlap with the calibrator's Brγ absorption. Thus, we fit a line profile to the calibrator's absorption 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. Fig. 3 shows the final GRAVITY mean spectrum for both Brγ and [Ca viii].
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. 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 since the continuum is much more compact.
Our choice of which line profile to use depends on the specific analysis we want to perform. For our BLR modelling (Sec. 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 of the narrow line. To produce the template, we fit the [O iii] line from previous X-Shooter spectra (Schnorr-Müller et al. 2016b) with 4 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 2 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. Fig. 4 plots the results of our line fitting and shows that this model very well reproduces Brγ profile. The narrow component was then subtracted from the spectrum to produce the profile we will use for our BLR analysis.

The 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 et al. 2018).
In contrast to IRAS 09149-6206 though, we are able to directly constrain the "continuum phase" contribution through our image reconstruction of the hot dust continuum (see Sec. 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.

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, where u is the uv coordinate of the baseline, x BLR,λ is the BLR coordinate of the photocentre at wavelength λ w.r.t. the photocentre of the continuum.  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 normalized by 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.  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 off, we still see a slight residual systematic offset of the photocentres from the origin. This is likely caused by our 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 "Sshape" signal in other baselines.
We further test the robustness of the velocity gradient by fitting for a single photocentre for all 5 blueshifted channels and all 8 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 ∼ 10 light-day time lags measured from re-  line from X-Shooter spectra (green) and a broad component consisting of 2 Gaussians (red). The residuals (purple) show that our fit very well matches the observed Brγ profile. verberation mapping (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.

BLR modelling
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 et al. (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 disk, θ 0 .
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 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 cloud's radial position and M BH . The distribution fol-lows an ellipse described by Equation 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 sky by an inclination angle, i, and position angle, PA and translated by an offset, (x 0 , y 0 ). Line-of-sight (LOS) 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 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 Equation 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 Observed wavelength ( m) Fig. 7: (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. (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. (c) Photocentres produced from our best fitting BLR model (see Section 3.2) which we calculated by first producing mock differential phase spectra and fit them in the same way as panel (a). Colours are the same as in panel (a).

Parameters
Best Fit Values Table 2: The best fit parameters and central 95% credible interval for the modelling of the BLR spectrum and differential phase of NGC 3783. Because f flow is a binary switch, we instead report P(inflow) which indicates the probability of inflow where f flow < 0.5. ∆v BLR is the difference between the velocity derived from the best-fit λ emit and the systemic velocity based on the redshift.
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 Our best fit BLR model well reproduces both the differential phase and Brγ line profile with a reduced chi-square, χ 2 r = 0.665. Fig. 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 Fig. 8a and Fig. 7c 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 sky away from the PA of the rotating disk.
With a best fit inclination of i = 23 • and θ 0 = 24 • , we find a relatively face-on and moderately thick disk describes the BLR of NGC 3783 well, as expected for a type 1 AGN and similar to both 3C 273 (Gravity Collaboration et al. 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. Fig. 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 0 indicating 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 green ellipse at the origin illustrates the uncertainty of the offset of the BLR centre. (b) The 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. (c) The 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).
the clouds on the outer surface of the BLR are not radiating more than the inner clouds. The motion of the clouds are split fairly 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. 4 Our finding of a significant fraction of inflowing clouds agrees with the simple velocity resolved time lag measurements from Bentz et al. (2020) which showed slightly shorter lags in the red wing of the line profile. 4 Instead of reporting the exact value of f flow which has no physical meaning we instead report the probability that f flow < 0.5. 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, R BLR = 71 µas which corresponds to 16 light-days and a minimum radius, 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 ionization lines such as He ii, C iv, and Si iv have been measured between 1-4 days (Reichert et al. 1994;Onken & Peterson 2002;Bentz et al. 2020) 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. 2016a;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 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 nearly 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 reverberation mapping and GRAVITY analysis of NGC 3783 that 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 reverberation mapping and discuss the SMBH mass and NGC 3783's place on the R-L relation.

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 reverberation mapping. However, in order to compare the model size from interferometry with reverberation mapping, 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 so differences between the two are not unexpected.
To test this, we generated a mock emission line light curve based off the real continuum light curve from Bentz et al. (2020). 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. Fig. 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 calculate the time lag between the two light curves using the cross correlation function (CCF). The results of the CCF are shown in Fig. 9c.
We find a peak time lag of 10.2 +4.2 −2.7 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 5100 = 42.93 erg s −1 , from Bentz et al. (2020) 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. (2020) 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 2 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 thus we only plot one point for them in Fig. 10.
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 NGC 3783's discrepancy with the R − L relation is different than 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 AGN 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, reverberation mapping time lags could be biased towards the inner radius.
The discrepancy we observe could also be related to the same trend we have observed with the dust sizes. (GC20b) showed that lower luminosity AGN 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 AGN significantly decreases with increasing Eddington ratio due to the AGN driving gas and dust out of the central few parsecs (Ricci et al. 2017). If the gas and dust is 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 is driven out primarily from the inner regions and reduces 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.

The 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 NIR continuum. Fig. 11 shows a general trend of decreasing V 2 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 thus 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.

Continuum Image Reconstruction
The high quality and large uv coverage of the FT data allows us to utilise image reconstruction codes to spatially map the hot dust. As we did for our imaging of NGC 1068 (Gravity Collaboration et al. 2020c), we used MiRA 5 (Multi-aperture Image Reconstruction Algorithm Thiébaut 2008) to reconstruct the Kband image of NGC 3783 based on the V 2 and closure phase data as shown in Fig. 11. The V 2 uncertainties shown have been multiplied by 10 to match the S/N of the closure phases and avoid overweighting V 2 . 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/NIR 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 hyperparame-5 publicly available at https://github.com/emmt/MiRA ters 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. Too small values of τ will produce a cartoon-like image with very smooth regions that are then sharply cutoff. Large values of τ instead lead to effectively 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 5 and 10 −3 respectively. Fig. 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.8x1.2 mas and PA ∼ −44 • (East of North) for the bright central source and 1.5x0.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.
The PA and extent of the individual sources follow the PA of the GRAVITY beam 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 2 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 2 , however, is only moderately well matched. In particular, the gradient of V 2 seems steeper for the image reconstruction which leads to a larger primary source. This leads to the factor of ∼ 2 difference between the image reconstructed size and the size found in GC20b (0.82 mas) from fitting a Gaussian model to only V 2 data. Therefore, as a further test of the detection and specific properties of each source, in the next section we apply the same model fitting.

Visibility model fitting
We fit the individual V 2 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 2 errors by a factor of 10 to match the Article number, page 11 of 28 A&A proofs: manuscript no. gravi-agn_ngc3783_astro-ph 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.
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 on the choice of regularisation especially since 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 next section.

The 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 be indicative 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 NIR luminosity of 10 42 erg s −1 and a bolometric luminosity of 10 42.7 erg s −1 using the NIR-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 of the black hole mass of 4 × 10 4 M and assuming it is accreting at 10% Eddington, similar to the primary, gives M BH = 4 × 10 5 M . As a first order calculation of the expected orbital velocity of the secondary, we can simply use v = √ GM/r with M = 10 8 M found from our BLR modelling since the secondary must be much less massive than the primary. At a radius of 0.6 pc, we would then expect v ≈ 900 Article number, page 12 of 28     (Lira et al. 2020) and the green line indicates the position angle from MIR interferometry (Hönig et al. 2013). km s −1 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 could not achieve any 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, the current data on NGC 3783 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 disk or potentially outflowing away from it. The dust emission then would 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 ionization 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 MIR 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 an nonphysical amount of dust would be needed to produce the NIR 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), 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 46 erg s −1 and α = −2.1. We adjust for NGC 3783's AGN luminosity (log L bol = 44.5) by scaling with L 1/2 bol 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 spectral energy distribution (SED) data point, we use a simple modified blackbody model to calculate the dust mass.
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 κ 0 = 167 m 2 kg −1 , β = 1.25, and ν 0 = 137 THz (2.19 µm) from Draine (2003). With T d = 604 K from Eq. 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. cloud directly measured from the image, the density of the cloud would be ∼ 10 4 cm −3 which is at the high end of the observed densities in the NLR of local AGN ). 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.

The 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 AGN also show lines from highly ionized atoms, the so-called "coronal lines", for their first observation in the solar corona. These lines all have ionization potentials ≥ 100 eV and are collisionally excited forbidden transitions with relatively high critical densities (10 7 − 10 10 cm −3 ). 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, 2011Lamperti et al. 2017) it is thought that the CLR lies between the BLR and NLR. A distinct increase in line FWHM with ionization potential for some AGN 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 photoionization as the primary ionization 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 obscur-ing structure (Pier & Voit 1995;Ferguson et al. 1997;Murayama & Taniguchi 1998). Arguments against this hypothesis though include the fact that coronal lines are observed equally in type 1 and type 2 AGN (e.g. Rodríguez-Ardila et al. 2011) and that significant coronal line emission is seen in even the most heavily Compton thick AGN 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;Reunanen et al. 2003;Rodríguez-Ardila et al. 2006;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 indicated 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.
The differential visibility amplitude across an emission line is described by the following equation, 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. Brγ, however, 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, 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 Eq. 4 and 5 we model all emitting regions as Gaussian sources. In the marginally resolved limit, the visibility of a Gaussian follows, where V 0 is the zero baseline visibility, r uv is the baseline length in units of mas −1 , and FWHM is the size of the source in mas.
While a single compact source should always have V 0 = 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 since the BLR is effectively unresolved. The measured size of 101 µas from spectroastrometry (see Section 3.2) results in V b = 0.995 at the longest baseline observed. Eq. 4 and 5 also show that the differential visibility amplitude depends on the normalised flux profile. While for the [Ca viii] line we could 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 (Pier & Voit 1995;Ferguson et al. 1997;Murayama & Taniguchi 1998). This however does not rule out an inner region origin for all coronal lines especially since [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 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 low velocity shift of the [Ca viii] and relatively low FWHM since at larger radii, the gravitational potential of the stellar bulge should slow down the clouds.
The [Ca viii] flux 6 of 4×10 −18 W m −2 observed from a region of radius 0.92 mas 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 f = 0.1, 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 2-level analysis, and a Ca abundance of 2.e-6 by number (Landi et al. 2004), n e ≈ 2 × 10 5 cm −3 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 narrow line region, see, e.g.,  for a recent assessment of NLR densities. The required CLR density scales with f −0.5 , i.e. the electron density in the [Ca viii] emitting region must be at least about 10 5 cm −3 but could be clearly higher if arising in very low filling factor clouds or filaments.

Connecting the nuclear and circumnuclear regions
With these new GRAVITY observations of NGC 3783, we have 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 with a 0.14 pc 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 2 (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 ionization potential (IP = 167 eV) and thus allows for tracing the ionized gas out to larger scales. The H 2 line traces hot molecular gas that is likely inflowing and feeding the AGN (e.g Hicks et al. 2009Hicks et al. , 2013Davies et al. 2014).
Each emission line was fit 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 2 ). Only spaxels which had at least one spectral channel with S/N> 3 were fit 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 2 (1-0) S(3) line is expected. Velocities were allowed to be ±1000 km s −1 and velocity dispersions were allowed to be 0-500 km s −1 . 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. Fig. 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 −1 ). In addition we also show example fits to two pixels in Fig. 17. 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.
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 disk 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 2 flux map of Fig. 16. The kinematics are more indicative of disk rotation with a kinematic major axis along a PA of ∼ −20 • (shown as a blue line in the H 2 velocity map of Fig. 16 and estimated visually from the gradient of the velocity field) matching previous SINFONI results Hicks et al. 2009;Müller-Sánchez et al. 2011). The LOS velocities are also quite low (±50 km s −1 ) and the dispersions relatively high (for a rotating disk) suggesting a low inclination, thick disk.
Interestingly, neither the [Si vi] nor the hot H 2 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 2 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 2 is blueshifted to the north but much more to the northwest. The hot molecular gas disk 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 kpc scale inclination of the host galaxy. Thus, there must be a warping of the gas disk as it flows from 50 pc down to sub-pc scales. Such a warping has been observed in other AGN, most notably NGC 1068 (Impellizzeri et al. 2019) which shows counter-rotation at larger scales compared to the sub-pc maser disk.
Our BLR orientation however is in relatively good agreement with the polar axis measured by Smith et al. (2002) and Smith et al. (2004). We measure a BLR polar axis of -65 • compared with a polar axis of -45 • measured through polarization. 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 disk+wind model of Hönig & Kishimoto (2017) lending more evidence in favour of the interpretation that the extended MIR 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 MIR component. The outflow then seems to shift northward first to −18 • and finally +10 • matching the host galaxy disk. 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 photoionization is producing the coronal lines and is expected to be fainter 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 photoionization alone but required additional shock excitation at radii larger than 100 pc. These shocks could be produced by a radio jet interacting with the interstellar medium of NGC 3783, however both VLA and VLBA observations show only an unresolved component at scales less than ∼20 pc (Schmitt et al. 2001;Orienti & Prieto 2010). More likely is that the shocks are produced by the AGN radiation pressure driven wind.
An interesting thing to note is 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 ionized gas from sub-pc to 10s of pc scale as seen in other AGN (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).

Summary
In this paper, we reported on 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 broad line region, 2) the hot dust, and 3) the coronal line region. Our main conclusions are as follows: -We detect and successfully model the BLR as a rotating, thick disk 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 reverberation mapping results. We find that the key parameters of interest for this Paper like M BH appear to be robustly determined but others, in particular R min , are not. In a future paper exploring joint modelling of GRAVITY and reverberation mapping data we will return to this issue. -We reconstruct an image of the hot dust which reveals 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 combine our results with past MIR interferometric and our VLT/SINFONI integral field unit data to establish a comprehensive view of the nuclear and circumnuclear region of NGC 3783 that includes an extended dusty outflow originating along the polar axis of the BLR. The AGN sits within a thick molecular gas disk that is feeding the AGN. Both the outflow or the molecular gas disk 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 into context with each other. We note the cartoon is not to scale so all inferred distances and sizes of components are not correct. We have further inferred 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 -100 pc around an AGN that was only achievable with the impressive capabilities of NIR interferometry and VLTI/GRAVITY. With GRAVITY, we aim to perform a similar analysis for the brightest and nearest AGN, but the upgrade to GRAVITY+ will open a window to fainter and higher redshift AGN to be able to trace the evolution of gas around AGN as a function of cosmic time.  Fig. 18: A cartoon of the nuclear and circumnuclear region as described in Sec. 6. Different coloured clouds within the ionization cone correspond to coronal and narrow emission line emitting clouds. Arrows for the BLR and molecular gas disk indicate the direction of rotation. The image is not to scale in order to be able to show all components together. Observed wavelength ( m) Differential phase ( )  Table 2 Article number, page 23 of 28      Table 4 Article number, page 26 of 28