Confronting a Thin Disk-wind Launching Mechanism of Broad-line Emission in Active Galactic Nuclei with GRAVITY Observations of Quasar 3C 273

Quasars are host to sets of emission lines that are single peaked and have widths on the order of thousands of kilometers per second, scales that are assumed to be set by the gravity of the central supermassive black hole (Peterson 2006). The region where these broad lines originate—aptly named the broad-line region, or BLR—is assumed to be at distances of ∼10³–10⁴ r_s (with the Schwarzschild radius ${r}_{s}=\tfrac{2{{GM}}_{\mathrm{BH}}}{{c}^{2}}$) from the central black hole, often inferred by measuring time lags between changes in the continuum and the BLR line profiles (see Zhang et al. 2019 for a recent study of this phenomenon in quasar 3C 273). The broad single peak in the line profile is usually explained via a cloud model, as in Gravity Collaboration et al. (2018, hereafter G18); but Chiang & Murray (1996, hereafter CM96) showed that an alternative "disk-wind" model could match this morphology in the line profile. If the ratio of outer to inner radius is sufficiently small, one would naively expect to observe a double peak in the line profile, which is often not observed (Jackson et al. 1991). The disk-wind model initially proposed by CM96 avoids this problem by adding in the effects of high-velocity gradients present within an optically thick, geometrically thin Keplerian disk, with the shears altering the escape probability of photons emitted at different locations in the disk via Sobolev theory (Sobolev 1957). The high central luminosity of quasars is usually assumed to originate from accretion onto their central supermassive black holes (Lynden-Bell 1969; Shakura & Sunyaev 1973; Rees 1984), and thus it is tempting to assume this simple disk-like geometry extends to the BLR, and there is significant observational evidence for the presence of winds in quasars (Hamann et al. 1993; Bottorff et al. 1997; Elvis 2000), which could be launched by the large velocity gradients in such a model. Figure 1 illustrates the general geometry of the Sobolev disk-wind type model considered in this work.

Zoom In Zoom Out Reset image size

If we can measure a characteristic size, R_BLR, and velocity scale, ΔV, for orbiting gas in the BLR, we can infer the mass of the central supermassive black hole via

$$\begin{eqnarray}&&{M}_{\mathrm{BH}}=f\displaystyle \frac{{R}_{\mathrm{BLR}}{({\rm{\Delta }}V)}^{2}}{G}.\end{eqnarray} \tag{ 1 }$$

Here, G is the gravitational constant and f is the "virial factor," which originates from assuming that the BLR is virialized and whose value is model dependent on both the geometry and kinematics of the BLR. The notation here matches what is given in Waters et al. (2016, hereafter Waters16). Reverberation-mapping techniques allow us to measure a characteristic time delay, t, between changes in the continuum of the source and changes in the line profile, giving a characteristic size for the BLR of R_BLR ≈ ct, assuming the change is propagated at the speed of light, c (Peterson 2006). GRAVITY spatially resolves the BLR in 3C 273, which provides another method for measuring R_BLR. There are significant model- and measurement-dependent uncertainties in both R_BLR and ΔV, and thus it is of critical importance to constrain what physical models best fit the BLR.

Clear evidence for ordered BLR rotation around 3C 273 (see Figure 1 of G18 or Figure 4) is consistent with both the cloud and the thin disk model. While the cloud model is fit satisfactorily to the data in G18, it has not been considered if a thin disk-wind launching model can also explain the data. If the clouds in the model G18 fit to the GRAVITY data are real, how they could survive sufficient time and create sufficient smoothness in the line profiles (Mathews & Capriotti 1985; Dietrich et al. 1999) is not well understood. A disk-wind model is thus easier to reconcile physically, and even if the clouds are not actual physical objects (and instead just a numerical convenience for fitting and describing the BLR kinematics) it is important to test whether other models can also explain the data.

There are many possible disk-wind model morphologies to choose from, but here we consider a simple, two-dimensional, optically thick but geometrically thin disk-wind model, which is essentially a combination of previous work first explored in CM96 and in a follow-up paper published a year later (Murray & Chiang 1997, hereafter MC97) with minor modifications to the velocity gradients and geometries originally considered in those works, a general schematic of which is shown in Figure 1. This model considers a fixed optical depth, τ, in the τ ≫ 1 limit—greatly simplifying the equations of radiative transfer—where an extended thin accretion disk with various hydrodynamical shears drives subtle differences in the escape probability for line photons at different locations in the disk. We include four such shears in our model (further discussed below) and allow their strengths to be artificially varied in fitting to determine what possible wind orientations might exist.

We note that this is not the only possible disk-wind model that could explain this data, as other authors have explored variations on this idea in a variety of capacities: Chajet & Hall (2013) presented a version of the CM96/MC97 model, where they adoped a hydromagnetic prescription for the wind—as opposed to purely hydrodynamical—following on the work of Emmering et al. (1992); Flohic et al. (2012) corrected a minor error in MC97 and added in the effects of relativity, sampling a wide parameter space to produce line profiles that could be compared to observations; and many groups have included variants on the calculation presented in MC97, where the wind has a substantial vertical extent (Waters16; Baskin & Laor 2018; Naddaf & Czerny 2022), with some applying Monte Carlo methods to better model the radiative transfer than we do here (Matthews et al. 2020). We choose the model we do for its flexibility, simplicity, and computational ease of fitting as a first test to GRAVITY data.

As in CM96 and MC97, we utilize the Sobolev (Sobolev 1957) approximation in modeling the line emission, as the macroscopic velocities of gas within the disk are much larger than the intrinsic line width. The Sobolev optical depth is proportional to the inverse of the line-of-sight velocity gradient, which in this work is modeled as various possible shears within the thin disk geometry. In the high optical depth limit within the Sobolev regime, the equation of radiative transfer is essentially dominated by the source function multiplied by the probability that a photon will escape the disk, which in the optically thick limit is directly proportional to the line-of-sight velocity gradient. Thus, any anisotropic terms present in the velocity gradients of the hydrodynamic thin disk can alter the morphology of the observed line profile, leading to a variety of interesting possible shapes outside of the expected double peak.

In the Sobolev approximation, we can calculate the frequency-dependent line profile from a steady continuum source in a thin disk, L_ν , as

$$\begin{eqnarray}&&{L}_{\nu }=\sin i\int {rdr}{\int }_{0}^{2\pi }d\phi k(r)\beta (r,\phi ,i)S(r)\delta [\nu -\tilde{\nu }(\phi ,r)].\end{eqnarray} \tag{ 2 }$$

Here, i is the viewing inclination angle of the system, with r and ϕ disk coordinates, as shown in Figure 1. k is the integrated line opacity (proportional to the density of the emitting gas), β the photon escape probability, and S the source function, which here is modeled as a power-law function of r alone, S(r) ∝ r^−α , with α = 1 in keeping with CM96. This relationship is essentially the familiar equation for the formation of a spectral line, but the usual line function is replaced with the Dirac delta function in keeping with the Sobolev approximation.

