In paper II we present new calibrations of the SE black hole mass estimators using the broad Hα, Hβ, MgII and CIV emission lines. The underlying assumptions is that Eqn. 1 holds for all broad emission lines and V_BLR can be estimated from the FWHM of the line in question (FWHM (line)). In this case: (3)MBHSE(FWHM,Lλ)=fFWHMG−1RBLRFWHM2 (e.g., Shen, 2013) where f_FWHM is the virial factor associated with the line FWHM and R_BLR is obtained from various RM studies (see e.g., Kaspi et al., 2000, 2005; Bentz et al., 2009, 2013, and references therein) and can be written as: (4)RBLR∝(Lλ)αline here L_λ stands for λL(λ) at various wavelengths: 6,200Å for the Hα line, 5,100Å for the Hβ line, 3,000Å for the Mg ii λ2798 line and 1,450Å for the C iv λ1549line.

It is important to note that M_BH estimates derived from high ionization lines such as C iv have a bias due the presence of an associated outflow blue shifted component. For objects with high Eddington ratios, the blue shifted component in C iv predominates over the component from the virialized region Sulentic, 2017.

In papers I and III we implemented an alternative method to estimate the black hole mass in type1-AGN based on fitting the SED of the accretion disk. We used the geometrically thin, optically thick accretion disc model from Slone and Netzer (2012) which are based on Shakura and Sunyaev (1973). Using this model we obtained successful fits in 37 out 39 objects in our sample. The model is fully determined by MBHAD, a_*, Ṁ, the inclination (i_LOS) with respect to the line-of-sight (LOS) and A_V. The procedure consisted of a Bayesian minimization over a grid of models covering a range in values for these parameters following the procedure described in paper III. Particularly, we assumed Gaussian priors centered around the observational estimations of MBHSE(Hα, L₆₂₀₀) and Ṁ. MBHSE(Hα, L₆₂₀₀) and Ṁ are calculated assuming a virial factor f_FWHM = 1. We also assume scatters in both quantities of 0.3 and 0.2 dex respectively. The minimization process uses 6 continuum windows and the posterior probability of the models is unaffected by the broad emission lines. The role of the priors is to penalize models which deviate significantly from the observational estimation of M_BH(Hα) and Ṁ_SE. In spite of this, the code can freely move in the parameter space either side of M_BH(Hα) and Ṁ_SE.

A major goal of the present paper is to set limits on M_BH under the assumption that thin accretion disks provide an accurate explanation for the observed continuum in AGNs. We therefore try to test the reliability of the results of our Bayesian procedures using different central values for our priors as well as different scatters around them (σ_MBH). In practice, this is done by considering a large range of virial factors (f_FWHM = {0.25, 0.4, 1, 2.5, 4}) to compute MBHSE(Hα, L₆₂₀₀) and Ṁ_SE and a range of σ_MBH for each. Each combination of f_FWHM and σ_MBH is used to find the best fit SED and calculate a range of posteriors exactly as done in papers I and III. We find that when the scatter around the central values is small (<0.4 dex), the posterior probability distributions of MBHAD and Ṁ depends on the chosen value of f_FWHM. However, when the scatter is large (>0.8 dex), we find that the posterior probability distributions of MBHAD and Ṁ are insensitive to the assumed f_FWHM and they all tend to converge to the posterior probability distributions of MBHAD and Ṁ when f_FWHM ~ 1 and the scatter is ~0.3 dex. This would indicate that f_FWHM ~ 1 is a close representation of the median virial factor of the AGN population within the M_BH range that our sample is covering (7.5 ≲ M_BH ≲ 10).

A major drawback of the above approach is the assumption that the simplified accretion disks SEDs used here, that neglect the effect of disk-winds, complex transfer in the disk atmosphere and other simplifications, provide an accurate description of the geometry and physics close to the central BH. In this approach, there is a degeneracy between the accretion rate and the inclination angle of the disk. For a given flux, larger inclinations will return larger intrinsic luminosities which in turn will return larger accretion rates. Fortunately, the derived black hole mass does not strongly depend on either inclination or accretion rate. As a consequence, the derived inclination does not bias the estimation of the black hole mass as it does in the virial SE method. From all the test described here, we can conclude AD mass estimations are reliable and of comparable accuracy to SE masses.

In this section we present different tests that we carried out to determine whether inclination is driving the correlation between f and FWHM_obs.

Hereafter when referring to log f_AD, MBHSE(FWHM, L_λ) and FWHM_obs we are meaning log f_AD(Hα), MBHSE (FWHM_obs(Hα)) and the observed value of FWHM_obs(Hα), unless otherwise specified. The reason to select the Hα line instead of the Hβ line for the following analysis is the better S/N and hence more accurate measurements of FWHM_obs(Hα) in our sample. As shown in earlier works, FWHM_obs in both Balmer lines are the same within uncertainties (Greene and Ho, 2005; Mejía-Restrepo et al., 2016).

