Spectroastrometric Reverberation Mapping of Broad-line Regions

Spectroastrometry measures source astrometry as a function of wavelength/velocity. Reverberations of spectroastrometric signals naturally arise in broad-line regions (BLRs) of active galactic nuclei (AGNs) as a result of the continuum variations that drive responses of the broad emission lines with time delays. Such signals provide a new diagnostic for mapping BLR kinematics and geometry, complementary to the traditional intensity reverberation mapping (RM) technique. We present the generic mathematical formalism for spectroastrometric RM and show that under realistic parameters of a phenomenological BLR model, the spectroastrometric reverberation signals vary on a level of several to tens of microarcseconds, depending on the BLR size, continuum variability, and angular-size distance. We also derive the analytical expressions of spectroastrometric RM for an inclined ring-like BLR. We develop a Bayesian framework with a sophisticated Monte Carlo sampling technique to analyze spectroastrometric data and infer the BLR properties, including the central black hole mass and angular-size distance. We demonstrate the potential of spectroastrometric RM in spatially resolving BLR kinematics and geometry through a suite of simulation tests. The application to realistic observation data of 3C~273 obtains tentative, but enlightening results, reinforcing the practical feasibility of conducting spectroastrometric RM experiments on bright AGNs with the operating Very Large Telescope Interferometer as well as possibly with the planned next-generation 30 m class telescopes.


