© The Authors 2024

1. Introduction

It is challenging to measure the mass of a supermassive black hole (BH) at the center of a galaxy because, ideally, one needs to resolve its sphere of influence (Kormendy & Ho 2013). This is true in the local Universe, where dynamical methods are preferred. At cosmic distances, one has to rely on scaling relations, except when the broad-line region (BLR) of an active galactic nucleus (AGN) is observable. The broad recombination lines with a typical full width at half maximum (FWHM) ≳1000 km s⁻¹ (Khachikian & Weedman 1974) are emitted by the ionized gas surrounding the accreting BH (Peterson 1997). The BH mass can be derived with the virial method once the size of the BLR has been measured, M_BH = f_virialRv²/G, where R is the BLR radius, v is a characteristic velocity of the BLR rotation, and f_virial is the virial factor that takes account of the geometry of the BLR. The structure and dynamics of the BLR strongly affect the virial factor and are critical to the BH mass measurement (Collin et al. 2006; Mejía-Restrepo et al. 2018).

The broad-line profile suggests that the BLR has a disk-like geometry (e.g., Wills & Browne 1986; Vestergaard et al. 2000; Kollatschny & Zetzl 2011; Shen & Ho 2014; Storchi-Bergmann et al. 2017). Its size is most often measured from the time lag between the AGN continuum and the broad emission line light curves by using the reverberation mapping (RM) technique (Blandford & McKee 1982; Peterson 2014). The characteristics of these time lags across different velocity channels have provided evidence of inflow and outflow motions in the BLR (e.g., Bentz et al. 2010a; Grier et al. 2013; Du et al. 2016). This has led to the development of comprehensive models that can constrain the BLR structure using high-quality RM data (Brewer et al. 2011; Pancoast et al. 2011, 2014a,b; Li et al. 2013). More recently, the BLR has been spatially resolved with spectro-astrometry, which is a powerful technique for measuring the BLR structure and BH mass (GRAVITY Collaboration 2018, 2020, 2021a, 2024), even out to the cosmic noon (Abuter et al. 2024). Attempts have been made to analyze spectro-astrometry and RM data jointly (hereafter, the SARM method) in order to measure the geometric distance of the BLR and better constrain the BLR structure and BH mass (Wang et al. 2020; GRAVITY Collaboration 2021b; Li et al. 2022).

The high-quality data needed for the detailed analyses described above are not widely available for large AGN samples, even with the ongoing large RM projects, such as SDSS-RM (Shen et al. 2015) and OzDES-RM (Malik et al. 2023). However, there is a wealth of AGN samples with good quality single-epoch spectra. To exploit these samples, Raimundo et al. (2019, 2020) modified the widely used BLR dynamical modeling code CARAMEL and used it to fit single-epoch line profiles. They were able to constrain some BLR model parameters, such as the inclination angle and disk thickness, and to estimate a BH mass by setting a prior on the BLR size based on the empirical size-luminosity relation (e.g., Bentz et al. 2013). We turn this idea around and focus on investigating the BLR structure and the virial factors derived from multiple broad lines, which are covered simultaneously in a spectrum across UV, optical and/or near-IR (NIR). By doing so, we can understand how the structure of the BLR changes between different lines within the same AGN. Previous RM observations found that higher ionization lines respond more promptly to continuum variations than lower ionization lines (e.g., Clavel et al. 1991; Gaskell 2009). The photoionization model, such as the “locally optimally emitting cloud” (LOC) model (Baldwin et al. 1995), can naturally produce such “radial ionization stratification”. Korista & Goad (2004) predicted decreasing time lags of Hα, Hβ, Hγ, He I, and He II using the LOC model, which were confirmed by the RM observation of nearby AGNs (Bentz et al. 2010b).

In this paper, we analyze the VLT/X-shooter (Vernet et al. 2011) spectrum of NGC 3783, which covers several strong, prominent broad emission lines of hydrogen and helium and for which the high spectral resolution enables a robust decomposition of the broad and narrow lines. These data were previously used in a study of BLR excitation and extinction in several AGN (Schnorr-Müller et al. 2016). For NGC 3783, Bentz et al. (2021a) reported the RM time lags of various lines, including Hβ and Hγ. Bentz et al. (2021b) performed dynamical modeling of Hβ and He II lines and derived a BLR size and BH mass consistent with the traditional method. GRAVITY Collaboration (2021a) reported the spectro-astrometry measurement of the broad Brγ line. This motivated a joint SARM analysis, which yielded consistent results (GRAVITY Collaboration 2021b). Here, we introduce the data and our method to decompose the broad-line profiles in Sect. 2. Then, we make a nonparametric characterization of them in Sect. 3. We discuss our modeling of the line profiles in Sect. 4, comparing our results with GRAVITY Collaboration (2021b) because the joint analysis provides the strongest constraint of the BLR model. We discuss the strengths and limitations of the current method in Sect. 5 and conclude in Sect. 6.

2. Data reprocessing of X-shooter spectra

2.1. Data reduction

The X-shooter data were acquired as part of the LLAMA project (Davies et al. 2015). A description of the observations and data reduction can be found in Schnorr-Müller et al. (2016) and Burtscher et al. (2021). We briefly summarize the key points as follows: NGC 3783 (11:39:01.7, −37:44:19.0) was observed with X-shooter at the Very Large Telescope in early 2014 using the IFU mode (program ID 092.B-0083). The spectral resolving power, R = λ/Δλ, is about 8400 (UVB), 13200 (VIS), and 8300 (NIR; Schnorr-Müller et al. 2016). The data were reduced with the ESO reflex pipeline (version 2.6.8) by using the Kepler GUI interface (Modigliani et al. 2010) and mostly the default configuration. The pipeline provides the data cubes of UVB, VIS, and NIR arms separately for each observation. Telluric and flux calibrator stars were also observed. Flux calibration was performed with a spectro-photometric standard from Moehler et al. (2014). From a comparison of the stars observed throughout the program, the spectrum was calibrated to an accuracy of about 2%. We extracted 1D spectra from each of the NGC 3783 data cubes using a rectangular slit with a width of 1.8″, and we applied minor scaling corrections to match the different spectral ranges.