This relationship was first shown by Rybicki & Hummer (1983), and here this is the same as CM96's Equation (8), where (as shown in Waters16) $\tilde{\nu }={\nu }_{0}\left(1+\tfrac{{v}_{l}}{c}\right)$ is the Doppler-shifted resonant frequency of the line as seen by the observer, i.e., ${v}_{l}=\sin i({v}_{r}\cos \phi -{v}_{\phi }\sin \phi )=-\sin {{iv}}_{\phi }\sin \phi$ (v_r = 0, and the l subscript here denotes that this quantity is along the observer's line of sight). Note that this is only the classical Doppler shift from circular motion, and does not include the effects of relativity or transverse motion. Both effects are radially dependent and act to make the line profile more skewed, with the relativistic effect being the more important of the two, with a scaling on the order of 1.5(r_s /r) × 10⁵ km s⁻¹ (Waters16). Fortunately, the characteristic size of the BLR is of order 10³ r_s , and thus we do not model this effect as it is not significant over the width of the line profile. Flohic et al. (2012) explored the effects of relativity on a similar disk-wind model and computed line profiles for a wide and robust set of parameters, which verify that for the parameter space surveyed in this work relativity is not a primary factor in the shape of the line profiles.

If we also consider continuum variability, there will be an associated time lag at each disk location in addition to the Doppler shift. Adding this in gives us the so-called "transfer function" (as shown in Waters16 and originally introduced by Blandford & McKee 1982):

$$\begin{eqnarray}&&{\rm{\Psi }}(t,\nu )=\iint {rk}(r)\beta (r,\phi ,i)S(r)\delta [\nu -\tilde{\nu }(\phi ,r)]\delta [t-\tilde{t}]{drd}\phi .\end{eqnarray} \tag{ 3 }$$

For a thin disk with Keplerian rotation, this resonance condition for the time delay, $\tilde{t}$, is given by $\tilde{t}=\tfrac{r}{c}(1-\cos \phi \sin i)$. The deprojected line profile is then just ${\rm{\Psi }}(\nu )={\int }_{0}^{\infty }{\rm{\Psi }}(t,\nu ){dt}$ and the so-called response function is ${\rm{\Psi }}(t)={\int }_{-\infty }^{\infty }{\rm{\Psi }}(t,\nu )d\nu$, which we illustrate with echo images of our best fit in the results section.

Using the standard Sobolev approach, we express the escape probability as $\beta =\tfrac{1-{e}^{-\tau }}{\tau }$, with τ here representing the line optical depth, given by $\tau =\left(\tfrac{k(r)c}{{\nu }_{0}}\right)| \tfrac{{{dv}}_{l}}{{dl}}{| }^{-1}$. Here, we consider optically thick line emission to simplify the equations of radiative transfer (i.e., τ ≫ 1), and in this limit the escape probability β reduces to simply $\beta \approx \tfrac{1}{\tau }$. In the Sobolev approximation τ is inversely proportional to the line-of-sight velocity gradient, giving $\beta \approx \tfrac{1}{\tau }=\left[\tfrac{{\nu }_{0}}{k(r)c}\right]| \tfrac{{{dv}}_{l}}{{dl}}|$. Thus, we can express the quantity $k(r)\beta (r)=\tfrac{{\nu }_{0}}{c}| \tfrac{{{dv}}_{l}}{{dl}}|$. The core assumptions of this approximation are that it is hard for line photons in general to escape the disk as the line is optically thick, with the probability for escape being shaped by the velocity gradients (shears) present within the disk. This allows for more complicated line morphologies than one might naively expect.

As discussed above, we choose a simple power-law dependence for the source function to match previous work, but note that, as shown in Waters16, one can think of the source function as going like $S(r)=A(r){F}_{X}^{\eta (r)}$, where F_X is the flux from the continuum and the power-law index $\eta (r)=\tfrac{\partial \mathrm{ln}{S}_{l}}{\partial \mathrm{ln}{F}_{X}}$ (where S_l is the source function along the line of sight; for further details, see Krolik et al. 1991). Typical η values from photoionization modeling are between 0 and 2, with A(r) = A₀ r^γ setting the overall response amplitude in the line. Equation (2) shows that the emissivity in the line j ∝ S(r), thus the responsivity of the line is $\tfrac{\partial j}{\partial {F}_{X}}\propto A(r)\eta (r){F}_{X}^{\eta (r)-1}$. In this formalism the values chosen for η and γ then set the overall radial scaling of the source function as well as how the line emission responds to the continuum, with higher values of γ enhancing the response at larger radii. Note that there are then many possible models of how the BLR responds to illuminating continuum flux on the disk for a given choice for S(r). Assuming illumination from a central source onto a thin disk BLR gives F_X ∝ r⁻³; thus, for the constraints on η, we can constrain γ to be between −1 (for η = 0) and 5 (for η = 2) to yield the source function S(r) ∝ r⁻¹ we use in this work.

Applying this prescription for the source function as well as the simplification for k(r)β(r) then reduces Equation (3) to

$$\begin{eqnarray}&&{\rm{\Psi }}(\nu ,t)={{\rm{\Psi }}}_{0}\displaystyle \frac{{\nu }_{0}}{c}{\int }_{{r}_{\min }}^{{r}_{\max }}{\int }_{0}^{2\pi }| \displaystyle \frac{{{dv}}_{l}}{{dl}}| \delta [\nu -\tilde{\nu }]{drd}\phi \delta [t-\tilde{t}],\end{eqnarray} \tag{ 4 }$$

where Ψ₀ is a normalization constant, whose value one would need to specify to compare to physical flux values but here we leave as an arbitrary constant as we are only attempting to match the shape of the line profile.

We will evaluate this integral numerically with the following general approach:

• 1.  
  First, we make a two-dimensional disk in log polar coordinates, where each cell has coordinates (r, ϕ) that correspond to an associated resonant Doppler shift $\tilde{\nu }$ and time delay $\tilde{t}$. The general geometry of the model is shown in Figure 1.
• 2.  
  We then calculate the intensity escaping toward the observer at each location within the disk, as given by Equation (4).
• 3.  
  Finally, we integrate to get the total line luminosity (${\rm{\Psi }}\nu ={\int }_{0}^{\infty }{\rm{\Psi }}(\nu ,t){dt}$), binning the disk according to $\tilde{\nu }$ and summing over all time delays to get the line profile and vice versa (${\rm{\Psi }}t={\int }_{-\infty }^{\infty }{\rm{\Psi }}(\nu ,t)d\nu$) to obtain the response function. This ensures each region of the disk only contributes to the total intensity integral at its corresponding resonant Doppler frequency as is required by the Dirac delta function in Equations (2)–(4).

We must now evaluate $| \tfrac{{{dv}}_{l}}{{dl}}|$. The line-of-sight velocity gradient can be found using the rate of strain tensor, $\hat{n}\cdot {\boldsymbol{\Lambda }}\cdot \hat{n}$, often denoted as $Q={\sum }_{i,j}\tfrac{1}{2}\left(\tfrac{\partial {v}_{i}}{\partial {r}_{j}}+\tfrac{\partial {v}_{j}}{\partial {r}_{i}}\right)$. Shears in the disk create the velocity gradients we seek to recover, which we will now derive. The $\hat{n}$ for this geometry at an observer of ϕ = 0 and inclination i is given by Waters16 (see there and/or CM96 and Flohic et al. 2012 for further discussion and derivation) as

$$\begin{eqnarray}\begin{array}{rcl}\hat{n} & = & (\sin \theta \cos \theta \sin i+\cos \theta \cos i)\hat{r}\\ & & +(\cos \theta \cos \phi \sin i-\sin \theta \cos i)\hat{\theta }\\ & & -(\sin \phi \sin i)\hat{\phi }.\end{array}\end{eqnarray} \tag{ 5 }$$