The anti-correlation between log f_AD and FWHM_obs could be explained by the inclination of the axis of symmetry of a planar BLR with respect to the LOS. If we consider the median LOS inclination, i_median, at which Type-1 AGN are typically observed, we can also define a median virial factor f_median at which the SE M_BH calibration represents an accurate black hole mass for objects observed at i_median. Objects with narrower than usual broad emission lines are more likely observed at i < i_median (face-on orientations) and objects with broader than usual emission are more likely observed at i > i_median (edge-on orientations). This will produce too large (too small) SE mass estimates for objects with very broad (very narrow) emission lines, and would translate into a virial factor that anti-correlates with the line FWHMs.

For a planar BLR with a thickness ratio H/R and inclination i with respect to the line-of-sight we will have FWHMobs= FWHMint×sin2(i)+(H/R)2. Thus, for an ensemble of randomly orientated BLRs the final distribution of FWHM_obs will depend on (1) the intrinsic FWHM_int distribution and (2) the range of possible random orientations at which the BLR can be observed, both of which are, a priori, not known.

To check the inclination hypothesis we first need to determine the distribution of FWHM_int that is consistent with the probability density distribution (PDF) of the observed FWHM_obs. We then need to test whether it is possible to recover the anti-correlation of f with FWHM_obs. In other words, we need to test whether a population of randomly generated inclinations and FWHM_int that satisfy the PDF of FWHM_obs, can also account for the bidimensional distribution of the parameter space given by f and FWHM_obs.

We first assumed a thin BLR by taking H/R = 0. We computed the PDF as the product of two independent random variables (Glen et al., 2004) and applied it to the special case where FWHM_obs = FWHM_int × sin (i) (Lopez and Jenkins, 2012). For the FWHM_int distribution, we assumed an underlying truncated normal distribution with certain mean (FWHM_mean) and dispersion (FWHM_std). Our normal distribution was truncated to allow FWHM_int to vary between 1,000 and 30,000 km s⁻¹. We also assumed that our sample is limited to objects with line-of-sight inclination angles between i_min = 0° and i_max = 70°, with i_max determined by the torus opening angle. For an optimal exploration of the parameter space we ran a Monte Carlo Markov Chain simulation using the python code EMCEE (Foreman-Mackey et al., 2013). For the simulation we used 20 independent walkers and 5,000 iterations that mapped a total of 10⁵ models.

In the Figure S1 we compare the observed cumulative PDF (FWHM_obs) and its uncertainty (magenta thin line and shadowed region, respectively) with the predicted cumulative PDF from the model with the highest posterior probability (orange line). The parameters of this model are: i_min = 19°, i_max = 45°, FWHM_mean = 8,500, FWHM_std = 2,150, FWHM_min = 4,200 and FWHM_max = 30,000. Our model successfully reproduces the observed cumulative PDF. However, a simple normal distribution (turquoise line) is also consistent with the data and cannot be rejected. We also determined the best fit model for a distribution with FWHM_std = 0, i.e., effectively a single velocity. This model (dashed blue-line) is able to reproduce the distribution at low values of FWHM_obs, but it is unable to account for the distribution at large velocity widths.

First, we tested whether our thin BLR model is successful in reproducing the f−FWHM_obs distribution seen in the data. In the right panel of Figure 1 we show the predicted bi-dimensional probability density distribution of the virial factor and the observed FWHM_obs(Hα) as predicted by the thin BLR model. The Figure includes contours showing 25, 50, 75, and 99% confidence limits contours (black-thin lines) centred around the maximum probability point. We also superposed the data from the right panel of Figure 2 (open-blue circles). The magenta line represents the derived relation f = (FWHM_obs(Hα)/4,000 km s⁻¹). The thick yellow line is the median of the f−FWHM_obs distributions derived using a quantile non-parametric spline regression (COBS, Ng and Maechler, 2007). Analogously, the blue-dashed lines represent the 25, 50, and 75% quantiles of the observational distribution. To obtained these quantiles, for each observed data we randomly generated 1,000 points following the error distributions in f_AD and FWHM_obs(Hα) and then applied the COBS method to characterize the resulting distribution. We can notice that the median (50%-quantile) of the theoretical and observational distributions are in very good agreement. The scattered open-blue circles also show excellent agreement with the the bi-dimensional probability density function from the best model. Explicitly, we find that from our 37 objects, 21% fall inside the central 25% confidence level region, 51% fall inside the 50% confidence level region, 74% fall inside the 75% confidence level region, and 85% fall inside the 99% confidence level region.

In order to test the effects introduced by a thick BLR (0 < H/R < 1), we assumed a single H/R for all objects and followed the same steps outlined for the case of a thin BLR. We found that a wide range in BLR thickness ratios (H/R < 0.5) is able to reproduce the cumulative FWHM_obs PDF. However, objects with large thickness ratios clearly fail to reproduce the bi-dimensional distributions of f−FWHM_obs as can be seen in the right panel of Figure 2. We generally find that only relatively thin BLRs, i.e., those with H/R < 0.1, are able to reproduce both the bi-dimensional distributions and the cumulative FWHM_obs(Hα) PDF. In particular, for a BLR with H/R → 0, we find that the derived f_AD values constrain the range of inclinations at which the BLR is observed in our sample to 15° ≲ i ≲ 50°. This upper limit is consistent with typical expectations of a central torus hiding the BLR. We also find that the median virial factor in our sample, f = 0.95, corresponds to a median orientation of i_median = 31°.

