Scatter Analysis along the Multidimensional Radius–Luminosity Relations for Reverberation-mapped Mg ii Sources

The full sample includes 68 objects with $42.8\,\lt \mathrm{log}{L}_{3000}\,\lt 46.8$ at $0.003\lt z\lt 1.89$. It includes all the objects with reverberation-mapped Mg ii measurements reported to date. The sample includes the 57 Mg ii time lags (henceforth SDSS-RM sample) from the recent SDSS-RM monitoring (Homayouni et al. 2020) and the six objects previously monitored by the same project (Shen et al. 2016). Only one object from the previous SDSS-RM monitoring is included in the recent one; both measurements are considered. We also include the high-luminosity objects: CTS 252 monitored by Lira et al. (2018), CTS C30.10 measured by Czerny et al. (2019a), HE 0413-4031 monitored by Zajaček et al. (2020), and the two old IUE measurements reported for NGC 4151 by Metzroth et al. (2006).

The full sample considered is relatively homogeneous, since $\sim 83 \%$ of the objects come from the SDSS-RM sample and $\sim 9 \%$ from their previous program. The rest of the objects with the lowest (NGC 4151) and the highest (CTS C30.10, HE 0413-4031, and CTS 252) luminosities are crucial for the detection of the trends in the R–L relation. For the recent SDSS-RM sample, the time delay is estimated using the JAVELIN method (Zu et al. 2011, 2013, 2016), while for the rest of the sources other methods are applied, specifically the interpolated cross-correlation function (ICCF; Gaskell & Peterson 1987; Peterson et al. 1998b), the discrete correlation function (DCF; Edelson & Krolik 1988), the z-transformed DCF (zDCF; Alexander 1997), the light-curve similarity estimators (Von Neumann or Bartels estimators; Chelouche et al. 2017), and the ${\chi }^{2}$ method (Czerny et al. 2013, 2019a; Zajaček et al. 2020), or an average of them. Recent statistical analyses (Li et al. 2019; Yu et al. 2020) point out that the JAVELIN method is more powerful than other traditional methods (ICCF and zDCF) to recover time delays. However, the analyses were performed on the mock sample where the continuum light curve is generated using the damped random-walk (DRW) process, and the JAVELIN method makes use of the DRW to recover the time delay. Therefore, there may be a bias since the stochastic, red-noise part of the variable continuum is omitted in the generated light curves from their power-density spectra. In particular, the model-independent, discrete methods such as zDCF and the von Neumann estimator may still be more suitable for the analysis of irregular and heterogeneous light curves of more distant quasars (Czerny et al. 2019a; Zajaček et al. 2019, 2020). In particular, for the highly accreting source HE 0413-4031 at z = 1.39, Zajaček et al. (2020) saw a difference between JAVELIN and ICCF on one hand and discrete methods (DCF, zDCF, and von Neumann) on the other hand, in the determination of the primary time delay peak, $\tau \sim 431$ and ∼303 days, respectively.

In addition, each research group assumes different time delay significance criteria, which may not be satisfied for the objects of other samples. For example, the recent SDSS-RM monitoring (Homayouni et al. 2020) considers that a time delay is statistically significant if (1) a primary time-lag peak includes at least 60% of the weighted time delay posterior samples and (2) the time delay is well detected at the 3σ level different from the zero time delay. Therefore, for completeness, we consider all the objects as we analyze a mixed sample, for which it is difficult to establish general statistically robust criteria.

For all the sources in the sample, we collected several measured parameters, which are summarized in Table B1 in Appendix B. For the SDSS-RM sample, the luminosity at 3000 Å, FWHM, and equivalent width (EW) of Mg ii and Fe ii (at 2250–2650 Å) were taken from the Shen et al. (2019) catalog, which assumes a flat ΛCDM cosmology with ${{\rm{\Omega }}}_{{\rm{\Lambda }}}$ = 0.7 and H₀ = 70 km s⁻¹ Mpc⁻¹. The time delay is taken from Homayouni et al. 2020. They estimated time delays of Mg ii by two methods: JAVELIN and CREAM. In this contribution we use the JAVELIN estimations since they are more reliable according to the authors. For the rest of the sample, the measurements were taken from the compilation done by Zajaček et al. (2020, see their Table 2).

Since in the optical range the strength of the Fe ii broad line shows a clear correlation with the accretion rate intensity (Boroson 2002; Negrete et al. 2018; Du & Wang 2019; Panda et al. 2019b), we explore a similar correlation considering the UV Fe ii. The relation between the UV Fe ii and the accretion parameters is scarcely studied (Dong et al. 2011), but since the UV and optical Fe ii are correlated, and in turn both show similar anticorrelations with the FWHM of Hβ (Dong et al. 2009, 2011; Kovačević-Dojčinović & Popović 2015; Śniegowska et al. 2020), it opens the possibility that the UV Fe ii is also related to the accretion rate intensity as is the optical Fe ii emission. In Appendix A, we include a discussion about this issue considering the parameters estimated below.

The UV Fe ii strength is estimated using the parameter R_{Fe ii} , which is defined as the ratio of the equivalent width of Fe ii pseudo-continuum measured at 2250–2650 Å to the EW of the Mg ii line:

$$\begin{eqnarray}&&{R}_{\mathrm{Fe}{\rm\small{ii}}}=\displaystyle \frac{\mathrm{EW}\ (\mathrm{UV}\ \mathrm{Fe}\,{\rm\small{ii}})}{\mathrm{EW}(\mathrm{Mg}\,{\rm\small{ii}})}.\end{eqnarray} \tag{ 1 }$$