We corrected the Galactic extinction based on A_V = 0.332 (Schlafly & Finkbeiner 2011) and the Cardelli et al. (1989) extinction model, where we specified R_V = 3.1 as the ratio of total to selective extinction, and we converted the spectrum to the rest frame, adopting a redshift of 0.00973 measured from the H I 21 cm line observations (Theureau et al. 1998). A final correction was made to the blue wing of the broad hydrogen lines in the UV arm due to absorption in the telluric star. We masked these narrow wavelength ranges (4286−4307 Å and 4800−4818 Å) when modeling the BLR profiles (Sect. 4). We also identified a few bad channels in the He Iλ5876 profile and some channels contaminated by absorption and emission lines of the sky in the Paβ profile. We masked these channels when we modeled the line profiles, noting that they are always much narrower than the broad-line profiles, so they will not affect the modeling. We present the optical and NIR parts of the spectrum that are relevant to this work in Fig. 1 together with the spectral decomposition that is discussed in Sect. 2.2.

Fig. 1. X-shooter spectrum of NGC 3783 simultaneously covering UV, optical, and NIR ranges. In each panel, the spectrum is shown with the full model that enables decomposition of the broad-line profiles overplotted, and the residual is shown underneath. The emission line components are also plotted separately. The gray-shaded regions in the residuals represent wavelength ranges of bad channels and features due to the poor telluric correction, which are masked in the fitting. Panel a: UV/optical spectrum with the Hα, Hβ, and Hγ lines. Panel b: part of the NIR spectrum used in this study, with the Paβ line. Panel c: [S II] doublet used as the narrow line template.

2.2. Spectral decomposition

We decomposed the broad-line components of Hα, Hβ, Hγ, and He Iλ5876 from the combined UVB and VIS spectra and the Paβ profile from the NIR spectrum. The final overall fit can be found in Fig. 1.

2.2.1. Narrow line template

The narrow line profiles of AGNs are usually more complex than a simple Gaussian profile, which argues for the use of a template based on isolated lines. The [S II] λλ6716,6731 and [O III] λλ4959,5007 lines are usually used for this purpose, and one can use multiple Gaussian components to generate a noiseless narrow line template (Ho et al. 1997). To generate the template, we fit each of the [S II] lines with three Gaussian profiles, tying their width and velocity shift from the laboratory wavelength between the pair but allowing the total scaling to vary. At the same time, we fit the local continuum with a third-order polynomial function. We found that neither more Gaussian components nor a higher-order polynomial function could improve the fitting. We adopted the [S II] doublet because the [O III] lines show stronger blue-shifted wind components. When trying a template based on [O III], we found it did not match the narrow line components of the H I lines well.

2.2.2. Decomposing the broad-line profiles

The continuum of an AGN optical spectrum consists of emission from the accretion disk, the host galaxy, and the pseudo-continuum of the Fe II lines (e.g., Barth et al. 2015). We adopted a power-law model for the AGN featureless continuum and found a host galaxy component based on Bruzual & Charlot (2003) stellar populations amounted to ≲10% of the continuum and was thus not needed. To fit the Fe II features, we incorporated the newly published high-quality template covering 4000−5600 Å (Park et al. 2022). We found the final normalized line profiles only change ≲0.05 if we adopt the Fe II template from Boroson & Green (1992), which has little effect on the BLR modeling results. We broadened and shifted the Fe II template as part of the fitting process. A wide wavelength range is useful for decomposing these components, so we opted to fit the entire optical spectrum over 4200−6800 Å and decompose the broad Hα, Hβ, Hγ, and He I profiles simultaneously (Fig. 1a). We describe fitting the more isolated Paβ line later in the text.

In the optical spectrum, the majority of narrow lines were fitted by the [S II] template with two free parameters, the amplitude and the velocity shift from the laboratory wavelength. We tied the velocity shifts of all the templates, so the fitting reflects a small deviation (≈ − 2 km s⁻¹) from the redshift measured by the atomic H I gas. Because the UVB arm has a lower spectral resolution than the VIS arm, we broadened the narrow line template with a Gaussian kernel (σ ≈ 35 km s⁻¹) for all lines at wavelengths shorter than 5600 Å. For the [O III] lines, which have a higher critical density and more contribution from blue-shifted components than [S II], we added additional Gaussian components: one for [O III] λ4363 and two (tied together and with a ratio of 2.98) for each line in the [O III] λλ4959,5007 doublet. The [N II] λ6550, 6585 doublets were fitted with the [S II] template with the amplitude ratio fixed to the theoretical value of 2.96. For completeness and to avoid influencing the continuum placement, we fitted several other narrow lines¹ in the spectrum, although they do not directly affect the broad line decomposition.

For the broad lines, we found three Gaussian components are sufficient to fit the Hα, Hβ, Hγ, and He I profiles. We fit the broad He IIλ4686 as well, although it is too faint to provide a robust line profile for our dynamical modeling.

Lastly, we multiplied the entire optical spectrum model by a fifth-order polynomial function to account for large-scale variations caused by the instrumental and calibration effects (Cappellari 2017). This method can moderately improve the fitting of the continuum in the line wings at a level of ∼10⁻¹⁶ erg s⁻¹ cm⁻² Å⁻¹.

In the NIR spectrum, we fitted the broad Paβ profile and continuum over a wavelength range of 12 400−13 300 Å. We included a power law for the continuum, the [S II] template (broadened to match the resolution of the NIR arm) for the narrow Paβ line, and three Gaussian components for the broad Paβ line. The velocity shift with respect to the theoretical wavelength of the narrow Paβ line was not tied to the optical narrow lines, but their velocity difference of ≲1.7 km s⁻¹ is consistent with the systematic uncertainty of X-shooter wavelength calibration over different arms². This high accuracy enabled us to tie the central wavelength offsets when we fit the broad-line profiles simultaneously (Sect. 4).

2.2.3. Resampling and uncertainties

While the high resolution of the X-shooter spectrum helps decompose the narrow- and broad-line profiles robustly, it is not necessary for the modeling, which is only constrained by the broad (≳1000 km s⁻¹) and smooth features. Therefore, we resampled the decomposed line profiles to an effective resolving power of 2000, which is high enough to retain the characteristic features of the profiles while reducing the computational time. We also verified that our conclusions do not change with a slightly lower (R = 1000) or higher (R = 4000) resolution.

We first convolved the decomposed broad-line profiles with Gaussian kernels to the required resolution. Then, we used the Python tool SpectRes (Carnall 2017) to resample them while preserving the integrated flux and propagating the uncertainties. We chose the channel width of the re-sampled profiles to be ∼75 km s⁻¹ so that the spectra are still Nyquist sampled.

To calculate the uncertainties of the line profiles, we summed two components in quadrature: (1) the uncertainty of the observed spectrum and the line profile decomposition, and (2) 5% of the line flux. To estimate the first component, we used a running root mean square (RMS) of the fitting residual of the original spectra, which includes the imperfectness of the fitting and other potential artifacts. The second component is important to avoid too much emphasis on the line center of the strongest lines (i.e., Hα and Hβ). Specifically, we adopted the following equation