Using this $\hat{n}$ and the rate of strain tensor terms in spherical coordinates from Batchelor (1968), we can write the line-of-sight velocity gradient as

$$\begin{eqnarray}\begin{array}{rcl}\hat{n}\cdot {\boldsymbol{\Lambda }}\cdot \hat{n} & = & {\sin }^{2}i\left[\displaystyle \frac{\partial {v}_{r}}{\partial r}{\cos }^{2}\phi -\left(\displaystyle \frac{\partial {v}_{\phi }}{\partial r}-\displaystyle \frac{{v}_{\phi }}{r}\right)\sin \phi \cos \phi \right.\\ & & \left.+\displaystyle \frac{{v}_{r}}{r}{\sin }^{2}\phi \right]\\ & & -\sin i\cos i\left[\left(\displaystyle \frac{1}{r}\displaystyle \frac{\partial {v}_{r}}{\partial \theta }+\displaystyle \frac{\partial {v}_{\theta }}{\partial r}-\displaystyle \frac{{v}_{\theta }}{r}\right)\cos \phi \right.\\ & & \left.-\displaystyle \frac{1}{r}\displaystyle \frac{\partial {v}_{\phi }}{\partial \theta }\sin \phi \right]\\ & & +{\cos }^{2}i\left[\displaystyle \frac{1}{r}\displaystyle \frac{\partial {v}_{\theta }}{\partial \theta }+\displaystyle \frac{{v}_{r}}{r}\right].\end{array}\end{eqnarray} \tag{ 6 }$$

Where in arriving at the form above we have assumed all of the $\tfrac{\partial }{\partial \phi }$ operator terms are 0 (the underlying disk is Keplerian) and the disk is in the equatorial plane ($\theta =\tfrac{\pi }{2}$), which allows us to significantly simplify the result. A more thorough derivation of Equation (6) is given in Appendix A. We keep the θ terms here to generalize the wind in both the "vertical" and radial directions. This is the same result as given in CM96/MC97 but with the ϕ convention for the observer given in Waters16, which differs from CM96/MC97 by $-\tfrac{\pi }{2}$. For the rest of this work, we will use the ϕ convention specified in Waters16 (Equation (6)).

In CM96/MC97, the authors assume that v_r ≈ 0 in the thin disk, but that there is an acceleration related to the escape velocity, i.e., $\tfrac{\partial {v}_{r}}{\partial r}\approx 3\sqrt{2}\tfrac{{v}_{\phi }}{r}$, where ${v}_{\phi }=\sqrt{\tfrac{{GM}}{r}}$ is the Keplerian v_ϕ , which gives us $\tfrac{\partial {v}_{\phi }}{\partial r}=\tfrac{-{v}_{\phi }}{2r}$. These accelerations are important, as in Sobolev theory the velocity gradients essentially give us the escape probability of photons resonating in the thick medium of the disk, and it is these escaping photons that we image (Sobolev 1957; Rybicki & Hummer 1983).

But what are the θ terms? Following in the footsteps of CM96/MC97, it makes sense to assume that on average v_θ ≈ 0, for the same reason v_r ≈ 0, but similarly we will assume a particle may be lifted by the wind and accelerated to the local escape velocity (but now in the $\hat{\theta }$ direction) such that $\tfrac{\partial {v}_{\theta }}{\partial \theta }\approx \tfrac{{v}_{{esc}}}{\left(H/R\right)}$ (where H/R ≪ 1 is the scale height of the disk) and $\tfrac{\partial {v}_{\theta }}{\partial r}\approx \tfrac{\partial {v}_{r}}{\partial r}$. Since v_ϕ (r) is a function of r alone, $\tfrac{\partial {v}_{\phi }}{\partial \theta }=0$, and we also set $\tfrac{\partial {v}_{r}}{\partial \theta }=0$ in keeping with the idea of a thin disk.

These approximations reduce Equation (6) to

$$\begin{eqnarray}\begin{array}{rcl}\hat{n}\cdot {\boldsymbol{\Lambda }}\cdot \hat{n} & = & 3\displaystyle \frac{{v}_{\phi }}{r}{\sin }^{2}i\cos \phi \left[\sqrt{2}\cos \phi +\displaystyle \frac{\sin \phi }{2}\right]\\ & & -\sin i\cos i\left[3\sqrt{2}\displaystyle \frac{{v}_{\phi }}{r}\cos \phi \right]\\ & & +{\cos }^{2}i\left[\displaystyle \frac{1}{r}\displaystyle \frac{{v}_{{esc}}}{(H/R)}\right].\end{array}\end{eqnarray} \tag{ 7 }$$

Rescaling v_ϕ into units of r_s gives us ${v}_{\phi }=\sqrt{\tfrac{1}{2r^{\prime} }}$ and ${v}_{{esc}}=\sqrt{\tfrac{1}{r^{\prime} }}$ (where $r^{\prime} =r/{r}_{s}$, so $r=r^{\prime} {r}_{s}$), which, after simplifying, gives us

$$\begin{eqnarray}\begin{array}{rcl}\hat{n}\cdot {\boldsymbol{\Lambda }}\cdot \hat{n} & = & \displaystyle \frac{1}{{r}_{s}}\sqrt{\displaystyle \frac{1}{2r{{\prime} }^{3}}}\left[3{\sin }^{2}i\cos \phi \left(\sqrt{2}\cos \phi +\displaystyle \frac{\sin \phi }{2}\right)\right.\\ & & \left.-3\sqrt{2}\sin i\cos i\cos \phi +\displaystyle \frac{\sqrt{2}}{(H/R)}{\cos }^{2}i\right].\end{array}\end{eqnarray} \tag{ 8 }$$

There are four different possible ϕ dependencies (within the angle brackets) as a result of these shears within the disk, which we qualitatively describe below:

• 1.  
  The first term (f₁ in Figure 1, $\propto {\cos }^{2}\phi$) largely describes radial shear from the wind, where the angular dependence allows photons with small Doppler shifts to escape more easily to the observer from regions with large radial shears at the near and far sides of the disk (ϕ ≈ 0 or π, where the line-of-sight projected velocities are small).
• 2.  
  The second term (f₂ in Figure 1, $\propto \cos \phi \sin \phi$) describes gradients caused by Keplerian shear in the disk, which by itself produces a double-peaked line profile (with peaks corresponding to the blue and red sides of the disk), replicating the "M profile" first shown in Rybicki & Hummer (1983, hereafter RH83) and further discussed in the context of double-peaked line profiles from cataclysmic variables by Horne & Marsh (1986).
• 3.  
  The third term(f₃ in Figure 1, $\propto \cos \phi$) represents the "lifting" shear as a function of radius, where again the angular dependence allows photons with small Doppler shifts to escape more easily (although less strongly than in the case of radial shear).
• 4.  
  The final term (f₄ in Figure 1) represents the "lifting" shear as a function of height off of the disk, and it interestingly has no ϕ dependence, meaning it represents a form of isotropic emission that by itself produces a doubly peaked line profile, albeit of a different shape than the profile given by just the Keplerian shear. For thin disks H ≪ R, and thus even a small velocity gradient in the θ direction will be amplified greatly by this term, so to keep its magnitude similar to the other terms we absorb this dependence into f₄ (i.e., ${f}_{4}=\tfrac{C}{(H/R)}$, where C is the unamplified wind contribution in this direction).

