Broad-line region geometry from multiple emission lines in a single-epoch spectrum

The broad-line region (BLR) of active galactic nuclei (AGNs) traces gas close to the central supermassive black hole (BH). Recent reverberation mapping (RM) and interferometric spectro-astrometry data have enabled detailed investigations of the BLR structure and dynamics, as well as estimates of the BH mass. These exciting developments motivate comparative investigations of BLR structures using different broad emission lines. In this work, we have developed a method to simultaneously model multiple broad lines of the BLR from a single-epoch spectrum. We apply this method to the five strongest broad emission lines (H$\alpha$, H$\beta$, H$\gamma$, Pa$\beta$, and He $I\;\lambda$5876) in the UV-to-NIR spectrum of NGC 3783, a nearby Type I AGN which has been well studied by RM and interferometric observations. Fixing the BH mass to the published value, we fit these line profiles simultaneously to constrain the BLR structure. We find that the differences between line profiles can be explained almost entirely as being due to different radial distributions of the line emission. We find that using multiple lines in this way also enables one to measure some important physical parameters, such as the inclination angle and virial factor of the BLR. The ratios of the derived BLR time lags are consistent with the expectation of theoretical model calculations and RM measurements.


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 typical full width at half maximum (FWHM) ≳ 1000 km s −1 (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 virial Rv 2 /G, where R is the BLR radius; v is a characteristic velocity of the BLR rotation; and f virial is the virial factor which 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, the reverberation mapping (RM) technique (Blandford & McKee 1982;Peterson 2014).The characteristics of these time lags across different velocity channels provide 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;Li et al. 2013;Pancoast et al. 2014a,b).More recently, the BLR has been spatially resolved with spectro-astrometry (SA), which is a powerful technique for measuring the BLR structure and BH mass (GRAVITY Collaboration et al. 2018, 2020, 2021b, 2023).Attempts have been made to analyze SA and RM data jointly (hereafter, the SARM method), to measure the geometric distance of the BLR and better constrain the BLR structure and BH mass (Wang et al. 2020;GRAVITY Collaboration et al. 2021a;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).But there is a wealth of AGN samples with good quality single-epoch spectra.To exploit these, Raimundo et al. (2019Raimundo et al. ( , 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 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 one UV/optical/NIR spectrum.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 is 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) report 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 et al. (2021b) report the SA measurement of the broad Brγ line.This motivated a joint SARM analysis, which yielded consistent results (GRAVITY Collaboration et al. 2021a).Here, we introduce the data and our method to decompose the broad-line profiles in Section 2. Then we make a nonparametric characterization of them in Section 3. We discuss our modeling of the line profiles in Section 4, comparing our results with GRAVITY Collaboration et al. (2021a) because the joint analysis provides the strongest constraint of the BLR model.We discuss the strengths and limitations of the current method in Section 5 and conclude in Section 6.

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 ESO reflex pipeline (version 2.6.8) with 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 is calibrated to an accuracy of about 2%.We extracted 1D spectra from each of the NGC 3783 datacubes using a rectangular slit with a width of 1.8 ′′ , and applied minor scaling corrections to match the different spectral ranges.
We corrected Galactic extinction based on A V = 0.332 (Schlafly & Finkbeiner 2011) and the Cardelli et al. (1989) extinction model where we specify R V = 3.1 as the ratio of total to selective extinction, and convert the spectrum to the rest frame adopting a redshift of 0.00973 measured from 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 modelling the BLR profiles (Section 4).We also identify 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 Figure 1, together with the spectral decomposition that will be discussed in Section 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 Figure 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 noise less narrow line template (Ho et al. 1997).To generate the template, we fit each of the [S II] lines with 3 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 3rd-order polynomial function.We found more Gaussian components or a higher-order polynomial function cannot improve the fitting.We have 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.

Decomposing the broad-line profiles
The continuum of an AGN optical spectrum consists of emission from the accretion disk, the host galaxy, as well as the pseudocontinuum 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 re- Rest Wavelength (Å)  sults.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 simultaneously (Figure 1a).Fitting the more isolated the Paβ line is described later.
In the optical spectrum, the majority of narrow lines are fitted by the [S II] template with two free parameters, the amplitude and the velocity shift from the laboratory wavelength.We tie the velocity shifts of all the templates, so it reflects a small deviation (≈ −2 km s −1 ) from the redshift measured by the atomic H I gas.Because the UVB arm has a lower spectral resolution than the VIS arm, we broaden the narrow line template with a Gaussian kernel (σ ≈ 35 km s −1 ) 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 add 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 are 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 nar-row lines 1 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 5th-order polynomial function to account for large-scale variations due to 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 −16 erg s −1 cm −2 Å −1 .
In the NIR spectrum, we fitted the broad Paβ profile and continuum over a wavelength range of 12400-13300 Å.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 is not tied to the optical narrow lines.But their velocity difference of ≲ 1.7 km s −1 is consistent with the systematic uncertainty of X-shooter wavelength calibration over different arms. 2 This high accuracy enables us to tie the central wavelength offsets when we fit the broad line profiles simultaneously (Section 4).

Resampling and uncertainties
While the high resolution of the X-shooter spectrum helps to 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 −1 ) 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 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 −1 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 orig-2 https://www.eso.org/sci/facilities/paranal/instruments/xshooter/doc/XS_wlc_shift_150615.pdfinal 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 adopt the following equation, where RMS run (res) is the running root-mean-squared (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,3 and F λ is the flux density.
As a final step, we normalize the line profiles and their uncertainties according to the peak of our multi-Gaussian model of the broad lines in Section 2.2.2.The resulting profiles are shown in Figure 2. The centroid and width of the line peaks are well 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 is fitted with the narrow line template plus a Gaussian component (Section 2.2.2).We opted to keep the [O III] λ4363 model simple to avoid biasing the Hγ profile.However, this results in a relatively large residual (0-2000 km s −1 ), and therefore uncertainty, of the line profile as shown in Figure 2(d).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.He I is the most susceptible to any artifact of the spectrum among the five broad lines studied in this work.

