11institutetext:  Max Planck Institute for Extraterrestrial Physics (MPE),
Giessenbachstr.1, 85748 Garching, Germany 22institutetext: LESIA, Observatoire
de Paris, Université PSL, CNRS, Sorbonne Université, Univ. Paris Diderot,
Sorbonne Paris Cité, 5 place Jules Janssen, 92195 Meudon, France
33institutetext: I. Institute of Physics, University of Cologne, Zülpicher
Straße 77, 50937 Cologne, Germany 44institutetext: Departments of Physics and
Astronomy, Le Conte Hall, University of California, Berkeley, CA 94720, USA
55institutetext: Department of Physics and Astronomy, University of
Southampton, Southampton, UK 66institutetext: Department of Physics, Kyoto
Sangyo University, Kita-ku, Japan 77institutetext: Université Côte d’Azur,
Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, Nice, France
88institutetext: School of Physics and Astronomy, Tel Aviv University, Tel
Aviv 69978, Israel 99institutetext: Univ. Grenoble Alpes, CNRS, IPAG, 38000
Grenoble, France 1010institutetext: Center for Computational Astrophysics,
Flatiron Institute, 162 5th Ave., New York, NY 10010, USA 1111institutetext:
European Southern Observatory, Casilla 19001, Santiago 19, Chile
1212institutetext: European Southern Observatory, Karl-Schwarzschild-Str. 2,
85748 Garching, Germany 1313institutetext: Sterrewacht Leiden, Leiden
University, Postbus 9513, 2300 RA Leiden, The Netherlands 1414institutetext:
Max Planck Institute for Radio Astronomy, Auf dem Hügel 69, 53121 Bonn,
Germany 1515institutetext: Universidade de Lisboa - Faculdade de Ciências,
Campo Grande, 1749-016 Lisboa, Portugal 1616institutetext: Faculdade de
Engenharia, Universidade do Porto, rua Dr. Roberto Frias, 4200-465 Porto,
Portugal 1717institutetext: CENTRA - Centro de Astrofísica e Gravitação, IST,
Universidade de Lisboa, 1049-001 Lisboa, Portugal 1818institutetext: Max
Planck Institute for Astronomy, Königstuhl 17, 69117, Heidelberg, Germany
1919institutetext: Instituto de Astrofísica de Canarias (IAC), E-38200 La
Laguna, Tenerife, Spain 2020institutetext: Department of Astrophysical &
Planetary Sciences, JILA, University of Colorado, Duane Physics Bldg., 2000
Colorado Ave, Boulder, CO 80309, USA 2121institutetext: Research School of
Astronomy and Astrophysics, Australian National University, Canberra, ACT
2611, Australia 2222institutetext: Department of Physics, Technical University
Munich, James-Franck-Straße 1, 85748 Garching, Germany 2323institutetext:
Department of Physics and Astronomy, Georgia State University, Atlanta, GA
30303, USA 2424institutetext: LCOGT, 6740 Cortona Drive, Suite 102, Goleta, CA
93117, USA

# A geometric distance to the supermassive black Hole of NGC 3783

GRAVITY Collaboration: A. Amorim 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.15151717 M. Bauböck 11 M. C. Bentz 2323
W. Brandner 1818 M. Bolzer 11 Y. Clénet 22 R. Davies 11 P. T. de Zeeuw 111313
J. Dexter 202011 A. Drescher 112222 A. Eckart 331414 F. Eisenhauer 11 N. M.
Förster Schreiber 11 P. J. V. Garcia 111116161717 R. Genzel 1144 S. Gillessen
11 D. Gratadour 222121 S. Hönig 55 D. Kaltenbrunner 11 M. Kishimoto 66 S.
Lacour 221212 D. Lutz 11 F. Millour 77 H. Netzer 88 C. A. Onken 2121 T. Ott 11
T. Paumard 22 K. Perraut 99 G. Perrin 22 P. O. Petrucci 99 O. Pfuhl 1212 M. A.
Prieto 1919 D. Rouan 22 J. Shangguan Corresponding author: J. Shangguan
e-mail<EMAIL_ADDRESS>T. Shimizu 11 J. Stadler 11 A. Sternberg 881010
O. Straub 11 C. Straubmeier 33 R. Street 2424 E. Sturm 11 L. J. Tacconi 11 K.
R. W. Tristram 1111 P. Vermot 22 S. von Fellenberg 11 F. Widmann 11 J. Woillez
1212

The angular size of the broad line region (BLR) of the nearby active galactic
nucleus (AGN) NGC 3783 has been spatially resolved by recent observations with
VLTI/GRAVITY. A reverberation mapping (RM) campaign has also recently obtained
high quality light curves and measured the linear size of the BLR in a way
that is complementary to the GRAVITY measurement. The size and kinematics of
the BLR can be better constrained by a joint analysis that combines both
GRAVITY and RM data. This, in turn, allows us to obtain the mass of the
supermassive black hole in NGC3783 with an accuracy that is about a factor of
two better than that inferred from GRAVITY data alone. We derive
$M_{\mathrm{BH}}=2.54_{-0.72}^{+0.90}\times 10^{7}\,M_{\odot}$. Finally, and
perhaps most notably, we are able to measure a geometric distance to NGC 3783
of $39.9^{+14.5}_{-11.9}$ Mpc. We are able to test the robustness of the BLR-
based geometric distance with measurements based on the Tully-Fisher relation
and other indirect methods. We find the geometric distance is consistent with
other methods within their scatter. We explore the potential of BLR-based
geometric distances to directly constrain the Hubble constant, $H_{0}$, and
identify differential phase uncertainties as the current dominant limitation
to the $H_{0}$ measurement precision for individual sources.

###### Key Words.:

galaxies: active – galaxies: nuclei – galaxies: Seyfert – (galaxies:) quasars:
individual: NGC 3783 – (cosmology:) distance scale

##  1 Introduction

Trigonometry is the basis of distance measurements. The parallax method uses
the motion of the Earth around the Sun to measure the angular displacement of
a nearby star (Bessel, 1838). From the Hipparcos satellite of the European
Space Agency (ESA; ESA 1997) to the recent Gaia mission (Gaia Collaboration et
al., 2016), the parallaxes of more than 1 billion stars in the Milky Way have
now been measured (Gaia Collaboration et al., 2018). The geometric method can
also be applied to any object, as long as both its physical size ($R$) and
angular size ($\Theta$) are measurable, as $D=R/\Theta$. For a distant
extragalactic target, the measured geometric distance is the angular diameter
distance ($D_{A}$), including cosmological expansion. In contrast to “standard
candles,” such as pulsating stars (e.g., Bhardwaj, 2020) and Type Ia
supernovae (SNe Ia; e.g., Phillips 1993; Riess et al. 1996), which measure the
luminosity distance ($D_{L}$), the geometric method does not rely on the
calibration based on the so-called distance ladder (e.g., Riess et al., 2009,
2021).

However, objects to which the geometric method can be applied are usually rare
as it is difficult to measure both $R$ and $\Theta$ for the same target. The
distance to the supermassive black hole (BH) at the center of the Milky Way is
measured to a $<0.5\%$ uncertainty level, based on the 27-year astrometric and
spectroscopic monitoring of the 16-year orbital motion of the star S2 (Gravity
Collaboration et al., 2019). For detached eclipsing binary stars, the linear
size of each component can be measured from photometric and spectroscopic
monitoring, while the angular size of each star can be derived from its
empirical relation with the color of the star (Lacy, 1977). This method has
been used to measure the distance to the Large Magellanic Cloud with high
accuracy (Pietrzyński et al., 2019). By monitoring the Keplerian motion of
water maser-emitting gas with very long baseline interferometry (VLBI)
observations, one can measure the geometric distance to nearby Type-2 Seyfert
galaxies that host an observable megamaser disk (Herrnstein et al., 1999;
Braatz et al., 2010). Baryonic acoustic oscillations can also provide $D_{A}$
using the clustering of galaxies at a certain redshift range (e.g., Eisenstein
et al., 2005; Anderson et al., 2014). The broad line region (BLR) of active
galactic nuclei (AGNs) has also been proposed as a probe of geometric distance
(Elvis & Karovska, 2002). The linear size of the BLR can be measured by the
reverberation mapping (RM) technique (Peterson, 2014), while the angular size
of the BLR can be measured by near-infrared (NIR) interferometry (Petrov et
al., 2001; Woillez et al., 2004). Unlike detached eclipsing binaries and
megamaser systems, which are difficult to discover, AGNs are luminous sources
that are commonly found locally as well as out to high redshifts, even at
$z>6$. A similar approach has been applied to NGC 4151 by resolving the hot
dust continuum emission instead of the BLR (Hönig et al., 2014).

Recently, GRAVITY, the second generation Very Large Telescope Interferometer
(VLTI) instrument, spatially resolved the BLR of 3C 273 for the first time
using the spectroastrometry (SA) technique (Bailey, 1998). GRAVITY combines
all four of the 8 m Unit Telescope (UT) beams to yield six simultaneous
baselines (Gravity Collaboration et al., 2017). Gravity Collaboration et al.
(2018) reported the mean radius of the BLR as $46\pm 10\,\mbox{$\mu$as}$ and a
BH mass of $(2.6\pm 1.1)\,\times 10^{8}\,M_{\odot}$, which is fully consistent
with that measured by Zhang et al. (2019) using 10 year RM data. Wang et al.
(2020) conducted the first joint analysis (hereafter, SARM, as the combination
of SA and RM) and derived an angular diameter distance for 3C 273 of
$D_{A}=551.5^{+97.3}_{-78.7}$ Mpc, corresponding to a Hubble constant
$H_{0}=71.5^{+11.9}_{-10.6}\,{\mbox{$\mathrm{km\,s^{-1}\leavevmode\nobreak\
Mpc^{-1}}$}}$.