The wavelength range selected for the EW of UV Fe ii is defined by the one reported in the Shen et al. (2019) catalog. In the case of quasars CTS C30.10 and HE 0413-4031, UV Fe ii was directly measured only in the 2700–2900 Å range and for a different UV Fe ii template, but for consistency we refitted these spectra using the UV Fe ii template of Vestergaard & Wilkes (2001) and rescaled the newly derived EW(UV Fe II) from 2700–2900 Å to 2250–2650 Å by the factor of 2.32 appropriate for the Vestergaard & Wilkes (2001) UV Fe ii template. For NGC 4151 and CTS 252, the EW of Mg ii and UV Fe ii is not reported in the required wavelength range, and hence we do not consider these objects in the analysis where R_{Fe ii} is used.

We also estimated the level of variability using the F_var parameter (Rodríguez-Pascual et al. 1997), defined by

$$\begin{eqnarray}&&{F}_{\mathrm{var}}=\displaystyle \frac{{\left({\sigma }^{2}-{{\rm{\Delta }}}^{2}\right)}^{1/2}}{\langle f\rangle },\end{eqnarray} \tag{ 2 }$$

where ${\sigma }^{2}$ is the variance of the flux, Δ is the mean square value of the uncertainties (${{\rm{\Delta }}}_{i}$) associated with each flux measurement (f_i ), and $\langle f\rangle$ is the mean flux. F_var or similar expressions of the variability level are also correlated with the accretion rate intensity (Wilhite et al. 2005; MacLeod et al. 2010; Sánchez-Sáez et al. 2018; Martínez-Aldama et al. 2019). For the SDSS-RM sample, we use the fractional rms variability provided by the Shen et al. (2019) catalog (see their Table 2) for estimating F_var, while the F_var values for the rest of the objects were taken from their originals works.

Furthermore, we consider the secondary parameters that can be derived from the basic measurements and in turn from the black hole mass: the dimensionless accretion rate ($\dot{{ \mathcal M }}$) and the Eddington ratio (${L}_{\mathrm{bol}}/{L}_{\mathrm{Edd}}$). Several results indicate that line dispersion (${\sigma }_{\mathrm{line}}$) is a better estimator of the velocity field and in turn a better entity for the estimation of the black hole mass, because it is less biased than the FWHM (e.g., Peterson et al. 2004; Denney et al. 2013; Dalla Bontà et al. 2020; Yu et al. 2020a). However, since ${\sigma }_{\mathrm{line}}$ is not reported for the SDSS-RM sample, which includes the majority of the sources used in this paper, our estimations are done considering the FWHM taken from the Shen et al. (2019) catalog. Then, to get the black hole mass (M_BH) estimation, we consider a virial factor anticorrelated with the FWHM of the emission line defined as ${f}_{{\rm{c}}}={\left({{\rm{FWHM}}}_{{\rm{Mg}}{\rm{II}}}/3200\pm 800\right)}^{-1.21\pm 0.24}$ (Mejía-Restrepo et al. 2018), the time delay (${\tau }_{\mathrm{obs}}$) of Mg ii reported in Table B1 (Appendix B), and the virial relation, M_BH = ${f}_{c}\,c\,{\tau }_{{obs}}\,{\mathrm{FWHM}}^{2}/G$, where c is the speed of light and G the gravitational constant.

The dimensionless accretion rate ($\dot{\,{ \mathcal M }}$) was introduced by Wang et al. (2014a) assuming a Shakura–Sunyaev (SS) disk. Originally, $\dot{{ \mathcal M }}$ is adjusted for the continuum at 5100 Å, and therefore it must be rescaled for other wavelengths. Following the standard SS accretion disk, we get

$$\begin{eqnarray}&&\dot{{ \mathcal M }}=\alpha \ {\left(\displaystyle \frac{\lambda }{5100}\right)}^{-1/2}{L}_{\lambda }^{3/2}\,{M}_{\mathrm{BH}}^{-2},\end{eqnarray} \tag{ 3 }$$

where λ is the wavelength of the continuum luminosity (${L}_{\lambda }$). Comparing with the formula given by Wang et al. (2014a), $\dot{\,{ \mathcal M }}$ at 3000 Å is given by

$$\begin{eqnarray}&&\dot{{ \mathcal M }}=26.2{\left(\displaystyle \frac{{L}_{44}}{\cos \theta }\right)}^{3/2}{m}_{7}^{-2},\end{eqnarray} \tag{ 4 }$$

where ${L}_{44}$ is the luminosity at 3000 Å in units of 10⁴⁴ erg s⁻¹, θ is the inclination angle of disk to the line of sight, and m₇ is the black hole mass in units of 10⁷ ${M}_{\odot }$. We considered cos ${\text{}}\theta$ = 0.75, which is the mean disk inclination for type 1 AGNs and is in agreement with the typical torus covering factor found (e.g., Lawrence & Elvis 2010; Ichikawa et al. 2015).

Previously, in Zajaček et al. (2020) we used the $\dot{\,{ \mathcal M }}$ definition adjusted for 5100 Å; therefore, the $\dot{\,{ \mathcal M }}$ values reported there must be rescaled simply multiplying by a factor of 1.3038. The corrected $\dot{\,{ \mathcal M }}$ values for the Zajaček et al. (2020) sample are already included in Table B1.

The Eddington ratio is defined as ${L}_{\mathrm{bol}}/{L}_{\mathrm{Edd}}$ = $\mathrm{bc}\,\cdot {L}_{3000}/{L}_{\mathrm{Edd}}$, where $\mathrm{bc}=5.62\pm 1.14$ (Richards et al. 2006) and ${L}_{\mathrm{Edd}}=1.5\times {10}^{38}\left(\tfrac{{M}_{\mathrm{BH}}}{{M}_{\odot }}\right)$. In Table B1, we report the dimensionless accretion rates and the Eddington ratios. We use both the dimensionless accretion rate and the Eddington ratio since it is not clear which of these two dimensionless values is determined more accurately. Eddington ratio requires the knowledge of the broadband spectrum or use of the bolometric correction, which neglects, for example, the spin issue in a quasar and is based on a specific spectral energy distribution (SED) shape. The $\dot{\,{ \mathcal M }}$ determination does not require the knowledge of the bolometric correction, but the dependence on the black hole mass determination is much stronger than for the Eddington ratio, and the spectral slope in the quasar might not be consistent with the theoretical power law expected from the Shakura & Sunyaev (1973) accretion disk model.