This is the crux of our model, as the emission intensity at each location in the disk is set by $| \tfrac{{{dv}}_{l}}{{dl}}|$ alone in the optically thick limit, as shown in Equation (4).

To better explore the parameter space, we want to examine how each of these terms affects the line profile shape. Equation (4) is directly proportional to the line-of-sight velocity gradient—a result of our assumption that τ ≫ 1 in the Sobolev approximation—and this makes it the critical component that shapes the line and response profiles. Discarding the normalization and primes, we can explore how this line-of-sight velocity gradient changes as a function of r, i, and ϕ alone:

$$\begin{eqnarray}&&\begin{array}{l}\hat{n}\cdot {\boldsymbol{\Lambda }}\cdot \hat{n}\approx \displaystyle \frac{{{dv}}_{l}}{{dl}}\approx \sqrt{\displaystyle \frac{1}{2{r}^{3}}}\left(3{\sin }^{2}i\cos \phi \left[\sqrt{2}{f}_{1}\cos \phi +{f}_{2}\displaystyle \frac{\sin \phi }{2}\right]\right.\\ \quad -\left.3{f}_{3}\sin i\cos i\cos \phi +\sqrt{2}{f}_{4}{\cos }^{2}i\right).\end{array}\end{eqnarray} \tag{ 9 }$$

Neglecting physical constants, this completes the mathematical description of our model, where f_1,2,3,4 are artificial constants that allow the fit to be more flexible, varying the strength of the different wind launching terms described above. In the high-inclination limit only the sin² i term is important, and our result reduces to the form considered in CM96 and is overall similar to that presented in MC97. At low inclinations the cos² i term is very important, and the addition of both terms at moderate inclinations make a significant difference when compared with the results shown in CM96 as demonstrated originally in MC97.

Figure 2 explores the general behavior of these terms and how they affect the shape of the line and phase profiles. Qualitatively, the terms have the following effects on the line profile:

• 1.  
  f₁ creates a symmetric bell-curve-like shape about the line center, as it is proportional to ${\cos }^{2}\phi$.
• 2.  
  f₂ in isolation creates an "M"-shaped line profile with two peaks, with the line center as a minimum. This replicates the result first shown in RH83's Figure 2. In conjunction with other terms, it acts to depress the region just to the left of the center of the line profile and to raise the region just to the right of the center as it is proportional to $\cos \phi \sin \phi$, providing a red–blue asymmetry to the line profile.
• 3.  
  f₃ in isolation creates a flat-topped line profile with the plateau centered on the line center, but in conjunction with other terms it acts to widen the line profile and alter the shape of the peak as it is proportional to $\cos \phi$ and in our definition the observer is at ϕ = 0, meaning there is no left/right asymmetry.
• 4.  
  f₄ creates a symmetric "twin-horned"-type feature, as it is an isotropic emission feature and thus only depends on the Doppler shift, i.e., the frequency shift is proportional only to $\sqrt{\tfrac{1}{2r}}\sin i\cos \phi$. This is analogous to the dashed–dotted profile shown in CM96's Figure 2. Its shape is thus entirely dependent on the delta-function Doppler resonance: gas on the bluer/redder sides of the disk resonates in blue/red wavelengths and thus only contributes on the corresponding side of the line profile, while the gas with no Doppler shift close to the front and back sides of the disk keeps the line profile from going to zero at line center.

Zoom In Zoom Out Reset image size

Regardless of whether the individual term would produce a single- or double-peaked line profile, all of the phase profiles display the standard "S" shape, but with different morphologies. This pattern in the phase profile indicates a red–blue asymmetry in the emission-line centroids, as expected for ordered rotation, with the differential phase being a measurement made possible through interferometry with the GRAVITY instrument on the Very Large Telescope Interferometer (GRAVITY Collaboration et al. 2017). As shown in the left panel of Figure 2, when all the terms of are of equal strength f₁ and f₄ dominate.

Our code uses this model prescription to generate a two-dimensional ray-traced image of an inclined disk assumed to represent the BLR, employing a polar coordinate grid with logarithmic radial spacing. The intensity at each grid cell in the disk is given by Equations (4) and (9), and we can calculate the line profile numerically as described above (with the summed intensities at each disk location weighted by their corresponding area element). In fitting, we allow the terms to vary independently with priors 0 ≤ f_1,2,3,4 ≤ 1, consistent with the disk launching outflows. The differential phase is calculated following the standard BLR photocenter and kinematic modeling prescription in the marginally resolved limit, where we keep only the first-order term in the expansion of the complex visibility such that ${\rm{\Delta }}{\rm{\Phi }}=-2\pi \left(u\cdot \bar{x}\right)\tfrac{f}{1\,+\,f}$ (Lachaume 2003; Waisberg et al. 2017). Here, u are the interferometric baselines, $\bar{x}$ the on-sky emission centroids, and f the normalized line flux such that $\tfrac{f}{1+f}$ represents the contrast between the line and the continuum. Note that this Δϕ does not refer to any physical ϕ in the disk and is only the differential phase angle.

We use data previously published by the GRAVITY collaboration on the quasar 3C 273 to fit our model; for a detailed description of the observations and reduction techniques used, see G18. 3C 273 is uniquely suited to test the model because the system's distance is close enough that GRAVITY obtains a spectro-astrometric differential phase signature across the broad-line profile, which we can fit for in addition to the shape of just the line profile (GRAVITY Collaboration et al. 2017; G18). 3C 273 is also oriented such that we observe the jet, allowing us to constrain what we believe the "true" inclination of the system to be (Kundt & Gopal-Krishna 1986). The fitting is done with flux and phase measurements along each wavelength channel (measurements taken in 40 channels between ∼2.13 and ∼2.22 μm over six baselines at four different epochs), using Markov Chain Monte Carlo (MCMC) methods. To ensure we sample a large region of the parameter space, we employ a parallel tempered MCMC method developed by Foreman-Mackey et al. (2013) and Vousden et al. (2016), using six different logarithmically spaced temperatures each with 24 walkers. The model as we fit to the data is fully described by the following 11 parameters:

• 1.  
  The inclination angle of the system, i, where i = 90° corresponds to an edge-on viewing angle and i = 0° a "face-on" viewing angle, as shown in Figure 1. Higher values of i lead the first term in Equation (6) to be the dominant drivers of the line and phase profiles, with lower values leading the last term to be most significant.
• 2.  
  The mass of the central supermassive black hole, M_BH. Increasing the mass of the black hole increases the amplitude of the phase profile.
• 3.  
  The mean radius of the BLR as weighted by the emissivity, $j(r)=k(r)\beta (r)S(r)\approx \tfrac{{\nu }_{0}}{c}| \tfrac{{{dv}}_{l}}{{dl}}| S(r)$, which for our scaling of S(r) gives j(r) ∝ r^−5/2 and thus

  $$\begin{eqnarray*}&&\bar{r}=\displaystyle \frac{{\int }_{{r}_{{\rm{\min }}}}^{{r}_{{\rm{\max }}}}{rj}(r){dr}}{{\int }_{{r}_{{\rm{\min }}}}^{{r}_{{\rm{\max }}}}j(r){dr}}=3\left(\displaystyle \frac{{r}_{{\rm{\max }}}^{-1/2}-{r}_{{\rm{\min }}}^{-1/2}}{{r}_{{\rm{\max }}}^{-3/2}-{r}_{{\rm{\min }}}^{-3/2}}\right).\end{eqnarray*}$$

  Higher values of $\bar{r}$ lead to the line/phase profiles being "squeezed" in wavelength space.