In summary, our results show that a population of randomly orientated, thin BLRs can successfully reproduce our observations. We can thus conclude that inclination is very likely the main reason for the observed f−FWHM_obs correlations.

We finally considered the possibility that radiation pressure perturbations to the BLR motions might cause the observed f_AD−FWHM_obs dependency. A recent model considers the effects of radiation pressure in a BLR composed of pressure confined clouds, hence allowing the gas density of individual clouds to decrease with distance to the central black hole (Netzer and Marziani, 2010). In this model the system is still bound by gravity and FWHM_obs becomes smaller with increasing λ_Edd. The reason for this trend is that as λ_Edd increases, the clouds spend more time at large distances from the black hole, therefore increasing the median R_BLR and decreasing the median BLR Keplerian velocities. To account for this effect, the authors of this model proposed a modified expression for R_BLR: (6)RBLR=RBLR0[a1Lλαline+a2(Lλ/MBH)] where a₁ and a₂ are constants. The first term accounts for the observational relation described in Equation (4) and the second term represents a radiation pressure perturbation quantified by L_λ/M_BH ∝ λ_Edd. When replaced into the virial mass equation (Equation 1) this relation leads to a simple quadratic equation on M_BH with solution: (7)MBHrad=a102LαlineFWHMobs2[1+1+4 a20Lλ1−2αlinea102 FWHMobs2] or equivalently: (8)frad∝[1+1+4 a20Lλ1−2αlinea102 FWHMobs2] where M_BH^rad and f_rad are the black hole mass and virial factor for a radiation pressure dominated BLR. a10=a1f0 RBLR0G−1, a20=a2f0 RBLR0G−1, and f₀ is a normalization constant. In the case when 4 a20Lλ1−2αline/a102 FWHMobs2≫1 this would result in a close agreement with the inverse proportionality between f_AD and FWHM_obs found in our data. Given that α_line is found to be ~0.6 for all lines (see Table S1), this would translate into an explicit dependency of f on L_λ. We would then expect that the scatter in the f_AD−FWHM_obs relation should be driven by L_λ. In the right panel of Figure 1 larger (smaller) values of L₅₁₀₀ are represented by redder (bluer) colors. We can see that there is no clear suggestion that the scatter in driven by L₅₁₀₀ in any of the lines. Note however that the relatively narrow range in L₅₁₀₀ covered by our sample (from L₅₁₀₀ = 2.0 × 10⁴⁴ to L₅₁₀₀ = 1.6 × 10⁴⁶ ergs/s, corresponding to a factor of 80), together with the uncertainties in our estimations of f, do not allow us to rule out this mechanism.

Testing this model further, we found the combination of parameters a₁, a₂, and f₀ that best reproduce our MBHAD measurements and the observed relation between f and FWHM_obs for the Hα line. To obtain dimensionless values for a₁ and a₂ we expressed M_BH, L_λ and FWHM in units of 10⁸ M_⊙, 10⁴⁴erg s⁻¹ and 1,000 km s⁻¹, respectively. Taking α_line = 0.63, as suggested by the observations (see Table S1), we carried out a Monte-Carlo Markov Chain exploration of the parameter space of the model and found that a₁ = 0.88, a₂ = 0.36 and f_o = 0.51 are able to reproduce our MBHAD measurements with a scatter of 0.12 dex, preserving the experimental dependence of R_BLR on L_λ as expressed in Equation 4 with a scatter of 0.05 dex. At the same time the results are able to reproduce the observed f−FWHM_obs relation with a scatter of 0.11 dex [see Figure 3, which presents our observations (black squares with error bars) together with the prescribed values for f as given by Equation 8 (colored circles without error bars)]. However, we also found that the residuals between the predicted values and the best fit to the correlation are heavily correlated with L₅₁₀₀ (r_s > 0.63, P_s < 2 × 10⁻⁵), as can be seen by the color gradient of our simulated points in the direction perpendicular to the correlation best fit in Figure 3. This bias is introduced by the explicit dependence of f_rad on L_λ which is not observed in our sample, although notice that the error bars of our derived f values are of the order of, if not larger, than the expected dependence (see Figure 3). Finally, the dependency on L₅₁₀₀ vanishes when α_line = 0.5. For this case, however, we were unable to reproduce any the observables. Extending our sample towards lower luminosities will allow us to determine a possible dependence of f_rad on L_λ. This should yield the final test to be able to confidently conclude whether this model can be the driving mechanism for the observed f−FWHM_obs correlation.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer, AN, and handling Editor declared their shared affiliation.