First, we use the full reverberation-mapped Mg ii sample and study the relation between the measured time delay of the Mg ii line and the continuum flux at 3000 Å. In Figure 1 (top left panel), a clear linear trend is visible, although the rms scatter, defined using the sample size N, observed values ${\tau }_{i}$, and predicted values ${\hat{\tau }}_{i}$ as

$$\begin{eqnarray}&&{\sigma }_{\mathrm{rms}}=\sqrt{\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}{({\tau }_{i}-{\hat{\tau }}_{i})}^{2}},\end{eqnarray} \tag{ 5 }$$

is considerable, with ${\sigma }_{\mathrm{rms}}=0.3014$ dex. The best-fit parameters are given in Table 1, where we also give the values of the rms scatter and the Pearson correlation coefficient (r). The correlation is moderately strong but clearly visible. The scatter is visibly larger than in Zajaček et al. (2020), where only 11 sources were available and the R–L relation for all the available sources at that time had a scatter of ${\sigma }_{\mathrm{rms}}=0.22$ dex and ${\sigma }_{\mathrm{rms}}=0.19$ dex when two sources with the largest offset were removed. The slope is now also much shallower than given by Zajaček et al. (2020) (0.42 ± 0.05 for all 11 sources and 0.58 ± 0.07 with two outliers removed; see their Figure 5, right panel) or obtained earlier by McLure & Jarvis (2002) at the basis of the Mg ii line shape. The R–L luminosity in the past was much better studied for the ${\rm{H}}\beta$ time delay with respect to the continuum at 5100 Å, and then the slope was close to 0.5 when the host contamination was carefully taken into account (e.g., Bentz et al. 2013). The correlation for the Mg ii line is partially (but not entirely) driven by the extreme points representing the lowest (NGC 4151) and the highest (CTS C30.10, HE 0413-4031 and CTS 252) luminosity sources.

Zoom In Zoom Out Reset image size

Since it is generally accepted now that the large dispersion is related to the range of dimensionless accretion rates (Du et al. 2016; Yu et al. 2020b; Zajaček et al. 2020), we mark the points in the plot with the color, corresponding to $\dot{{ \mathcal M }}$ values. A certain degree of dependence is visible, but the trend is not at all clear—yellower points occupy mostly the lower part of the diagram, but they also concentrate more toward higher values of the luminosity. Zajaček et al. (2020) achieve the considerable reduction of the scatter when the observed time delay was corrected for the trend with $\dot{\,{ \mathcal M }}$ or ${L}_{\mathrm{bol}}/{L}_{\mathrm{Edd}}$, particularly if the virial factor f_c depending on FWHM (Mejía-Restrepo et al. 2018; Martínez-Aldama et al. 2019) instead of a constant value has been used. In their work, Lusso & Risaliti (2017) use the additional dependence on FWHM to reduce the scatter in their quasar relation between the X-ray and the UV flux, and Du & Wang (2019) and Yu et al. (2020b) applied the additional dependence on the $\dot{\,{ \mathcal M }}$ and R_{Fe ii} to reduce the scatter in the R–L relation for Hβ. Thus, we can expect that including more parameters in the fit will lead to the reduction in the scatter to a smaller or larger extent.

The need to reduce the large scatter along the R–L relation motivates us to search for extended, multidimensional radius–luminosity relations that involve linear combinations of the logarithms of monochromatic luminosity (L₄₄) with additional parameters, which can generally be written as $\mathrm{log}{\tau }_{\mathrm{obs}}\,={K}_{1}\mathrm{log}{L}_{44}+{\sum }_{i=2}^{n}{K}_{i}\mathrm{log}{Q}_{i}+{K}_{n+1}$, where the parameters Q_i are typically related to the accretion rate and we use $\dot{\,{ \mathcal M }}$, ${L}_{\mathrm{bol}}/{L}_{\mathrm{Edd}}$, F_var, and R_{Fe ii} for this purpose. The total number of quantities is typically n = 2; in a few cases we have n = 3. FWHM is either included or in some cases omitted from the linear combination. An overview of all studied cases is in Table 1 (third column).

We use the Python packages sklearn and statsmodels to perform a multivariate linear regression and to obtain regression coefficients K₁, K_i , and ${K}_{n+1}$, including their standard errors, and the correlation coefficients (r), which are listed in Table 1 for each linear combination. In addition, we also include the rms scatter calculated using Equation (5).

We list all the relations in Table 1, and graphically they are shown in Figures 1 and 2. In these figures, for an easier comparison of the slopes ${\tau }_{\mathrm{obs}}\propto {L}_{44}^{\alpha }$, we plot $\mathrm{log}{\tau }_{\mathrm{obs}}$ versus $\mathrm{log}{L}_{44}+{\sum }_{i=2}^{n}{K}_{i}/{K}_{1}\mathrm{log}{Q}_{i}+{K}_{n+1}/{K}_{1}$, where K₁ represents the slope α.

Table 1. Results of the Parameter Inference Applied to the Logarithm of the Observed Time Delay (${\tau }_{\mathrm{obs}}$) Expressed as a Linear Combination of the Logarithm of the Monochromatic Luminosity (Expressed as L₃₀₀₀ or L₄₄) and Other Quantities