• 4.  
  A size scaling factor, r_fac, which in conjunction with $\bar{r}$ gives the minimum and maximum radii of the BLR via ${r}_{{\rm{\min }}}=\tfrac{\bar{r}}{3}\left(\tfrac{{r}_{{\rm{fac}}}^{-3/2}-1}{{r}_{{\rm{fac}}}^{-1/2}-1}\right)$ and ${r}_{{\rm{\max }}}={r}_{{\rm{fac}}}{r}_{{\rm{\min }}}$. Increasing r_fac slightly "stretches" the line/phase profiles in wavelength space, and also steepens the slope in the S-curve of the phase profile connecting the negative and positive peaks. This also increases the total flux of the line, but this does not affect the fit as we are only seeking to match the characteristic line shapes.
• 5.  
  The proportional strength of the radial shear wind term, f₁.
• 6.  
  The proportional strength of the Keplerian shear term, f₂.
• 7.  
  The proportional strength of the radial lifting shear term, f₃.
• 8.  
  The proportional strength of the height lifting shear term, f₄.
• 9.  
  The rotation of the model with respect to the orientation of the baselines in the data, θ_PA, reported in the standard convention in reference to the orientation of the jet (90° offset from the disk).
• 10.  
  A parameter, n, that can vary the normalization of the line profile with respect to the data slightly as the data points may not be exactly at the peak of the line, n ≥ 1, where n = 1 corresponds to scaling the model exactly to the maximum flux measurement in the data.
• 11.  
  A parameter, Δλ_c , that varies the line center, thus slightly shifting the models left and right in λ space. Here, we model Paα line emission, which has a known center near λ_c ≈ 2.172 μm for 3C 273 at a redshift of ∼0.16.

After fitting the model as described in the previous section, we reached convergence after ∼10,000 iterations in each walker at each temperature. We consider the fit converged when the maximum autocorrelation time of any parameter is 1% of the total number of steps taken. We show only the lowest temperature in our results presented below as the upper temperatures are designed to explore the parameter space and "trickle" down to the lowest temperature for further refinement (Vousden et al. 2016). Most importantly, the fit generally prefers higher inclinations and thus lower black hole mass. Table 1 presents the mean values in the fit with 1σ confidence intervals.

Table 1. MCMC Fit Parameters

i (°)	MBH (107 M⊙)	$\bar{r}\ ({r}_{s})$	rfac	f1	f2	f3	f4	θPA (°)	n	Δλc (μm)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
${78}_{-15}^{+4.9}$	${8.4}_{-1.7}^{+0.93}$	${6400}_{-1100}^{+320}$	${47}_{-7.6}^{+3.3}$	${0.69}_{-0.32}^{+0.14}$	${0.80}_{-0.53}^{+0.11}$	${0.65}_{-0.51}^{+0.17}$	${0.44}_{-0.40}^{+0.23}$	${240}_{-8.1}^{+3.9}$	${1.0}_{-0.008}^{+0.005}$	${0.00032}_{-0.0}^{+0.0}$

Note. Means with 1σ percentile confidence intervals (to two significant figures) on each of our 11 fit parameters, each of which is fully described in the previous section (the ordering 1–11 matches the ordering of the table). f₃ and f₄ are particularly poorly constrained, a result of the higher inclination preference of the sampler.

Download table as:  ASCIITypeset image

Figure 3 shows our best fit to the line and phase profiles, with the fainter red lines showing the distribution of the samples represented in Table 1. Appendix C shows a corner plot (Figure 9) of all of our parameters with their associated one-dimensional histograms. The fit is good, with a reduced χ² value of ∼1.36 for the mean parameters (and ∼1.35 for the best fit). The fit prefers higher inclinations, leading to only the f_1,2 terms being significant in the fit. The uncertainty on the importance of the wind terms with respect to each other is large and thus it is difficult to draw conclusions on the importance of any aspect over another, aside from the larger importance in f_1,2 that is largely driven by the inclination dependence. However, this does appear to match observational evidence that there may be a large radial-velocity component in any disk-wind outflows (Vestergaard et al. 2000).

Zoom In Zoom Out Reset image size

The black hole mass and on-sky position angle are the physical parameters best constrained by our model, but one should note that the black hole mass is strongly correlated with inclination. If we restrict the sampler to inclinations of less than 45°, we find a best low-inclination fit at $i={37}_{-10}^{+4.6}$ deg, which then prefers a higher black hole mass of ${M}_{\mathrm{BH}}\,={2.9}_{-1.3}^{+1.0}\times {10}^{8}{M}_{\odot }$. In either case the inferred mean BLR size is of order ∼20 μas. Taking the same distance and observed luminosity of 3C 273 as in G18 mean that the black hole mass as presented in Table 1 implies the system is super-Eddington by a factor of a few (Greenstein & Schmidt 1964), but this can be rectified if we restrict the prior to low inclinations, as shown above. Figure 4 shows how our model centroids compare to those in the data, in good agreement with ordered rotation around the jet.

Zoom In Zoom Out Reset image size

Radio observations of the jet indicate that the true inclination angle of 3C 273 is ∼20° ± 10° (Kundt & Gopal-Krishna 1986). This disfavors the specific disk-wind model considered here, indicating the cloud model previously presented in G18 remains the best-fit model to the BLR in 3C 273. However, we do not rule out disk-wind models completely, both for other active galactic nuclei (AGNs) and/or in the case of different morphologies/kinematics, such as those considered by Flohic et al. (2012), Chajet & Hall (2013), Waters16, Baskin & Laor (2018), Matthews et al. (2020), Naddaf & Czerny (2022), and others. It remains to be explored whether other kinds of disk-wind models such as these can accurately fit GRAVITY data and reconcile the inclination discrepancy.

If the disk-wind model presented here correctly described the BLR physics, our model would predict a black hole mass lower by a factor of ∼5 compared with results previously published in G18, while implying an angular size of the BLR that is smaller by a factor of ∼2. The reduced χ² presented for the cloud model fit in G18 is given as ∼1.3, which to two significant figures is slightly better than our reduced χ² ≈ 1.35. Still, both models appear to fit the data with roughly the same quality, and without external knowledge of the inclination angle this would imply possible systematic errors in estimating both the size of the BLR and the mass of the central black hole from interferometry data by factors of ∼2 and ∼5, respectively.

If we use a prior that allows only for low-inclination models, we find a low-inclination fit with a black hole mass that increases to ∼3 × 10⁸ M_⊙, in better agreement with the cloud model mass, but the on-sky size remains ∼20 μas. This implies that the black hole mass is roughly model independent and is instead simply strongly correlated with the inclination angle, a widely noted and expected correlation in astronomy. The BLR size is a model-dependent systematic uncertainty, however, indicating that the models produce different values for the virial factor f in Equation (1). This low-inclination best fit has a reduced χ² ≈ 1.37, which is slightly worse than the best fit at higher inclinations described but still acceptable. The average low-inclination fit is much poorer, with a reduced χ² ≈ 2.3. This lower-inclination best fit suffers from a problem of fine-tuning, however, as it essentially forces the f₄ and f₃ terms to 0 in order to fit the data (in the low-inclination case, the average fit results for these terms are—with 1σ confidence intervals—${f}_{4}={0.011}_{-0.01}^{+0.009}$ and ${f}_{3}={0.06}_{-0.038}^{+0.013}$, while ${f}_{1}={0.63}_{-0.22}^{+0.091}$ and ${f}_{2}={0.87}_{-0.50}^{+0.048}$). As Figure 5 shows, the f₁ term must dominate the model in order to produce the broad single peak in the line profile; any deviation from a Keplerian thin disk with a radial-velocity outflow thus breaks the fit at lower inclinations. The low-inclination fit to the data also fits the position angle signficantly worse, preferring ${\theta }_{\mathrm{PA}}={257}_{-13}^{+1.9}$°, in disagreement with the measured value of roughly 220°. In comparing to G18, however, this result is interesting because the kinematics they consider in their cloud model allow for only Keplerian rotation/shears, thus the kinematics that are required by our low-fit model (radial and Keplerian shears only) are not entirely at odds with those required by the cloud model fit.