Nonparametric properties of the broad line profiles
Because broad-line profiles cannot be described by a simple analytical function, we characterize them nonparametrically and later use 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 Section 2.1.All of the broad lines peak near the systemic velocity.We calculate the Full Width at 25% (W 25 ), 50% (FWHM, W 50 ), and 75% (W 75 ) Maximum.We also calculate the Asymmetry Index (A.I.) and Kurtosis Index (K.I.) of each line profile as defined in Marziani et al. (1996), 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 calculate the A.I. at 50% (instead of 25%) of the peak flux because our line profiles become more symmetric towards the wings.A positive A.I. indicates that the line profiles are skewed towards 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 75 /W 25 .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 calculate the second moment of the line profiles, σ line , following the definition of the Equation (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 perturb each profile using the uncertainties from Section 2.2.3, re-measure these quantities 500 times, and calculate the standard deviation of the results.
We compare these parameters in Figure 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 25 , while Hβ shows the highest W 25 among the H I lines.The line wing of Hβ is significantly broader than the other H I lines (see also Figure 2).Although the Hβ region is complicated, the line wing cannot be biased by the spectral decomposition.As shown in Figure 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 75 ) but becomes much broader towards the wing (W 50 and W 25 ).The σ line of all lines are similar, and in between their W 50 and W 75 , respectively.We do not plot them in Figure 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.