Sample	Inference	$\mathrm{log}{\tau }_{\mathrm{obs}}=$	K1	K2	K3	K4 or f	${\sigma }_{\mathrm{rms}}$ [dex]	r
All	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}$	0.298 ± 0.047	1.670 ± 0.053	⋯	⋯	0.3014	0.62
Low	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}$	0.520 ± 0.078	1.732 ± 0.056	⋯	⋯	0.2815	0.76
High	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}$	0.414 ± 0.058	1.382 ± 0.085	⋯	⋯	0.2012	0.78
All	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}{\mathrm{logFWHM}}_{3}+{K}_{3}$	0.30 ± 0.05	0.29 ± 0.28	1.49 ± 0.18	⋯	0.2990	0.63
Low	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}{\mathrm{logFWHM}}_{3}+{K}_{3}$	0.54 ± 0.08	−0.59 ± 0.40	2.12 ± 0.27	⋯	0.2722	0.78
High	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}{\mathrm{logFWHM}}_{3}+{K}_{3}$	0.42 ± 0.06	−0.12 ± 0.32	1.44 ± 0.17	⋯	0.2007	0.78
All	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}\mathrm{log}\dot{{ \mathcal M }}+{K}_{3}$	0.694 ± 0.022	−0.432 ± 0.018	1.442 ± 0.019	⋯	0.0947	0.97
All	MCMC	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}\mathrm{log}\dot{{ \mathcal M }}+{K}_{3}$	${0.6965}_{-0.0098}^{+0.0102}$	$-{0.4550}_{-0.0054}^{+0.0054}$	${1.4618}_{-0.0128}^{+0.0124}$	⋯	0.0985	0.97
All	MCMC (f)	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}\mathrm{log}\dot{{ \mathcal M }}+{K}_{3}$	${0.6984}_{-0.0234}^{+0.0236}$	$-{0.4558}_{-0.0125}^{+0.0127}$	${1.4600}_{-0.0293}^{+0.0297}$	$f={2.3730}_{-0.1916}^{+0.2292}$	0.0985	0.97
All	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}\mathrm{log}L/{L}_{{Edd}}+{K}_{3}$	0.910 ± 0.029	−0.863 ± 0.035	0.380 ± 0.056	⋯	0.0947	0.97
All	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}\mathrm{log}{R}_{\mathrm{Fe}{\rm\small{}}\mathrm{ii}}+{K}_{3}$	0.21 ± 0.06	0.16 ± 0.16	1.77 ± 0.07	⋯	0.2862	0.51
Low	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}\mathrm{log}{R}_{\mathrm{Fe}{\rm\small{}}\mathrm{ii}}+{K}_{3}$	0.27 ± 0.11	0.35 ± 0.18	1.92 ± 0.08	⋯	0.2521	0.68
High	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}\mathrm{log}{R}_{\mathrm{Fe}{\rm\small{}}\mathrm{ii}}+{K}_{3}$	0.47 ± 0.06	1.04 ± 0.33	1.25 ± 0.08	⋯	0.1718	0.84
High	MCMC	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}\mathrm{log}{R}_{\mathrm{Fe}{\rm\small{}}\mathrm{ii}}+{K}_{3}$	${0.4749}_{-0.0178}^{+0.0177}$	${1.0647}_{-0.1146}^{+0.1120}$	${1.2678}_{-0.0328}^{+0.0327}$	⋯	0.1743	0.85
High	MCMC (f)	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}\mathrm{log}{R}_{\mathrm{Fe}{\rm\small{}}\mathrm{ii}}+{K}_{3}$	${0.4761}_{-0.0397}^{+0.0404}$	${1.0735}_{-0.2725}^{+0.2663}$	${1.2659}_{-0.0751}^{+0.0718}$	$f={2.3420}_{-0.2671}^{+0.3458}$	0.1743	0.85
All	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}\mathrm{log}{F}_{\mathrm{var}}+{K}_{3}$	0.307 ± 0.056	0.048 ± 0.178	1.712 ± 0.165	⋯	0.3012	0.62
Low	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}\mathrm{log}{F}_{\mathrm{var}}+{K}_{3}$	0.550 ± 0.087	0.173 ± 0.223	1.879 ± 0.198	⋯	0.2788	0.77
High	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}\mathrm{log}{F}_{\mathrm{var}}+{K}_{3}$	0.370 ± 0.058	-0.406 ± 0.179	0.980 ± 0.194	⋯	0.1863	0.82
All	OLS	${K}_{1}\mathrm{log}{L}_{44}+{K}_{2}{\mathrm{logFWHM}}_{3}+{K}_{3}\mathrm{log}{F}_{\mathrm{var}}+{K}_{4}$	0.306 ± 0.056	0.280 ± 0.282	0.032 ± 0.178	${K}_{4}=1.522\pm 0.252$	0.2989	0.63

Note. We analyze parameter values and the scatter for the whole Mg ii sample (denoted as "All") or its low- or high-accretion subsamples (denoted as "Low" or "High," respectively), which is specified in the first column. The parameter inference type—ordinary least squares (OLS) or Markov Chain Monte Carlo (MCMC)—is specified in the second column. The third column lists the analyzed parameter combination. The other columns contain coefficient values with standard errors within 1σ, rms scatter (${\sigma }_{\mathrm{rms}}$) along the linear relations in dex, and the Pearson correlation coefficient (r). For the two cases with the smallest rms scatter, namely, combinations with $\mathrm{log}{L}_{44}$, $\mathrm{log}\,\dot{\,{ \mathcal M }}$ and $\mathrm{log}{L}_{44}$, $\mathrm{log}{R}_{\mathrm{Fe}{\rm\small{}}\mathrm{ii}}$, we perform MCMC fitting to cross-check our linear regression results with and without the underestimation factor f, which, if included, is listed in the same column as the K₄ coefficient.

Download table as:  ASCIITypeset image