INTRODUCTION
The so-called broad-line regions (BLRs) of active galactic nuclei (AGNs), responsible for prominent broad emission lines in AGN electromagnetic spectra, are of central importance to the unveiling of gaseous environments surrounding the central supermassive black holes (SMBHs) in general and to the direct measurement of SMBH masses in particular. The characteristic sizes of BLRs range from light-days to light-months, yet too compact at cosmic distances to be spatially resolved by existing instruments. The reverberation mapping (RM) technique stemming from the early works of Bahcall et al. (1972) and Blandford & McKee (1982) provides an effective tool to resolve BLRs by swapping spatial resolution for time resolution (Peterson 1993). Its basic ideas are straightforward: the BLR reprocesses the incident ionizing continuum emitted by the accretion disk into broadline emissions with time delays due to the light-travel time from the accretion disk to the BLR. Different parts of the BLR have different time delays and line-of-sight (LOS) velocities; therefore, analyzing reverberation properties of the broad emission line with respect to the continuum variations deliver information about the BLR geometry and kinematics.
Such information is encoded in the two-dimensional transfer function in velocity and time-delay plane (also called the velocity-delay map), which can be inferred from spectroscopic monitoring data. However, strictly speaking, RM probes the BLR structure along two dimensions, namely, the time delay and LOS velocity dimensions. On the isodelay and/or isovelocity surfaces, the BLR structures are indeed degenerated from the viewpoint of RM. In this sense, RM bears its own limitations in probing the full dimensions of the BLR structures and further improvements and/or alternative approaches are still warranted for a more thorough understanding of BLR geometry and kinematics.
The spectroastrometry (SA) technique provides such an alternative in the sense that SA probes new dimensions perpendicular to the LOS (Beckers 1982;Bailey 1998a). In particular, measuring a source's SA, namely, the photocenter as a function of wavelength, can achieve a higher positioning accuracy than the angular resolution of the image by a factor of ∼ 1/ dN ph /dν, where dN ph /dν is the number of photons received in a spectral pixel bin (Beckers 1982;Bailey 1998a). This means that the SA technique can deliver spatial information about a source once it is sufficiently bright. Indeed, the SA technique has long seen application in radio and millimeter observations, where making the so-called velocity channel maps is a quite standard procedure, from which source positions can be routinely derived. In the optical/infared, however, early full-fledged applications came around the turn of the 21st century and were mainly concentrated on detecting close binary stars (e.g., Bailey 1998b; Baines et al. 2004;Porter et al. 2004) and spatially resolving circumstellar environments and stellar surface structures (e.g., Takami et al. 2001Takami et al. , 2003Whelan et al. 2004) at the level of milliarcseconds. Subsequently, Gnerucci et al. (2010) explored the potential of using the SA of narrow emission lines from the rotating gaseous disk in galactic nuclei to constrain the kinematics of the disk as well as to measure the mass of the central SMBH. They soon after applied this method to two nearby bright galaxies (Gnerucci et al. 2011(Gnerucci et al. , 2013.
Because of the compact sizes of BLRs, until recently, applications to BLRs had been made possible with the successful observation of the infrared Paα line of 3C 273 by the GRAVITY instrument on board the Very Large Telescope Interferometer (VLTI; GRAVITY Collaboration et al. 2017). GRAVITY can achieve an angular astrometric resolution down to ∼10 µas (GRAVITY Collaboration et al. 2018), sufficient to resolve the BLRs of bright AGNs. This achievement, together with the subsequent successful observations of two other AGNs (GRAVITY Collaboration et al. 2020, 2021a, ushered in a new pathway towards BLR physics. As in RM, arising from BLR reprocesses, the SA of the BLR also reverberates to the continuum variations. Because of different time delays and LOS velocities at different BLR parts, the reverberation of the spectroastrometric signals with respect to the continuum can be used to map BLR geometry and kinematics, complementary to the traditional RM. Such an idea had previously been outlined by Shen (2012), which numerically showcased astrometric signals of BLR models with simple geometry. In this work, we further carry forward this idea by developing a generic mathematical framework for spectroastrometric RM and deriving all invoked formulae. As such, we can calculate the spectroastrometric RM signals given an arbitrary BLR model as well as construct a Bayesian approach to explore the BLR model parameters.
Another remarkable capability of spectroastrometric RM is that it can directly measure the geometric distance of the source (see also Wang et al. 2020). This is because SA observations yield both the intensity and photocenter of the broad emission line with wavelength. The intensity reverberation provides information on the physical size of the BLR as in the traditional RM, while the photocenter reverberation provides information on the angular size of the BLR. A combination of these two lines of information naturally constitutes an elegant probe of the geometric distance using AGNs, which therefore has great potential for cosmology (Elvis & Karovska 2002). Previously, Wang et al. (2020) made the first effort in this direction. They integrated the SA observations of the infrared Paα line and RM observations of the optical Hβ line in 3C 273 and inferred the angular-size distance to 3C 273 (see also Li et al. 2022). A subsequent work from the GRAVITY Collaboration et al. (2021b) made an application to another nearby bright AGN NGC 3783. Those applications, however, all neglected the SA reverberation and resorted to combining observations of different emission lines, which might correspond to BLRs different in sizes and/or kinematics (see the discussion in Li et al. 2022) so some systematics will likely arise. The spectroastrometric RM proposed in this work is more straightforward and surmounts those issues as it measures a single emission line. This paper is organized as follows. Section 2 develops the generic framework for BLR spectroastrometric RM, including the basic equations, realistic implementation, and a crosscorrelation analysis. Section 3 presents illustrative examples in which there exists simplified or analytical expressions for spectroastrometric RM and Section 4 presents spectroastrometric RM signals for a generic BLR model. In Section 5, we construct a Bayesian approach to infer BLR model parameters and demonstrate the validity of our approach through a suite of simulation tests. With this approach, we also study the dependence of the Bayesian inferences on the sampling rate and measurement errors of simulated SA data. In Section 6, we then apply our approach to the observations of 3C 273 and show the tentatively obtained black hole mass and angular-size distance of 3C 273. A discussion and conclusions are given in Sections 7 and 8, respectively.
Throughout the paper, we use intensity RM to refer to the traditional RM (only involving line intensity) so as to distinguish it from spectroastrometric RM.

SPECTROASTROMETRIC RM 2.1. Basic Equations
The emission lines stemming from BLRs respond to the continuum variations with time delays. The spatially extended distribution of BLR gases leads to a distribution of response amplitude over time delays and velocities. At time t, the emission line flux is given by where F c is the continuum flux, δ(x) is the Dirac δ function, n denotes the unit vector of the LOS (which points from the BLR to the observer), (r) denotes the spatial distribution of the response coefficient, f (r, w) denotes the distribution of velocity w at the position r, and v is the LOS velocity related to the observed wavelength as where λ 0 is the rest-frame wavelength of the emission line under consideration. By defining an intensity transfer function as Equation (1) can be simplified into (e.g., Blandford & McKee 1982;Peterson 1993) The intensity transfer function Ψ(v, τ ) is determined by the emissivity and velocity distributions of BLR gases (Blandford & McKee 1982). Here, we assume that the BLR remains dynamically stable during the period of observations so that the time dependence of (r) and f (r, w) is negligible. This assumption is reasonable provided the period of observations is shorter than the dynamical time of BLRs. Integrating Equation (4) over velocity yields the integrated line fluxF whereΨ(τ ) is the velocity integral of the intensity transfer function, i.e.,Ψ Hereafter, we use the tilde symbol over a variable to denote its velocity integral. Due to the responses of BLRs to continuum variations, the SA may also vary with time, namely, giving rise to spectroastrometric reverberations. We will show below that SA can be measured by either a spectrometer or an interferometer. Similar to Equation (1), we define the moment of the BLR photons as where extracts the component of r perpendicular to n, namely, the projection of r onto the sky. It is worth stressing that the moment defined above is not normalized by the source's total flux, somehow different from the normal convention. As shown below, a such defined moment can be expressed in the framework of RM and facilitates theoretical analysis and model calculations. Similarly, by defining a spectroastrometric transfer function as Equation (7) can be simplified into By noting that Equations (4) and (10) have the same integral form, we therefore call such astrometric responses of BLR emissions to continuum variations as spectroastrometric RM.
With the above defined moment, the photocenter of the BLR is calculated as Integrating Equation (10) over velocity v results iñ whereΠ(τ ) is the velocity integral of the spectroastrometric transfer function, i.e., Meanwhile, we define the velocity integral of the photocenter Θ(v, t) weighted by the line profile as As a result, we obtain the relation among the photocenter, moment, and line flux We can also calculate the delay integral of Ψ(v, t) and Π(v, t) aŝ where hereafter we use the hat symbol over a variable to denote its delay integral. It is easy to prove that the mean line profileF BLR (v) is proportional toΨ(v) and the mean mo-mentM BLR (v) is proportional toΠ(v). Here,F BLR (v) and M BLR (v) are defined as time averages over a time duration T ,F whereF c is the mean continuum flux over the time duration. Similarly, if defining the mean photocenter as the time average of Θ(v, t) weighted by the line profile, namely, we have the relation Delay integral of the intensity transfer function Velocity integral of the spectroastrometric transfer function Π(v) Delay integral of the spectroastrometric transfer function Fc(t) Driving continuum flux density F c (t) Continuum flux density underlying the emission line Time averaging of emission line flux density Velocity integral of the photocenter Θ(v) Time averaging of the photocenter ∆Θ(v, t) Differential photocenter ∆Θ(t) Velocity integral of the differential photocenter ∆Θ(v) Time averaging of the differential photocenter In a nutshell, the essence of spectroastrometric RM can be expressed in a concise equation 1 where ⊗ represents a convolution. This equation implies that variations of the emission line F BLR and moment M BLR can be regarded as blurred echoes of the variations of the continuum F c . In Table 1, we summarize the major notations and their meanings. There had been a number of methods developed to perform deconvolution for intensity RM, among which include the maximum entropy technique (Horne 1994), the regularized linear inverse method (Krolik & Done 1995;Anderson et al. 2021), a non-parameteric Bayesian method (Li et al. 2016), and the Pixon-based method (Li et al. 2021). In addition, a Bayesian forward approach, dynamical modeling of BLRs, had also been proposed for RM analysis (Pancoast et al. 2011(Pancoast et al. , 2014Li et al. 2013Li et al. , 2018. This approach starts with a flexible dynamical model of BLRs, from which the intensity transfer function can be directly determined, and then employs a Bayesian framework to constrain the model parameters and hence BLR geometry and kinematics. These methods can also be applied to spectroastrometric RM.

Observational Perspectives
In practice, we can observe the SA of BLRs using a spectrometer (e.g., Stern et al. 2015;Bosco et al. 2021) or interferometer (e.g., GRAVITY Collaboration et al. 2018; see Section 7.1 below for a brief discussion of the observational challenges of SA). A spectrometer yields the angular photocenters of BLRs with velocity/wavelength along a specific direction, related to physical photocenters through the cosmic distance of the object, namely, where θ(v, t) is the observed angular photocenters, j is the spatial direction of the spectrometer's slit, and D A is the angular-size distance. An interferometer measures phases of BLRs with velocity/wavelength, related to physical photocenters through the baselines and the cosmic distance of the object, namely, where λ is the wavelength, φ(v, t) is the observed phases, and B is the baseline of the interferometer. In realistic observations, the lights admitted to a telescope always consist of two sources: one from the continuum and the other from the BLR. As a result, the continuum emission also contributes to the observed photocenters, which are now written as where F c (t) is the flux and Θ c (v, t) is the photocenter of the continuum underlying the emission line. Note that usually F c (t) is not the same as the driving continuum F c (t) in Equation (21). However, AGN continuum across UV/optical and infrared bands are well correlated (e.g., Edelson et al. 2019;Minezaki et al. 2019 and references therein). Therefore, F c (t) can be regarded as an echo of F c (t) with a time delay and possibly time blurring (if there is one). Unless stated otherwise, below we neglect this time delay and time blurring and directly use F c (t) to replace F c (t) for the sake of simplicity. It is generally reasonable to assume that the continuum's photocenter does not change with wavelength/velocity, i.e., Θ c (v, t) = Θ c (t). Therefore, we define differential photocenters to simplify the analysis where Here, the continuum's photocenter Θ c (t) can be measured in wavelength regions without the presence of emission lines. Below, by default, we neglect Θ c (t) and simply adopt Θ c (t) = 0.
It is worth pointing out that for SA observed with a spectrometer, one can alternatively first fit and subtract the continuum underlying the emission line in the recorded 2dimensional spectrum so as to directly measure pure photocenters of the emission line. The procedure of continuum subtraction effectively adds extra noises to the measured pure photocenters (e.g., see Whelan & Garcia 2008). Throughout the paper, we by default use the differential photocenters defined by Equation (25), in which there is a scaling factor of the flux ratio F BLR /(F c + F BLR ). Below we also simply use "photocenter" to refer to "differential photocenter".

Cross-correlation Functions
It is known that the cross-correlation function (CCF) between the continuum and emission line is related to the intensity transfer function as (e.g., Welsh 1999;Li et al. 2013) where the definition of CCF is given in Appendix A and ACF c represents the auto-correlation function of the continuum itself, namely, Similarly, the CCF between the continuum and moments is related to the spectroastrometric transfer function as (29) The CCF between the continuum and photocenters is not straightforward because photocenters are not linearly dependent on continuum. Nevertheless, we can implement the following approximations. For small continuum variations (δF c /F c 1), we have As a result, the CCF between the continuum and photocenters can be written where σ(x) represents the standard deviation of the time series x. We can deduce from the above equation that the CCF of the photocenters is proportional to the difference between the CCFs of the moment and line profile with weights.

ILLUSTRATIVE EXAMPLES 3.1. Two Extreme Cases
Equations (4) and (10) specify how the spectral flux and photocenter of the BLR reverberate to the continuum variations. Under a generic framework, it is not straightforward to calculate the integrals analytically in Equations (4) and (10). However, in the following two cases, there exist simple expressions.
Single-pulse Continuum -If the continuum only consists of a single pulses at a time t 0 and vanishes at other times, i.e., (4) and (10) we have This implies that the variations of line flux and SA directly reflect the corresponding transfer functions.
Constant Continuum -Conversely, if the continuum remains constant, i.e., F c (t) = F c,0 , the BLR's spectral flux, moment of photons, and photocenter all do not vary with time, and These equations are analogous in form to Equations (17),  (17), (18), and (20) refer to time averaging so that any time-dependent information is eliminated. When there is only one-epoch observation, one can use the above equations as the first order of approximation.

An Inclined Planar Ring
For an inclined planar ring, there exist analytical expressions for the intensity and spectroastrometric transfer functions. We create a Cartesian coordinate frame X Y Z with its origin located at the center of the ring and Z -axis aligned with its rotating axis. We then rotate the X Y Z frame around the Y -axis by an angle of π −i to create a new Cartesian coordinate frame XY Z. Here, i is the inclination angle of the ring. We set the LOS along the X-axis so that the Y Z plane defines the observer's sky plane. For simplicity, we assume that the emissivity is isotropic and constant along the ring. Figure 1 shows a schematic of the coordinates and the planar ring. In Appendix B, we derive the analytical expressions of spectroastrometric RM for an inclined planar ring.
In the three top panels of Figure 2, we plot the intensity transfer function Ψ(v, τ ), the yand z-components of the spectroastrometric transfer function Π y (v, τ ) and Π z (v, τ ) for a planar ring with an inclination angle of 60 • . The two bottom right panels of Figure 2 show the photocenters Θ y (v, t) and Θ z (v, t) in a case where the continuum pulses at t = 0. The bottom leftmost panel shows the photocenter along the direction j = (0, cos π/4, sin π/4). All these quantities are nonzero only along the ellipse v 2 /V 2 + (cτ /R − 1) 2 = sin 2 i. It is easy to show that for a BLR composed of a series of coplanar rings, the above plots are just a superimposition of corresponding ellipses with different velocities V and radius R. In Figure 3, we illustrate velocity integrals and delay integrals of the transfer functions.
In the left bottom panel, we plot the change in the photo-centerΘ z (t) with time in a case where the continuum pulses at t = 0. In the right bottom panel, we plot the photocen-terΘ y (v) with velocity in a case where the continuum is a constant.    (2) construct a BLR model. We assume that the continuum variations follow the damped random walk (DRW) model and generate mock light curves using the procedure detailed in Appendix C. Regarding the BLR model, without losing the generality, we simply assume a disk-like axisymmetric geometry and Keplerian rotation. The BLR consists of a large number of discrete point-like clouds, which rotate coherently around the central SMBH and reprocess the ionizing continuum. The clouds' emissions are isotropic and the shadowing among clouds is neglected for simplicity. By using cylindrical coordinates, the clouds follow a gamma distribution in the radial r-direction, parametrized by the mean BLR radius R BLR , the inner edge parameter F , and the shape pa-rameter β (see Pancoast et al. 2014 andLi et al. 2018 for details); in the vertical θ direction, the clouds subtend an opening angle θ opn and have a uniform distribution in terms of θ; in the azimuthal ϕ direction, the clouds are also distributed uniformly. The BLR is viewed at an inclination angle of θ inc . In Table 2, we summarize the BLR model parameters and list their values used in our calculations. The flux ratio of the line peak to the continuum is set to about 0.8 (see Equation 25). The spectral broadening is set to 235 km s −1 to account for instrumental broadening effects. These values are consistent with the inferences from the spectroastrometric observations of 3C 273 by the GRAVITY Collaboration et al. (2018). In Appendix D, we demonstrate how to calculate the transfer functions (Equations 3 and 9) given a BLR model.
In Figure 4, we show the obtained intensity transfer function Ψ(v, τ ) and spectroastrometric transfer functions Π y (v, τ ) and Π z (v, τ ), as well as their corresponding de-  Table 2. In each panel, the right and bottom subpanels show the velocity and delay integrals of the transfer functions, respectively. All these transfer functions are in arbitrary units, but the relative scaling of the units is the same for Πy(v, τ ) and Πz(v, τ ). lay and velocity integrals. The intensity transfer function has a bell-like shape, with the significant responses concentrated around 50 days and a long tail extending to several hundred days. This is because the radial distribution of the BLR clouds has a steep decay with radius for the shape parameter β = 1.4. The y-component of the spectroastrometric transfer function Π y has an S-shape along the velocity axis, arising from the Keplerian rotation that causes the redshifting and blueshifting clouds to offset oppositely on the sky. The velocity integralΠ y (τ ) goes to zero because the red and blue parts exactly cancel out. The z-component Π z displays a different response pattern withΠ y . The nearside of the BLR responds earlier and has a negative z-coordinate (see Figure 1); therefore,Π z is negative at short time delays. At long time delays, the farsize of the BLR starts responding andΠ z turns positive. This also leads the delay integralΠ z (v) to vanish because the parts with short and long delays cancel out in the integral.
In Figure 5, we plot a randomly generated continuum light curve and its driven time series of the emission line and photocenters. As expected, the variations in the emission line flux are delayed with respect to those of the continuum. Regarding the photocenters, first of all, we note that while the intensity transfer function is always positive, the spectroastrometric transfer function can be either positive or negative because, as mentioned above, the defined photocenters of BLR clouds can either be positive or negative (see Figure 1). This is also seen in the velocity integrals of the line flux and photocenters shown in Figure 6. The line fluxF BLR (t) varies as a delayed and blurred echo of the continuum and a positive CCF peak appears around 50 days (because the intensity transfer function peaks at about 50 days; see Figure 4). The y-component ∆Θ y (t) is almost zero since the photocenter has positive values at the blue wavelength and negative values at the red wavelength so that their velocity integral cancel out. The z-component ∆Θ z (t) displays an inverse variation pattern compared to that of the continuum due to the strong negative response ofΠ z (t) around 50 days (see Figure 4), which leads to a negative CCF as shown in the right bottom panel of Figure 6. However, we note that such negative CCFs are not intrinsic and just caused by the definition of the photocenter axes shown in Figure 1.
From Figure 5, we can also find that the overall variation amplitude of the y-component photocenters ∆Θ y around the line core is at a level of several microarcseconds, depending on the continuum variability, BLR size, and angular-size distance as well. Generally speaking, larger continuum variability, a larger BLR, or a smaller angular-size distance will yield larger variations in photocenters. However, a larger BLR usually corresponds to a more luminous AGN and thereby a longer variation timescale.

BAYESIAN INFERENCES 5.1. A Bayesian Framework
We now develop a generic Bayesian framework to infer BLR parameters from spectroastrometric data. Given a DRW model and BLR dynamical model, we can reconstruct the continuum light curve from the observed continuum light curve and calculate the spectroastrometric signals using the procedures described in the preceding sections. The calculated spectroastrometric signals are then compared against the observed data, namely, the flux (D line ) and SA (D sa ) data of the emission line. By assuming that the data errors are Gaussian and uncorrelated, the likelihood probability is given by where θ represents the model parameter set, and where the superscript "m" represents the corresponding quantities calculated from the BLR model and i and j represent the epoch and wavelength bin. The posterior probability is then given by where P (θ) is the prior probability of model parameters and P (D line , D sa ) is the Bayesian evidence. We employ the Markov Chain Monte Carlo (MCMC) technique to optimize the posterior probability and the diffusive nested sampling algorithm ) to generate Markov chains. We implement the above procedures based on our previously developed package BRAINS for BLR dynamical modeling (Li et al. 2013(Li et al. , 2018(Li et al. , 2022, which is publicly available at https://github.com/LiyrAstroph/BRAINS. This package is written in C language and uses the diffusive nested sampling library CDNest (Li 2020a) based on the original work of the diffusive nested sampling algorithm by Brewer et al. (2011). This package supports the standardized message-passing interface and therefore can run on a wide range of supercomputer clusters without any reliance on special features of proprietary compilers.

Simulation Configurations
We generate simulated datasets by injecting Gaussian noise to mimic realistic observations and then run our Bayesian package to test its validity. As mentioned in Section 2.2, there are two approaches to measure the SA of BLRs, namely, using a spectrometer (e.g., Bosco et al. 2021) or an interferometer (e.g., GRAVITY Collaboration et al. 2018). Their respective observables, the photocenter and phase, are indeed related through Equations (22) and (23). Hereafter, unless stated otherwise, we only use photocenters for our following simulation tests.
As illustrated in the middle panel of Figure 7, we generate spectroastrometric data under three slit orientations rotated by 60 • from one another at each epoch (see, e.g., Pontoppidan et al. 2008). Without loss of generality, we set the three orientations to j = (1, 0), (cos π/3, sin π/3), and (cos 2π/3, sin 2π/3). We neglect the possible errors in slit positioning and assume that the orientations remain stable for all epochs. We adopt the values of the model parameters listed in Table 2 and show the generated distributions of  Figure 7. In the right panel of Figure 7, we plot the line profile and photocenters as a function of wavelength at the three orientations for a randomly selected epoch. To account for measurement errors, we add Gaussian noises with a standard deviation of 0.01 (arbitrary unit) for the line profiles and 0.1 µas for the photocenters. As a reference, the peak flux is about 0.8 and the maximum shift of the true photocenters is about 5 µas (see Figure 7), which correspond to a signal-to-noise ratio (S/N) of 80 and 50, respectively. We use 42 equally spaced velocity bins over a range between -6000 and 6000 km s −1 (corresponding to a spectral resolution of λ/∆λ ≈ 1000). The spectral instrumental broadening is set to 235 km s −1 .
The continuum light curve is generated using the DRW model with a typical timescale of 300 days and a variation amplitude parameter that results in an overall relative variability of about 40%. The mean continuum flux density is set to be unity, resulting in a flux ratio of F BLR /F c ≈ 0.8 at the line peak. Again, Gaussian noise with a standard deviation of 1% is injected into the continuum light curve. The time   Figure 9. The inferred posterior distributions of model parameters using spectroastrometric RM data (blue) and only intensity RM data (yellow). Green lines represent the input values. The contours are at the 1σ, 1.5σ, and 2σ levels. The numbers above each diagonal panel mark the media value and 68.3% confidence intervals from the spectroastrometric RM data. span is set to 1500 days and the sampling cadence is set to 3 days apart. Such a cadence is feasible considering the fact that one can synthesize data from different monitoring campaigns as well as from public time-domain surveys, such as the All-Sky Automated Survey for Supernovae (Kochanek et al. 2017) and the Zwicky Transient Facility (Graham et al. 2019). The seasonal gaps are not included for the sake of simplicity.
We stress that the above configurations are designed only for illustration purposes. In particular, the injected S/N of the photocenters is somehow idealized considering the astrometric accuracy achievable at current facilities (see Section 7.1 below). Besides, there are several additional factors that are not included in simulations. (1) We use uniform errors for all epochs. In reality, the errors might vary among epochs because of different observing conditions; (2) We do not include the narrow-line component superimposed on the broad emission line, which might affect the observed photocenters, depending on its flux ratio compared to the broad component. A possible economic solution is masking out the wavelength range with the narrow line in the BLR modeling; (3) AGNs show a wide range of variability, not all of which are conducive to doing spectroastrometric RM analysis. This issue can be resolved by the preselection of targets based on variability. Detailed investigations into these factors are quite beyond the scope of this work. Below we will only show how the two key configurations, namely, sampling rate and photocenter errors, affect the Bayesian inference.

A Test Case
As a test case, we generate 50 equally spaced epochs of SA over a time span of 1500 days. Figure 8 shows the generated mock data of the continuum (F c ), line profile (F BLR ), and photocenters (Θ) at three slit orientations (illustrated in Figure 7). Figure 9 plots the posterior distributions of the model parameters. As can be seen, all the parameters are well consistent with the input values at a level of 1σ confidence. In particular, both the black hole mass and angularsize distance are well constrained. For the sake of comparison, in the bottom panels of Figure 8, we plot the recovered line profile and photocenters, which are again well consistent with the simulated data. In Figure 9, we also superimpose the posterior distributions only using the simulated spectral data (namely, intensity RM). The obtained parameters have relatively broader contours, in particular, for the parameters R BLR and β, reflecting the potential of spectroastrometric RM in constraining BLR geometry and kinematics. In addition, if we use a more stringent prior for the angular-size distance from the standard cosmology, all parameter inferences can be further improved.

Dependence on the Sampling Rate
We randomly discard a fraction of epochs of the line profile and photocenter data generated in the preceding section and obtain a set of new data with line epochs ranging from 2 to 50. Note that the sampling of the continuum light curve remains unchanged. We rerun our package and summarize the recovered parameter values and uncertainties for different numbers of epochs as shown in the left panel of Figure 10. As expected, the uncertainties gradually decrease as the epochs increase. It is worth mentioning that even in the case of two epochs, we can still reasonably constrain the BLR size, black hole mass, and angular-size distance, albeit with relatively large uncertainties. Figure 11 compares the posterior distributions of the model parameters for the cases of two epochs and nine epochs, from which, as expected, we can find stronger degeneracy between the BLR size and angular-size distance in the former case. The reasons that the BLR model parameters are reasonably constrained for a few line epochs are twofold. (1) The high fidelity of the continuum light curve ensures a meaningful detection of time delays.
(2) The same BLR model is used for input and outputs so that there are no systematic errors arising from a possible BLR model mismatch. It is worth further investigating the second point in a future work.

Dependence on the Measurement Errors of Photocenters
Considering that the errors of the continuum and line profile are reasonable with existing telescopes, we only concentrate on the errors of the photocenters. We change the input errors of the photocenters from 0.1 to 10 µas (while keeping the errors of the line profiles unchanged) and generate a set of new simulated data. The right panel of Figure 10 illustrates the recovered parameter values and uncertainties for different input errors. Figure 12 shows the simulated data and the recovery for the case of input errors of 5 µas, which are comparable to the maximum shift of the true photocenters (see Figure 7). Because the errors of the line profiles  Figure 11. A comparison of the inferred posterior distributions of model parameters for nine (blue) and two (yellow) randomly selected epochs of spectroastrometric RM data. Green lines represent the input values. The contours are at the 1σ, 1.5σ, and 2σ levels. The numbers above each diagonal panel mark the media value and 68.3% confidence intervals from the case of nine epochs. Note that the sampling of continuum light curves are the same for both cases of two and nine epochs.
are not changed, large input photocenter errors mean that the likelihood from intensity RM becomes dominant over that from spectroastrometric RM. As a result, the BLR model parameters can still be well constrained, except for angular-size distance and P.A., of which the uncertainties increase rapidly as the input photocenter errors are comparable to or larger than 5 µas. In Figure 13, we compare the posterior distributions of the model parameters between the input photocenter errors of 1 and 10 µas. In the latter case, the photocenter error is 2 times the maximum shift of the true photocenters (5 µas). The angular distance and P.A. have broad distributions, but from which reasonable inferences can be made.  Figure 14 and Appendix E.
In running the analysis code, we use the same parameter priors as listed in Table 2 and set the redshift of 3C 273 to be z = 0.158. Because the GRAVITY observations used the adaptive optics systems, there was difficulty in calibrating the absolute fluxes and the measured Paα line fluxes were normalized by the underlying continuum. We assume that the continuum variations underlying the Paα follow those of the V -band light curve (Li et al. 2020b) but with a time delay of about 450 days (Sobrino Figaredo et al. 2020). As such, we determine the continuum fluxes underlying the Paα line by interpolating the the V -band light curve after correcting the time delay and thereby obtaining the absolute fluxes of the Paα line, for which the uncertainties from the interpolation are also included. We finally use the V -band light curve as a proxy for the driving continuum and take into account the delay of 450 days in the spectroastrometric RM analysis.
In the left panels of Figure 14, we plot the reconstructions to the continuum light curve, Paα line profiles, and the differential phase curves of one baseline UT4-UT3. The full fits to the differential phase curves of all six baselines are shown in Appendix E. The right panels of Figure 14 plot the obtained posterior distributions of the BLR size (R BLR ), black hole mass (M • ), and angular-size distance (D A ). There appears strong degeneracy among the three parameters, implying that they are not well constrained and have large uncertainties. This is not surprising considering the noisy differential phase curves and only three epochs of data. While M • tends to approach the lower prior limit of 10 7 M , R BLR and D A are peaked within their prior ranges, which are set to (10, 1300) light-days and (10, 10 4 ) Mpc, respectively. For the sake of comparison, we superimpose the best inferences on the BLR size and black hole mass from GRAVITY Collaboration et al. (2018). Our results are comparable to these inferences within a confidence level of 3σ. Note that the analysis of GRAVITY Collaboration et al. (2018) did not include the spectroastrometric RM, therefore, they needed to preset the angular-size distance (550 Mpc). Remarkably, our obtained angular-size distance log(D A /Mpc) = 2.40 +0.28 −0.34 is marginally consistent within uncertainties with this fiducial value.
We stress that our present application to 3C 273 is somehow tentative in light of the data quality. In particular, the obtained black hole mass reaches its lower prior limit, indicating that to some extent, the adopted priors play a role in constraining the model parameters. Nevertheless, the results are still enlightening. In the future, more epochs of observations with improved data quality for 3C 273 would be highly worthwhile to reinforce our present analysis and most importantly, to reliably measure the black hole mass and geometric distance of 3C 273. This is practically feasible considering the forthcoming upgraded instrument GRAVITY+ (see GRAVITY+ Collaboration et al. 2022).
For the GRAVITY interferometer instrument, the baselines are on the order of 100 m, corresponding to σ s ∼ 3 µas at the K band. For the 30m class single-aperture telescopes, the expected accuracy increases up to σ s ∼ 10 µas. Those estimates are based on the instruments being in ideal conditions (excluding all other sources of errors) and a presumption of 10 6 photons per spectral bin. However, in practice, there are a variety of subtle statistical and systematic error sources that limit the achievable astrometric accuracy, such as atmospheric differential tilt jitter, anisoplanatism of AO systems, and atmospheric differential chromatic refraction, etc (e.g., Cameron et al. 2009;van Belle 2009;Trippe et al. 2010;Rodeghiero et al. 2021). All these error sources can dilute the final achievable astrometric accuracy at different levels. An advantage here is that spectroastrometric RM only requires differential astrometry across small spatial scales so that some of the issues related to absolute astrometry are no longer important. By taking into account all practical sources of errors, previous studies have shown that the forthcoming 30m class single-aperture telescopes are likely to achieve a resolution at the level of several tens to approximately 100 microarcseconds at the K band (e.g., Trippe et al. 2010;Stern et al. 2015;Bosco et al. 2021;Rodeghiero et al. 2021). The observations with the GRAVITY instrument have demonstrated that it can achieve an angular resolution down to ∼10 µas at      Figure 13. A comparison of the inferred posterior distributions of the model parameters for the relative photocenter errors of 20% (blue) and 200% (yellow). Here, the relative photocenter error means the photocenter error relative to the maximum shift of the true photocenters (∼5 µas). Green lines represent the input values. The contours are at the 1σ, 1.5σ, and 2σ levels. The numbers above each diagonal panel mark the media value and 68.3% confidence intervals from the case of the photocenter error of 20%. the K band given sufficiently bright targets as well as adequate exposure times (usually several hours; GRAVITY Collaboration et al. 2017GRAVITY+ Collaboration et al. 2022). These factors, together with our application to 3C 273 in Section 6, indicate that it is viable to conduct spectroastrometric RM experiments on bright AGNs. Such experiments would span several months or years, depending on the BLR sizes, so as to capture the reverberation signals. Simulation tests in Section 5.2 imply that several epochs of SA observations already yield meaningful constraints on BLR model parameters. Nevertheless, we bear in mind the difficulties in performing SA observations. Currently, there are only a few AGNs with SA observation results published (GRAVITY Collaboration et al. 2018, 2021aBosco et al. 2021). This situation might change with the upgraded GRAVITY+ and next-generation 30 m class telescopes.

Angular Sizes of BLRs from the Sloan Digital Sky Survey
In this section, we estimate angular-size distributions of the quasar catalog from the 14th data release (DR) of the Sloan Digital Sky Survey (SDSS; Pâris et al. 2018) to illustrate that there are available potential candidate targets for future spectroastrometric RM experiments. Rakshit et al. (2020) performed detailed spectral measurements for the quasars with the continuum S/N >3 per pixel through multicomponent spectral decompositions, which included host-galaxy subtraction so that the central AGN luminosities could be obtained directly. We convert the luminosities to BLR sizes (or time delays) using the BLR size-luminosity relationship of Bentz et al. (2013) (see also Du & Wang 2019). This relationship relies on the 5100Å luminosity and only corresponds to the Hβ BLR size. For simplicity, we assume that all broad emission lines appropriate for infrared SA observations have the same BLR sizes as the Hβ line. For those high-redshift AGNs (z > 0.8) without available 5100Å luminosities, we make estimates from the given 3000 or 1350Å luminosities using the bolometric correction factors in Richards et al.  Equation (25) demonstrates that the observed photocenters scale with the normalized line fluxes relative to the underlying continuum flux, namely, the flux ratio F BLR (v, t)/F c (t). The peak flux ratio generally differs among AGN emission lines, e.g., the typical value is ∼0.06 for the Brγ line, ∼0.3 for Paβ, 0.6 for Paα, 0.12 for Paγ, and >1 for Hα (Landt et al. 2008;Rakshit et al. 2015). Using the suitable emission line for the K band at different redshifts, we simply assume a peak flux ratio of 0.06 for z < 0.08, 0.6 for 0.08 z < 0.4, 0.3 for 0.4 z < 0.87, 0.12 for 0.87 z < 1.26, and 1.0 for z 1.26. All the calculations assume a ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7. Figure 15 shows the distributions of the maximum shifts of the photocenters, observed time delays, and K-band (Vega) magnitudes. Here, the K-band magnitudes are compiled by Pâris et al. (2018) from the Two Micron All Sky Survey data. The yellow points in Figure 15 represent those AGNs with the maximum shift of photocenters > 5 µas, observed time delays < 1000 days, as well as magnitudes K < 15, which can serve as candidates for future spectroastrometric RM experiments within a reasonable time span. The present dynamical modeling approach requires presuming a BLR model. A naturally arising issue is how the results depend on the presumed BLR model. As in our previous works of Li et al. (2018) and Li et al. (2022), an appropriate way to address this issue is testing a suite of BLR models and singling out the most probable one using statistical methods (such as Bayesian model selection). In some cases, there exist other independent observational measurements/constraints, with which we can further prove the validity of the selected model.

Constraining the BLR Models
For the sake of illustration, in Figure 16, we showcase the transfer functions for three exemplary BLR models: the inflow and outflow prescription of Pancoast et al. (2014), and the disk wind model of Higginbottom et al. (2013). We stress that these three phenomenological BLR models are still restrictive and do not represent the whole story of BLRs, considering the complicated BLR kinematics (e.g., Baskin et al. 2014;Czerny et al. 2017;Wang et al. 2017). Notwithstanding, as can be seen, there are distinctive patterns in both the intensity and spectroastrometric transfer functions among the three models. We expect that appropriate analysis with spectroastrometric RM data can better constrain and prove the validity of different BLR models, compared to the traditional intensity RM.

Comparison with the Previous Joint Analysis of SA and RM
As mentioned above, Wang et al. (2020) proposed jointly analyzing SA and intensity RM (hereafter SARM) data to probe the BLR geometry and kinematics and measure the geometric distance. They made the first such application to 3C 273, which had been observed with infrared SA of the Paα (GRAVITY Collaboration et al. 2018) line and optical intensity RM of the Hβ line (Zhang et al. 2019). Currently, SA is only feasible in infrared whereas almost all intensity RM campaigns are undertaken in the optical, therefore, their respectively observed emission lines are indeed different. As a result, in such joint SARM analysis, one needs to presume that the two different lines share the same BLR. However, Li et al. (2022) showed that the profiles of the Hβ and Paα lines in 3C 273 differ regarding both widths and shapes. Such differences appear to be common in AGNs (Landt et al. 2008; see also Figure 2 of Li et al. 2022), implying that the BLRs must be different in some respects. Li et al. (2022) proposed using velocity-resolved intensity RM data and treating the respective BLRs separately, but leting them share only the inclination angle and SMBH mass. As such, the differences between the two lines are naturally taken into account. The disadvantage of this approach is a longer parameter list, which likely results in relatively large parameter uncertainties.
The spectroastrometric RM circumvents the above issue since it observes the intensity and SA of the same line; therefore, it gets rid of possible systematic errors arising from different lines in previous joint SARM analyses. Currently, the challenges of SA observations still restrict the targets of spectroastrometric RM to a few bright AGNs; however, as mentioned above, the situation might significantly change with the forthcoming upgraded instrument GRAVITY+, which is planned to make performance improvements in several key respects, such as the interferometric fringe tracking and sensitivity magnitude (GRAVITY+ Collaboration et al. 2022).

CONCLUSIONS
We propose that spectroastrometric signals of BLRs in AGNs reverberate to the continuum variations and the responses of different BLR parts show different time delays as a result of spatial distributions of BLR gas (see also Shen 2012). Considering that SA resolves the BLR structure perpendicular to the LOS, spectroastrometric RM, therefore, provides a new diagnostic for BLR geometry and kinematics, complementary to the traditional intensity RM technique. We present the basic mathematical framework for spectroastrometric RM, which indeed can be regarded as a deconvolution problem (see Equation 21) so that a variety of wellestablished mathematical methods are applicable. The underlying essence of spectroastrometric RM is to determine the intensity and spectroastrometric transfer functions, which encode the full information about the BLR. We derive analyt-ical expressions for the case of an inclined ring-like BLR and also show that in extreme cases where the continuum has a pulsing variation or remains constant, there exist simple expressions. For a generic BLR with realistic parameter values, the spectroastrometric signals vary on a level of several to tens of microarcseconds, mainly depending on the BLR size, continuum variability, and cosmic distance.
We developed a forward Bayesian dynamical modeling approach to analyze spectroastrometric RM data and infer BLR properties, in which the posterior probability is explored with the MCMC technique. We constructed a suite of simulation tests to demonstrate the validity of our approach and show the potential of spectroastrometric RM in resolving BLR geometry and kinematics and most importantly measuring the SMBH mass and angular-size distance. An application to the spectroastrometric data of 3C 273 yields tentative, but enlightening constraints on its BLR size, central SMBH mass, and angular-size distance, although there are large uncertainties (see Figure 14). Despite the challenges remaining in SA observations, these results remarkably imply the feasibility of conducting pilot spectroastrometric RM experiments on nearby bright AGNs, in particular, considering the forthcoming upgraded GRAVITY+/VLTI and the panned nextgeneration 30 m class telescopes.
B. DERIVING THE EQUATIONS FOR A RING-LIKE BLR In the XY Z frame, the LOS is n = (1, 0, 0) and a point P (see Figure 1) in the ring has a coordinate of R = R(sin i cos θ, cos θ, − cos i cos θ) and a velocity of w = V (− sin i sin θ, cos θ, cos i sin θ), where R and V are the the radius and rotating velocity of the ring, respectively. As a result, the corresponding time delay of the point P with respect to the origin O is τ = R − R · n c = R c (1 − cos θ sin i).
The LOS velocity is v = −w · n = V sin i sin θ, (B5) and the projected location in the observer's sky is r ⊥ = R − (R · n)n = R(0, sin θ, − cos i cos θ). (B6) The intensity transfer function is given by where R(1 − sin i)/c τ R(1 + sin i)/c, otherwise Ψ(v, τ ) = 0, and The spectroastrometric transfer function has a zero x-component, and its yand z-components are given by and The velocity integrals of the above transfer functions areΨ andΠ y (τ ) = 0, Π z (τ ) = − cos i sin i 2c (R − cτ ) ρ .
The delay integrals of the above transfer functions arê andΠ C. SIMULATING LIGHT CURVES USING THE DRW MODEL We generate mock light curves using the DRW model as follows. Given with the covariance matrix C of a DRW model, its Cholesky decomposition is written as C = M M T , where M is a lower triangular matrix (Press et al. 1992, Chapter 2.9). A mock light curve is obtained with u = M r, where r is a series of Gaussian random numbers with a zero mean and unity deviation. It is easy to show that such a light curve u has a covariance matrix of uu T = M rr T M T = C. Here, the covariance matrix of a DRW model is given by where t i and t j are the times of ith and jth points of the light curve, respectively, and σ d and τ d are parameters that represent the long-term standard variation and typical damping time scale of the DRW process.
D. CALCULATING SPECTROASTROMETRIC SIGNALS GIVEN A BLR MODEL After generating BLR clouds' velocities and positions according to the given BLR model (see, e.g., Li et al. 2022), the transfer functions defined in Equations (3) and (9) are calculated as where i , v i , τ i , y i , and z i are the response coefficient, LOS velocity, time delay, and y and z coordinates of the i cloud, respectively. The flux and SA of the emission line are then calculated using Equations (4) and (10). For simplicity, we assume that all clouds have a uniform response coefficient i and the possible nonlinear response of the line emission to the continuum is neglected (Li et al. 2013).
Wavelength ( Wavelength ( Figure 17, we show the full fits to the differential phase curves of the six baselines observed by the GRAVITY/VLTI for the Paα line of 3C 273 (GRAVITY Collaboration et al. 2018). See Section 6 for the details of the model fitting.