$$\begin{array}{c}\hfill F_{\mathrm{err}}=\sqrt{{({\mathrm{RMS}}_{\mathrm{run}}(\mathrm{res})\times \sqrt{N_{\mathrm{org}}/N_{\mathrm{dwn}}})}^2+{(0.05F_{\lambda })}^2},\end{array}$$(1)

where RMS_run(res) is the running RMS over the spectrum residuals with a window size of 30, N_org/N_dwn is the ratio of the number of channels in the original spectrum over the downgraded spectrum³, and F_λ is the flux density.

As a final step, we normalized the line profiles and their uncertainties according to the peak of our multi-Gaussian model of the broad lines in Sect. 2.2.2. The resulting profiles are shown in Fig. 2. The centroid and width of the line peaks are consistent, while the width of the line wings varies for different lines. We note that the red wing of the broad Hγ overlaps with the [O III] λ4363 line, which we fit with the narrow line template plus a Gaussian component (Sect. 2.2.2). We opted to keep the [O III] λ4363 model simple to avoid biasing the Hγ profile. However, this resulted in a relatively large residual (0−2000 km s⁻¹) and therefore uncertainty of the line profile, as shown in Fig. 2d. The Hγ line shows the strongest asymmetry due to the “shoulder” on the red wing, which is very difficult to model. We believe it can be at least partly explained in terms of the decomposition. In addition, the entire broad He I profile shows relatively large uncertainties because this line is very weak compared to the H I lines, and its uncertainties are dominated by the RMS term. Among the five broad lines studied in this work, He I is the most susceptible to any artifact of the spectrum.

Fig. 2. Normalized broad emission line profiles of Hα, Hβ, Hγ, Paβ, and He I. These were extracted as described in Sect. 2.2. Panel a shows the profiles superimposed for an easier comparison, while panels b–f display the line profiles individually. The uncertainties are shown in gray. The masked data have been replaced by the multi-Gaussian model, as shown in Fig. 1 for clarity.

3. Nonparametric properties of the broad-line profiles

Because broad-line profiles cannot be described by a simple analytical function, we characterized them nonparametrically and later used these quantities (Table 1) to assess the validity of our model. The peak velocity, v_peak, is the deviation of the peak wavelength of the line from its expected wavelength after shifting to the rest frame, as described in Sect. 2.1. All of the broad lines peak near the systemic velocity. We calculated the full width at 25% (W₂₅), 50% (FWHM; W₅₀), and 75% (W₇₅) maximum. We also calculated the asymmetry index (A.I.) and Kurtosis Index (K.I.) of each line profile as defined in Marziani et al. (1996),

$$\begin{array}{cc}\hfill \mathrm{A}.\mathrm{I}.& =\frac{v_R(1/2)+v_B(1/2)-2v_{\mathrm{peak}}}{v_R(1/2)-v_B(1/2)},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \mathrm{K}.\mathrm{I}.& =\frac{v_R(3/4)-v_B(3/4)}{v_R(1/4)-v_B(1/4)},\hfill \end{array}$$(3)

where v_R(x) > 0 and v_B(x) < 0 (x = 1/4, 1/2, or 3/4) are the velocity of line profiles at the corresponding fraction of the line peak on the red and blue wings, respectively. The values of A.I. indicate the direction and degree of asymmetry in the line profile shape. Differing slightly from the definition of Marziani et al. (1996), we calculated the A.I. at 50% (instead of 25%) of the peak flux because our line profiles become more symmetric toward the wings. A positive A.I. indicates that the line profiles are skewed toward the red side relative to the profile center. This suggests an excess of line emission or broader velocity distribution on the red side compared to the blue side. K.I., on the other hand, is essentially W₇₅/W₂₅. A Gaussian profile has K.I. ≈0.46. A smaller K.I. indicates that the line profile has a broader wing than a Gaussian profile, and vice versa. We also calculated the second moment of the line profiles, σ_line, following the definition of the Eq. (3) of Dalla Bontà et al. (2020) for completeness because σ_line is widely used in RM works to derive the BH mass (Peterson et al. 2004; Wang et al. 2019). To determine the uncertainty for each reported value, we randomly perturbed each profile using the uncertainties from Sect. 2.2.3, remeasured these quantities 500 times, and calculated the standard deviation of the results.

We compare these parameters in Fig. 3. The H I lines show almost the same widths at 75% and 50% of the line peaks. The Hα line shows the lowest W₂₅, while Hβ shows the highest W₂₅ among the H I lines. The line wing of Hβ is significantly broader than the other H I lines (see also Fig. 2). Although the Hβ region is complicated, the line wing cannot be biased by the spectral decomposition. As shown in Fig. 1, the broad Hβ line is much stronger than the (pseudo-)continuum and far enough from the [O III] lines. The He I line shows a comparable line width to the H I lines at its peak (W₇₅) but becomes much broader toward the wing (W₅₀ and W₂₅). The σ_line of all lines are similar and in between their W₅₀ and W₇₅, respectively. We do not plot the σ_line in Fig. 3 for clarity. All the lines show measurable and positive A.I., and although the Balmer lines show differing values, there are large uncertainties. The K.I. values are similar among the lines and are smaller than that expected for a single Gaussian profile, indicating that the profiles are peaky with broad wings.

Fig. 3. Comparison of the nonparametric properties of the line profiles (denoted by different colors and with 1σ uncertainties). Panel a: line widths (circles, triangles, and squares represent W25, W50, and W75, respectively). Panel b: asymmetry index (A.I.). Panel c: kurtosis index (K.I.). In the panels, the gray lines correspond to equivalent measurements on the model line profiles for the separate (dotted) and combined (solid) fitting results.

4. Modeling the broad-line profiles

In this section, we first introduce our BLR dynamical model, its limitations, and our inference strategy (Sect. 4.1). We then model the broad-line profiles in two steps: (1) We fit the line profiles separately and study the consistency of the model parameters (Sect. 4.2). (2) We fit the line profiles with almost all the BLR parameters tied and only allow the radial distribution of different line emissions to vary freely (Sect. 4.3).

4.1. Broad-line region dynamical model and the inference