The smallest scatter of ${\sigma }_{\mathrm{rms}}\sim 0.1\,\mathrm{dex}$ is for combinations that include L₄₄ and either $\dot{\,{ \mathcal M }}$ or ${L}_{\mathrm{bol}}/{L}_{\mathrm{Edd}}$ (see Figure 2). When, in addition to $\dot{\,{ \mathcal M }}$ or ${L}_{\mathrm{bol}}/{L}_{\mathrm{Edd}}$, FWHM is added to the combination, the scatter is of the order of ${\sigma }_{\mathrm{rms}}\sim {10}^{-4}\,\mathrm{dex}$ only. However, in this case, the correlation is artificially enhanced (the correlation coefficient is essentially 1.00), as both $\dot{{ \mathcal M }}$ and ${L}_{\mathrm{bol}}/{L}_{\mathrm{Edd}}$ depend on ${\tau }_{\mathrm{obs}}$ and FWHM via the black hole mass. The use of such interdependent quantities is not welcome, as that may easily create an apparent correlation, and the subsequent error determination when such a relation is used has to take that into account. For this reason, we do not include these two cases in the overview of linear combinations in Table 1.

Zoom In Zoom Out Reset image size

Considering independent quantities like FWHM, R_{Fe ii} , and F_var decreases the scatter, but the effect is relatively small and is comparable to the original R–L relation for Mg ii, ${\sigma }_{\mathrm{rms}}\sim 0.3\,\mathrm{dex}$ (see Figure 1). The scatter drops by $0.8 \%$ when FWHM is added to L₄₄ and only by $0.07 \%$ when F_var is added to ${L}_{44}$. This is not a significant improvement in terms of the scatter, but the relations define planes instead of a line and thus connect three or four independent observables, which can be relevant in terms of understanding mutual relations among them.

An additional dependence on R_{Fe ii} gives a better result: the scatter drops from 0.3 dex down to 0.286 dex (see Table 1), but the correlation coefficient also drops, so apparently R_{Fe ii} adds considerably to the scatter. Also, the slope of the relation is shallower. Although R_{Fe ii} and F_var are physically related to the Eddington ratio, neither of these quantities leads to the considerable improvement when the whole sample is considered.

Combining even more quantities does not provide an improvement either: a combination of L₄₄, FWHM, F_var, and R_{Fe ii} still leads to almost the same scatter and the same value of the correlation coefficient as in the basic R–L relation. Other combinations also do not provide an improvement.

To make use of the strong dependence on $\dot{\,{ \mathcal M }}$ in the reduction of the scatter in Section 3.2, we divided the full sample into two subsamples: low and high accretors. Although limits for low and high accretors have previously been discussed (e.g., Marziani et al. 2003; Du et al. 2015), in this work we consider a division that gives us the subsamples containing comparable numbers of objects. Thus, as a reference, we consider the median $\dot{\,{ \mathcal M }}$ value (log $\dot{\,{ \mathcal M }}$ = 0.2167) to get an equal number of sources (34 objects) in each subsample. In Figure 3, the $\dot{\,{ \mathcal M }}$ distribution is shown for low and high accretors. In Section 4.3, we include a discussion of the accretion rates observed in our sample in a general context considering samples like the DR7 (Shen et al. 2011) and DR14 (Rakshit et al. 2020), which also support the $\dot{\,{ \mathcal M }}$ division considered in this work.

Zoom In Zoom Out Reset image size

The division of the sample into high and low accretors results in significantly reducing the scatter, particularly for the highly accreting subsample (see Table 1 and Figure 4). It is also interesting to note that the slopes in the R–L relation in both subsamples are steeper than for the whole sample and much closer to the theoretically expected value of 0.5. Comparing both cases, for low accretors, the slope is steeper, closer to the canonical 0.5 value than for the highly accreting subsample, but the difference is within the quoted slope errors.

Zoom In Zoom Out Reset image size

The Pearson correlation coefficient increased significantly for both subsamples, supporting the view that mixing sources with different accretion rates spoils the pattern. The dispersion in both subsamples decreased in comparison with the whole sample, but here the effect in the two subsamples is clearly different. In the low-$\dot{\,{ \mathcal M }}$ subsample, the reduction in the dispersion is not that strong, from 0.3 dex down to 0.28 dex. However, the drop in the dispersion for the highly accreting sources is spectacular, from 0.3 dex to 0.2 dex. We should stress here again that this division has been set at $\dot{\,{ \mathcal M }}$ = 1.65, and the definition of $\dot{\,{ \mathcal M }}$ does not include the accretion efficiency, so for a standard efficiency of 10% this corresponds to mild accretion rates above 0.165 in dimensionless units. Hence, a considerable fraction of quasars in the SDSS-RM objects (Homayouni et al. 2020) belong to this category (see Section 4.3).

For the low accretion rate subsample, our dispersion in the R–L relation for Mg ii is still higher than the dispersion of 0.19 dex obtained by Bentz et al. (2013) for R–L in Hβ, without the removal of outliers. However, our sources are, on average, brighter than those in Bentz et al. (2013), where most of the sources are at $\mathrm{log}{L}_{5100}\sim 43.5$, and we study a different emission line. For high accretion rate sources, our dispersion is comparable to Bentz et al. (2013), and we did not remove any outliers in our analysis. Removing outliers (e.g., by 3σ clipping) would clearly tighten the correlation, but we use the data from the literature and we do not think we can reliably eliminate some of the available measurements. This shows that for Mg ii quasar population the highly accreting sources are much more attractive for cosmological applications.

We also studied the subsamples allowing for additional parametric dependencies: FWHM, R_{Fe ii} , and F_var. The inclusion of FWHM in the fit gave some further decrease of the dispersion, but the effect is not strong (see bottom panels of Figure 4). Concerning the fractional variability F_var, for the low-accretion sources, the improvement was marginal in comparison with the base R–L relation for these sources; see Figure 5 (bottom left panel). However, moving toward the highly accreting subsample, the reduction in scatter is significant, down to 0.19 dex (see Figure 5, bottom right panel), and the correlation coefficient also increased. An even larger reduction in the scatter was achieved when we included R_{Fe ii} (see Figure 5, top panels). The effect was clearly visible for both subsamples, and for the high accretion rate subsample, our dispersion reduced to only 0.17 dex, again without the removal of any outliers. This was the smallest scatter achieved in our study without the use of interdependent quantities.