Modeling the broad-line profiles
In this section, we first introduce our BLR dynamical model, its limitations, and our inference strategy (Section 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 (Section 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 (Section 4.3).

BLR 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 parameterizing the geometry, kinematics, and photoionization physics flexibly, enabling interpretation of observations, in particular recent high-quality RM and interferometric data (Korista & Goad 2004;Pancoast et al. 2011Pancoast et al. , 2014a;;Li et al. 2013Li et al. , 2018;;Williams & Treu 2022;GRAVITY Collaboration et al. 2018, 2020, 2021b).However, the physics of the cloud model may be oversimplified.For example, there is ongoing debate regarding the confinement of gas within highdensity (≳ 10 9 cm −3 ) 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), stateof-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 2-5, likely due to the failure of escape probability formalism in photoionization codes like 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 understand the photoionization physics of the BLR (Waters et al. 2016;Matthews et al. 2016Matthews et al. , 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 employ 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 et al. 2020Collaboration et al. , 2021a,b),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.DyBEL 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 et al. ( 2020), we provide a brief introduction to the key parameters, summarized in Table 2.
The model comprises a large number4 of non-interacting mass-less point particles orbiting the central BH with mass M BH , 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 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 circ = √ GM 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 esc = √ 2v 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 exclude them because they are generally unconstrained and do not influence our fitting (GRAVITY Collaboration et al. 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 need two nuisance parameters, the central wavelength (λ c ) and the peak flux ( f peak ) of the line.
We use 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 model the normalized line profiles.This approach is adopted by the recent works of GRAVITY Collaboration et al. (2018, 2020, 2021b, 2023) when they model 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 as 'mirrors' reflecting the continuum emission, the limitations of which are nicely summarized by Raimundo et al. (2020, in their Section 2.2).Their arguments make it clear that photoionization physics is needed in the application of RM modeling.Indeed, there is recent progress in addressing this problem (Williams & Treu 2022;Rosborough et al. 2023).Nevertheless, the recent study of NGC 3783 shows that the modeling using GRAVITY and RM data is remarkably consistent (Bentz et al. 2021b;GRAVITY Collaboration et al. 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 Section 5.3.We use 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.> 20) parameters and a potentially multimodal posterior distribution.By design, the nested sampling algorithm (Skilling 2004) can estimate the Bayes evidence, which enables us to compare different models.We use 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 use 1200 live points for the baseline run and add 500 points for each of the 10 re-samplings.We adopt the random walk algorithm (rwalk) to sample the prior space and use the multi-ellipse method (multi) to create the nest boundaries.We adopt the default values for all the remaining options of dynesty.
The metric we use to define the goodness of fit is the likelihood function, 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 and fl,i (Θ 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).
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 include physical prior information by fixing M BH = 10 7.4 M ⊙ (GRAVITY Collaboration et al. 2021a).The exact value of the M BH does not affect the derived BLR model parameters except the BLR radius.We will 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 et al. 2021a,b), using line profiles alone cannot distinguish between inflow and outflow (either in the model or via the Bayes evidence) because there is no spatial information, and the specific choices of f flow do not affect our results.Therefore, we adopt an inflow model setting f flow < 0.5.
When we fit the line profiles simultaneously with DyBEL, we can choose to tie additional parameters besides fixing the same BH mass.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 leave µ, F, and β free to vary for each line, while the remaining model parameters are tied.We also tie 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 adopt a Gaussian prior centered at 0 with a small standard deviation of 0.01 for ϵ because λ c is expected to be close to λ air .Similarly, the normalized line peaks are expected to be close to 1, so we only adopt a single nuisance parameter, f peak , in the fitting with a Gaussian prior centered at 1 with a standard deviation of 0.1.The priors of the remaining parameters are adopted from GRAVITY Collaboration et al. (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 calculate several parameters derived from the fit: the minimum radius, R min ≡ µF; the weighted mean radius, R mean , of the BLR clouds; and the virial factor, f virial (Equation 5).

Fitting the line profiles separately
In this section, we investigate what can be learned from fitting the five line profiles separately.We compare the data and the fitting results in Figure 4. We generate 500 model line profiles with parameters randomly selected from the posterior samples of each line, and plot 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.Note.The upper part of the table consists of the parameters tied for the five lines, and the lower part presents the model parameters (β, F, and µ) that are sampled freely for individual lines, and the corresponding derived quantities (R min , R mean , and f virial ).The median and 68% confidence interval values are reported.
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 tighter 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 Sec.2.2.3.The He I model profile shows some deviation from the data too, although always within the 1-σ uncertainty level.
The nonparametric parameters (line widths, A.I., and K.I.) for the model profiles are shown as vertical dotted lines in Figure 3.The model line widths follow the corresponding lines re-markably 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, although the model posteriors tend to be higher.
There are 12 free parameters for each profile.The posterior distributions are displayed in Figure 5.Following Raimundo et al. (2020), we consider a model parameter constrained if its 68 per cent confidence range is less than half of its prior range.Consistent with these authors, we find the geometric parameters, i, θ o , µ, and β, of most of the lines can be constrained, while the In particular, the large asymmetry of Hγ leads to a higher probability of a high inclination angle, so we see 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 have 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.

Fitting the line profiles simultaneously
We perform a fit in which we tie all of the parameters of the BLR model for the five lines, except those defining the radial distribution of the clouds, namely, µ, β, and F. In this approach, we can 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: 9 tied between all the line profiles and 3 left separate for each of the 5 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 Figure 6.The panel (a) shows that the model profiles of Hα, Hγ, and Paβ are very similar, while Hβ and He I are wider.This is similar to the results of the separate fitting.Because most of the model parameters are tied now, we can conclude that the line profile differences can be explained by the radial distributions of the line emission.Moreover, the combined fit results show tighter 68% confidence intervals than the separate fit.The improvement is the most obvious in 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.
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 ξ) are tied.We verified that if these are left free, the model yields different values of ξ to fit the individual asymmetries of the lines better.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 the model profiles yield consistent σ line to 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 so 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 et al. (2021a) tested fixing the R min = 4 ld in the SARM analysis and found a reasonable fitting result with β ≈ 1, interestingly, close to our β.However, GRAVITY Collaboration et al. (2021a) caution that additional restrictions bias their geometric distance measurement by about 30%.
The parameter that differs most from the SARM result is γ, where a large value is preferred, indicating that the vertical distribution of line emission is more concentrated towards the surface of the BLR.However, there is a long tail towards low values in the posterior distribution of γ.The number of clouds in the midplane 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 enable it to better reproduce the asymmetric features in broad lines of NGC 3783.

BLR 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. (2019Raimundo et al. ( , 2020)).Because we fix 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 distributions are somewhat different, as expected (see below).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).We now compare the BLR radius (R mean ) of the combined fit to the five lines (Figure 8a) with published measurements based on the Hβ and Brγ lines (Bentz et al. 2021a,b;GRAVITY Collaboration et al. 2021a,b).We focus on the joint SARM analysis because it obtains the tightest constraints of the model parameters using the spectro-astrometry and RM data simultaneously, and also because we have adopted its BH mass here.The R mean from the separate fitting show 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 con-trast, R mean derived from the combined fitting show much smaller uncertainties, with the Hβ R mean consistent with the SARM result as well.The He I R mean is the smallest among the five lines.We note that the posterior distributions of R mean are strongly correlated between different lines (Figure A.1).While such a correlation increases the uncertainties of the mean radii (≳ 30%, Figure 8a), as discussed below, it also means that the ratios of the mean radii have smaller uncertainties (≲ 10%).
As shown in Table 4, we measure 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 Hγ and He I lines are larger than the RM results.We also calculate the time lag ratios based on the radiation pressure confined (RPC) BLR model as 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. 2009Bentz et al. , 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), 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 line emissivity at all locations, are obtained naturally from the assumption that the clouds' column densities are large, and they are in total gas and radiation pressure equilibrium.The mean emissivity radius for each line is then computed from the model and then translated to 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 are also 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 (Figure 1a).Therefore, the small but systematic residual of the decomposition may lead to a too large Hγ BLR size compared with its actual dimensions.The clear systematic deviations between the data and model in Figure 6d support this point to some extent.This issue highlights the importance of high-quality line profiles to reveal 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 will make 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 [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  Article number, page 11 of 16 A&A proofs: manuscript no.ms The BLR size ratios derived from our simultaneous fitting.We randomly selected 500 sets of model parameters from the posterior samples and calculated the median R mean ratios and the 68% confidence intervals.Column (3): The time lag ratios measured by Bentz et al. (2010b).The Paβ time lag is not measured.Column (4): The time lag ratios derived by the radiation pressure confined (RPC) BLR model as described by Netzer (2020).We explored a range of model parameters that are in agreement with most RM measurements of various hydrogen and helium lines (e.g.Bentz et al. 2009Bentz et al. , 2010b)).The Hβ time lag is normalized to unity in all of the results.
spectro-astrometry and RM data of different lines in a joint analysis by sharing most of the BLR model parameters of the two lines.

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 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 have fixed the BH mass of NGC 3783, the value of f virial that we derive does not depend on M BH because the BLR radius is a free parameter in the fit.The reason, explained in Section 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 Figure 8b, the derived values from both separate fit and combined fit are consistent with that of the SARM joint analysis ( f virial = 2.52 +0.62 −0.53 ).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.