The nature of the BLR is still an open question, and many models have been proposed to explain its various aspects (Peterson 2006; Czerny 2019, and references therein). One major class of models assumes that the BLR consists of many discrete clouds (e.g., Rees et al. 1989; Baldwin et al. 1995; Czerny & Hryniewicz 2011; Baskin & Laor 2018; Rosborough et al. 2023). The cloud model offers advantages in flexibly parameterizing the geometry, kinematics, and photoionization physics, enabling interpretation of observations, particularly recent high-quality RM and interferometric data (Korista & Goad 2004; Pancoast et al. 2011, 2014a; Li et al. 2013, 2018; Williams & Treu 2022; GRAVITY Collaboration 2018, 2020, 2021a). However, the physics of the cloud model may be oversimplified. For example, there is ongoing debate regarding the confinement of gas within high-density (≳10⁹ cm⁻³) clouds (Mathews 1986; Rees 1987; Krolik 1988; Baskin et al. 2014; Proga et al. 2014; Proga & Waters 2015). Moreover, as discussed in detail by Netzer (2020), state-of-the-art cloud photoionization models tend to underproduce the luminosity of the Balmer lines (and other non-resonant H I lines) by a factor of two to five, likely due to the failure of escape probability formalism in photoionization codes such as CLOUDY (Ferland et al. 1998) for high densities and optical depths in the BLR. Radiation hydrodynamic simulations of the disk wind, coupled with radiative transfer calculations, have been deemed crucial to understanding the photoionization physics of the BLR (Waters et al. 2016; Matthews et al. 2016, 2020; Mangham et al. 2017). However, the high computational expense impedes the development of more comprehensive models for detailed data interpretation.

In this work, we employed our self-implemented BLR dynamical model to characterize the distribution and kinematics of the BLR line emission (or “emissivity”). The model parameterization was initially developed in Pancoast et al. (2014a). Our model has been utilized to fit the normalized line profile and differential phase signal, tracking the spatially resolved kinematics of recent interferometric observations of BLRs (GRAVITY Collaboration 2020, 2021a,b). Where our implementation differs from Pancoast et al. (2014a) is in using a Monte Carlo cloud model to depict the line emission at the moment of the observation, excluding variable continuum light curves and reverberation mapping physics. Thus, our modeling approach circumvents the challenges associated with the aforementioned photoionization physics of the BLR by modeling the line emission distribution instead of the physical clouds and photoionization physics. In making this statement, we emphasize that our use of the term “cloud” should be taken to mean “line emitting entity”. It does not refer to physical “gas clouds”, nor does it indicate a preference for the cloud model over the so-called disk-wind model (e.g., Chiang & Murray 1996; Matthews et al. 2016; Long et al. 2023). Here, we adapted the BLR model into a Python package, DyBEL, allowing the fitting of line profiles exclusively. The DyBEL package can be used to fit either a single line profile or multiple lines of an AGN simultaneously. As the detailed BLR model is presented in GRAVITY Collaboration (2020), we provide a brief introduction to the key parameters, which we summarize in Table 2.

The model comprises a large number⁴ of non-interacting mass-less point particles orbiting the central BH with mass M_BH and forming a disk-like structure. The radial distribution of the clouds follows a shifted Gamma function governed by three parameters: μ as the mean radius, F = R_min/μ with R_min as the minimum cloud radius, and β for the shape of the profile (Gaussian: 0 < β < 1; exponential: β = 1; and heavy-tailed: 1 < β < 2). The angular thickness of the disk is θ_o, and the vertical distribution of the clouds is governed by γ, with a higher value (γ > 1) corresponding to more clouds concentrating on the disk surface. The structure is viewed at an inclination angle i (with 0° corresponding to a face-on view). Each cloud is randomly assigned to be on a quasi-circular orbit with radial and tangential velocities (v_r, v_ϕ) around (0, $v_{\mathrm{circ}}=\sqrt{G{M}_{\mathrm{BH}}/r}$), or a quasi-radial orbit. The fraction on quasi-radial orbits is controlled by f_ellip (with f_ellip = 1 meaning all clouds are on such orbits). A binary parameter, f_flow, governs the direction of the cloud radial motion. These clouds are inflowing if f_flow < 0.5 or outflowing if f_flow > 0.5. Their radial and tangential velocities are controlled by an angular parameter θ_e such that when θ_e = 0, the velocity vector is ($v_{\mathrm{esc}}=\sqrt{2}{v}_{\mathrm{circ}}$, 0), and when θ_e = π/2, it is (0, v_circ). While Pancoast et al. (2014a) had additional parameters defining how the cloud velocities are dispersed around these points, we excluded them because they are generally unconstrained and do not influence our fitting (GRAVITY Collaboration 2020, 2021a,b). Finally, the weight of each cloud is controlled by −0.5 < κ < 0.5, reflecting the anisotropy of the cloud illumination. Clouds closer to the observer have higher weights if κ > 0, and vice versa. The ratio of clouds below and above the midplane is controlled by ξ, reflecting the “midplane obscuration” of the BLR. There are equal amounts of clouds between the midplane if ξ = 1, while there is no cloud below the midplane if ξ = 0. To fit the line profiles, we needed two nuisance parameters: the central wavelength (λ_c) and the peak flux (f_peak) of the line.

We used the above BLR model to describe the line emission distribution of the BLR without including any photoionization physics. We cannot predict the physical line luminosity with this model, so the line peak flux is a free parameter in the fitting, and we only modeled the normalized line profiles. This approach was adopted in the recent works GRAVITY Collaboration (2018, 2020, 2021a, 2024) when they modeled the line profile and differential phase data (equivalent to spatially resolved kinematic data). In contrast, the original application of this model to the RM data assumes the point particles to be “mirrors” reflecting the continuum emission, the limitations of which are nicely summarized by Raimundo et al. (2020, in their Sect. 2.2). Their arguments make it clear that photoionization physics are needed in the application of RM modeling. Indeed, there has been recent progress in addressing this problem (Williams & Treu 2022; Rosborough et al. 2023). Nevertheless, recent studies of NGC 3783 have shown that the modeling using GRAVITY and RM data is remarkably consistent (Bentz et al. 2021b; GRAVITY Collaboration 2021a,b), which supports the application of this simple model, at least for the particular case of NGC 3783. We discuss the caveats of applying this model to the single-epoch spectra, which is the main goal of this work, in Sect. 5.3.

We used the Python package dynesty (Speagle 2020) to fit the data with a nested sampling algorithm, which is more powerful than the typical Markov chain Monte Carlo method for complex models with many (e.g., more than 20) parameters and a potentially multimodal posterior distribution. By design, the nested sampling algorithm (Skilling 2004) can estimate the Bayes evidence, which enabled us to compare different models. We used the dynamic nested sampler (DynamicNestedSampler), which better estimates the likelihood function by re-sampling the posterior function a few more times after the “baseline run”. We used 1200 live points for the baseline run and added 500 points for each of the 10 re-samplings. We adopted the random walk algorithm (rwalk) to sample the prior space and used the multi-ellipse method (multi) to create the nest boundaries. We adopted the default values for all the remaining options of dynesty.