Zoom In Zoom Out Reset image size

For all the studied relations $\mathrm{log}{\tau }_{\mathrm{obs}}={K}_{1}\mathrm{log}{L}_{44}\,+{\sum }_{i=2}^{n}{K}_{i}\mathrm{log}{Q}_{i}+{K}_{n+1}$ listed in Table 1, we used the higher-dimensional ordinary least squares (OLS) method (2D or 3D including the constant factor) as implemented in Python packages sklearn and statsmodels. This allowed us to quickly infer the relevant parameters including standard errors and compare the rms scatter and the correlation coefficient for as many as 15 combinations of relevant quantities. However, hidden correlations between quantities are not apparent when using the multidimensional least squares. Therefore, for the relations with the smallest scatter below 0.2 dex, specifically $\mathrm{log}{\tau }_{\mathrm{obs}}={K}_{1}\mathrm{log}{L}_{44}+{K}_{2}\mathrm{log}\,\dot{\,{ \mathcal M }}\,+{K}_{3}$ (whole sample) and $\mathrm{log}{\tau }_{\mathrm{obs}}={K}_{1}\mathrm{log}{L}_{44}+{K}_{2}\mathrm{log}{R}_{\mathrm{Fe}{\rm\small{}}\mathrm{ii}}+{K}_{3}$ (highly accreting sources), we apply the Markov Chain Monte Carlo (MCMC) inference of parameters, using the Python sampler emcee. The two inference techniques—OLS and MCMC—are specified in the second column of Table 1. The MCMC is a robust Bayesian technique that takes into account measurement errors while inferring the parameters with the maximum likelihood, which in our case is defined as

$$\begin{eqnarray}&&{ \mathcal L }=-\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}{\left(\displaystyle \frac{{\tau }_{i}-{\hat{\tau }}_{i}}{{\sigma }_{i}^{\tau }}\right)}^{2},\end{eqnarray} \tag{ 6 }$$

where ${\tau }_{i}$ are individual time delay measurements, ${\sigma }_{i}^{\tau }$ are the measurement errors of the time delay, and ${\hat{\tau }}_{i}$ are the predicted time delay values according to the inferred model. In addition, MCMC allows us to construct 2D histograms to see potential degeneracies among different parameters.

The results for the two above-mentioned relations are shown in Figure 6 for the combination including $\dot{\,{ \mathcal M }}$ in the left panels and for the combination including R_{Fe ii} in the right panels. From 2D histograms we can see that the combination with $\dot{\,{ \mathcal M }}$ exhibits a degeneracy between the constant K₃ coefficient and K₁ coefficient, as well as between the K₃ and K₂ coefficients. This can be understood in terms of $\dot{\,{ \mathcal M }}$ being intrinsically correlated with ${\tau }_{\mathrm{obs}}$. The degeneracy between the K₃ and K₂ parameters is lifted for the combination with R_{Fe ii} , as R_{Fe ii} and ${\tau }_{\mathrm{obs}}$ are intrinsically not correlated.

Zoom In Zoom Out Reset image size

As can be seen from the bottom panels of Figure 6, the parameters with the maximum likelihood are consistent with the best-fit parameters from the least-squares technique. The rms scatter and correlation coefficients are also comparable within uncertainties; see Table 1. However, the MCMC uncertainties based on the 16th and 84th percentiles (1σ) are smaller by a factor of about two than the 1σ errors inferred using the OLS for the combination involving $\dot{\,{ \mathcal M }}$; see Table 1. Similarly, for the combination with R_{Fe ii} , the reduction in uncertainty is by almost a factor of three. It may seem that the MCMC inference using the maximization of the likelihood function defined by Equation (6) is more robust for constraining individual parameters even for general priors of uniformly distributed coefficients in the studied linear combinations. However, the likelihood function defined by Equation (6) considers only the measurement errors of ${\tau }_{i}$. When we take into account that the measurement errors ${\sigma }_{i}^{\tau }$ could generally be underestimated by a factor f, then the likelihood function takes the following form:

$$\begin{eqnarray}&&{{ \mathcal L }}_{{\rm{f}}}=-\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}\left[\displaystyle \frac{{({\tau }_{i}-{\hat{\tau }}_{i})}^{2}}{{s}_{i}^{2}}+\mathrm{log}{s}_{i}^{2}+\mathrm{log}(2\pi )\right],\end{eqnarray} \tag{ 7 }$$

where ${s}_{i}=f{\sigma }_{i}^{\tau }$. Maximizing the functional ${{ \mathcal L }}_{{\rm{f}}}$ that depends on the additional underestimation parameter f yields the posterior 2D and 1D marginalized distributions as shown in Figure 7. The parameter uncertainties are now comparable to those inferred using the OLS method; see Table 1, where we list the parameters inferred using the likelihood function ${{ \mathcal L }}_{{\rm{f}}}$ under the abbreviation MCMC(f). The difference between uncertainties inferred using either ${ \mathcal L }$ or ${{ \mathcal L }}_{{\rm{f}}}$ shows the importance of the defined likelihood function for a given problem. In our case, the OLS and MCMC methods are in better agreement in terms of the uncertainties when the factor f is included, which suggests that the measurement time delay uncertainties are underestimated by a factor of about two or three.

Zoom In Zoom Out Reset image size

Since OLS and MCMC can be considered as independent inference techniques, their overall consistency confirms the statistical robustness of our results, mainly concerning the low scatter for highly accreting sources.

The R–L relation is of key importance for black hole mass measurements in the sources that were not studied through reverberation mapping. It is also a promising tool to be applied in cosmology if we have a large enough sample showing a small scatter around the best-fit relation. In the present paper, we studied the R–L relation based on the Mg ii line, using all available data, and we focused on the selection of additional parameters/methods that help to reduce the observed scatter.