The Hubble constant has been measured to an accuracy of a few percent in the
low-$z$ Universe with various methods. Using the distance ladder, the SH0ES
(supernovae H0 for the equation of state) project recently derived
$H_{0}=73.5\pm 1.4\,{\mbox{$\mathrm{km\,s^{-1}\leavevmode\nobreak\
Mpc^{-1}}$}}$ (Reid et al., 2019). This value is consistent with other
measurements independent of the distance ladder, such as megamaser
observations (e.g., Pesce et al., 2020) and observations of gravitationally
lensed quasars (e.g., Wong et al., 2020) based on the so-called time-delay
distance. In comparison, $H_{0}$ inferred from conditions in the early
Universe tends to be substantially smaller than the low-$z$ values (4–6
$\sigma$ significance level, Riess 2020; however, see Freedman et al. 2019).
For instance, Planck Collaboration et al. (2020b) derive $67.4\pm
0.5\,{\mbox{$\mathrm{km\,s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$}}$ based on
cosmic microwave background observations. The high luminosity and remarkably
uniform spectral properties of AGNs have motivated many attempts to use them
as standard candles (Baldwin, 1977; Collier et al., 1999; Elvis & Karovska,
2002; Watson et al., 2011; Wang et al., 2013; Hönig, 2014; La Franca et al.,
2014; Risaliti & Lusso, 2015). In particular, the good correlation between the
BLR radius and continuum luminosity (e.g., Kaspi et al., 2000; Bentz et al.,
2013) makes it promising to probe the $D_{L}$ based on the RM measurements
(e.g., Watson et al., 2011; Wang et al., 2013). However, recent studies reveal
significant deviations from the R–L relation, primarily driven by the increase
in the Eddington ratio (Du et al. 2016, 2018; Du & Wang 2019; Martínez-Aldama
et al. 2019; Dalla Bontà et al. 2020; however, see Grier et al. 2017; Fonseca
Alvarez et al. 2020). It is therefore of great importance to explore the power
of the SARM method, namely to study the BLR structure in detail and to measure
the BLR-based geometric distance with the goal to independently test the
$H_{0}$ tension in the future.

With $z\approx 0.009730$ (Theureau et al., 1998) and a much shorter time lag
than that of 3C 273, NGC 3783 provides the best opportunity to compare the
geometric distance derived with SARM to other independent measurements. A new
RM campaign of NGC 3783 was reported by Bentz et al. (2021) with a time lag of
about 10 days, consistent with previous measurements (Onken & Peterson, 2002).
Gravity Collaboration et al. (2021, hereafter, GC21) reported a BLR mean
angular radius of about 70 $\mu$as, and the RM-measured time lag can be
reproduced with the measured continuum light curve and the best-fit BLR model
inferred only from GRAVITY data at an assumed $D_{A}=38.5$ Mpc. In this paper,
we construct a Bayesian model to fit the GRAVITY and RM data simultaneously
(Section 3). The inferred BLR model is entirely consistent with the results of
GC21 (Section 4). We find that the inferred $D_{A}$ of NGC 3783 is fully
consistent with distances measured using the Tully-Fisher relation and other
indirect methods (Section 5). The application of the SARM method to $H_{0}$ is
discussed in Section 5.4, together with improvements that may come with the
ongoing upgrade of the GRAVITY instrument. This work adopts the following
fiducial parameters for a $\Lambda$CDM cosmology: $\Omega_{m}=0.3$,
$\Omega_{\Lambda}=0.7$, and $H_{0}=70$ km s-1 Mpc-1, unless otherwise
specified.

##  2 GRAVITY interferometer and reverberation mapping observations

GRAVITY interferometric data were collected over 3 years from 2018 to 2020
through a series of Open Time programs and an ESO Large
Programme111Observations were made using the ESO Telescopes at the La Silla
Paranal Observatory, program IDs 0100.B-0582, 0101.B-0255, 0102.B-0667,
2102.B-5053, and 1103.B-0626. with the aim to measure the size of the BLR and
the mass of the central BH. Details of the data reduction and analysis can be
found in GC21. We only briefly summarize the main points here. Phase
referenced with its hot dust continuum emission, NGC 3783 was observed with
MEDIUM spectral resolution ($R\approx 500$) in the science channel and
combined polarization. After the pipeline reduction, we exclude data with
poorer performance in terms of phase reference, keeping those with fringe
tracking ratio $>80\%$. We find that the Br$\gamma$ profiles in the GRAVITY
spectra are consistent considering the calibration uncertainty, which
indicates that the BLR did not change significantly between observations. We
therefore average the differential phase of all of the data into three epochs
according to their $uv$ coordinates. The differential phase signal of the BLR
is dominated by the so-called continuum phase, indicating that the center of
the BLR is offset from the photocenter of the hot dust continuum emission. The
reconstructed image of the NGC 3783 continuum emission displays a secondary
component, which is the main cause of the continuum phase signal, although the
asymmetry of the primary continuum emission can also contribute to the
continuum phase (Gravity Collaboration et al. 2020, hereafter, GC20).
Following GC21, for our analysis we adopt the differential phase data after
subtracting the continuum phase that is calculated from the reconstructed
continuum image. There could be still some residual continuum phase due to the
uncertainty of the spatial origin of the image reconstruction, so we still
allow the BLR to be offset from the continuum reference center in the modeling
(Section 3.1). The amplitude of the differential phase due to BLR rotation is
about 0.2∘ (Figure 8(b) of GC21).

Following GC21, for the model inference we adopt the Br$\gamma$ profile from
our $R\approx 4000$ adaptive optics observations obtained on 20 April 2019
using SINFONI (Eisenhauer et al., 2003; Bonnet et al., 2004). The nuclear
spectrum is extracted using a circular aperture centered on the peak of the
continuum and a diameter of about 0.1 arcsec, which is larger than the field
of view of GRAVITY $\sim 60$ milliarcsec. This results in higher narrow
Br$\gamma$ emission in the SINFONI spectrum; meanwhile, the broad Br$\gamma$
profile from the SINFONI spectrum is consistent with those from GRAVITY
observations. The high spectral resolution of the former guarantees a robust
decomposition of the narrow-line component, while the latter suffer from low
spectral resolution and systematic uncertainties due to the calibration.
Therefore, we prefer to use the SINFONI data in the BLR modeling.

The RM data were collected in the first half of 2020. The details of the
observations and analyses are reported in Bentz et al. (2021). Briefly, the
photometric and spectroscopic monitoring was conducted with the Las Cumbres
Observatory global telescope (LCOGT) network. A total of 209 $V$-band images
and 50 spectra were obtained. The [O III] $\lambda\lambda$4959,5007 doublet
region was used to calibrate the relative flux of the spectra, resulting in a
$\sim 2\%$ accuracy for relative spectrophotometry of H$\beta$ throughout the
monitoring. The emission line light curves were derived from the calibrated
spectra by integrating the emission lines with the local continuum fitted and
subtracted. The continuum light curve at rest-frame 5100 Å is also measured
from the spectra and it is used to cross-calibrate the $V$-band light curve.
Both of them are merged into the final continuum light curve. The RMS spectrum
shows that H$\beta$, He II $\lambda$4686, H$\gamma$, and H$\delta$ are
variable. A time lag of about 10 days is derived from the well calibrated
light curves of the continuum and H$\beta$ line.

Previous studies have found that the higher order Balmer lines tend to show
shorter time lags (e.g., Kaspi et al., 2000; Bentz et al., 2010). High-
ionization lines such as He II $\lambda$4686 also show smaller lags than that
of low-ionization lines such as Balmer lines (e.g., Clavel et al., 1991;
Peterson & Wandel, 1999; Bentz et al., 2010; Grier et al., 2013; Fausnaugh et
al., 2017; Williams et al., 2020). Photoionization models (Netzer, 1975; Rees
et al., 1989; Baldwin et al., 1995; Korista & Goad, 2004) provide physical
explanations for the radial stratification of the BLR. The higher order Balmer
lines, with lower optical depth, are expected to be more efficiently emitted
by gas with higher density and closer to the central BH, resulting in a
shorter time lag and higher “responsivity” (see Korista & Goad, 2004).
However, the optical depth of H$\beta$ is $\gtrsim 10^{3}$ based on typical
BLR models (e.g., Netzer, 2020). It is very challenging for photoionization
models such as CLOUDY (Ferland et al., 2017) to theoretically calculate the
Hydrogen line emission.

Our joint analysis assumes that RM and GRAVITY probe the same regions of the
BLR, which is encouraged by the consistent size measured from the H$\beta$
time lag and GRAVITY measurements of the Br$\gamma$ line. Wang et al. (2020)
conducted a joint analysis using the H$\beta$ light curve and GRAVITY
measurements of the Pa$\alpha$ line. The choice is mainly limited by the
available measurements, although Zhang et al. (2019) show that the time lags
of H$\beta$ and H$\gamma$ for 3C 273 are consistent. It is still unclear
however how well the BLR sizes of Hydrogen lines in the optical and NIR are
consistent with each other. Unfortunately, the H$\gamma$ and H$\delta$ light
curves of NGC 3783 are not robustly calibrated, because they are far from the
reference [O III] lines. Their time lags are 2–4 times smaller than that of
H$\beta$, while their FWHMs are also smaller than that of H$\beta$. This
strongly suggests that the lags of H$\gamma$ and H$\delta$ are not physically
robust. Bentz et al. (2021) also report a $\sim 2$-$\sigma$ detection of the
lag of He II about 2 days; however, the discrepancy of BLR sizes between
H$\beta$ and He II lines are not unexpected. We therefore only adopt the
H$\beta$ light curve together with the continuum light curve in our joint
analysis. One of our main goals in this work is to test whether H$\beta$ and
Br$\gamma$ BLR radii are consistent by comparing our geometric distance with
other independent distance measurements.

##  3 SARM joint analysis

###  3.1 BLR model and spectroastrometry

Our BLR model has been introduced in detail by Gravity Collaboration et al.
(2018) and GC20. Here, we only provide a brief description of the model that
is necessary for this work. The model was first developed by Pancoast et al.
(2014a) with the original purpose to model velocity resolved RM data. The BLR
is assumed to consist of a large number of non-interacting clouds, whose
motion is governed only by the gravity of a central BH with mass
$M_{\mathrm{BH}}$. A shifted gamma distribution,
$r=R_{\mathrm{S}}+F\,R_{\mathrm{BLR}}+g\,(1-F)\,\beta^{2}\,R_{\mathrm{BLR}}$,
is used to describe the radial distribution of the clouds, where
$R_{\mathrm{S}}$ is the Schwarzschild radius, $g=p(x|1/\beta^{2},1)$ is drawn
randomly from a Gamma distribution,
$p(x|a,b)=x^{a-1}e^{-x/b}/(\Gamma(a)\,b^{a})$, and $\Gamma(a)$ is the gamma
function. The shape parameter, $\beta$, controls the radial profile to be
Gaussian ($0<\beta<1$), exponential ($\beta=1$), or heavy-tailed
($1<\beta<2$). The weighted mean cloud radius, $R_{\mathrm{BLR}}$, and
fractional inner radius, $F=R_{\mathrm{in}}/R_{\mathrm{BLR}}$, are fitted.
Clouds are then randomly distributed in a disk with the angular thickness
$\theta_{0}$, which ranges from $0^{\circ}$ (thin disk) to $90^{\circ}$
(sphere). $\gamma$ controls the concentration of the cloud distribution toward
the edge of the disk (Equation (4) of GC20). The weight of the emission from
the clouds is controlled by $\kappa$. With respect to the observer, the near
side clouds have higher weight if $\kappa>0$, while the far side clouds have
higher weight otherwise (Equation (5) of GC20). The fractional difference in
the number of clouds above and below the mid-plane is controlled by $\xi$.
Setting $\xi=0$ means there are equal numbers of clouds each side of the mid-
plane, while there are no clouds below the mid-plane if $\xi=1$.

We define the velocity of the clouds according to their position and the BH
mass. We draw cloud velocities from the parameter space distribution of radial
and tangential velocities centered around either the circular orbits for bound
clouds or around the escape velocity for inflowing or outflowing clouds
(Equation (6) of GC20). The fraction of clouds in bound elliptical orbits is
controlled by $f_{\mathrm{ellip}}$. Clouds in radial orbits (inflowing or
outflowing) are allowed to be mostly bound, mainly controlled by $\theta_{e}$.
A random distribution is drawn with Gaussian dispersion along and
perpendicular to the ellipse connecting circular orbit velocity and radial
escape velocity. Whether the radial motion of the cloud ensemble is inflowing
or outflowing is controlled by a single parameter $f_{\mathrm{flow}}$, $<0.5$
for inflowing and $>0.5$ for outflowing. The model further considers a line-
of-sight velocity dispersion to model the macroturbulence. As in GC21, we find
the dispersion parameters of the Gaussian distributions in the phase space are
not crucial to fit the data. We fix them to zero, so that the BLR model is the
same as that adopted in GC21.

The BLR model is rotated with inclination angle $i$, ranging from $0^{\circ}$
(face-on) to $90^{\circ}$ (edge-on), and position angle PA. Line-of-sight
velocities of the clouds account for the full relativistic Doppler effect and
gravitational redshift. The flux of each spectral channel is calculated by
summing the weights of clouds in each velocity bin. The model line profile is
scaled according to the maximum, $f_{\mathrm{peak}}$. The photocenter of each
channel is the weighted average position of all clouds in each bin. The
differential phase at wavelength $\lambda$ is calculated as

$\Delta\phi_{\lambda}=-2\pi\frac{f_{\lambda}}{1+f_{\lambda}}\@vec{u}\cdot(\@vec{x}_{\mathrm{BLR,\lambda}}-\@vec{x}_{c}),$
(1)

where $f_{\lambda}$ is the line flux at wavelength $\lambda$ to a continuum
level of unity, $\@vec{u}$ is the $uv$ coordinate of the baseline, and
$\@vec{x}_{c}$ is the offset of BLR center from the reference center of the
continuum emission, which can introduce the continuum phase (GC20).

###  3.2 Light curve modeling

In order to generate the model H$\beta$ light curve, we need to calculate the
continuum flux reverberated by clouds with different time lag at a given
observed time. We need, therefore, to interpolate and extrapolate the
continuum light curve taking into account the measurement uncertainty. The
variability of an AGN can be described by a damped random walk model (Kelly et
al., 2009) for which the covariance function between any two times $t_{1}$ and
$t_{2}$ is

$S(t_{1},t_{2})=\sigma_{\mathrm{d}}^{2}\exp\left(-\frac{|t_{1}-t_{2}|}{\tau_{\mathrm{d}}}\right),$
(2)

where $\sigma_{d}$ is the long-term standard deviation and $\tau_{d}$ is the
typical correlated timescale of the continuum light curve. We use a Gaussian
process to model the continuum light curve (e.g., Pancoast et al. 2011, 2014a;
Li et al. 2018). In the fitting, we adopt $\tau_{d}$ and
$\hat{\sigma}_{\mathrm{d}}=\sigma_{\mathrm{d}}/\sqrt{\tau_{\mathrm{d}}}$ in
order to relax the correlation between the two parameters.

We can calculate the model H$\beta$ light curve by convolving the model
continuum light curve $l_{c}(t)$ (Appendix A) with a so-called transfer
function $\Psi(\tau)$,

$\tilde{l}_{\mathrm{H\beta}}(t)=A\int\Psi(\tau)\,l_{c}(t-\tau)\,d\tau,$ (3)

where $A$ is a scaling factor. The transfer function is the normalized
distribution of the time lag taking into account the weight of clouds of the
BLR.

###  3.3 Bayesian inference

Following GC20, we fit the observed Br$\gamma$ profile and differential phase
with the likelihood function of GRAVITY spectroastrometry,

$\displaystyle\mathcal{L}_{\mathrm{SA}}$
$\displaystyle=\prod_{i=1}\frac{1}{\sqrt{2\pi}\sigma_{f,i}}\exp\left(-\frac{(f_{i}-\tilde{f}_{i})^{2}}{2\sigma_{f,i}^{2}}\right)$
$\displaystyle\times\prod_{i=1}\frac{1}{\sqrt{2\pi}\sigma_{\phi,i}}\exp\left(-\frac{(\phi_{i}-\tilde{\phi}_{i})^{2}}{2\sigma_{\phi,i}^{2}}\right),$
(4)

where $f_{i}$ and $\tilde{f}_{i}$ are the observed and model fluxes of the
Br$\gamma$ profile; $\phi_{i}$ and $\tilde{\phi}_{i}$ are the observed and
model differential phases; $\sigma_{f,i}$ and $\sigma_{\phi,i}$ are the
measurement uncertainties of $f_{i}$ and $\phi_{i}$ respectively; and $i$
denotes the $i$th channel. The measured H$\beta$ light curve is fitted with
the likelihood function of RM,

$\mathcal{L}_{\mathrm{RM}}=\prod_{i=1}\frac{1}{\sqrt{2\pi}\sigma_{\mathrm{H\beta},i}}\exp\left(-\frac{(l_{\mathrm{H\beta},i}-\tilde{l}_{\mathrm{H\beta},i})^{2}}{2\sigma_{\mathrm{H\beta},i}^{2}}\right),$
(5)

where $l_{\mathrm{H\beta},i}$ and $\sigma_{\mathrm{H\beta},i}$ are the $i$th
measurements of H$\beta$ flux and uncertainty. Therefore, the joint likelihood
function is the multiplication of Equation (3.3) and Equation (5),
$\mathcal{L}_{\mathrm{SARM}}=\mathcal{L}_{\mathrm{SA}}\times\mathcal{L}_{\mathrm{RM}}$.

The physical parameters of the BLR model are the primary parameters to be
inferred from the joint analysis. Parameters of the continuum light curve
model ($\@vec{x_{s}}$, $x_{q}$, $\tau_{\mathrm{d}}$, and
$\hat{\sigma}_{\mathrm{d}}$) are also involved in the fitting. $\@vec{x_{s}}$
and $x_{q}$ are the deviation of the continuum light curve fluxes and their
long-term average from the maximum a posteriori conditioned by the observed
light curve using the Gaussian process regression (Appendix A). Moreover,
$\@vec{x_{s}}$ consists of 200 parameters because it is important to densely
sample the model continuum light curve. It is however worth emphasizing that
$\@vec{x_{s}}$ is well constrained by the prior information of the observed
continuum data and it only enters the likelihood function via Equation (3). We
tested the fitting using 100–300 points for $\@vec{x_{s}}$. While there is no
significant difference in the fitting results, we find that the reconstructed
continuum light curve with 100 points does not capture some features of the
observed data in some densely sampled region. The goodness of the fitting does
not improve when we adopt 300 points. Therefore, for the sake of computation
power, we adopted 200 points in our analysis. With 221 free parameters in
total, it is very challenging to sample the parameter space. We utilized the
diffusive nested sampling code CDNest (Li, 2018) to do so. Diffusive nested
sampling (Brewer et al., 2011) has been shown to be effective for fitting RM
data with a BLR model (e.g., Pancoast et al. 2014a, b; Li et al. 2018), as
well as for the joint analysis of 3C 273 (Wang et al., 2020).

We find that $\hat{\sigma}_{\mathrm{d}}$ and $\tau_{\mathrm{d}}$ are easily
biased in the joint fitting if no informative prior is used. The covariance
model (Equation (2)) should be able to describe the continuum light curve well
independently. Therefore, we opted to optimize the likelihood function of the
continuum light curve data given the covariance model to constrain
$\hat{\sigma}_{\mathrm{d}}$ and $\tau_{\mathrm{d}}$ in advance of the joint
fitting,

$\mathcal{L}_{c}=\frac{1}{\sqrt{(2\pi)^{m}|\@vec{C_{11}}|}}\times\exp\left(-\frac{(\@vec{l_{c}}-\@vec{E}q)^{T}\@vec{C_{11}}^{-1}(\@vec{l_{c}}-\@vec{E}q)}{2}\right),$
(6)

where
$\@vec{C_{11}}=S(\@vec{t_{1}},\@vec{t_{1}})+\@vec{\sigma_{c}}^{2}\@vec{I}$ is
the model covariance matrix of the measured continuum light curve
($\@vec{l_{c}}$) and uncertainties ($\@vec{\sigma_{c}}$); and the scalar $q$
is the long-term average of the light curve (see Appendix A for more details).
The flat priors with $(-2,2)$ and $(-1,3)$ for
$\log\,(\tau_{\mathrm{d}}/\mathrm{day})$ and
$\log\,(\hat{\sigma}_{\mathrm{d}}/\sqrt{\mathrm{day}})$ are wide enough for
this purpose. We are able to obtain good constraints of
$\log\,(\tau_{\mathrm{d}}/\mathrm{day})=1.75^{+0.87}_{-0.37}$ and
$\log\,(\hat{\sigma}_{\mathrm{d}}/\sqrt{\mathrm{day}})=-0.75^{+0.03}_{-0.03}$.
Therefore, we adopted the priors $\mathrm{Norm(1.75,0.5^{2})}$ and
$\mathrm{Norm(-0.75,0.03^{2})}$ for $\tau_{\mathrm{d}}$ and
$\hat{\sigma}_{\mathrm{d}}$, respectively in logarithmic scale in the joint
fitting. Unlike 3C 273 (Wang et al., 2020; Li et al., 2020), we do not find a
long-term trend in the light curve of NGC 3783. Therefore, a simple constant
$q$ is enough to describe the long-term average of the light curve. We
emphasize that the parameters that we are mainly interested in are the
parameters of the BLR model as well as $D_{A}$, while the parameters of the
light curves are considered as nuisance parameters.

##  4 BLR model inference

Parameters | Model Inference
---|---
$\Theta_{\mathrm{BLR}}\,(\mu\mathrm{as})$ | $71_{-14}^{+22}$
$R_{\mathrm{BLR}}\,(\mathrm{ld})$ | $16.2_{-1.8}^{+2.8}$
$D_{A}$ (Mpc) | $39.9_{-11.9}^{+14.5}$
$\log\,({\mbox{$M_{\mathrm{BH}}$}}/M_{\odot})$ | $7.40_{-0.15}^{+0.13}$
Table 1: Inferred median of posterior sample and central 68% credible interval
for the key parameters of the analysis. $\Theta_{\mathrm{BLR}}$ is the angular
radius of the BLR, which is mainly constrained by the differential phase of
GRAVITY data. $R_{\mathrm{BLR}}$ is the physical radius of the BLR, which is
primarily constrained by RM data. $D_{A}$ is the angular diameter distance.
$M_{\mathrm{BH}}$ is the BH mass.

The inferred model parameters from the joint analysis are entirely consistent
with what we have obtained from fitting GRAVITY data only (GC21). Table 1
presents the key parameters of this work, while the full model parameters are
reported in Table 3. All of the parameters are consistent within their
2-$\sigma$ uncertainty levels. One difference in our approach here is that,
instead of calculating the maximum a posteriori as GC21, we simply report the
median of the marginalized posterior distribution, because the dimension of
the parameter space ($>200$) is prohibitively high to robustly calculate the
maximum a posteriori. Nevertheless, we confirm that the median and maximum a
posteriori of the posterior sampling using only GRAVITY data always show
difference within 1-$\sigma$. The fitting results are discussed in detail in
Appendix B. We would like to highlight that the uncertainties of several key
parameters, in particular $M_{\mathrm{BH}}$, are reduced by a factor of $\sim
2$ after adding the RM light curves in the fitting, which indicates that RM
data help to constrain the model in a consistent manner with GRAVITY data (see
also Wang et al. 2020). Our inferred BH mass, $2.54_{-0.72}^{+0.90}\times
10^{7}\,M_{\odot}$, is fully consistent with the value based on RM-only data
by Bentz et al. (2021). Our uncertainty is about a factor of two larger than
their statistical uncertainty; but one should bear in mind that this does not
include the $\sim 0.3$ dex uncertainty due to the virial factor, which is the
primary uncertainty of integrated RM (e.g., Ho & Kim, 2014; Batiste et al.,
2017). Our best-fit model favors moderate inflow. This is expected because
both the GRAVITY data and RM data indicate a preference for gas inflow
individually. Bentz et al. (2021) derived the time lags of H$\beta$ in 5
velocity bins. The time lag profile across the line is slightly asymmetric
with the longest wavelength showing the shortest time lag, indicating that the
gas motion of the BLR is a combination of rotation and inflow (however, see
Mangham et al. 2019).

Figure 1: Best-fit line profile and differential phase signal of the BLR are
displayed with red curves. The Br$\gamma$ profile from SINFONI spectrum is in
black. The averaged differential phase data of UT4–UT2, UT4–UT1, and UT3-UT1
(see Figure 6) with the continuum phase subtracted are displayed in blue.
Figure 2: Fitting results of (a) continuum and (b) H$\beta$ light curves. We
randomly selected 50 reconstructed continuum light curves and the reverberated
H$\beta$ light curves from the posterior sample and display them in gray. The
median light curves are in red. The data points are plotted in blue. Figure 3:
Posterior probability distribution of the BLR radius, BH mass, and angular
diameter distance of the joint analysis. The dashed lines indicate the 16%,
50%, and 84% percentiles of the of the posterior distributions. The contours
indicate $1\sigma$, $1.5\sigma$, and $2\sigma$.

The line profile and differential phase signal based on the inferred median
parameters are plotted in Figure 1. In this figure, the differential phase is
averaged for UT4–UT2, UT4–UT1, and UT3–UT1. The continuum phase (Section 2) is
subtracted based on the median offset from the data of each baseline before
the average. The data averaged phase displays moderate asymmetry, the excess
in data points with respect to the best-fit model seen on the blue side of the
“S” shape, which can be explained by the uncertainty of those offsets. The
reconstructed continuum light curves and the reverberated H$\beta$ light
curves are shown with the observed data in Figure 2. The detailed features of
the continuum light curve are captured by the model, while the H$\beta$ light
curve is reasonably well fitted. We infer a correlation timescale
$\tau_{\mathrm{d}}\approx 87$ day, slightly longer than the center of the
prior distribution. The mean time lag of the BLR model, based on the median
inferred parameters, is $14.1_{-1.9}^{+2.3}$ day, which is close to the light-
weighted radius of the BLR, $16.2_{-3.5}^{+5.4}$ ld (Figure 3). However, the
centroid of the cross correlation function (CCF) of the model continuum and
H$\beta$ light curves is $\tau_{\mathrm{cent}}=8.2$ day. Following Bentz et
al. (2021), we derive $\tau_{\mathrm{cent}}$ as the first moment of the CCF
above 0.8 of the peak of the CCF (Koratkar & Gaskell, 1991). Our
$\tau_{\mathrm{cent}}$ is slightly lower than, but within 2-$\sigma$ to, that
reported by Bentz et al. (2021). This is consistent with the conclusion in
GC21: The measured time lag from the CCF underestimates the BLR radius of NGC
3783 because the BLR radial distribution is strongly heavy-tailed (large
$\beta$; Figure 4). Nevertheless, detailed comparisons of BLR structures of
Br$\gamma$ and H$\beta$ lines, including the velocity-resolved RM, are needed
in the future. We discuss the fitting results considering an outer truncation
of the BLR as well as the effect of nonlinear response of the continuum light
curve in Appendix B. The primary model introduced above is preferred to
interpret the data. While no disagreement is found statistically significant
by introducing more physical constraints to the model, we caution that the
inferred distance may be biased by the BLR model assumptions.

In order to assess the reliability of the uncertainties on the derived
parameters, we considered the fitting process itself as well as exploring how
sensitive they are to the light curve. It is possible to adjust the
temperature ($T\geq 1$), by which the logarithmic likelihood is divided, with
CDNest in order to enlarge the uncertainty of the data when the model is not
flexible enough to fit the data (Li, 2018; Brewer et al., 2011). We are always
able to find the clear peak of the posterior weights as a function of the
prior volume with $T=1$, indicating that the fitting has properly converged.
The fitting results with $T>1$ are consistent with those with $T=1$, although
the uncertainties increase. Therefore, we always report the fitting results
obtained using $T=1$. The measured uncertainty of the H$\beta$ light curve is
typically $\sim 1\%$ of the flux, which is lower than the typical uncertainty
based on the intercalibration using the [O III] doublet (Bentz et al., 2021).
We fit the data with the H$\beta$ light curve uncertainties increased to be 2%
of the flux if they are smaller than the latter, and find that the results are
fully consistent with those using the measured H$\beta$ uncertainties. The
uncertainty of $M_{\mathrm{BH}}$ reaches 0.2 dex, while that of $D_{A}$ barely
increases as it is dominated by the phase data. In order to further test
whether the joint analysis is sensitive to the light curve measurement, we
measured the H$\beta$ light curve by decomposing the spectra with a relatively
wide wavelength range (see Appendix B.3 for more details). Using the
decomposition method, we measured the light curve with an approximatively 20%
lower variation amplitude. We are able to obtain consistent fitting results
once the nonlinear responsivity is taken into account in the joint analysis.

Figure 4: Transfer function (Equation (3)) of the best-fit model peaks at the
time lag close to $\tau_{\mathrm{cent}}$, measured from the CCF method. But
the mean lag ($\tau_{\mathrm{mean}}$) is longer because the transfer function
has a long tail, which extends beyond the limit of the plot. The inset on the
right displays the cloud distribution of the best-fit BLR model. The color
code is the light-of-sight velocity.

##  5 Distance and peculiar velocity of NGC 3783

Our joint analysis infers that the angular diameter distance of NGC 3783 is
$39.9_{-11.9}^{+14.5}$ Mpc. NGC 3783 has an observed redshift of 0.009730 and
heliocentric velocity $v_{h}\approx 2917\,{\mbox{$\mathrm{km\,s^{-1}}$}}$
(Theureau et al., 1998). In this section, we derive the peculiar velocity of
NGC 3783, compare our measured distance with other direct and indirect
methods, and discuss the uncertainty of our measurements in detail.

###  5.1 Peculiar velocity of NGC 3783

To correct for the motion of the Sun in the Milky Way and the peculiar motion
of the Milky Way, we calculate the velocity of NGC 3783 in the frame of the
Local Sheet (Tully et al., 2008; Kourkchi et al., 2020),

$\displaystyle v_{\mathrm{ls}}$ $\displaystyle=v_{h}-26\cos l\cos b+317\sin
l\cos b-8\sin b$ (7) $\displaystyle=2628\,{\mbox{$\mathrm{km\,s^{-1}}$}},$

where $(l,b)$ are the Galactic longitude and latitude. The Local Sheet
reference frame is a variant of the Local Group rest frame (Yahil et al.,
1977; Karachentsev & Makarov, 1996). Tully et al. (2008) advocated that the
Local Sheet is preferable to the Local Group because it is more stable. For
comparison, NGC 3783 has a Local Group velocity 2627 $\mathrm{km\,s^{-1}}$,
according to NASA/IPAC Extragalactic Database (NED) velocity calculator. The
redshift corrected for the peculiar motion of the observing frame is
$z_{\mathrm{ls}}=v_{\mathrm{ls}}/c\approx 0.008767$, which may still deviate
from the cosmological redshift in the Hubble flow due to the peculiar motion
of NGC 3783. Assuming $(\Omega_{m},\Omega_{\Lambda},H_{0})=(0.3,0.7,70)$, we
can infer a cosmological redshift $\bar{z}=0.0094_{-0.0028}^{+0.0035}$
according to the measured $D_{A}$. The peculiar velocity can be calculated
(Davis & Scrimgeour, 2014),

$v_{\mathrm{pec}}=c\frac{z_{\mathrm{ls}}-\bar{z}}{1+\bar{z}}=-195_{-1029}^{+843}\,{\mbox{$\mathrm{km\,s^{-1}}$}}.$
(8)

This is very close to, albeit with large uncertainties, the peculiar velocity
estimated with the Cosmicflows-3 distance–velocity calculator (see below),
$-182\,{\mbox{$\mathrm{km\,s^{-1}}$}}$. We also estimate the peculiar velocity
using the 6dF galaxy redshift survey peculiar velocity map (Springob et al.,
2014). We find 11 galaxies in the 6dF peculiar velocity catalog within
$8h^{-1}$ Mpc centered on the position of NGC 3783. We averaged their peculiar
velocities and uncertainties with equal weight. The estimated peculiar
velocity of NGC 3783, $-158\pm 43\,{\mbox{$\mathrm{km\,s^{-1}}$}}$, is
consistent with the other estimates above.222We adopt the recession velocity,
3197 $\mathrm{km\,s^{-1}}$ in the CMB frame, of the galaxy group that
comprises NGC 3783 (Kourkchi & Tully, 2017) in order to obtain the 3-D
supergalactic coordinate of NGC 3783. We adopt the typical smoothing scale,
$8h^{-1}$ Mpc, to find the galaxies sharing the same large-scale structure
with NGC 3783. We find it makes little difference if we weight the peculiar
velocity according to their separation or not, so we average them with equal
weight for simplicity. The peculiar velocity has been converted from the CMB
frame (Planck Collaboration et al., 2020a) to the Local Sheet frame in order
to be compared with other estimates.

Distance (Mpc) | Method description
---|---
$39.9_{-11.9}^{+14.5}$ | SARM (this work)
$49.8\pm 19.6$ | Tully-Fisher relation (Robinson et al., 2021)
$38.5\pm 14.2$ | Tully-Fisher relation (Tully & Fisher, 1988)
35.1 | Cosmicflows-3 (EDD; Kourkchi et al., 2020)
37.9 | NED (Virgo + GA + Shapley)
42.1 | Galaxy group (Kourkchi & Tully, 2017)
Table 2: Distance of NGC 3783 measured with different methods. The SARM and
Tully-Fisher relation are the direct methods, while the other methods are
indirect. The distance uncertainties of the indirect methods are discussed in
the text whenever they can be estimated.

###  5.2 Comparison with other distance measures

Besides our measured $D_{A}$, the distance of NGC 3783 can be measured
“directly” with the Tully-Fisher relation (Tully & Fisher, 1977) on the one
hand and estimated “indirectly” with various methods relying on the large-
scale structure and velocity field of the local universe on the other.
Strictly speaking, these methods, as discussed below, provide the luminosity
distance. Since the redshift of NGC 3783 is too low for there to be a
significant difference between angular diameter distance and luminosity
distance ($<2\%$), we do not distinguish the two in the following discussion.
The direct distance measurements and indirect estimates that we discuss below
are summarized in Table 2.

The distance of NGC 3783 has been measured as $38.5\pm 14.2$ Mpc by Tully &
Fisher (1988) using the Tully-Fisher relation. This measurement was based on
the magnitude of the galaxy without removing the nuclear emission from the BH
accretion. This will bias the distance toward a lower value. Robinson et al.
(2021) recently derived the distance of NGC 3783 to be $49.8\pm 19.6$ Mpc.
They derived the maximum rotational velocity measured from the newly observed
H I spectrum and measured the multiband optical/NIR magnitudes of the galaxy
after carefully decomposing the nuclear emission. Although they adopt a 20%
uncertainty of the distance throughout the entire sample, the best-estimate
distance of NGC 3783 is based on the HST photometry as the value quoted above.
Nevertheless, Robinson et al. (2021) emphasized that the distance of NGC 3783
is quite uncertain mainly because the galaxy is rather face-on, leading to a
large uncertainty of the maximum rotational velocity. The strong bar of NGC
3783 may also influence the rotation velocity of the galaxy (e.g.,
Randriamampandry et al., 2015) and affect the distance measured from the
Tully-Fisher relation.

Based on the Cosmicflows-3 catalog (CF3, Tully et al., 2016) of the distance
of $\sim 18000$ galaxies, Graziani et al. (2019) reconstruct the smoothed
peculiar velocity field, for the first time, up to $z\approx 0.05$ using a
linear density field model. They assume a fiducial cosmology
$(\Omega_{m},\Omega_{\Lambda},H_{0})=(0.3,0.7,75)$ in the modeling while
considering deviations of $H_{0}$ from the fiducial value, so that the
reconstructed velocity field does not depend on the assumed $H_{0}$. They find
the cosmic expansion is consistent with their fiducial $H_{0}$ (see also Tully
et al., 2016). We use the CF3 distance–velocity calculator provided by the
Extragalactic Distance Database (EDD)333http://edd.ifa.hawaii.edu/. to
estimate the luminosity distance of NGC 3783 based on its $v_{\mathrm{ls}}$
corrected for cosmological effects (Davis & Scrimgeour, 2014; Kourkchi et al.,
2020),

$\displaystyle v^{c}_{\mathrm{ls}}$
$\displaystyle=v_{\mathrm{ls}}\left[1+\frac{1}{2}(1-q_{0})z_{\mathrm{ls}}-\frac{1}{6}(2-q_{0}-3q_{0}^{2})z_{\mathrm{ls}}^{2}\right],$
(9) $\displaystyle q_{0}$
$\displaystyle=\frac{1}{2}(\Omega_{m}-2\Omega_{\Lambda}).$ (10)

Adopting $\Omega_{m}=0.3$ and $\Omega_{\Lambda}=0.7$, we find $q_{0}=-0.55$
and NGC 3783 has $v^{c}_{\mathrm{ls}}=2646\leavevmode\nobreak\
{\mbox{$\mathrm{km\,s^{-1}}$}}$.444It is worth noting that
$v^{c}_{\mathrm{ls}}$ does not depend on $H_{0}$ (Davis & Scrimgeour, 2014).
And the correction with Equation (9) is $<1\%$ as $z_{\mathrm{ls}}$ is very
small. The corresponding distance of NGC 3783 is 35.1 Mpc according to the CF3
calculator. It is not straightforward to estimate the uncertainty associated
with this method. We expect that the principal source of uncertainty comes
from the peculiar velocity of the galaxy due to the nonlinear effects that
cannot be described by the linear model of Graziani et al. (2019). They
approximate the nonlinear effects in the model with a single parameter of
nonlinear velocity dispersion $\sigma_{\mathrm{NL}}$. Taking the face value of
$\sigma_{\mathrm{NL}}\approx 280\,{\mbox{$\mathrm{km\,s^{-1}}$}}$ from their
Bayesian inference, the linear $v^{c}_{\mathrm{ls}}$ of NGC 3783 is in the
range 2366–2926 $\mathrm{km\,s^{-1}}$. The CF3 calculator yields distances
ranging from 30.1 to 41.8 Mpc. Therefore, we estimate the uncertainty of the
distance from the CF3 calculator is at least 6 Mpc.555Another calculator, NAM,
based on nonlinear model of galaxies within 38 Mpc (or 2850
$\mathrm{km\,s^{-1}}$) is provided by EDD. However, NGC 3783 is just on the
upper boundary of the NAM calculator. It yields $D_{L}=37.5$ Mpc fully
consistent with that of the CF3 calculator, although the NAM result is likely
much more uncertain as the nonlinear model is poorly constrained at the edge
of application (Kourkchi et al., 2020). We therefore prefer the result from
the CF3 calculator.

For completeness, NED provides the estimates of galaxy distances based on the
multiattractor model by Mould et al. (2000). We obtain a Hubble flow velocity
2638 $\mathrm{km\,s^{-1}}$ after correcting for the infall due to Virgo
cluster, Great Attractor, and Shapley supercluster. Taking our fiducial
$H_{0}=70\,{\mbox{$\mathrm{km\,s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$}}$, the
distance of NGC 3783 is 37.9 Mpc, very close to the result of CF3. NGC 3783 is
found in a group of 9 galaxies according to an updated nearby galaxy catalog
by Kourkchi & Tully (2017). The weighted distance using the 2 galaxies with
measured distance from Cosmicflows-3 is $42.1\pm 5.9$ Mpc. Although with
considerably large uncertainty, this estimate provides an indication on the
absolute uncertainty of the indirect approach.

As displayed in Figure 5, our newly measured distance is fully consistent with
the results based on other direct and indirect methods, among which the RMS
scatter is $\sim 5$ Mpc. The uncertainty of our measurement is comparable to
those of the Tully-Fisher relation.

###  5.3 Statistical and systematic uncertainties

The spectral resolution of GRAVITY is moderate and noise of adjacent channels
is likely correlated. This may lead to an underestimate of the inferred
uncertainty. To gauge the effect of the correlated noise, we perform the
analysis in two ways, (1) using only half of GRAVITY differential phase
spectra (every second channel) and (2) rebinning across every two channels in
the differential phase spectra. The inferred parameter uncertainties are
typically $\lesssim 10\%$ higher than those reported in Table 3. Therefore, we
conclude that the effect of the correlated phase noise is moderate.666In
contrast to the 2-$\sigma$ credible interval used in GC20 and GC21, we report
the 1-$\sigma$ uncertainty level throughout this paper for the simplicity of
making fair comparisons to the distance measurements with other methods. The
posterior distributions of most of the key parameters have profiles close to
Gaussian, so the choice of the credible interval does not affect our
conclusions.

The relative uncertainty of our measured distance is about 33%, which is the
combination of the uncertainty of the linear radius ($R_{\mathrm{BLR}}$) and
the angular radius ($\Theta_{\mathrm{BLR}}$) of the BLR. We find
$R_{\mathrm{BLR}}$ has a relatively small uncertainty of about 14%, and is
mainly constrained by the RM data. In contrast, $\Theta_{\mathrm{BLR}}$ has
about 25% uncertainty, mainly due to the relatively large uncertainty of the
differential phase. The differential phase measurements have about 0.1∘
uncertainty per baseline, so the relative uncertainty of the $\sim
0.2{\mbox{${}^{\circ}$}}$ phase signal on three baselines is about 30%.
Therefore, we find
$\delta_{D_{A}}\approx\sqrt{\delta_{R_{\mathrm{BLR}}}^{2}+\delta_{\phi}^{2}}$,
where $\delta_{D_{A}}$, $\delta_{R_{\mathrm{BLR}}}$, and $\delta_{\phi}$ are
the fractional uncertainties of the distance, linear radius of the BLR, and
the differential phase combining 3 baselines. This is overall consistent with
the conclusion of Songsheng et al. (2021) based on tests with mock data.

A primary concern in terms of systematic uncertainty is that the BLR structure
of different lines measured by RM and GRAVITY may be different. A detailed
comparison of H$\beta$ and Br$\gamma$ BLR structures for NGC 3783 is out of
the scope of this work and will be studied in a separate paper where the
velocity-resolved RM modeling will be used. Alternatively, the broad line
profile contains the information of the BLR structure. Following Wang et al.
(2020), we can estimate the relative size difference between the H$\beta$ and
Br$\gamma$ BLRs using $(q-1)/(q+1)$, where
$q=\left(\mathrm{FWHM_{H\beta}/FWHM_{Br\gamma}}\right)^{2}$. We measure
$\mathrm{FWHM}_{\mathrm{Br\gamma}}\approx
3680\,{\mbox{$\mathrm{km\,s^{-1}}$}}$ from our SINFONI spectrum (GC21). As
introduced in more detail in Appendix B.3, we fit globally the H$\beta$
complex of 50 RM spectra individually, including the narrow and broad
components of H$\beta$ and He II $\lambda$4686, [O III]
$\lambda\lambda$4959,5007, and Fe II emission, as well as the power-law
continuum (e.g., Barth et al., 2015; Hu et al., 2015). Accounting for the
instrumental broadening (Bentz et al., 2021), we find the FWHM of H$\beta$ is
stable over the RM campaign and $\mathrm{FWHM}_{\mathrm{H\beta}}=4341\pm
197\,{\mbox{$\mathrm{km\,s^{-1}}$}}$. Our result is slightly smaller than the
$\mathrm{FWHM}_{\mathrm{H\beta}}$ of the mean spectrum reported by Bentz et
al. (2021) mainly due to our different approaches to fit the continuum.
Therefore, the size difference between the H$\beta$ and Br$\gamma$ BLRs can be
about $\sim 15\%$, which is comparable to the 13% RMS of different distance
measurements (Section 5.2). Calculations with photoionization models will also
provide useful insights (Zhang et al., 2021), although special attention needs
to be paid to the difficulty of reproducing the observed flux ratios of
Hydrogen lines in CLOUDY (Netzer, 2020). RM observations of the same line as
that observed with GRAVITY can avoid this problem.

The short time lag of NGC 3783 makes it much more straightforward to compare
the BLR of different lines from RM and GRAVITY than that of 3C 273. The
$\sim$100 day time lag of 3C 273 requires the campaign to be at least several
years. The continuum light curve suffers from a further complication caused by
an overall long-term trend (Li et al., 2020). In fact, the newly measured lag
is about 2 times smaller than the result of early RM campaigns (Kaspi et al.,
2000). The dynamical timescale of the BLR (Peterson, 1993) in NGC 3783 is
$t_{\mathrm{dyn}}\approx R_{\mathrm{BLR}}/v_{\mathrm{FWHM}}\approx 3$ years.
The BLR structure likely varies on a timescale $\gtrsim t_{\mathrm{dyn}}$
(e.g., Peterson 1993; Peterson et al. 2004; Lu et al. 2016). This explains why
we find that the Br$\gamma$ spectra in GRAVITY measurements from 2018 to 2020
show overall consistent line profiles, and that these are also consistent with
the profile from the SINFONI observation in 2019. Although it is not possible
to entirely exclude it, we do not expect significant size variation in the BLR
between GRAVITY and RM observations. Meanwhile, we emphasize that quasi-
simultaneous observation of GRAVITY and RM within $t_{\mathrm{dyn}}$ is
necessary to avoid issues associated with the variability of the BLR
structure.

Figure 5: Comparison of the distances of NGC 3783 measured with the methods
listed in Table 2. The 1-$\sigma$ uncertainties of direct measurements are
plotted, while those of indirect measurements are not likely smaller than the
former. We discuss the uncertainties of the indirect measurements in the text
whenever they can be estimated.

###  5.4 Toward a new estimate of the Hubble constant

With a measurement of $D_{A}$, one natural step forward is to estimate the
Hubble constant. We derive
$H_{0}=v^{c}_{\mathrm{ls}}/\left[D_{A}(1+z_{\mathrm{ls}})^{2}\right]=65_{-17}^{+28}\,{\mbox{$\mathrm{km\,s^{-1}}$}}\,\mathrm{Mpc}^{-1}$
from the joint analysis of NGC 3783. The peculiar velocity of NGC 3783 is not
included in this estimate. Our estimate (Sec. 5.1) suggests the additional
systematic uncertainty introduced by the peculiar velocity is $\lesssim 10\%$.
At this level, the differential phase errors dominate the uncertainty of
$H_{0}$. With upgrades to GRAVITY including adaptive optics with laser guide
star and wide angle off-axis phase referencing, GRAVITY+ will allow us to
observe at least several hundred suitable Type 1 AGN targets up to redshift
2–3, with Pa$\alpha$, Pa$\beta$, Pa$\gamma$, H$\alpha$, and H$\beta$
redshifted into $K$ band. These lines are much stronger than Br$\gamma$ with
respect to the continuum. Therefore, we expect much less statistical
uncertainty (e.g., $\lesssim 10\%$) on $D_{A}$ with future observations.
Assuming the same sensitivity level as the current GRAVITY performance,
Songsheng et al. (2021) show that it might in principle be possible to obtain
a precision measurement of $H_{0}$ using a large sample of differential phase
and RM measurements. In practice, more effort would be necessary to address
the potential bias due to the properties of different hydrogen lines and the
assumptions of the BLR model.

##  6 Conclusions

NGC 3783 is the second AGN observed by both GRAVITY interferometry and RM
campaigns. We fitted both data sets simultaneously with a SARM joint analysis.
The inferred model parameters are fully consistent with the previous study
using only GRAVITY data, and the uncertainties of key model parameters are
significantly reduced. In particular, the inferred BH mass is
$2.54_{-0.72}^{+0.90}\times 10^{7}\,M_{\odot}$ (68% credible interval). For
this parameter, the uncertainty is a factor of two smaller than both GRAVITY-
only inference and RM measurements using the integrated H$\beta$ light curve
(dominated by the 0.3 dex calibration uncertainty of the virial factor). We
also confirm the previous finding of GC21 that the BLR is highly concentrated
with an extended tail of clouds out to large radius, which leads to an
apparent discrepancy between the mean radius of the BLR and the measured time
lag.

With the joint analysis, we are able to constrain the angular diameter
distance of NGC 3783 to be $39.9^{+14.5}_{-11.9}$ Mpc. Because NGC 3783 is in
the nearby Universe, this distance can be compared to other independent direct
and indirect measurements based on the Tully-Fisher relation, galaxy flow
models, and its galaxy group. The BLR-based geometric distance is fully
consistent with these other results. The dominant uncertainty for the distance
comes from the differential phase, which has a relative uncertainty of about
30%, while the relative uncertainty from RM data is about 14%.

With the ongoing upgrade of GRAVITY, we expect to observe many AGNs with
improved sensitivity. Our analysis indicates that the uncertainties of BLR
model parameters, such as $\Theta_{\mathrm{BLR}}$ and $M_{\mathrm{BH}}$, will
be reduced as the phase uncertainty decreases. Being able to substantially
reduce the statistical uncertainty of $H_{0}$ looks promising by both
improving the precision of individual measurements and averaging the
measurements over many targets. In this context, it will be important to
better understand the systematic uncertainties associated with the SARM
method. Further, measuring the same broad line in $K$ band with both RM and
GRAVITY will be important for future study.

###### Acknowledgements.

We thank the referees for their careful reading of the manuscript and their
suggestions that have helped to improve the clarity of it. We thank Bradley
Peterson for his helpful comments on this paper. J.S. would like to thank the
important help from Yan-Rong Li and Yu-Yang Songsheng to implement and test
the joint analysis model. He also thanks the useful discussion from Jian-Min
Wang and the IHEP group. J.D. was supported in part by NSF grant AST 1909711
and an Alfred P. Sloan Research Fellowship. M.C.B. gratefully acknowledges
support from the NSF through grant AST-2009230 to Georgia State University.
A.A. and P.G. were supported by Fundação para a Ciência e a Tecnologia, with
grants reference UIDB/00099/2020, PTDC/FIS-AST/7002/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. P.O.P. acknowledges
financial support from the CNRS « Programme national des hautes énergies » and
from the french space agency CNES. 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,777http://www.astropy.org a community-developed core Python package
for Astronomy (Astropy Collaboration et al., 2013, 2018), numpy (Van Der Walt
et al., 2011), scipy (Jones et al., 2001), and matplotlib (Hunter, 2007)

## References

  * Anderson et al. (2014) Anderson, L., Aubourg, É., Bailey, S., et al. 2014, MNRAS, 441, 24
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Bailey (1998) Bailey, J. A. 1998, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 3355, Optical Astronomical Instrumentation, ed. S. D’Odorico, 932–939
  * Baldwin et al. (1995) Baldwin, J., Ferland, G., Korista, K., & Verner, D. 1995, ApJ, 455, L119
  * Baldwin (1977) Baldwin, J. A. 1977, ApJ, 214, 679
  * Barth et al. (2015) Barth, A. J., Bennert, V. N., Canalizo, G., et al. 2015, ApJS, 217, 26
  * Batiste et al. (2017) Batiste, M., Bentz, M. C., Raimundo, S. I., Vestergaard, M., & Onken, C. A. 2017, ApJ, 838, L10
  * Bentz et al. (2013) Bentz, M. C., Denney, K. D., Grier, C. J., et al. 2013, ApJ, 767, 149
  * Bentz et al. (2021) Bentz, M. C., Street, R., Onken, C. A., & Valluri, M. 2021, ApJ, 906, 50
  * Bentz et al. (2010) Bentz, M. C., Walsh, J. L., Barth, A. J., et al. 2010, ApJ, 716, 993
  * Bessel (1838) Bessel, F. W. 1838, MNRAS, 4, 152
  * Bhardwaj (2020) Bhardwaj, A. 2020, Journal of Astrophysics and Astronomy, 41, 23
  * Bonnet et al. (2004) Bonnet, H., Conzelmann, R., Delabre, B., et al. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5490, Advancements in Adaptive Optics, ed. D. Bonaccini Calia, B. L. Ellerbroek, & R. Ragazzoni, 130–138
  * Boroson & Green (1992) Boroson, T. A. & Green, R. F. 1992, ApJS, 80, 109
  * Braatz et al. (2010) Braatz, J. A., Reid, M. J., Humphreys, E. M. L., et al. 2010, ApJ, 718, 657
  * Brewer et al. (2011) Brewer, B. J., Pártay, L. B., & Csányi, G. 2011, Statistics and Computing, 21, 649
  * Clavel et al. (1991) Clavel, J., Reichert, G. A., Alloin, D., et al. 1991, ApJ, 366, 64
  * Collier et al. (1999) Collier, S., Horne, K., Wanders, I., & Peterson, B. M. 1999, MNRAS, 302, L24
  * Dalla Bontà et al. (2020) Dalla Bontà, E., Peterson, B. M., Bentz, M. C., et al. 2020, ApJ, 903, 112
  * Davis & Scrimgeour (2014) Davis, T. M. & Scrimgeour, M. I. 2014, MNRAS, 442, 1117
  * Du et al. (2016) Du, P., Lu, K.-X., Zhang, Z.-X., et al. 2016, ApJ, 825, 126
  * Du & Wang (2019) Du, P. & Wang, J.-M. 2019, ApJ, 886, 42
  * Du et al. (2018) Du, P., Zhang, Z.-X., Wang, K., et al. 2018, ApJ, 856, 6
  * Eisenhauer et al. (2003) Eisenhauer, F., Abuter, R., Bickert, K., et al. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, ed. M. Iye & A. F. M. Moorwood, 1548–1561
  * Eisenstein et al. (2005) Eisenstein, D. J., Zehavi, I., Hogg, D. W., et al. 2005, ApJ, 633, 560
  * Elvis & Karovska (2002) Elvis, M. & Karovska, M. 2002, ApJ, 581, L67
  * ESA (1997) ESA. 1997, in ESA Special Publication, Vol. 1200, ESA Special Publication
  * Fausnaugh et al. (2017) Fausnaugh, M. M., Grier, C. J., Bentz, M. C., et al. 2017, ApJ, 840, 97
  * Ferland et al. (2017) Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mexicana Astron. Astrofis., 53, 385
  * Fonseca Alvarez et al. (2020) Fonseca Alvarez, G., Trump, J. R., Homayouni, Y., et al. 2020, ApJ, 899, 73
  * Freedman et al. (2019) Freedman, W. L., Madore, B. F., Hatt, D., et al. 2019, ApJ, 882, 34
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Gaskell & Sparke (1986) Gaskell, C. M. & Sparke, L. S. 1986, ApJ, 305, 175
  * Glass (1992) Glass, I. S. 1992, MNRAS, 256, 23P
  * Goad & Korista (2014) Goad, M. R. & Korista, K. T. 2014, MNRAS, 444, 43
  * Gravity Collaboration et al. (2017) Gravity Collaboration, Abuter, R., Accardo, M., et al. 2017, A&A, 602, A94
  * Gravity Collaboration et al. (2019) Gravity Collaboration, Abuter, R., Amorim, A., et al. 2019, A&A, 625, L10
  * Gravity Collaboration et al. (2021) Gravity Collaboration, Amorim, A., Bauböck, M., et al. 2021, A&A, 648, A117 (GC21)
  * Gravity Collaboration et al. (2020) Gravity Collaboration, Amorim, A., Bauböck, M., et al. 2020, A&A, 643, A154 (GC20)
  * Gravity Collaboration et al. (2018) Gravity Collaboration, Sturm, E., Dexter, J., et al. 2018, Nature, 563, 657
  * Graziani et al. (2019) Graziani, R., Courtois, H. M., Lavaux, G., et al. 2019, MNRAS, 488, 5438
  * Grier et al. (2013) Grier, C. J., Peterson, B. M., Horne, K., et al. 2013, ApJ, 764, 47
  * Grier et al. (2012) Grier, C. J., Peterson, B. M., Pogge, R. W., et al. 2012, ApJ, 755, 60
  * Grier et al. (2017) Grier, C. J., Trump, J. R., Shen, Y., et al. 2017, ApJ, 851, 21
  * Herrnstein et al. (1999) Herrnstein, J. R., Moran, J. M., Greenhill, L. J., et al. 1999, Nature, 400, 539
  * Ho & Kim (2014) Ho, L. C. & Kim, M. 2014, ApJ, 789, 17
  * Hönig (2014) Hönig, S. F. 2014, ApJ, 784, L4
  * Hönig et al. (2014) Hönig, S. F., Watson, D., Kishimoto, M., & Hjorth, J. 2014, Nature, 515, 528
  * Hu et al. (2015) Hu, C., Du, P., Lu, K.-X., et al. 2015, ApJ, 804, 138
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python
  * Karachentsev & Makarov (1996) Karachentsev, I. D. & Makarov, D. A. 1996, AJ, 111, 794
  * Kaspi et al. (2000) Kaspi, S., Smith, P. S., Netzer, H., et al. 2000, ApJ, 533, 631
  * Kelly et al. (2009) Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895
  * Koratkar & Gaskell (1991) Koratkar, A. P. & Gaskell, C. M. 1991, ApJS, 75, 719
  * Korista & Goad (2004) Korista, K. T. & Goad, M. R. 2004, ApJ, 606, 749
  * Kourkchi et al. (2020) Kourkchi, E., Courtois, H. M., Graziani, R., et al. 2020, AJ, 159, 67
  * Kourkchi & Tully (2017) Kourkchi, E. & Tully, R. B. 2017, ApJ, 843, 16
  * La Franca et al. (2014) La Franca, F., Bianchi, S., Ponti, G., Branchini, E., & Matt, G. 2014, ApJ, 787, L12
  * Lacy (1977) Lacy, C. H. 1977, ApJ, 213, 458
  * Li (2018) Li, Y.-R. 2018, CDNest: A diffusive nested sampling code in C
  * Li et al. (2018) Li, Y.-R., Songsheng, Y.-Y., Qiu, J., et al. 2018, ApJ, 869, 137
  * Li et al. (2013) Li, Y.-R., Wang, J.-M., Ho, L. C., Du, P., & Bai, J.-M. 2013, ApJ, 779, 110
  * Li et al. (2020) Li, Y.-R., Zhang, Z.-X., Jin, C., et al. 2020, ApJ, 897, 18
  * Lira et al. (2011) Lira, P., Arévalo, P., Uttley, P., McHardy, I., & Breedt, E. 2011, MNRAS, 415, 1290
  * Lu et al. (2016) Lu, K.-X., Du, P., Hu, C., et al. 2016, ApJ, 827, 118
  * Mangham et al. (2019) Mangham, S. W., Knigge, C., Williams, P., et al. 2019, MNRAS, 488, 2780
  * Martínez-Aldama et al. (2019) Martínez-Aldama, M. L., Czerny, B., Kawka, D., et al. 2019, ApJ, 883, 170
  * Mould et al. (2000) Mould, J. R., Huchra, J. P., Freedman, W. L., et al. 2000, ApJ, 529, 786
  * Netzer (1975) Netzer, H. 1975, MNRAS, 171, 395
  * Netzer (2020) Netzer, H. 2020, MNRAS, 494, 1611
  * Onken & Peterson (2002) Onken, C. A. & Peterson, B. M. 2002, ApJ, 572, 746
  * Pancoast et al. (2011) Pancoast, A., Brewer, B. J., & Treu, T. 2011, ApJ, 730, 139
  * Pancoast et al. (2014a) Pancoast, A., Brewer, B. J., & Treu, T. 2014a, MNRAS, 445, 3055
  * Pancoast et al. (2014b) Pancoast, A., Brewer, B. J., Treu, T., et al. 2014b, MNRAS, 445, 3073
  * Pesce et al. (2020) Pesce, D. W., Braatz, J. A., Reid, M. J., et al. 2020, ApJ, 891, L1
  * Peterson (1993) Peterson, B. M. 1993, PASP, 105, 247
  * Peterson (2014) Peterson, B. M. 2014, Space Sci. Rev., 183, 253
  * Peterson et al. (2004) Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004, ApJ, 613, 682
  * Peterson & Wandel (1999) Peterson, B. M. & Wandel, A. 1999, ApJ, 521, L95
  * Peterson et al. (1998) Peterson, B. M., Wanders, I., Bertram, R., et al. 1998, ApJ, 501, 82
  * Petrov et al. (2001) Petrov, R. G., Malbet, F., Richichi, A., et al. 2001, Comptes Rendus Physique, 2, 67
  * Phillips (1993) Phillips, M. M. 1993, ApJ, 413, L105
  * Pietrzyński et al. (2019) Pietrzyński, G., Graczyk, D., Gallenne, A., et al. 2019, Nature, 567, 200
  * Planck Collaboration et al. (2020a) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020a, A&A, 641, A1
  * Planck Collaboration et al. (2020b) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020b, A&A, 641, A6
  * Randriamampandry et al. (2015) Randriamampandry, T. H., Combes, F., Carignan, C., & Deg, N. 2015, MNRAS, 454, 3743
  * Rees et al. (1989) Rees, M. J., Netzer, H., & Ferland, G. J. 1989, ApJ, 347, 640
  * Reid et al. (2019) Reid, M. J., Pesce, D. W., & Riess, A. G. 2019, ApJ, 886, L27
  * Riess (2020) Riess, A. G. 2020, Nature Reviews Physics, 2, 10
  * Riess et al. (2021) Riess, A. G., Casertano, S., Yuan, W., et al. 2021, ApJ, 908, L6
  * Riess et al. (2009) Riess, A. G., Macri, L., Casertano, S., et al. 2009, ApJ, 699, 539
  * Riess et al. (1996) Riess, A. G., Press, W. H., & Kirshner, R. P. 1996, ApJ, 473, 88
  * Risaliti & Lusso (2015) Risaliti, G. & Lusso, E. 2015, ApJ, 815, 33
  * Robinson et al. (2021) Robinson, J. H., Bentz, M. C., Courtois, H. M., et al. 2021, ApJ, 912, 160
  * Rybicki & Press (1992) Rybicki, G. B. & Press, W. H. 1992, ApJ, 398, 169
  * Shen et al. (2016) Shen, Y., Horne, K., Grier, C. J., et al. 2016, ApJ, 818, 30
  * Songsheng et al. (2021) Songsheng, Y.-Y., Li, Y.-R., Du, P., & Wang, J.-M. 2021, ApJS, 253, 57
  * Springob et al. (2014) Springob, C. M., Magoulas, C., Colless, M., et al. 2014, MNRAS, 445, 2677
  * Theureau et al. (1998) Theureau, G., Bottinelli, L., Coudreau-Durand, N., et al. 1998, A&AS, 130, 333
  * Tully et al. (2016) Tully, R. B., Courtois, H. M., & Sorce, J. G. 2016, AJ, 152, 50
  * Tully & Fisher (1977) Tully, R. B. & Fisher, J. R. 1977, A&A, 500, 105
  * Tully & Fisher (1988) Tully, R. B. & Fisher, J. R. 1988, Catalog of Nearby Galaxies
  * Tully et al. (2008) Tully, R. B., Shaya, E. J., Karachentsev, I. D., et al. 2008, ApJ, 676, 184
  * van der Marel & Franx (1993) van der Marel, R. P. & Franx, M. 1993, ApJ, 407, 525
  * Van Der Walt et al. (2011) Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  * Wang et al. (2013) Wang, J.-M., Du, P., Valls-Gabaud, D., Hu, C., & Netzer, H. 2013, Phys. Rev. Lett., 110, 081301
  * Wang et al. (2020) Wang, J.-M., Songsheng, Y.-Y., Li, Y.-R., Du, P., & Zhang, Z.-X. 2020, Nature Astronomy, 4, 517
  * Watson et al. (2011) Watson, D., Denney, K. D., Vestergaard, M., & Davis, T. M. 2011, ApJ, 740, L49
  * Williams et al. (2020) Williams, P. R., Pancoast, A., Treu, T., et al. 2020, ApJ, 902, 74
  * Woillez et al. (2004) Woillez, J. M., Sol, H., Perrin, G. S., & Lai, O. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5491, New Frontiers in Stellar Interferometry, ed. W. A. Traub, 97
  * Wong et al. (2020) Wong, K. C., Suyu, S. H., Chen, G. C. F., et al. 2020, MNRAS, 498, 1420
  * Yahil et al. (1977) Yahil, A., Tammann, G. A., & Sandage, A. 1977, ApJ, 217, 903
  * Zhang et al. (2021) Zhang, X., He, Z., Wang, T., & Guo, H. 2021, ApJ, 914, 143
  * Zhang et al. (2019) Zhang, Z.-X., Du, P., Smith, P. S., et al. 2019, ApJ, 876, 49

##  Appendix A The model of the continuum light curve

The observed continuum light curve can be described by (Rybicki & Press 1992;
Li et al. 2018),

$\@vec{l_{c}}=\@vec{s}+\@vec{E}q+\@vec{n_{c}},$ (11)

where the vector $\@vec{s}$ is the variable component of the light curve with
mean zero, $\@vec{E}q$ represents the mean flux of the light curve, and
$\@vec{n_{c}}$ is the noise of the observation. $\@vec{E}$ is a vector of
unity with the same length of the observed light curve and $q$ is the mean of
the light curve. Assuming that the light curve is a multivariate Gaussian with
covariance described by Equation (2) and the measurement uncertainties are
Gaussian and uncorrelated, we can interpolate and extrapolate the light curve
with a uniformly and densely sampled continuum light curve model. The observed
light curve is sampled in a sparse uneven time series, $\@vec{t_{1}}$,
meanwhile the model light curve is densely sampled in a uniform time series,
$\@vec{t_{2}}$. The posterior mean and covariance of the variable component of
the model light curve can be calculated,

$\displaystyle\hat{\@vec{s}}$
$\displaystyle=(\@vec{C_{11}}^{-1}\@vec{C_{12}})^{T}(\@vec{l_{c}}-\@vec{E}q),$
(12) $\displaystyle\@vec{C_{s}}$
$\displaystyle=\@vec{C_{22}}-(\@vec{C_{11}}^{-1}\@vec{C_{12}})^{T}\@vec{C_{12}},$
(13)

where
$\@vec{C_{11}}=S(\@vec{t_{1}},\@vec{t_{1}})+\@vec{\sigma_{c}}^{2}\@vec{I}$,
$\@vec{C_{12}}=S(\@vec{t_{1}},\@vec{t_{2}})$, and
$\@vec{C_{22}}=S(\@vec{t_{2}},\@vec{t_{2}})$ are matrices calculated based on
Equation (2), and superscript “T” denotes transposition. We denote the
observation uncertainties as $\@vec{\sigma_{c}}$ and the identity matrix as
$\@vec{I}$. The mean and variance of the long-term average of the model light
curve are,

$\displaystyle\hat{q}$
$\displaystyle=C_{q}\@vec{E}^{T}\@vec{C_{11}}^{-1}\@vec{l_{c}},$ (14)
$\displaystyle C_{q}$
$\displaystyle=(\@vec{E}^{T}\@vec{C_{11}}^{-1}\@vec{E})^{-1}.$ (15)

Therefore, the model continuum light curve is,

$\@vec{\tilde{l}_{c}}=(\@vec{L}\@vec{x_{s}}+\hat{\@vec{s}})+\@vec{E}(\sqrt{C_{q}}x_{q}+\hat{q}),$
(16)

where $\@vec{L}$ is the lower triangular matrix of $\@vec{C_{s}}$, so that
$\@vec{C_{s}}=\@vec{L}\@vec{L}^{T}$, $\@vec{x_{s}}$ describes the deviation of
the light curve variable component from its posterior mean values, and $x_{q}$
describes the deviation of the long-term mean from its posterior mean value.
We fit $\@vec{x_{s}}$ and $x_{q}$ in the joint analysis as free parameters for
the continuum light curve model.

##  Appendix B Model fitting results and further tests

Figure 6: 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 best-fit residual continuum phases are shown with dotted lines. Figure 7: Posterior distribution of physical parameters from the Bayesian joint analysis. Parameters | Model Inference
---|---
$\Theta_{\mathrm{BLR}}\,(\mu\mathrm{as})$ | $71_{-13.8}^{+21.9}$
$R_{\mathrm{BLR}}\,(\mathrm{ld})$ | $16.2_{-1.8}^{+2.8}$
$\Theta_{\mathrm{min}}\,(\mu\mathrm{as})$ | $31_{-7}^{+11}$
$\beta$ | $1.95_{-0.08}^{+0.04}$
$\theta_{\mathrm{o}}\,(^{\circ})$ | $38_{-5}^{+5}$
$i\,(^{\circ})$ | $32_{-4}^{+4}$
PA (∘ E of N) | $302_{-19}^{+20}$
$\kappa$ | $-0.03_{-0.14}^{+0.19}$
$\gamma$ | $1.7_{-0.5}^{+1.1}$
$\xi$ | $0.5_{-0.3}^{+0.4}$
Offset ($\mu$as) | $\left(-1.1_{-8.4}^{+7.2},-21.6_{-6.9}^{+7.4}\right)$
$\log\,({\mbox{$M_{\mathrm{BH}}$}}/M_{\odot})$ | $7.40_{-0.15}^{+0.13}$
$f_{\mathrm{ellip}}$ | $0.65_{-0.23}^{+0.19}$
$P$(inflow) | 0.79
$\theta_{\mathrm{e}}\,(^{\circ})$ | $32_{-21}^{+28}$
$\Delta v_{\mathrm{BLR}}\,(\mathrm{km\,s^{-1}})$ | $328_{-97}^{+81}$
$D_{A}$ (Mpc) | $39.9_{-11.9}^{+14.5}$
$\log\,(\tau_{\mathrm{d}}/\mathrm{day})$ | $1.94_{-0.35}^{+0.42}$
$\log\,(\hat{\sigma}_{\mathrm{d}}/\sqrt{\mathrm{day}})$ | $-0.67_{-0.03}^{+0.03}$
$\chi_{\mathrm{r}}^{2}$ | 0.672
Table 3: Posterior median and 68% credible intervals for the BLR model
parameters. The central offset of the BLR is in (R. A., Dec.). $P$(inflow) is
the probability of inflow where $f_{\mathrm{flow}}<0.5$ in the posterior
sample. $\Delta v_{\mathrm{BLR}}$ is the difference between the velocity
derived from the best-fit $\lambda_{\mathrm{emit}}$ and the systemic velocity
based on the redshift.

As discussed in Section 2, we use the differential phase data after
subtracting the continuum phase primarily due to the secondary continuum
emission component according to the reconstructed image (GC21). In order to
display the goodness of the fit, we display the differential phases of
UT4–UT2, UT4–UT1, and UT3–UT1, averaged over the three epochs, in Figure 6
together with the best-fit models. The phase signal is dominated by the BLR
component while the residual continuum phase is moderate. The BLR phase signal
is not significantly detected in the other three short baselines.

The posterior probability distributions of all of the physical parameters are
shown in Figure 7. As discussed in Section 4, the posterior distributions of
the BLR model of the joint analysis show very good consistency with those
derived using GRAVITY data alone (GC21). The uncertainties of several key
parameters, such as $\beta$, $\theta_{\mathrm{o}}$, $i$, and
$M_{\mathrm{BH}}$, are reduced by a factor of $\sim 2$, while those of the
other parameters stay the same. The reduction of uncertainties of the inferred
parameters is also found in the joint analysis of 3C 273 (Wang et al. 2020).
This indicates that RM data help to constrain the model in a consistent manner
with GRAVITY data. The $\beta$ parameter is closer to the upper boundary than
GRAVITY only inference, reinforcing our previous finding that the BLR is
heavy-tailed. With GRAVITY data only, we found $\xi$ peaks at $\sim 1$ with a
long tail to $\sim 0$, indicating that $\xi$ is not strongly constrained. With
the joint analysis, $\xi$ is even less constrained with a nearly flat
posterior distribution (see Figure 7). This means that the fitting is not
sensitive to the midplane obscuration, mainly because most of the BLR clouds
($65\%$) are in elliptical orbits rather than radial orbits. The best-fit
model favors moderate inflow, which is consistent with Bentz et al. (2021).

###  B.1 Constraining the BLR outer radius

We further test whether the fitting results are significantly different when
we constrain the BLR inner and outer radii based on previous RM of high-
ionization lines and NIR continuum. Following the discussion of GC21, we cast
additional prior information that $R_{\mathrm{min}}=4$ ld and
$R_{\mathrm{out}}=80$ ld. The key results are shown in Table 4, where $D_{A}$
is reduced by about 30%. The angular size of the BLR is almost the same as our
primary model without additional constraints, as it is directly constrained by
the GRAVITY phase data. At the same time, $R_{\mathrm{BLR}}$ is reduced by
about 30% due to the hard boundary $R_{\mathrm{out}}=80$ ld of the cloud
distribution. Given the observed angular size of the BLR, a smaller $D_{A}$ is
required in order to incorporate a smaller BLR as required by the prior. In
addition, the inferred $\beta$ is around 1, meaning that the cloud
distribution is now exponential rather than heavy-tailed. This is because
$R_{\mathrm{out}}$ limits the extent of the cloud distribution. The peak of
the model differential phase is correspondingly reduced, although the model
fits the data almost as well as our primary model given the data uncertainty.
All of the other inferred model parameters are consistent with our primary
model results. Nevertheless, the radius-constrained model is not preferred
over our primary model because the inferred $D_{A}$ would be the smallest of
the measurements discussed in Section 5. Quantitatively, it is about
1.5-$\sigma$ lower than the other measurements, which cautions against the
additional constraints to the model. While $R_{\mathrm{out}}=80$ ld is a
reasonable choice based on NIR continuum RM analyses (Glass 1992; Lira et al.
2011), the contribution of NIR emission from the accretion disk may bias the
time lag to a small value (Lira et al. 2011). Nevertheless, the joint analysis
may also be biased by the assumptions of the BLR model. More comprehensive
comparisons will be performed in future paper that considers different model
assumptions, e.g., different cloud radial distributions and BLR size
constraints.

Parameters | Model Inference with
---|---
(1) Boundary | (2) Nonlinearity | (3) New light curve
$\Theta_{\mathrm{BLR}}\,(\mu\mathrm{as})$ | $71_{-14}^{+19}$ | $70_{-16}^{+27}$ | $74_{-16}^{+21}$
$R_{\mathrm{BLR}}\,(\mathrm{ld})$ | $12.7_{-1.6}^{+1.8}$ | $18.9_{-4.9}^{+4.6}$ | $16.4_{-5.6}^{+6.5}$
$D_{A}$ (Mpc) | $30.9_{-7.0}^{+8.6}$ | $43.6_{-17.8}^{+24.7}$ | $36.1_{-12.6}^{+21.6}$
$\log\,({\mbox{$M_{\mathrm{BH}}$}}/M_{\odot})$ | $7.51_{-0.17}^{+0.31}$ | $7.40_{-0.19}^{+0.14}$ | $7.50_{-0.19}^{+0.23}$
Table 4: Same as Table 1, but showing the fitting results: (1) when
$R_{\mathrm{min}}=4$ ld and $R_{\mathrm{max}}=80$ ld are required (Appendix
B.1); (2) when the nonlinear response is considered (Appendix B.2); (3) when
we use the light curve measured by the decomposition method (Appendix B.3).

###  B.2 The nonlinear response of RM

The BLR reverberates the Hydrogen ionizing (UV) photon variability from the
accretion disk. The optical (i.e., $V$ band) light curve is usually used as an
approximation to the UV light curve. However, their correlation is not
necessarily linear. Following (Li et al. 2018, see also Li et al. 2013), we
revise Equation (3) into

$\tilde{l}_{\mathrm{H\beta}}(t)=A\int\Psi(\tau)\,l_{c}^{1+\delta_{c}}(t-\tau)\,d\tau,$
(17)

where $\delta_{c}$ captures the nonlinearity of the response. The model
inference is conducted in the same way for the other parameters as discussed
in Section 3, with the uniform prior of $\delta_{c}$ between -1 and 3 (Li et
al. 2018). We find the inferred model parameters are entirely consistent with
our primary fitting results (Table 3). Some of the key parameters are listed
in Table 4. Meanwhile, the uncertainty of $R_{\mathrm{BLR}}$ increases by
about 70%, resulting in an increase to the uncertainty of $D_{A}$ by the same
amount. This mainly originates from the degeneracy between $\delta_{c}$ and
the scaling factor $A$ (Equation (17)). Optimizing the sampling method will
reduce the uncertainties of some parameters due to the degeneracy. However,
this is beyond the scope of the current paper. On the other hand, the fitting
of the light curve is not substantially improved. Therefore, we opt not to
include the nonlinear effect of the RM response in the modeling.

###  B.3 Measuring the light curve using the decomposition method

Here we test whether the joint analysis is sensitive to the method measurement
the light curve. Bentz et al. (2021) measured the H$\beta$ flux by fitting a
local linear continuum underneath the H$\beta$ line. This method has been
widely used and has been shown to be effective in producing accurate light
curves for strong emission lines, such as H$\beta$ (e.g., Peterson et al.
1998; Kaspi et al. 2000; Grier et al. 2012). Recent studies have developed
methods to decompose the emission lines by simultaneously fitting various
physical components to the spectra over a wide wavelength range (Barth et al.
2015; Hu et al. 2015; Shen et al. 2016). We therefore conducted a joint fit
using the H$\beta$ light curve measured with the decomposition method in order
to test whether the joint analysis is sensitive to the measurement of the
light curve. Following the approach of Barth et al. (2015), our model consists
of the AGN power-law continuum, [O III] doublet, broad and narrow H$\beta$ and
He II lines, and the Fe II blended lines. We find the stellar continuum,
reddening, and He I lines at 4471, 4922, and 5016 Å are not necessary to
achieve a good fit so, to improve the stability of the decomposition of
different spectra, we do not include these components. The broad and narrow
emission lines are all fitted with a fourth-order Gauss–Hermite function (van
der Marel & Franx 1993), whose centroid, dispersion, $h_{3}$, $h_{4}$, and
amplitude are the parameters in the fitting. The profiles of the other narrow
lines (the dispersion, $h_{3}$, and $h_{4}$) are tied to that of [O III]
$\lambda$5007\. The ratio of [O III] $\lambda\lambda$4959,5007 is fixed to
1/3. The rest-frame wavelength separation of the two lines is fixed to 47.9 Å.
The Fe II emission lines are modeled with the template from Boroson & Green
(1992). We adjust the redshift, dispersion, and amplitude of the template in
the fitting. We fit the spectra at rest-frame 4500–5500 Å, masking the small
wavelength range of [N I] lines at 5199 and 5201 Å. We subtract all of the
best-fit components except the narrow and broad H$\beta$ line, and the
H$\beta$ flux is integrated over rest-frame 4800–4940 Å, in order to be
consistent with Bentz et al. (2021). To estimate the uncertainty of the light
curve, we perturb the data according to the measurement uncertainty of the
spectra and redo the measurement. The final uncertainty is the quadrature sum
of the statistical uncertainty from the bootstrapping and 1% of the integrated
flux to make sure that the resulting uncertainty is not too small.

The light curve measured by the decomposition method has a very similar shape
to that from Bentz et al. (2021) as shown in Figure 8. The former is about 5%
higher than the latter, probably because the best-fit continuum from the
decomposition method is slightly lower than that from the local continuum
underneath the H$\beta$ line. Using the method introduced in Section 3, we are
able to infer model parameters consistent with those reported in Table 3
within 1-$\sigma$. However, the peak-to-valley variation amplitude of the new
light curve is reduced by about 20%. Correspondingly, the new inference
prefers about 30–50% larger BLR radius, inclination angle, and $D_{A}$ in
order to reduce the variation amplitude of H$\beta$ light curve. The model
inference is slightly biased in this way. To compensate for this, we need to
take the nonlinearity of the H$\beta$ response (Equation 17) into account.
This yields inferred parameters almost exactly the same as those reported in
Table 3 with the additional parameter $\delta_{c}\approx-0.2$. The key
parameters are listed in Table 4. The responsivity of H$\beta$ is physically
expected to be $\lesssim 1$ (i.e., $\delta_{c}<0$) based on photoionization
models (e.g., Gaskell & Sparke 1986; Korista & Goad 2004) as well as
observations of e.g., NGC 5548 (Goad & Korista 2014). We emphasize that the
observational effects, such as using the $V$-band continuum as a proxy for the
ionization continuum and spectral decomposition over a wide wavelength range,
may also introduce a nonlinear effect. The change of the light curve variation
amplitude does not influence much the time lag of the RM. However, our test
shows that the joint analysis may be biased if the nonlinear effect is not
properly taken into account.

Figure 8: Comparison of the H$\beta$ light curves measured by Bentz et al.
(2021, in blue) and the decomposition method (in red, see the text).