Zoom In Zoom Out Reset image size

We also calculate an echo image of our best fit (shown in Figure 6), to compare to reverberation-mapping techniques and illustrate a mapping of the resonance conditions within our disk that produce the line and response profiles. Here, we plot the change in frequency, Δv, in units of megameters per second, i.e., ${\rm{\Delta }}v=c{\rm{\Delta }}\nu =-{v}_{\phi }\sin \phi \sin i$. The mean light travel time to our best-fit value for ${\bar{r}}_{{BLR}}$ is ≲50 days, in both the low- or high-inclination fit cases.

Zoom In Zoom Out Reset image size

The most recent results, published by Zhang et al. (2019), have measured time delays of ∼145 ± 10 days in the rest frame of 3C 273 using the Hβ and Hγ, while older results favor longer time delays of ∼300 days (Kaspi et al. 2000; Peterson et al. 2004). Previously published results in G18 produced a characteristic time delay of ∼145 ± 35 days at their value of $\bar{r}$ (corresponding to ∼50 μas on the sky), a size roughly twice as large as ours. While the time delay to our value of $\bar{r}$ is lower than this, this does not necessarily indicate that reverberation-mapping data are in conflict with this model, as this value is weighted by the intensity and not the responsivity of the disk. To do this would require a complete photoionization model of the BLR, and is complicated by many factors as described in Goad et al. (1993) and Goad & Korista (2014) and papers cited therein. A better estimate could be made by convolving our model transfer function, Ψ(t) (as shown at right in Figure 6), with the continuum lightcurve of 3C 273, with the peak of the resulting cross-correlation function indicating the characteristic delay, but we leave this for future work. We include the delay number discussed here as a comparison to the previously published work in G18, and note that it is roughly consistent with reverberation-mapping data for such a rough estimate. In future work it may be possible to jointly fit reverberation-mapping and GRAVITY data, in which case one would essentially fit for the transfer function in addition to the line profile and phase data.

Our model is simple, including only four possible wind launching terms simplified under the Sobolev approximation to rely on strong velocity gradients motivated by the local escape speed in the optically thick regime. It may be possible that further extensions/modifications of the model we use could help to mitigate the problems discussed here, but we leave this for future work. As it stands, it appears that the simple disk-wind model considered in this work is in tension with the assumption that type 1 AGNs are generally viewed at low inclinations, as we only expect to see single-peaked lines at higher inclinations (or the disk-wind model must be exceedingly well tuned at lower inclinations). Thus, we plan to test different kinds of disk-wind models to see if this tension can be resolved, additionally including further resolved GRAVITY sources as well as reverberation-mapping data in our fits.

We have fit a simple disk-wind launching model to GRAVITY data, showing that such a model can fit the line profiles, phase profiles, and emission centroids observed by GRAVITY. Our fit results in a smaller black hole mass for 3C 273 (≈8 × 10⁷ M_⊙) than other modeling results, which would make the system likely super-Eddington given its observed luminosity (Greenstein & Schmidt 1964), but this smaller inferred mass is a result of our fit preferring higher inclinations, which observations of the radio jet do not support. Our fit also prefers a smaller size for R_BLR, which is independent of inclination and thus a model-dependent uncertainty. Combining these two uncertainties would lead to larger systematic errors in inferring the black hole mass (see Equation (1)), but we emphasize that (at least in the case of 3C 273) this disk-wind model is disfavored.

In order to create a single peak in the line profile our fit requires the f₁ term be dominant over f₃ and f₄, which can be obtained by either forcing the viewing inclination of the system to be high or by fine-tuning the wind at low inclinations, such that only a Keplerian thin disk with radial outflows is present. Overall, our fit prefers a higher inclination, and our best-fit value is much higher than what is inferred for 3C 273 from observations of the jet. Our fitting results prefer an inclination angle of ∼75°, while the jet for 3C 273 indicates an inclination angle of ∼20° ± 10° (Kundt & Gopal-Krishna 1986). Assuming the jet and the disk are not significantly misaligned, the cloud model presented by G18 appears to thus better match the data. The phase and line profile fits in that work are at least as good as in our model but the cloud model produces a fit at much lower inclinations (∼12°), as expected given the observational constraints on the jet orientation of 3C 273. We show that if we restrict the sampler to lower inclinations we obtain a fit that is plausible, albeit slightly statistically worse than the higher-inclination fit, but that this lower-inclination fit suffers from a fine-tuning problem in the wind terms.

Our model is simple, and there may be extensions that further improve the fit which we have not considered in this work, such as those considered by Flohic et al. (2012), Chajet & Hall (2013), Waters16, Baskin & Laor (2018), Matthews et al. (2020), and Naddaf & Czerny (2022). Furthermore, it may be possible that not all quasar BLRs are governed by the same physics, and while 3C 273 may not be governed by disk-wind launching dynamics other quasars may still be. We hope to extend this work to include other disk-wind morphologies and kinematics, as well as more robustly test these models with reverberation-mapping data in addition to GRAVITY data. Based on these initial results, however, it seems difficult to fit a disk-wind model such as this to any type 1 AGNs as they are preferentially viewed at low inclinations. If the evidence continues to favor the cloud model as this result does, it is of increasing importance to try to better understand the physical processes that can result in the BLR being best modeled as a distribution of cold and dense puffed-up clouds of atomic gas, as well as how this picture can be connected to the strong observational evidence for outflows.

This work was supported in part by NSF grant No. AST-1909711 and an Alfred P. Sloan Research Fellowship (J.D.). We are grateful to the computing resources made available to us by Research Computing at CU Boulder, as this work utilized the Summit supercomputer, which is supported by the National Science Foundation (awards ACI-1532235 and ACI-1532236), the University of Colorado Boulder, and Colorado State University. The Summit supercomputer is a joint effort of the University of Colorado Boulder and Colorado State University. K.L. is especially grateful to Sajal Gupta for many helpful discussions over the course of the project, and to Marcel Corchado-Abelo for discussions on the rate of strain tensor. We are grateful to the anonymous referee for their thorough report, which greatly improved the quality of the paper. The code used in this work is available free and open-source on GitHub; ³ a copy of this material has been deposited to Zenodo: doi:10.5281/zenodo.8117776.

Facility: VLTI(GRAVITY)(GRAVITY Collaboration et al. 2017) - The Very Large Telescope Interferometer.

Software: Julia, ptemcee, python.