Considering a virialized and photoionized gas, it is expected that log ${\tau }_{\mathrm{obs}}\propto \alpha \mathrm{log}L$, where the slope of the luminosity is given by $\alpha =0.5$. In the optical range, the Hβ reverberation mapping results give a slope for the optical luminosity at 5100 Å of $\alpha ={0.533}_{-0.033}^{+0.035}$ in the most accepted R–L relation (Bentz et al. 2013). This slope was nicely consistent with simple predictions of the BLR location based on a fixed ionization parameter (see Czerny et al. 2019b and the references therein) or the failed radiatively accelerated dusty outflow (FRADO) model (Czerny & Hryniewicz 2011).

The slope of the R–L relation based on the Mg ii line shows a large diversity (Vestergaard & Osmer 2009; Trakhtenbrot & Netzer 2012; Homayouni et al. 2020; Zajaček et al. 2020). The most recent Mg ii monitoring from the SDSS-RM project (Homayouni et al. 2020), which increased significantly the number of Mg ii time lags, provides an R–L relation with a shallower slope than that seen for the Hβ line by Bentz et al. (2013). Based on their time delay significance criteria, the SDSS-RM sample (Homayouni et al. 2020) is divided into two subsamples (significant and gold sample), where the slope of the luminosity in both cases is 0.22 ± 0.06 and ${0.31}_{-0.10}^{+0.09}$, respectively.

In this work we combine the data of Homayouni et al. (2020) with the low- and high-luminosity sources collected by Zajaček et al. (2020). When all the measurements of Mg ii delay up to now are included, we get a slope of 0.298 ± 0.047, which is in agreement with the one inferred by Homayouni et al. (2020). However, when the sample is divided based on the median $\dot{\,{ \mathcal M }}$ intensity, the slope becomes much steeper. In the case of the low-accretion subsample, the slope is very close (0.520 ± 0.078) to the expected value and the one provided by Hβ results. As we mentioned previously, most of the first Hβ reverberation-mapped objects were selected based on their high variability (high F_var) and strong [O iii] $\lambda \lambda 5007,4959$ emission, which indicates a low accretion rate and is expected to show an agreement with the standard slope. In the case of the high-$\dot{\,{ \mathcal M }}$ subsample, the slope is 0.414 ± 0.058, which, within uncertainties, is very close to the expected value. A large sample of highly accreting objects is needed to confirm a real deviation of the slope of the luminosity for this kind of object.

When additional independent observational parameters are considered in linear combinations with the luminosity at 3000 Å, the slope shows different values. Since the typical slope of 0.5 is estimated considering only the luminosity, we cannot claim a deviation from it in the other studied cases where R_{Fe ii} , F_var, and FWHM are used, because new theoretical models taking into account these parameters should be considered.

Newer, larger samples for the Hβ line brought an additional scatter in comparison with the sample of Bentz et al. (2013), and also the scatter in the Mg ii full sample is relatively large, ${\sigma }_{\mathrm{rms}}\sim 0.3$ dex. However, a more advanced approach presented in this paper helped to reduce this scatter considerably. The smallest scatter has been achieved when interdependent quantities ($\dot{\,{ \mathcal M }}$ and ${L}_{\mathrm{bol}}/{L}_{\mathrm{Edd}}$) are used. These results, thus, have to be treated with care since they are intrinsically correlated with the time delay. On one hand, it is an attractive hypothesis that the scatter in the original R–L relation (i.e., log ${\tau }_{\mathrm{obs}}$ vs. log L₄₄) is due to the spread in the accretion rate intensity. The underlying mechanism for shortening the time delay with an increase in the accretion rate has already been suggested by Wang et al. (2014c), who argued that the self-shielding of the accretion disk also leads to the selective shielding of the BLR and its division into two distinctly different regions. In addition, within the FRADO model combined with the shielding effect, such a trend is expected (Naddaf et al. 2020).

The use of the interdependent quantities has a drawback, that, despite the small scatter in the final plot, the recovery of the values predicted by the relation comes with a large error. If we want to use the linear combination with $\dot{{ \mathcal M }}$ (left panel in Figure 2; see also Table 1, rows 7 and 8) to predict the value of the time delay for a given source, if we measure the FWHM and the L₄₄ assuming a known cosmology, the error of this prediction will become larger than the scatter visible in the plot: the minimum value of the error of $\mathrm{log}{\tau }_{\mathrm{obs}}$ in this case would be $\delta {K}_{3}/0.136=0.147$, where $\delta {K}_{3}$ is the error of the coefficient K₃ in Table 1, row 7. The same will happen if we use the linear combination with $\dot{{ \mathcal M }}$ (left panel of Figure 2) to obtain the absolute luminosity from the measured time delay and FWHM—the minimum error of the predicted $\mathrm{log}$ L₄₄ would then be $\delta {K}_{3}/0.046=0.43$. Thus, for the black hole mass measurements or for the cosmological applications, the use of the relations based on independent quantities still gives much better results. An attractive possibility would be to use the accretion rate or the Eddington ratio as independent parameters. This would require an independent measurement of the black hole mass, e.g., from the broadband SED fitting (Capellupo et al. 2015).

The use of the MCMC fitting method with the interdependent $\dot{\,{ \mathcal M }}$ parameter reduces the error for all the coefficients, so the numbers mentioned above will be formally up to 3 times lower. However, the dispersion around the fit is not affected by the inference method, so it is not clear whether the use of MCMC indeed reduces the errors when predictions are made. As we also showed in Section 3.4, the uncertainty intervals of the posterior parameter distributions also depend on whether the underestimation parameter is included in the likelihood function or not. If it is included, then the uncertainties are consistent with the OLS method, which suggests that the time delay uncertainties are generally underestimated by about a factor of two.