Caveats
In this work, we investigate the BLR structure by modeling the normalized single-epoch line profiles.We model the broadline emission and the associated kinematics with a Monte Carlo model of points without a physical size.This model is intended to avoid considering 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 fitted simultaneously, can provide most of the BLR model parameters consistently with the fitting including 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 estimate 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 adopt the ratio of Hβ line and hydrogen recombination coefficients of α eff Hβ /α B ≈ 1/8.5 and the UV photon flux ∞ ν 0 L ν hν dν ≈ 1.6-3.7 ×10 53 s −1 (ν 0 = 13.6 eV), which is estimated with the measured λL λ (5100 Å) ≈ 4.1 × 10 42 erg s −1 and the assumed AGN SED with low and intermediate Eddington ratios (Ferland et al. 2020;Jin et al. 2012) to enclose the range of typical Seyfert galaxy SEDs.We derive L Hβ ≈ 0.4-0.9×10 41 erg s −1 .The measured Hβ luminosity5 , ∼ 1.1 × 10 41 erg s −1 , 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 assume the maximum absorption of the UV photons with the model BLR geometry, and 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 Section 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.

Conclusions
We investigate the BLR structure of NGC 3783 using multiple broad lines in a high-resolution single-epoch spectrum obtained with VLT/X-Shooter.We decompose the strongest five broad lines (Hα, Hβ, Hγ, Paβ, and He I λ5876), and model 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 opt to fix the BH mass to a value reported in the literature and focus on the BLR structure and emissivity that can be derived from the line profiles.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.Our main results are, 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 develop 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 some 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 derive 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.