The metric we used to define the goodness of fit is the likelihood function,

$$\begin{array}{c}\hfill ln\mathcal{L}=-\frac{1}{2}{\displaystyle \sum _l^N\sum _i^n\left(\frac{{(f_{l,i}-{\stackrel{\sim }{f}}_{l,i}({\mathbf{\Theta }}_l))}^2}{{\sigma }_{l,i}^{2}T}\right),}\end{array}$$(4)

where the first summation over l is for different lines and the second summation over i is for different channels of a line profile; f_l, i and σ_l, i are the line profile flux and uncertainty; ${\stackrel{\sim }{f}}_{l,i}({\mathbf{\Theta }}_l)$ is the corresponding line model with the set of parameters Θ_l; and a temperature parameter T > 1 is included to effectively scale up the uncertainties of the data. The temperature makes the likelihood function less peaky, which facilitates proper estimation of the posterior distribution. We found that T = 16 is a suitable setting, and our fitting results are not sensitive to the specific choice of temperature (e.g., T = 8, 16, or 32).

We find it is important to bear in mind that when fitting only the line profiles, M_BH and μ are fully degenerate. This is because the cloud velocities always scale with v_circ, which depends on M_BH/r, and it means that one needs to fix the BH mass in order to investigate the BLR sizes. Therefore, we included physical prior information by fixing M_BH = 10^7.4 M_⊙ (GRAVITY Collaboration 2021b). The exact value of the M_BH does not affect the derived BLR model parameters except for the BLR radius. We discuss this point in more detail in the following sections. Next, we set f_flow, a binary flag, to decide the direction of the radial velocity of the clouds. While the previous modeling efforts all indicate radial inflow (Bentz et al. 2021b; GRAVITY Collaboration 2021a,b), by using line profiles alone, one cannot distinguish between inflow and outflow (either in the model or via the Bayes evidence) because there is no spatial information. We note that the specific choices of f_flow do not affect our results. Therefore, we adopted an inflow model setting f_flow < 0.5.

When fitting the line profiles simultaneously with DyBEL, we could choose to tie other parameters in addition to fixing the same BH mass. As our aim is to assess whether the difference in the line profiles can be attributed to radial differences in the line emission distribution for an otherwise fixed BLR geometry, we therefore left μ, F, and β free to vary for each line, while the remaining model parameters were tied. We also tied the central wavelength offsets for each line, allowing them to shift together by a small amount ϵ = λ_c/λ_air − 1, where λ_air is the nominal wavelength in air. We adopted a Gaussian prior centered at zero with a small standard deviation of 0.01 for ϵ because λ_c is expected to be close to λ_air. Similarly, we expected the normalized line peaks to be close to one, so we only adopted a single nuisance parameter, f_peak, in the fitting with a Gaussian prior centered at one with a standard deviation of 0.1. The priors of the remaining parameters were adopted from GRAVITY Collaboration (2020).

For comparison, we first fit each profile separately, then mainly discuss the results fitting all of the five lines simultaneously. In each case, we calculated several parameters derived from the fit, including the minimum radius, R_min ≡ μF; the weighted mean radius, R_mean, of the BLR clouds; and the virial factor, f_virial (Eq. (5)).

4.2. Fitting the line profiles separately

In this section, we investigate what can be learned from fitting the five line profiles separately. We show a comparison of the data and the fitting results in Fig. 4. We generated 500 model line profiles with parameters randomly selected from the posterior samples of each line and plotted the median line profiles in panel a and the 68% confidence interval of each line in panels b–f. Qualitatively, the BLR model can fit the line profiles reasonably well. The median model profiles reflect the differences of the lines. The 68% confidence profiles largely enclose the data of all lines. The Hα and Hβ model profiles are more tightly constrained than the other three lines thanks to their smaller uncertainties. The most obvious mismatch comes from the Hγ line, due to the “shoulder” on the red wing, as noted in Sect. 2.2.3. The He I model profile shows some deviation from the data too, although always within the 1-σ uncertainty level.

Fig. 4. Median and 68% confidence intervals of the model line profiles when the lines are fitted separately. The model profiles are overplotted in panel a, and the 68% confidence intervals of (b) Hα, (c) Hβ, (d) Hγ, (e) Paβ, and (f) He I are compared with the data (black). The masked data are not plotted, so gaps are visible in some profiles.

The nonparametric parameters (line widths, A.I., and K.I.) for the model profiles are shown as vertical dotted lines in Fig. 3. The model line widths follow the corresponding lines remarkably well in all cases. The A.I. values show clear differences to the data, although largely within the uncertainties. In particular, the A.I. of Hγ is much higher than the model value. The K.I. values are consistent with the data, though the model posteriors tend to be higher.

There are 12 free parameters for each profile. The posterior distributions are displayed in Fig. 5. Following Raimundo et al. (2020), we considered a model parameter to be constrained if its 68% confidence range is less than half of its prior range. Consistent with these authors, we found the geometric parameters, i, θ_o, μ, and β, of most of the lines can be constrained, while the remaining parameters, f_ellip, κ, γ, ξ, and θ_e, mostly could not be constrained. The i and θ_o of Hα and Hβ agree well with the inclination derived by the SARM joint analysis (GRAVITY Collaboration 2021b), although these two parameters were less well constrained for Hγ, Paβ, and He I. In particular, the large asymmetry of Hγ led to a higher probability of a high inclination angle, so we observed tentative double-peaked posterior distributions for i and θ_o. In contrast, the distributions for Paβ and He I are broad and smooth, likely because the uncertainties of these two profiles are relatively large. The fits yielded generally small values for β (except Paβ), compared to the SARM result (β ≈ 1.95). This differs from the high value found in the SARM joint analysis but is consistent with the RM result for Hβ reported by Bentz et al. (2021a). Similarly, the minimum cloud radii (R_min) for all lines except Paβ prefer smaller values than the SARM result. Nevertheless, the posteriors of the BLR sizes (R_mean) are largely consistent with the SARM results.

Fig. 5. Posterior probability distributions of the five lines fitted separately shown in colored histograms. For comparison, the dashed vertical line and the shaded region indicate the model inference best-fit results and 1-σ intervals from the SARM joint analysis (GRAVITY Collaboration 2021b). Panels a–i present the physical model parameters that are directly sampled in the fitting. The circles and crosses in these panels indicate whether these parameters are constrained by the data for each line profile (by comparing the width of their posterior distributions with their prior range). The remaining panels present parameters that are derived from the posterior samples. The y-axis tick labels are unimportant, so we removed them for clarity. Details are discussed in Sect. 4.2.