The dominant role of the dimensionless accretion rate or the Eddington ratio is supported by the fact that the division of the sample into two parts representing low and high accretion rates also reduced the scatter in the R–L relation considerably, particularly for the case of the highly accreting subsample. In combination with the measurement of R_{Fe ii} , this reduced the scatter in the R–L relation down to ∼0.17 dex. The scatter for the lower accretion rate subsample remained at ${\sigma }_{\mathrm{rms}}\,\sim 0.25$ dex. This level of scatter is most likely related to the red-noise character of AGN variability in the optical band (e.g., Czerny et al. 1999; Kelly et al. 2009; Kozłowski et al. 2010; Kozłowski 2016). A relatively short monitoring allows one to determine the time delay of the lines with respect to the continuum, but years of monitoring are needed to determine the mean luminosity level, instead of a part of the light curve catching the source in a relatively high or a relatively low state. As was shown by Ai et al. (2010) for SDSS Stripe 82 AGNs, higher Eddington ratio sources vary less in the optical/UV bands. This variability, in particular for a low-accretion subsample, may lead to an irreducible scatter in the R–L relation, as is clearly seen for our results (Figures 4 and 5). The same scatter was discussed by Risaliti & Lusso (2019) in the context of the broadband UV–X-ray relation, where they argue that the variability is relatively unimportant for high-redshift quasars, leading to the scatter of 0.04 dex. However, their selection of predominantly blue quasars contributed to the reduction of this scatter. During the 16 yr quasar monitoring, the variability varied from 0.04 to 0.1 dex, depending on the quasar absolute luminosity (Hook et al. 1994). In our sample, no preselection of objects based on the UV slope has been made.

Alternatively, the reduction in the scatter in the subsample with high $\dot{\,{ \mathcal M }}$ could be related to the fact that sources radiating close to their Eddington limit saturate toward a limiting value, which leads to the stabilizing of the ratio between the luminosity and black hole mass (which is basically $\dot{\,{ \mathcal M }}$ or ${L}_{\mathrm{bol}}/{L}_{\mathrm{Edd}}$), making the sources steadier (Marziani & Sulentic 2014). Other explanations of the additional scatter include the spin effect and the possibility of a retrograde accretion (Wang et al. 2014b; Czerny et al. 2019b).

We stress the fact that R_{Fe ii} and F_var show the smallest scatter (0.17 and 0.19 dex, respectively) when the sample is divided into the two considering the $\dot{\,{ \mathcal M }}$ intensity. Both observational properties are correlated with $\dot{\,{ \mathcal M }}$ and ${L}_{\mathrm{bol}}/{L}_{\mathrm{Edd}}$ (Marziani et al. 2003; Wilhite et al. 2008; Dong et al. 2009, 2011; MacLeod et al. 2010; Sánchez-Sáez et al. 2018; Martínez-Aldama et al. 2019; Du & Wang 2019; Yu et al. 2020b), but they are independent. This result suggests that the accretion rate drives the scatter in the R–L relation.

The division based on the $\dot{\,{ \mathcal M }}$ for the current sample could be affected by inclusion of newer sources or reanalyses of existing ones. As a check for completeness, we compare the distribution of $\dot{\,{ \mathcal M }}$ estimated for two large SDSS quasar catalogs—for the DR7 release (Shen et al. 2011, hereafter S11) and for a more recent DR14 release (Rakshit et al. 2020, hereafter R20), with various spectral parameters estimated for 105,783 and 526,265 sources, respectively. The $\dot{\,{ \mathcal M }}$ formalism used in this paper (see Equation (4)) is a function of the black hole mass and the monochromatic luminosity at 3000 Å (the associated inclination term, cos ${\text{}}\theta$, is set to 0.75, which is the mean disk inclination for type 1 AGNs). We filter the catalogs first by limiting to the values that are reported to be positive and nonzero. Additionally, in the latter case (R20), the authors also provide quality flags for selected parameters including the monochromatic luminosity at 3000 Å (L₃₀₀₀). For the DR7 QSO catalog, no such quality flags were provided; thus, we use the full sample in this case. For the black holes masses, we use three variants common to the two catalogs: (a) from Vestergaard & Osmer (2009, hereafter VO09), (b) from S11, and (c) from the fiducial virial black hole mass values calculated based on (a) the Hβ line (for z > 0.8) using the calibration of Vestergaard & Peterson (2006), (b) the Mg ii line (for 0.8 $\leqslant$ z >1.9) using the calibration provided by VO09, and (c) the C iv line (for z $\geqslant$ 1.9) using the Vestergaard & Peterson (2006) calibration. In the case of the black hole masses, the DR14 QSO catalog provides the quality flag only for the fiducial masses and is unavailable for the VO09 and S11 mass estimates. We thus use the quality flags to control the sample for the fiducial mass estimates in this case. Figure 8 demonstrates the $\dot{\,{ \mathcal M }}$ distributions computed for the two catalogs. The three panels are synonymous to the three cases of black hole mass estimates incorporated to estimate the $\dot{\,{ \mathcal M }}$ values.

Zoom In Zoom Out Reset image size

Due to the quality control and filtering, the source sample drops to ∼79% (DR7) and to ∼67% (DR14) of their respective original source counts. The effective numbers of sources per case of the black hole mass estimates remain almost alike. ³ To predict the variation with respect to our small sample of Mg ii RM-reported sources, we extract the mean (μ) and the standard deviation (σ) from each $\dot{\,{ \mathcal M }}$ distribution (in log scale) shown in Figure 8 using simple Gaussian fits. The ($\mu \pm \sigma$) values for each panels are (VO09) 0.29 ± 0.72 (DR7) and –0.14 ± 0.95 (DR14), (S11) –0.02 ± 0.71 (DR7) and –0.17 ± 0.93 (DR14), and (fiducial) –0.01 ± 0.74 (DR7) and 0.05 ± 0.81 (DR14). The median value for our sample, log $\dot{\,{ \mathcal M }}$ = 0.2167, is well within 1σ limits regardless of the distribution taken from the larger catalogs, and hence it will not have significant effects on the correlations quoted in this paper with the inclusion of more sources in the future.