As discussed in the text, the line-of-sight velocity gradient $\tfrac{{{dv}}_{l}}{{dl}}$ can be approximated with the rate of strain tensor. In computing $\tfrac{{{dv}}_{l}}{{dl}}\approx \hat{n}\cdot {\boldsymbol{\Lambda }}\cdot \hat{n}$, as given in Equation (6), we need the components of the rate of strain tensor Λ_ij in spherical geometry, which we obtain from Batchelor (1968). As in the text (see Equation (5)), we use

$$\begin{eqnarray}&&\hat{n}=(\sin \theta \cos \phi \sin i+\cos \theta \cos i)\hat{r}+(\cos \theta \cos \phi \sin i-\sin \theta \cos i)\hat{\theta }-(\sin \phi \sin i)\hat{\phi }.\end{eqnarray} \tag{ A1 }$$

Using this, we evaluate $\hat{n}\cdot {\boldsymbol{\Lambda }}\cdot \hat{n}$ as

$$\begin{eqnarray}\begin{array}{rcl}\hat{n}\cdot {\boldsymbol{\Lambda }}\cdot \hat{n} & = & {{\rm{\Lambda }}}_{{rr}}{\left(\sin \theta \cos \phi \sin i+\cos \theta \cos i\right)}^{2}+{{\rm{\Lambda }}}_{\theta \theta }{\left(\cos \theta \cos \phi \sin i-\sin \theta \cos i\right)}^{2}+{{\rm{\Lambda }}}_{\phi \phi }{\left(\sin \phi \sin i\right)}^{2}\\ & & +2{{\rm{\Lambda }}}_{r\theta }\left(\sin \theta \cos \phi \sin i+\cos \theta \cos i\right)\left(\cos \theta \cos \phi \sin i-\sin \theta \cos i\right)\\ & & -2{{\rm{\Lambda }}}_{\theta \phi }\left(\cos \theta \cos \phi \sin i-\sin \theta \cos i\right)\left(\sin \theta \sin i\right)-2{{\rm{\Lambda }}}_{r\phi }\left(\sin \theta \cos \phi \sin i+\cos \theta \cos i\right)\left(\sin \phi \sin i\right).\end{array}\end{eqnarray} \tag{ A2 }$$

Applying the approximation that the disk is very thin and at the midplane, we set $\theta =\tfrac{\pi }{2}$, reducing Equation (A2) to

$$\begin{eqnarray}\begin{array}{rcl}\hat{n}\cdot {\boldsymbol{\Lambda }}\cdot \hat{n} & = & {{\rm{\Lambda }}}_{{rr}}{\left(\cos \phi \sin i\right)}^{2}-{{\rm{\Lambda }}}_{\theta \theta }{\left(\cos i\right)}^{2}+{{\rm{\Lambda }}}_{\phi \phi }{\left(\sin \phi \sin i\right)}^{2}-2{{\rm{\Lambda }}}_{r\theta }\left(\cos \phi \sin i\right)\left(\cos i\right)\\ & & +2{{\rm{\Lambda }}}_{\theta \phi }\left(\cos i\right)\left(\sin i\right)-2{{\rm{\Lambda }}}_{r\phi }\left(\cos \phi \sin i\right)\left(\sin \phi \sin i\right).\end{array}\end{eqnarray} \tag{ A3 }$$

Simplifying, this leaves us with

$$\begin{eqnarray}&&\hat{n}\cdot {\boldsymbol{\Lambda }}\cdot \hat{n}={\sin }^{2}i\left({{\rm{\Lambda }}}_{{rr}}{\cos }^{2}\phi +{{\rm{\Lambda }}}_{\phi \phi }{\sin }^{2}\phi -2{{\rm{\Lambda }}}_{r\phi }\sin \phi \cos \phi \right)-\sin i\cos i\left(2{{\rm{\Lambda }}}_{r\theta }\cos \phi -2{{\rm{\Lambda }}}_{\theta \phi }\sin \phi \right)+{{\rm{\Lambda }}}_{\theta \theta }{\cos }^{2}i.\end{eqnarray} \tag{ A4 }$$

From Batchelor (1968), the Λ_ij terms in spherical coordinates are

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Lambda }}}_{{rr}}\,=\,\displaystyle \frac{\partial {v}_{r}}{\partial r}{\rm{;}}{{\rm{\Lambda }}}_{\theta \theta }=\displaystyle \frac{1}{r}\displaystyle \frac{\partial {v}_{\theta }}{\partial \theta }+\displaystyle \frac{{v}_{r}}{r}{\rm{;}}{{\rm{\Lambda }}}_{\phi \phi }\,=\,\displaystyle \frac{1}{r\sin \theta }\displaystyle \frac{\partial {v}_{\phi }}{\partial \phi }+\displaystyle \frac{{v}_{r}}{r}+\displaystyle \frac{{v}_{\theta }\cot \theta }{r}{\rm{;}}\\ {{\rm{\Lambda }}}_{r\theta }=\displaystyle \frac{1}{2}\left(r\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{v}_{\theta }}{r}\right)+\displaystyle \frac{1}{r}\displaystyle \frac{\partial {v}_{r}}{\partial \theta }\right){\rm{;}}{{\rm{\Lambda }}}_{\theta \phi }\,=\displaystyle \frac{1}{2}\left(\displaystyle \frac{\sin \theta }{r}\displaystyle \frac{\partial }{\partial \theta }\left(\displaystyle \frac{{v}_{\phi }}{\sin \theta }\right)+\displaystyle \frac{1}{r\sin \theta }\displaystyle \frac{\partial {v}_{\theta }}{\partial \phi }\right){\rm{;}}\\ {{\rm{\Lambda }}}_{r\phi }=\displaystyle \frac{1}{2}\left(\displaystyle \frac{1}{r\sin \theta }\displaystyle \frac{\partial {v}_{r}}{\partial \phi }+r\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{v}_{\phi }}{r}\right)\right).\end{array}\end{eqnarray} \tag{ A5 }$$

Noting again that $\theta =\tfrac{\pi }{2}$, and assuming that the disk is axisymmetric in ϕ (thus all of the $\tfrac{\partial }{\partial \phi }$ terms are 0) gives

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Lambda }}}_{{rr}}\,=\,\displaystyle \frac{\partial {v}_{r}}{\partial r}{\rm{;}}{{\rm{\Lambda }}}_{\theta \theta }\,=\,\displaystyle \frac{1}{r}\displaystyle \frac{\partial {v}_{\theta }}{\partial \theta }+\displaystyle \frac{{v}_{r}}{r}{\rm{;}}{{\rm{\Lambda }}}_{\phi \phi }\,=\,\displaystyle \frac{{v}_{r}}{r}{\rm{;}}{{\rm{\Lambda }}}_{r\theta }=\displaystyle \frac{1}{2}\left(r\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{v}_{\theta }}{r}\right)+\displaystyle \frac{1}{r}\displaystyle \frac{\partial {v}_{r}}{\partial \theta }\right){\rm{;}}\\ {{\rm{\Lambda }}}_{\theta \phi }\,=\,\displaystyle \frac{1}{2}\left(\displaystyle \frac{\sin \theta }{r}\displaystyle \frac{\partial }{\partial \theta }\left(\displaystyle \frac{{v}_{\phi }}{\sin \theta }\right)\right){\rm{;}}{{\rm{\Lambda }}}_{r\phi }\,=\,\displaystyle \frac{1}{2}\left(r\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{v}_{\phi }}{r}\right)\right).\end{array}\end{eqnarray} \tag{ A6 }$$