4.3. Fitting the line profiles simultaneously

We performed a fit in which we tied all of the parameters of the BLR model for the five lines except those defining the radial distribution of the clouds, namely, μ, β, and F. With this approach, we could test whether the difference in the line profiles can be solely explained by the radial stratification of the BLR. There are 24 free parameters in the fit: nine tied between all the line profiles and three left separate for each of the five lines. We report the median and 68% confidence interval values of the combined fit posterior samples in Table 3. The model profiles and the 68% confidence intervals are shown in Fig. 6. Panel a of the figure shows that the model profiles of Hα, Hγ, and Paβ are very similar, while those of Hβ and He I are wider. This is similar to the results of the separate fitting. Because most of the model parameters were tied in this approach, we could conclude that the line profile differences can be explained by the radial distributions of the line emission. Moreover, the combined fit results showed tighter 68% confidence intervals than the separate fit. The improvement was most obvious in the Hγ, Paβ, and He I lines, whose uncertainties are relatively large. As a result, the deviation between the model and data of Hγ is more obvious.

Fig. 6. Median and 68% confidence intervals of the model line profiles when the lines are fitted simultaneously with most of the parameters tied. The panels and symbols are the same as those of Fig. 4.

The nonparametric values for these model profiles are similar to those from the separate fits, while the uncertainties of the simultaneous fits are smaller than those of the separate fits. Again, they reproduce the line widths of the data well. The A.I. values are similar between the lines because the model parameters controlling asymmetry in the model (κ and ξ) were tied. We verified that if these are left free, the model yields different values of ξ to better fit the individual asymmetries of the lines. Meanwhile, the other parameters do not change substantially, and in neither case is the high asymmetry of Hγ reached. The K.I. values of the tied model profiles are more consistent with the data than those of the separate fits. We note that both model profiles yield σ_line consistent with the data profiles of all lines when fitting the data both separately and simultaneously.

Figure 7 shows that almost all of the tied model parameters are properly constrained and consistent with the SARM joint analysis. In particular, the distributions of β are similar to those from the separate fits except that Paβ now too prefers a lower value and thus matches the other lines. The lines also have similar R_min ∼ 4 ld; hence, the line emission is concentrated in an inner ring and extends to large radii. In passing, we tested to further tie R_min in the fit and found the results stay almost entirely the same with the tied R_min ≈ 3.5 ld, confirming that the simultaneous fit favors the different lines sharing the same R_min. GRAVITY Collaboration (2021b) tested fixing the R_min = 4 ld in the SARM analysis and found a reasonable fitting result with β ≈ 1, which, interestingly, is close to our β. However, GRAVITY Collaboration (2021b) cautioned that additional restrictions bias their geometric distance measurement by about 30%.

Fig. 7. Posterior probability distribution of the five lines fitted simultaneously with tied geometry parameters. The panels and symbols are similar to those of Fig. 5. For comparison, as before, the dashed vertical line and the shaded region indicate the model inference results from the joint SARM fit (GRAVITY Collaboration 2021b).

The parameter that differs the most from the SARM result is γ, where a large value is preferred, indicating that the vertical distribution of line emission is more concentrated toward the surface of the BLR. However, there is a long tail toward low values in the posterior distribution of γ. The number of clouds in the mid-plane increases by less than 20% when γ decreases to 1.7 (the value derived in the SARM joint analysis), so the model does not change substantially. Nevertheless, the edge-concentrated structure with high γ may indicate that a biconical wind-like structure (e.g., Matthews et al. 2016; Waters et al. 2021) could fit the line profiles. The flexibility enabled by including such a capability in the model may allow it to better reproduce the asymmetric features in broad lines of NGC 3783.

5. Broad-line region geometry and virial factor from single-epoch line profiles

5.1. Broad-line region model parameters

In this section, we discuss the model parameters of the BLR that can be measured from single-epoch line profiles. For individual line profiles, we can constrain the inclination (i) and disk thickness (θ_o), a finding consistent with Raimundo et al. (2019, 2020). Because we fixed the BH mass, the radial distribution of the line emission can be constrained too. More interestingly, we find the BLR model can simultaneously fit multiple line profiles using the same geometry and kinematics. Almost all of the model parameters can be constrained, and only the radial distributions change for each line, indicating that the radial emissivity distributions are somewhat different, as expected (see the next paragraph). The H I lines tend to favor heavy-tailed distributions with β > 1, while the He I line is close to exponential (β ≈ 1). While it is beyond the scope of this work, we note that β ≳ 1 is close to a truncated power-law distribution that the photoionization model may produce (e.g., Netzer 2020).

Next, we compare the BLR radius (R_mean) of the combined fit of the five lines (Fig. 8a) with published measurements based on the Hβ and Brγ lines (Bentz et al. 2021a,b; GRAVITY Collaboration 2021a,b). We focused on the joint SARM analysis because it obtains the tightest constraints of the model parameters using the spectro-astrometry and RM data simultaneously and because we adopted its BH mass. The R_mean from the separate fitting shows relatively large uncertainties. The R_mean of Hβ, which happens to have the smallest uncertainty, is consistent with the R_mean derived by the SARM analysis. In contrast, R_mean derived from the combined fitting shows much smaller uncertainties, with the HβR_mean consistent with the SARM result as well. The He IR_mean is the smallest among the five lines. We note that the posterior distributions of R_mean are strongly correlated between different lines (Fig. A.1). While such a correlation increases the uncertainties of the mean radii (≳30%, Fig. 8a), as discussed in the next paragraph, it also means that the ratios of the mean radii have smaller uncertainties (≲10%).

Fig. 8. Comparison of derived parameters. Values of (a) mean radius and (b) virial factor are shown for the separate (dotted lines) and combined (solid lines) fits. The three horizontal bars on each line correspond to the 16th, 50th, and 84th percentiles of the posterior distribution. The dashed lines and gray regions represent the results from the SARM joint analysis (GRAVITY Collaboration 2021b).