Fig. 1 .
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, with the residual underneath.The emission line components are also plotted separately.The shaded grey regions in the residuals represent wavelength ranges of bad channels and features due to the poor telluric correction, which are masked in the fitting.(a) the UV/optical spectrum with the Hα, Hβ, and Hγ lines.(b) the part of the NIR spectrum used in this study, with the Paβ line.(c) the [S II] doublet used as the narrow line template.

Fig. 2 .
Fig. 2. Normalized broad emission line profiles of Hα, Hβ, Hγ, Paβ, and He I.These were extracted as described in Sec.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 are replaced by the multi-Gaussian model as shown in Figure 1 for clarity.

Note.Fig. 3 .
Fig. 3. Comparison of the nonparametric properties of the line profiles (denoted by different colors, and with 1σ uncertainties).(a) line widths (circles, triangles, and squares represent W 25 , W 50 , and W 75 , respectively); (b) asymmetry index (A.I.); (c) kurtosis index (K.I.).The gray lines correspond to equivalent measurements on the model line profiles for the separate (dotted) and combined (solid) fitting results.

Fig. 4 .
Fig. 4. The 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.Table 3. Summary of the combined fitting results

Fig. 5 .
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 et al. 2021a).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 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 remove them for clarity.Details are discussed in Section 4.2.

Fig. 6 .
Fig.6.The 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 that of Figure4.

Fig. 7 .
Fig. 7. Posterior probability distribution of the five lines fitted simultaneously with tied geometry parameters.The panels and symbols are similar as that of Figure 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 et al. 2021a).

Fig. 8 .
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 grey regions represent the results from the SARM joint analysis (GRAVITY Collaboration et al. 2021a).

Table 1 .
Nonparametric properties of the broad line profiles

Table 2 .
Summary of the BLR model parameters, including a short explanation and prior range

Table 4 .
BLR size ratios ratios