We now apply a final approximation to the disk in which we assume that v_r ≈ v_θ ≈ 0, in keeping with the standard assumptions for thin disks. Plugging this result into Equation (13) then gives us

$$\begin{eqnarray}\begin{array}{rcl}\hat{n}\cdot {\boldsymbol{\Lambda }}\cdot \hat{n} & = & {\sin }^{2}i\left(\displaystyle \frac{\partial {v}_{r}}{\partial r}{\cos }^{2}\phi -\left(r\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{v}_{\phi }}{r}\right)\right)\sin \phi \cos \phi \right)\\ & = & -\sin i\cos i\left(\left(r\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{v}_{\phi }}{r}\right)+\displaystyle \frac{1}{r}\displaystyle \frac{\partial {v}_{r}}{\partial \theta }\right)\cos \phi \right.\left.\,-\,\left(\displaystyle \frac{\sin \theta }{r}\displaystyle \frac{\partial }{\partial \theta }\left(\displaystyle \frac{{v}_{\phi }}{\sin \theta }\right)\right)\sin \phi \right)\,=\,+\displaystyle \frac{1}{r}\displaystyle \frac{\partial {v}_{\theta }}{\partial \theta }{\cos }^{2}i.\end{array}\end{eqnarray} \tag{ A7 }$$

We now turn our attention to the derivatives. We assume Keplerian orbits such that ${v}_{\phi }=\sqrt{\tfrac{{GM}}{r}}$ and thus $\tfrac{\partial {v}_{\phi }}{\partial r}=\tfrac{-{v}_{\phi }}{2r}$. Applying this in conjunction with the chain rule, we can write $r\tfrac{\partial }{\partial r}\left(\tfrac{{v}_{\phi }}{r}\right)=\tfrac{\partial {v}_{\phi }}{\partial r}-\tfrac{{v}_{\phi }}{r}=-\tfrac{3{v}_{\phi }}{2r}$. Similarly, $r\tfrac{\partial }{\partial r}\left(\tfrac{{v}_{\theta }}{r}\right)=\tfrac{\partial {v}_{\theta }}{\partial r}-\tfrac{{v}_{\theta }}{r}$ and $\tfrac{\sin \theta }{r}\tfrac{\partial }{\partial \theta }\left(\tfrac{{v}_{\phi }}{\sin \theta }\right)=\tfrac{1}{r}\tfrac{\partial {v}_{\phi }}{\partial \theta }-\tfrac{{v}_{\phi }}{r}\cot \theta =\tfrac{1}{r}\tfrac{\partial {v}_{\phi }}{\partial \theta }$.

This reduces Equation (A7) to

$$\begin{eqnarray}\begin{array}{lcl}\hat{n}\cdot {\boldsymbol{\Lambda }}\cdot \hat{n} & = & {\sin }^{2}i\left(\displaystyle \frac{\partial {v}_{r}}{\partial r}{\cos }^{2}\phi +\displaystyle \frac{3{v}_{\phi }}{2r}\sin \phi \cos \phi \right)\\ & & -\sin i\cos i\left(\left(\displaystyle \frac{\partial {v}_{\theta }}{\partial r}-\displaystyle \frac{{v}_{\theta }}{r}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial {v}_{r}}{\partial \theta }\right)\cos \phi -\left(\displaystyle \frac{1}{r}\displaystyle \frac{\partial {v}_{\phi }}{\partial \theta }\right)\sin \phi \right)\,+\displaystyle \frac{1}{r}\displaystyle \frac{\partial {v}_{\theta }}{\partial \theta }{\cos }^{2}i.\end{array}\end{eqnarray} \tag{ A8 }$$

We now must apply some kind of physical prescriptions to the remaining derivatives, which we do in keeping with the Sobolev approximation that there must be large velocity gradients present. Using this, we approximate the gradients in terms of local escape velocities, i.e., $\tfrac{\partial {v}_{r}}{\partial r}\approx 3\sqrt{2}\tfrac{{v}_{\phi }}{r}$, $\tfrac{\partial {v}_{\theta }}{\partial \theta }\approx \tfrac{{v}_{{esc}}}{\left(H/R\right)}$ (with H/R ≪ 1 being the scale height of the disk), and $\tfrac{\partial {v}_{\theta }}{\partial r}\approx \tfrac{\partial {v}_{r}}{\partial r}$. Here, we assume the local escape velocity is ${v}_{e}=\sqrt{\tfrac{2{GM}}{r}}=\sqrt{2}{v}_{\phi }$. This would seem to imply that $\tfrac{\partial {v}_{r}}{\partial r}=\sqrt{2}\tfrac{{v}_{\phi }}{r}$, but CM96 adopted an arbitrary extra factor of 3 (seemingly to assume the wind launching regions generate an outflow that travels at a substantial velocity with respect to the source) that we keep to better compare with their results. Since v_ϕ is function of r alone $\tfrac{\partial {v}_{\phi }}{\partial \theta }=0$, and finally we also set $\tfrac{\partial {v}_{r}}{\partial \theta }=0$ in keeping with the idea of a geometrically thin disk. This allows us to arrive at the form presented in Equation (7) in the text:

$$\begin{eqnarray}&&\hat{n}\cdot {\boldsymbol{\Lambda }}\cdot \hat{n}\approx \displaystyle \frac{{{dv}}_{l}}{{dl}}=3\displaystyle \frac{{v}_{\phi }}{r}{\sin }^{2}i\cos \phi \left[\sqrt{2}\cos \phi +\displaystyle \frac{\sin \phi }{2}\right]-\sin i\cos i\left[3\sqrt{2}\displaystyle \frac{{v}_{\phi }}{r}\cos \phi \right]+{\cos }^{2}i\left[\displaystyle \frac{1}{r}\displaystyle \frac{{v}_{{esc}}}{(H/R)}\right].\end{eqnarray} \tag{ A9 }$$

Note that there is a sign difference between the term proportional to $\sin i\cos i$ than in similar work done by Flohic et al. (2012), which is the result of different assumptions for the underlying kinematics/velocity gradient fields, but as noted therein this has minimal effect on the shape of the line profiles.

Here we show the phase data with average model phases (from Table 1) from each GRAVITY baseline at each observation, and illustrate why only a subset of these data are averaged to show the phase profile presented in Figures 3 and 5. The UV coverage of the observations and their orientation relative to the jet axis in 3C 273 is shown in Figure 7, where the black points correspond to the bolded "off" (the jet) axis phase measurements shown in Figure 8, while the lighter grey points correspond to the "on" axis measurements.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 9 we show the parameter space explored by the sampler and the correlations between parameters in the resulting fit, with the units the same as presented in Table 1. The non-Gaussian shape of several of the histograms illustrates the importance of using multiple temperatures in the MCMC fitting to ensure the sampler does not get stuck in a local peak. The best fit solution (shown with blue lines) is essentially the high inclination model originally considered by Chiang and Murray in CM96, but note the difference in the reduced χ² between the best fit and the average parameters is only ∼0.01. While the various wind launching terms are poorly constrained, it is interesting to note that it appears (from the 1D histograms) that the sampler prefers f₂ approach 1 and f₄ 0. f₂ should be 1 in an ideal thin disk, as the Keplerian shear is non-negotiable. f₄ represents a form of isotropic emission which creates a double-horned profile, which means the contributions from it need to be small in order to preserve the observed single peak.

Zoom In Zoom Out Reset image size