As shown in Table 4, we measured the BLR mean radius ratios of Hα:Hβ:Hγ:Paβ:He I from the simultaneous fitting to be 1.47:1.0:1.22:1.36:0.72. They are largely consistent with the RM observation results reported by Bentz et al. (2010b). Our relative lags of the Hγ and He I lines are larger than the RM results. We also calculated the time lag ratios based on the radiation pressure confined (RPC) BLR model described by Netzer (2020). The theoretical model calculation explores a range of parameters that are in agreement with most RM measurements of various hydrogen and helium lines (e.g., Bentz et al. 2009, 2010b). The most important parameters of the RPC model are the radial dependence of the covering factor of the clouds, R_min which is somewhat arbitrary, and the BLR outer radius, which is determined by graphite dust sublimation radius, gas metallicity, and turbulent velocity within individual clouds. The level of ionization and hence the line emissivity at all locations were obtained naturally from the assumption that the clouds’ column densities are large and that they are in total gas and radiation pressure equilibrium. The mean emissivity radius for each line was then computed from the model and then translated into time lags, which depend on the light curve of the driving continuum. All such models predict τ_Hα > τ_Hβ > τ_Hγ > τ_He I, and the range is illustrated in Table 4. Similar tendencies have also been predicted by the very different LOC model computed by Korista & Goad (2004).

We suspect our large Hγ radius is due to the error of the Hγ line profile. The Hγ line, the weakest among the three Balmer lines, is blended with [O III] λ4363 line (Fig. 1a). Therefore, the small but systematic residual of the decomposition may lead to a too large Hγ BLR size compared to its actual dimensions. The clear systematic deviations between the data and model in Fig. 6d support this point to some extent. This issue highlights the importance of high-quality line profiles in revealing robust BLR properties. Another effect that may influence the observed line emission distribution is the polar dust around the BLR (Hönig et al. 2013). Higher extinction in the center makes the BLR look more extended. Since extinction is larger at shorter wavelengths (Li 2007), this effect might enlarge the observed Hγ size more than Hα and Hβ. It is worth mentioning that Bentz et al. (2021a) reported a very small time lag of 3.7 ld for Hγ in NGC 3783, as small as that of He II. As discussed by the authors, the Hγ time lag was likely underestimated due to imperfect internal flux calibration using the [O III] λ5007 line flux. Observations of more targets would be useful to disentangle these effects on the measured BLR size across different lines.

We remind the readers that the models shown here do not include time-dependent variations of the ionizing source, and we caution that the flux- weighted radius (R_mean) may be different from the time lag, depending on the structure of the BLR and the ionization continuum light curve (Netzer & Maoz 1990). However, we expect the R_mean ratios from our simultaneous fitting to be close to the ratios of the time lags because the different lines are assumed to share the same BLR structure. To conclude, our analysis shows that different broad emission lines appear to share most of the geometry and kinematics of the BLR. The relative BLR sizes are largely consistent with the theoretical expectation of photoionization models. This suggests that one can combine spectro-astrometry and RM data of different lines in a joint analysis by tying most of the BLR model parameters of the two lines.

5.2. The virial factor

The virial factor encapsulates the geometry of the BLR by linking the BH mass to the measurable properties of BLR size and line width. In this paper, we define it as

$$\begin{array}{c}\hfill f_{\mathrm{virial}}=\frac{G{M}_{\mathrm{BH}}}{{\sigma }_{\mathrm{line}}^2{R}_{\mathrm{mean}}},\end{array}$$(5)

where σ_line is the second moment of the model line profile, R_mean is the mean radius of the BLR, and M_BH is the fixed BH mass. A key property is that f_virial scales with M_BH/r, where r indicates the BLR size. Although we fixed the BH mass of NGC 3783, the value of f_virial that we derived does not depend on M_BH because the BLR radius is a free parameter in the fit. The reason, as explained in Sect. 4.1, is that M_BH and μ (or R_mean) are fully degenerate in our model. These parameters scale together without changing the line profile. Therefore, we can expect to derive meaningful virial factors from the fit. To confirm this point, we fit the data with a fixed M_BH = 10^6.4 M_⊙ (10 times smaller) and got the same f_virial results.

As shown in Fig. 8b, the derived values from both the separate fit and the combined fit are consistent with that of the SARM joint analysis ($f_{\mathrm{virial}}=2.{52}_{-0.53}^{+0.62}$). The combined fit shows smaller uncertainties than the separate fit. This is an encouraging result for investigating the BLR dynamics and measuring the BH mass. Our method could enable one to constrain the individual virial factor for each AGN by modeling the broad emission line(s) without the need for dynamical modeling of RM data (e.g., Villafaña et al. 2023). As a practical approach, one can fix the BH mass according to the single-epoch estimate (based on an averaged virial factor; McLure & Dunlop 2001; Ho & Kim 2015; Dalla Bontà et al. 2020) and model the broad-line profile(s) to derive the f_virial for individual AGNs. The f_virial can then be used to refine the BH mass estimate. As such, it could reduce the uncertainty in the BH mass that is otherwise introduced by adopting an average virial factor (Collin et al. 2006; Shen & Ho 2014). We caution that the derived virial factor may be biased by the oversimplified BLR model, the influence of which will be investigated with many more sources in the future.

5.3. Caveats

In this work, we have investigated the BLR structure by modeling the normalized single-epoch line profiles. We modeled the broad-line emission and the associated kinematics with a Monte Carlo model of points without a physical size. This model is intended to avoid consideration of the details of the photoionization physics and cannot predict the line strength physically according to the AGN luminosity. Without the data spatially resolving the BLR structure (e.g., GRAVITY differential phase), the model parameters may be degenerate when only fitted with the line profiles. Interestingly, for NGC 3783, we find the single-epoch line profiles, especially when fitted simultaneously, can provide most of the BLR model parameters in a manner consistent with the fitting that includes the size measurements. We caution, however, that more studies on different BLRs are needed to understand whether the conclusion holds widely.

To test whether our BLR model is quantitatively plausible with photoionization physics, we estimated the Hβ luminosity assuming the Case B recombination (Osterbrock & Ferland 2006) and the geometric covering factor based on our model fitting result. With θ_o ≈ 28.8° (Table 3), the geometric covering factor is sin θ_o ≈ 0.5. We adopted the ratio of the Hβ line and hydrogen recombination coefficients of ${\alpha }_{\mathrm{H}\beta }^{\mathrm{eff}}/{\alpha }_B\approx 1/8.5$ and the UV photon flux ${\int }_{{\nu }_0}^{\infty }\frac{L_{\nu }}{h\nu }\mathrm{d}\nu \approx 1.6-3.7\times {10}^{53}$ s⁻¹ (ν₀ = 13.6 eV), which is estimated with the measured λL_λ(5100 Å) ≈ 4.1 × 10⁴² erg s⁻¹ and the assumed AGN spectral energy distribution (SED) with low and intermediate Eddington ratios (Ferland et al. 2020; Jin et al. 2012) in order to enclose the range of typical Seyfert galaxy SEDs. We derive L_Hβ ≈ 0.4 − 0.9 × 10⁴¹ erg s⁻¹. The measured Hβ luminosity⁵, ∼1.1 × 10⁴¹ erg s⁻¹, is comparable to the estimated range, if not slightly higher, indicating that our BLR model is plausible in terms of line luminosity. While this estimate is admittedly oversimplified, it illustrates that the shape of the SED may easily influence the line luminosity by a factor of a few. It is worth noting that we assumed the maximum absorption of the UV photons with the model BLR geometry, and the Case B recombination does not consider the self-absorption of the Hβ photons. More detailed photoionization calculations with CLOUDY may only provide weaker line emission, which reflects the aforementioned problem in Sect. 4.1. Although this problem is beyond the scope of this work, our method provides a new approach to address it with single-epoch spectra.

6. Conclusions

We investigated the BLR structure of NGC 3783 using multiple broad lines in a high-resolution single-epoch spectrum obtained with VLT/X-shooter. We decomposed the five strongest broad lines (Hα, Hβ, Hγ, Paβ, and He Iλ5876), and modeled their profiles using the newly developed tool DyBEL, which allows one to tie parameters of the dynamical model between the lines. Since the BH mass and the BLR radius are fully degenerate, we opted to fix the BH mass to a value reported in the literature and focused on the BLR structure and emissivity that can be derived from the line profiles. Our main results are the following:

1. All lines analyzed here show broader wings than a Gaussian profile and are asymmetric (skewed to the red side). The He Iλ5876 profile is broader than the hydrogen profiles studied here.

2. We developed a fitting tool to model the line profiles with a dynamical BLR mode. Fitting multiple lines simultaneously by tying many of their parameters together yields a solution that is better constrained than when fitting them individually. In particular, it yields useful constraints on parameters such as inclination, BLR size, and the virial factor.

3. The difference in line profiles can be explained almost entirely in terms of differing radial distributions of the line emission. The derived relative BLR time lags are mostly consistent with the RM observation and with theoretical model calculations. Our results support that it is possible to combine spectro-astrometry and RM data in a joint analysis.

4. The virial factor we derived is nearly the same for the five lines and is independent of the adopted BH mass. We argue that by enabling one to constrain the virial factor for an individual AGN using a single-epoch spectrum, this method can reduce the uncertainty in BH masses derived from single-epoch spectra.

In the future, more comprehensive analyses will be useful to explore the potential of this method with more broad-line AGNs, different spectral resolutions and signal-to-noise ratios, and different BLR dynamical models.

Acknowledgments

L.K. thanks the financial support from the Mitacs Globalink Research Award. A.W.S.M. and L.K. acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) through grant reference number RGPIN-2021-03046. J.S. dedicates this paper in memory of his grandmother, G.X. (1930–2023). J.S. thanks Leonard Burtscher for sharing the X-shooter data.

References

Appendix A: Simultaneous fitting results with all parameters tied

The corner plot of the simultaneous fitting is shown in Figure A.1. The separate fittings show similar but less constrained results. We opted not to show all of the corner plots of the separate fitting for simplicity because all of the useful information has been shown in Figure 5.

Fig. A.1. Corner plot of the simultaneous fitting. The first nine parameters are tied in the fitting, while the remaining parameters are fitted for individual lines.

All Tables

All Figures

Fig. 1. X-shooter spectrum of NGC 3783 simultaneously covering UV, optical, and NIR ranges. In each panel, the spectrum is shown with the full model that enables decomposition of the broad-line profiles overplotted, and the residual is shown underneath. The emission line components are also plotted separately. The gray-shaded regions in the residuals represent wavelength ranges of bad channels and features due to the poor telluric correction, which are masked in the fitting. Panel a: UV/optical spectrum with the Hα, Hβ, and Hγ lines. Panel b: part of the NIR spectrum used in this study, with the Paβ line. Panel c: [S II] doublet used as the narrow line template.

In the text

	Fig. 2. Normalized broad emission line profiles of Hα, Hβ, Hγ, Paβ, and He I. These were extracted as described in Sect. 2.2. Panel a shows the profiles superimposed for an easier comparison, while panels b–f display the line profiles individually. The uncertainties are shown in gray. The masked data have been replaced by the multi-Gaussian model, as shown in Fig. 1 for clarity.
In the text

Fig. 3. Comparison of the nonparametric properties of the line profiles (denoted by different colors and with 1σ uncertainties). Panel a: line widths (circles, triangles, and squares represent W25, W50, and W75, respectively). Panel b: asymmetry index (A.I.). Panel c: kurtosis index (K.I.). In the panels, the gray lines correspond to equivalent measurements on the model line profiles for the separate (dotted) and combined (solid) fitting results.

In the text

	Fig. 4. Median and 68% confidence intervals of the model line profiles when the lines are fitted separately. The model profiles are overplotted in panel a, and the 68% confidence intervals of (b) Hα, (c) Hβ, (d) Hγ, (e) Paβ, and (f) He I are compared with the data (black). The masked data are not plotted, so gaps are visible in some profiles.
In the text

Fig. 5. Posterior probability distributions of the five lines fitted separately shown in colored histograms. For comparison, the dashed vertical line and the shaded region indicate the model inference best-fit results and 1-σ intervals from the SARM joint analysis (GRAVITY Collaboration 2021b). Panels a–i present the physical model parameters that are directly sampled in the fitting. The circles and crosses in these panels indicate whether these parameters are constrained by the data for each line profile (by comparing the width of their posterior distributions with their prior range). The remaining panels present parameters that are derived from the posterior samples. The y-axis tick labels are unimportant, so we removed them for clarity. Details are discussed in Sect. 4.2.

In the text

	Fig. 6. Median and 68% confidence intervals of the model line profiles when the lines are fitted simultaneously with most of the parameters tied. The panels and symbols are the same as those of Fig. 4.
In the text

	Fig. 7. Posterior probability distribution of the five lines fitted simultaneously with tied geometry parameters. The panels and symbols are similar to those of Fig. 5. For comparison, as before, the dashed vertical line and the shaded region indicate the model inference results from the joint SARM fit (GRAVITY Collaboration 2021b).
In the text

Fig. 8. Comparison of derived parameters. Values of (a) mean radius and (b) virial factor are shown for the separate (dotted lines) and combined (solid lines) fits. The three horizontal bars on each line correspond to the 16th, 50th, and 84th percentiles of the posterior distribution. The dashed lines and gray regions represent the results from the SARM joint analysis (GRAVITY Collaboration 2021b).

In the text

	Fig. A.1. Corner plot of the simultaneous fitting. The first nine parameters are tied in the fitting, while the remaining parameters are fitted for individual lines.
In the text