Current and future applications of Reverberation-mapped quasars in Cosmology

Reverberation mapping technique is an important milestone in AGN demographics, kinematics and the structure of the Broad Line Region (BLR) based on the time-delay response between the continuum and emission line. The time delay is directly related to the size of BLR which is related to the continuum luminosity of the source, producing the well-known Radius-Luminosity (RL) relation. The majority of sources have been monitored for their H$\beta$ emission line in low redshift sources ($z<0.1$), while there are some attempts using the MgII line for higher redshifts. We present a recent MgII monitoring for the quasar CTS C30.10 ($z=0.90$) observed with the 10-meter SALT, for which RL scaling based on MgII holds within measurement and time-delay uncertainties. One of the most important advantages of reverberation mapping technique is the independent determination to the distant source, and considering the large range of redshifts and luminosities found in AGNs, their application to cosmology is promising. However, recently it has been found that highly accreting sources show the time delays shorter than expected from RL relation. We propose correction for this effect using a sample of 117 H$\beta$ reverberating-mapped AGNs with $0.02<z<0.9$, recovering the low scatter along with the relation. We are able to determine the cosmological constants, $\Omega_m$ and $\Omega_\Lambda$ within 2$\sigma$ confidence level. We present the first steps in modelling of light curves for H$\beta$ and MgII and discuss the quasar selection in the context of photometric reverberation mapping with Large Synoptic Survey Telescope (LSST). With the onset of the LSST era, we expect a huge rise in the overall quasar counts and redshift range covered ($z\lesssim 7.0$), which will provide a better constraint of AGN properties with cosmological purposes.


INTRODUCTION
Reverberation mapping (hereafter RM) is a robust observational technique to determine the time-lag τ 0 between the variability of the ionizing continuum of an active galactic nucleus (AGN) and the line emission associated with the broad-line region (BLR, Blandford and McKee, 1982). This is only possible since the emission-line flux densities are highly correlated with the continuum flux, which implies that the AGN continuum related to the power generated by the accretion disk is the main source of photoionization. The main source of recombination emission lines is the BLR material, that is optically thick with respect to the ionizing UV/optical continuum.
The second application combines τ 0 and the line width of the broad emission line. It can be considered that the large line width of the BLR gas ∆v FWHM 2 of several thousand km/s, is due to the Doppler broadening of clouds that are gravitationally bound to the supermassive black hole (SMBH). The 3D Keplerian velocity of the gas v BLR = f vir ∆v FWHM can be inferred, though it depends on the rather uncertain virial factor f vir . By combining these two pieces of information, one can estimate the virial black hole mass, where the virial factor f vir is of the order of unity and depends on the overall geometrical distribution of the BLR clouds, their kinematics as well as their line of sight emission properties, c is the velocity of light and G is the Gravitational constant. The virial factor mainly depends on the inclination i of the BLR plane with respect to the observer and the thickness of the BLR, H BLR /R BLR , and can be expressed as (Collin et al., 2006;Mejía-Restrepo et al., 2018;Panda et al., 2019), Therefore, the virial factor dependence on the viewing angle of the BLR and its geometrical properties, introduces an overall uncertainty in the virial mass determination. Fixing the virial factor to a constant value, which is frequently applied to single-epoch measurements (Woo et al., 2015), can lead to a virial-mass difference of a factor of 2-3 (Mejía- . Comparing the black hole mass estimations determined from the spectral energy distribution (SED) fitting of the AGN continuum, with the viral black hole mass, Mejía-Restrepo et al. (2018) found that the virial factor is inversely proportional to the FWHM of the broad lines, f vir ∝ FWHM −1 line . This could be interpreted as the effect of the BLR viewing angle (inclination) or as the effect of the radiation pressure on the BLR distribution. The strong effect of the BLR inclination directly results from the plane-like geometry of the BLR, in which the BLR consists of cloudlets that have an overall flattened geometry (Gaskell, 2009). This "nest-like" model of the BLR is also consistent with the first direct velocity-resolved observation of the BLR in 3C273 (Gravity Collaboration et al., 2018) using the Very Large Telescope Interferometer GRAVITY at the European Southern Observatory (Gravity Collaboration et al., 2017). The flattened geometry of the BLR that follows the disc geometry could be explained by the formation of the BLR clouds from the disc material beyond the dust sublimation radius (Czerny and Hryniewicz, 2011). The radiation pressure acting on the dust in the BLR clouds then leads to their lift-off from the disc plane and subsequent fall-back when the dust evaporates. This is one possible scenario for the origin of the low-ionization line component of the BLR (LIL), so-called failed radiatively accelerated disc outflow (FRADO, Czerny et al., 2017), while the high-ionization line (HIL) component seems to be associated with nuclear outflows (Collin-Souffrin et al., 1988).
Finally, the third application of the RM is the power-law radius-luminosity relation R BLR ∝ L α AGN between the BLR radius (or time-delay) and the AGN monochromatic optical luminosity. Initially, Kaspi et al. (2005) found the slope α = 0.67 ± 0.05 between the optical monochromatic luminosity and the broad Hβ line, which implied a deviation from simple theoretical models predicting a slope of 0.5. The relation was further modified for lower-redshift sources using the Hβ RM and the luminosity at 5100Å (Bentz et al., 2006(Bentz et al., , 2009(Bentz et al., , 2013 taking into account the host galaxy starlight. After the removal of the host-galaxy starlight from the total luminosity, the power-law slope was determined to be α = 0.533 +0.035 −0.033 (Bentz et al., 2013), which is consistent with the theoretical expectation from simple photoionization models. The slope of α = 1/2 simply follows from the ionization parameter of a BLR cloud, U = Q(H)/(4πR 2 cn H ), where R is the cloud distance from an ionizing continuum, c is the speed of light, n H is the number density, and Q(H) = ∞ ν i L ν /hν dν is the ionizing photon flux. In general, one can assume that the ionization conditions as well as the gas density in the BLR are comparable for all AGN, i.e. the product U n H is constant, from which follows R BLR ∝ Q(H) 1/2 ∝ L 1/2 (Netzer, 2013).
The importance of the R Hβ − L 5100 relationship lies in the fact that it is the basis of secondary black-hole mass estimates, especially for high-z AGN (Laor, 1998;Wandel et al., 1999;McLure and Jarvis, 2002;Vestergaard and Peterson, 2006). In principle, by knowing the optical luminosity of the source, one can infer the black-hole mass from a single-epoch spectrum (by extracting the broad-line FWHM), which makes it especially useful for large statistical surveys of sources throughout the cosmic history. Given the importance of R Hβ − L 5100 relation, it is necessary to verify it using different emission lines from the LIL BLR region and higher-redshift sources, which we will analyze in subsequent section for the case of MgII line and the redshift of 0.9 (a complete and updated MgII based RM sample is reported in Zajaček et al. 2019, submitted to ApJ).
RM of quasars, using both the BLR reverberation technique and the accretion disc radiation reprocessing, was suggested to be utilized in cosmological studies (Collier et al., 1999;Elvis and Karovska, 2002;Horne et al., 2003). In particular, the R Hβ − L 5100 relation directly enables to turn reverberation-mapped quasars into "standardizable" luminosity candles (Haas et al., 2011;Watson et al., 2011;Czerny et al., 2013;Bentz et al., 2013). Under the assumption that the R Hβ − L 5100 power-law relationship applies not only to Hβ line but also other LIL lines, one could apply it to intermediate-and high-redshift quasars whose co-moving time-delay was determined to infer their absolute luminosities. From measured monochromatic flux densities of the sources, one can infer the luminosity distance D L = L 5100 /4πF and the Hubble diagram of reverberation-mapped quasars. The critical point here is to reduce the scatter along the radiusluminosity relation, both by improving the RM time-lag precision and by correctly subtracting the host starlight (Haas et al., 2011). In both panels, colors indicate the variation in dimensionless accretion rate,Ṁ c .
The paper is structured as follows. In Section 2, we report on the detection of MgII line time-delay for the intermediate-redshift quasar CTS C30.10 and its consistency with the radius-luminosity relation previously derived for Hβ emission in lower redshift sources. We continue with Section 3 to analyze in detail how accretion rate affects the position of sources along the R Hβ − L 5100 relation. We summarize the initial cosmological constraints as derived from the current reverberation-mapped quasar sample in Section 4. First steps of modelling stochastic lightcurves in the context of reverberation mapping with the Large Synoptic Survey Telescope (LSST) are analyzed highlighting a novel filtering algorithm to select higher quality photometric data in Section 5. We discuss a few open problems in Section 6 and finally conclude with Section 7. extend the redshift range to intermediate values. MgII line with the rest wavelength ∼ 2798Å is suitable for the reverberation-mapping studies in the redshift range z = 0.4 − 1.5 for ground based telescopes, which is also utilized for the sample of seven objects monitored by the South African Large Telescope (SALT) telescope (Czerny et al., 2013).
Between December 6, 2012 and December 10, 2018, we monitored the quasar CTS C30.10 using the 10-meter SALT telescope in order to detect the time-lag of the LIL emission-line of MgII with respect to the quasar continuum. The hypothesis was to probe the consistency of the radius-monochromatic luminosity relation towards higher redshifts and higher absolute luminosities. The quasar CTS C30.10 was discovered as a bright source in the Calan-Tololo survey (Maza et al., 1993)  The MgII emission-line light curve was constructed based on the long-slit spectroscopy mode of the SALT telescope with the slit width of 2" (for details, see Modzelewska et al., 2014;Czerny et al., 2019). The MgII line was extracted from the spectral energy distribution of CTS C30.10 in the wavelength range 2700 − 2900Å by subtracting the two spectral components, namely the power-law component due to the accretion disc thermal emission, and FeII -line pseudo-continuum. The MgII broad line was found to consist of two kinematic components -redshifted and blueshifted Lorentzian profiles 3 -each of which further consists of a doublet at 2796.35Å and 2803.53Å with the doublet ratio of 1.6. The equivalent width (EW) of the MgII line, EW∼ F line /F cont , was calculated by numerical integration of all four spectral components -2 kinematic and 2 doublet components. Because of the two kinematic components of MgII line, CTS C30.10 belongs to type B quasars (Sulentic et al., 2007;Modzelewska et al., 2014) characterized by lower Eddington ratios in the range log (L/L Edd ) = 0.01 − 0.2 (Sulentic et al., 2011), while type A sources exhibit a single component MgII line of a Lorentzian shape (Laor et al., 1997;Véron-Cetty et al., 2001;Sulentic et al., 2002;Zamfir et al., 2010;Sulentic et al., 2009Sulentic et al., , 2011Shapovalova et al., 2012). The origin of two components is still quite uncertain although the presence of the second, asymmetric component and the overall FWHM is consistent with our source belonging to Pop. B. It could either imply the presence of second emission region due to absorption or scattering or it could hint at the origin of LIL lines close to the disc plane, resembling thus the disc kinematics, which would naturally lead to two components when viewed off-axis.
The continuum light curve was obtained from the OGLE-IV survey in the V -band and SALTICAM in the g-band. The SALTICAM observations were shifted freely to match the overlapping V-band values. Since the SALT observations are not spectrophotometric, MgII flux density was obtained using its EW and the continuum flux density.
The normalized continuum dispersion was 6.0%, while the line dispersion was 5.2%, which is slightly lower but comparable within uncertainties. This is consistent with the simple reprocessing scenario, where the central disc provides all the ionizing UV photons which are absorbed and scattered by BLR clouds at larger distances. This allowed us to use both the continuum and the line-emission light curves to infer the time-delay τ MgII of the MgII broad line. Cross-correlation centroid distribution. Right panel: Cross-correlation peak distribution.

Time-lag determination of MgII line using cross-correlation function
First, we used the standard interpolated cross-correlation function (ICCF) to infer the time-delay between the MgII line-emission and the V-band continuum, which corresponds to the continuum around the redshifted MgII line. The ICCT requires regularly sampled datasets, while realistic light curves are unevenly sampled. The regular timestep is achieved by interpolating the continuum light curve to timeshifted emission-line light curve or vice versa. Typically, both interpolations are averaged to obtain the symmetric ICCF. Given the two light curves x i and y i sampled at discrete time intervals t i (i = 1, ..., N ) with the regular time-step ∆t = t i+1 − t i , we can define the cross-correlation function (CCF) as, where x and y are light curve mean values, respectively, and τ k = k∆t (k = 0, ..., N − 1) is the time-shift of the second light curve with respect to the first one at which the CCF is evaluated.
We make use of the PYTHON code PYCCF (Sun et al., 2018), which is the implementation of earlier CCF analysis by Gaskell and Peterson (1987) and Peterson et al. (1998). We calculate CCF separately for interpolated continuum, interpolated emission-line curves as well as the symmetric case. Using thousand Monte Carlo runs of the combined random subset selection (RSS) and flux randomization (FR), we obtain distributions of the CCF centroid (CCFC) and the CCF peak (CCFP), from which we can calculate the corresponding uncertainties, see Fig. 2.
In Table 1, we summarize the time-delay results for all considered cases -the time-delay values are expressed with respect to the observer's frame (τ observer = (1 + z) × τ rest ). For the symmetric interpolation, the centroid time-delay is τ ICCF = 1060 +262 −457 days. In general, the cross-correlation centroid and peak distributions have three peaks, see Fig. 2, while the peak at ∼ 1060 is the most prominent. At this value, ICCF also reaches the largest value for the symmetric and interpolated-continuum cases, see Fig. 2 (left panel), while the interpolated-line case peaks at larger time-delays. This offset may be due to the fact that the line-emission dataset is not so densely covered as the continuum light curve, which leads to artefacts when interpolating it to denser continuum light curve.
Second, we applied the discrete correlation function (DCF), which does not interpolate between the light curves and is thus better suited for unevenly sampled datasets with known measurement errors. The DCF was described by Edelson and Krolik (1988) and applied routinely to search for light curve correlations and time-lags. The first step in the DCF analysis is to look for data pairs (x i , y j ) that fall into the time-delay bin of τ − δτ /2 ≤ δt ji ≤ τ + δτ /2, where τ is the given time-delay, δτ is the time-delay bin and δt ji = t j − t i . Given M such data pairs, one calculates a corresponding number of unbinned discrete correlation coefficients U DCF ij in the following way, where x, y are the light curve means for a given time-delay bin. Other parameters s x and s y stand for the variances, and σ 2 x and σ 2 y are the mean measurement errors for a given time-delay bin. Subsequently, the DCF coefficient can be calculated for a given time-lag bin by averaging in total M U DCF M ij values, where, M is the number of light curve pairs that fall into a particular time-lag bin. The uncertainty can be estimated by the relation, We applied the PYTHON code of Damien Robertson (see Robertson et al., 2015, for the code implementation), which calculates the DCF with the possibility of Gaussian weighting for the matching pairs of light curve points. The code allows us to select a different size for equal time-bins as well as the searched interval for the time-delay. We expand the possibilities of the code by adding the bootstrap technique to assess the significance of individual DCF peaks and to better estimate their uncertainties. The DCF uses equal time-step binning and it is thus quite sensitive with respect to the time-bin size. We test this using three different timesteps -50, 70, and 100 days -which leads to the decrease in the mean DCF values for the time-delay at ∼ 1275 days, which is the most prominent for the smaller time-bins but significantly gets smaller for larger time-bins, see Fig. 3 (left panel). For the time-bin size of 100 days, the largest DCF of 0.62 is for 1050 days, followed by the peak at 550 days. Moreover, we verify the significance of individual peaks by running 1000 bootstrap realizations where we randomly create subsamples from both light curves simultaneously. We obtain the distribution of time-delay peaks (for which the DCF value is the largest for each run) in Fig. 3 (right panel), from which we obtain the peak time-delay with upper and lower 1σ uncertainties, τ DCF = 1030 +106 −140 days, which we also include in Table 1.
Finally, several biases and problems of the cross-correlation function which primarily stand out when dealing with sparse and heterogeneous light curves were taken into account in z-transformed DCF (Alexander, 1997). In comparison with the discrete correlation function (Edelson and Krolik, 1988), which bins data pairs into equal time-bins, z-transformed DCF applies equal population binning. It is, therefore, a more suitable and robust method for sparse, irregularly sampled, and heterogeneous pairs of light curves with as few as 12 points per population bin.
We show the calculated zDCF as a function of time delay in the observer frame in Fig. 4. The largest zDCF is for the time-delay at τ = 1050 days (zDCF= 0.76), followed by the smaller peak at τ = 571 days (zDCF= 0.68). We verify the significance of the peak using the maximum likelihood, which also enables us to estimate its uncertainty. We obtain the peak value of τ zDCF = 1050 +32 −27 days, which is also listed among other methods in Table 1.
We summarize the time-delay value for CTS C30.10 by calculating the average value for the most prominent peak at ∼ 1050 days in the observer frame using Table 1. We get the mean value of τ obs = 1048 +94 −156 days in the observer frame, which translates into τ source = 551 +49 −82 days in the source frame.

Radius-luminosity relation for MgII
As for Hβ broad line, we construct the radius-luminosity relation for MgII line taking into account our detection of time-lag for MgII in CTS C30.10 as well as the measurements of other sources that also have RM data of MgII: 6 sources from Shen et al. (2016), CTS252 from Lira et al. (2018), and NGC4151 from Metzroth et al. (2006). For an overview of the sources and their characteristics, see Table 3 in Czerny et al. (2019).
width, CTS252 is a higher Eddington-ratio source in comparison with CTS C30.10 , but the trend needs to be confirmed for a higher number of MgII measurements.
Excluding CTS252, we fit the radius-luminosity dataset with the general relation log(R BLR /1lt-day) = K + α log (λL λ /10 44 erg s −1 ). We obtain the best-fit parameters of K = 1.47 ± 0.06 and α = 0.59 ± 0.06. In Fig. 5, we see that the best-fit radius-luminosity relation passes nearly through the CTS C30.10 point in comparison with the Bentz relation, which is due to a slighly higher mean slope -α = 0.59 instead of α = 0.53. However, the slopes do not differ within uncertainties and hence the radius-luminosity relation with the expected scaling of R ∼ L 1/2 so far holds for MgII time-lag measurements.

ACCRETION RATE EFFECT OVER THE RADIUS-LUMINOSITY RELATION
The R Hβ − L 5100 relation offers the possibility of estimate the luminosity distance (D L ) independently of the redshift, which is suitable for cosmological applications (Martínez-Aldama et al., 2019a, and references therein). Also, it shows a low scatter (0.13dex, Bentz et al., 2013) and although objects, like NGC5548, show a large variation, the scatter is within uncertainties. However, the recent inclusion of AGN radiating close to the Eddington limit (super-Eddington sources) has increased the scatter significantly. This kind of sources are located on the right-bottom side of the relation (See Figure 1), which implies that the time delay of super-Eddington sources is shorter than the predicted by the R Hβ − L 5100 relation.
Super-Eddington sources tend to show an extreme behavior respect to the general AGN population: large densities (n H ∼ 10 13 cm −3 ), low-ionization parameters (log U < −2), high intensity of lowest ionization emission lines like FeII and strong outflows in the high-ionization like CIVλ1549  and references therein). These features are probably explained by an optically and geometrically thick slim disk (Wang et al., 2014b), which shields the line emitting region gas from the most intense UV radiation (Marziani et al., 2018). This peculiar behavior is also reflected in the R Hβ − L 5100 relation by the super-Eddington sources of SEAMBH sample, but also by some objects of SDSS-RM sample. One-third of the SDSS-RM shows an Eddington ratio higher than the average Eddington ratio in Bentz sample, then they also a departure from the expected R Hβ − L 5100 relation. Therefore, a correction is needed not just for super-Eddington sources, but all sources have to be rescaled.
To estimate this correction to the measured time delay as a function of the accretion rate, we estimate the black hole mass (Equation 1) considering a virial factor anti-correlated with the FWHM of the line recently proposed by Mejía-Restrepo et al. (2018), which in some sense corrects for the orientation effect. For the accretion rate, we will consider the dimensionless accretion rate (Ṁ c ) introduced by Du  , which is the difference between the time delay observed and estimated from the R Hβ − L 5100 relation. In the bottom panel of Fig. 1 is shown the behavior of ∆R Hβ as a function of L 5100 . Also, in colors, it is shown the variation of the dimensionless accretion rate, where it is clearly observed that largest departures correspond to highesṫ M c values. The relation between ∆R Hβ andṀ c can be described by the linear relation: ∆R Hβ,Ṁ c = (−0.283 ± 0.017) logṀ c + (−0.228 ± 0.016) .
With this relation, the observed time delay can be corrected by the dimensionless accretion rate effect, decrease the scatter and recover the time delay predicted by the R Hβ − L 5100 relation. However, this relation is strongly dependent on the virial factor selected for the black hole mass and the accretion rate estimations. The virial factor is still an open problem and although many formalisms have been proposed anyone can be applied to the general AGN population. f c BLR was modeling with a scarcity of narrow profiles (FWHM<2000 km s −1 ), which represent 30% of our full sample. The inclusion of narrow profiles in its modeling could modify the exponent of the anti-correlation and modify our results.

COSMOLOGY WITH REVERBERATION-MAPPED SOURCES
Recovering the low scatter of the R Hβ − L 5100 relation after the correction by the dimensionless accretion rate, we were able to build a Hubble diagram (Fig. 6) which relates the distance to the sources with the redshift (z) or the velocity recession. The slope of the relation between these two parameters is the value of Looking for the best fit using a minimization method (χ 2 ), we get the best values for Ω m and Ω Λ . Our results are in agreement with the standard model within 2σ confidence level, however, it is not yet ready to provide new for cosmological constraints, considering that error associated with CMB, Cepheids stars and SNIa are less than 2%.
There are still some systematic errors related that have to be corrected. In the previous section, we proposed a correction by the accretion rate effect, but this correction has to be improved since the virial factor has associated large uncertainties. Then, we have to work for clarifying the dynamics of the BLR and the orientation effect. Another important source of error is the method employed in the determination of the time delay. The delay measurements used here are a compilation of the results from the literature, and each group uses different criteria and methods, which also introduces an error when different samples are compared. In particular, Czerny et al. (2019) and Zajaček et al. (2019) applied six different methods to determine the time-lag between the continuum and MgII line response. These included cross-correlation based techniques ICCF, DCF, zDCF, χ 2 -based method (Czerny et al., 2013, the measure of data regularity (von Neumann's estimator; Chelouche et al., 2017), and the JAVELIN code (Just Another Vehicle for Estimating Lags In Nuclei, formerly SPEAR; Zu et al., 2011Zu et al., , 2013Zu et al., , 2016. They showed that for the intermediateredshift quasar CTS C30.10 (see also Section 2.1), the uncertainty in the time-lag determination depends on the approach. If a narrow surrounding of the main peak in time delay measure is taken into consideration, the uncertainty can be of the order of a few days only for the total observing run of 1000 days and more. The small uncertainty of the order of a few days only was obtained by JAVELIN. On the other hand, if the total range of possible time delays is considered, that is between zero and half of the duration of the observation, the uncertainty can be of the order of 10 and even 100 days depending on the quality of the data. This uncertainty further propagates into the radius-luminosity relation and the Hubble diagram of reverberation mapped quasars. In Fig. 7, we illustrate this in terms of constraining cosmological parameters from monitoring of a single source CTS C30.10 mentioned above (dark matter and dark energy fractions) for two different uncertainties -2 days as provided by the JAVELIN code and 22 days that is more representative when other methods are taken into account.
Although the error provided by the JAVELIN is underestimated and can be corrected by performing many bootstrap runs to recover the overall time-lag distribution, by combining several robust methods one can achieve the time-lag uncertainty of 5%. Under this assumption, approximately 625 sources are needed to achieve 2% uncertainty for all cosmological parameters including w 0 parameter. However, special attention is needed for the process of a systematic time-lag analysis for the whole sample of reverberation-mapped quasars since currently, the time-lag analysis is largely heterogeneous with different authors using different time-delay methods and their associated uncertainties. This is a source of large systematic uncertainties that will need to be addressed in future RM samples of quasars.
Quasars are complex objects and many corrections are still needed to consider these objects like standard candles suitable for cosmology. However, Cepheids star and SNIa have been also corrected for the shape of the lightcurve, extinction and the host galaxy, and sophisticated observational methods have been proposed for getting better observations (e.g. Scolnic et al., 2018). All this effort is reflected in the low uncertainties associated with the cosmological results. A similar situation should be happening with quasars taking advantage of their broad range of redshift and luminosity.

AGN Science with LSST
LSST will be a public optical/NIR survey of ∼half the sky in the ugrizy bands upto r∼27.5 based on ∼ 820 visits across all the six photometric bands over a 10-year period (Ivezić et al., 2019). The 8.4m telescope with a state-of-the-art 3.2 Gigapixel flat-focal array camera will allow performing rapid scans of the sky with 15 seconds exposure and thus providing a moving array of color images of objects that change. The whole observable sky is planned to be scanned every ∼4 nights. It is expected that upon the completion of the main-survey period, LSST will have mapped ∼20 billion galaxies and ∼ 17 billion stars using these six photometric bands 4 .
The LSST AGN survey will produce a high-purity sample of at least 10 million well-defined, opticallyselected AGNs (see Section 5.2). Utilizing the large sky coverage, depth, the six filters extending to 1µm, and the valuable temporal information of LSST, this AGN survey will supersede the largest current AGN samples by more than an order of magnitude. Each region of the LSST sky will receive ∼200 visits in each band in the decade long monitoring, allowing variability to be explored on timescales from minutes to a decade (see Chapter 10 of LSST Science Collaboration et al. (2009) for more details). LSST will conduct more intense observations of at least 4 Deep Drilling Fields or DDFs, i.e. COSMOS, XMM-LSS, W-CDF-S and ELAIS-S1, each of which cover ∼10 deg 2 (more details in Brandt et al. 2018). Brandt et al. (2018) estimates an increase in the depth in the DDFs by over an order of magnitude relative to the main-survey. This will be made possible owing to a ∼18-20X increase in the number of visits per photometric band.
LSST will be assisted by an array of campaigns (Brandt and Vito, 2017) during its proposed run. Moreover, there is a need for follow-up spectroscopic campaigns i.e. SDSS-V Black Hole Mapper (Kollmeier et al., 2017), which aim to derive BLR properties and reliable SMBH masses for distant AGNs with expected observed-frame reverberation lags of 10-1000 days. The SDSS-V Black Hole Mapper survey will also perform reverberation mapping campaigns in three out of the four LSST DDFs. LSST will provide additional high quality photometry for these sources that will substantially improve these estimations. The proposed cadence (∼2 nights) of the observations of AGNs and the expected re-visits during the 10-year run will allow to (1) ensure a high lag recovery fraction for the relatively short accretion-disk lags expected (Yu et al., 2018); and (2) allow investigations of secular evolution of these lags that test the underlying model. LSST will also perform quality accretion disk reverberation mapping for ∼3000 AGNs in the DDFs.

Suitability of quasar monitoring
Every night, LSST will monitor ∼75 million AGNs and is estimated to detect ∼300+ million AGNs in the ∼18000 deg 2 main-survey area (Luo et al., 2017). Taking into account factors like obscuration and host-galaxy contamination that will hinder this AGN selection, optimistic estimates predict around ∼20 million (from LSST alone) and ∼50+ million (LSST+ others) AGNs to be reliably monitored.   (Zheng et al., 1997) and the dashed box in purple highlights the importance of quasar selection. The spectrum is shown at a redshift, z = 2.56. The composite spectrum is downloaded from https://archive.stsci.edu/prepds/composite quasar/

Type-1 quasar counts
In the context of the reverberation mapping, one intends to focus on the broad emission lines that are a characteristics of Type-1 AGNs (Netzer and Peterson, 1997;Peterson et al., 2004;Gravity Collaboration et al., 2018). With LSST, we would want to constrain the quasar counts in terms of only Type 1 AGNs with respect to the observed broad emission lines such as Hβ, MgII, CIV. We perform a simple filtering to estimate these numbers in terms of total predicted number of quasars to be observed over the full duration of LSST.
Three entities affect our estimation of the "good" Type 1 quasars, namely (1)  that will be acquired from the pipeline, is already convolved with the filters and to recover the true intensity one needs to scale by 1 T , where T is the throughput (or transmittance) of the corresponding band which is estimated from the band function. We adopt a throughput value of 0.2 to exclude the gaps between the consecutive bands (see Figure 8). To get an idea of the impact of the line FWHM on the quasar number counts, we adopted the mean values for the broad emission lines -Hβ, MgII, from the SDSS DR7 Quasar catalogue (Shen et al., 2011). The average values for the FWHM are: 4662.26 ± 2181.53 km s −1 (MgII) and 4903.63 ± 3119.26 km s −1 (Hβ). For simplicity, we adopt the values 4500 km s −1 (MgII) and 5000 km s −1 (Hβ) for this preliminary analysis. Since, the emission lines in consideration are "broad", we estimate the equivalent width of MgII and Hβ to account for the broadness in these line profiles. We estimate this additional width from the mean FWHMs and the values obtained are ∼42Å and ∼80Å for MgII and Hβ respectively. Next, we retrieve the wavelength windows where these photometric bands overlap (see Table 2). As we have assumed the throughput value of 0.2, we exclude the u-band since the peak transmittance in this band is ∼0.14. This is the reason why the first entry in Table 2 is in the g-band (for MgII; rest-frame wavelength = 2799.12Å) and the wavelength is recovered corresponds to the intersection with the throughput value, T=0.2 on the lower wavelength side of the band. In the case for Hβ, the starting value z=0 is inclusive (rest-frame wavelength for Hβ = 4861.36Å). The remaining entries in the Table 2 are identical. We estimate that ∼26% (for MgII) and ∼36% (for Hβ) of the supposed Type 1 quasars will end up being observed in these overlapping windows taking into account the respective emission lines broadening and our naive assumption of the throughput cutoff. Figure 9 shows an illustrative view of this analysis in a more general context.
The expected number of Type 2 quasars is estimated to be ∼50-70% of the total population (Schmitt et al., 2001;Hao et al., 2005;Elitzur, 2012). So, if we start with a total population, for example, 1 million, we will have about 300,000-500,000 Type 1 AGNs based on this filtering. Next, the actual fraction of the "good" quasars will be ∼22-37% (MgII) and ∼19-32% (Hβ) after applying our filter due to the band overlap. It is expected that the filter wheel to be installed on LSST will have five slots (the filter wheel will prioritize griz and the y and u filters will be alternated for the fifth slot) complementary to the SDSS filter design (Fukugita et al., 1996). This means that at any given time, the image will be taken with one of the filters in-line with the camera. The filter replacement takes about 2 minutes to complete and this filter replacement will be performed 2-6 times per night and the total number of filter changes through the survey is 14,194 (LSST Science Collaboration et al., 2017). The remaining filter will be substituted for one of the existing filters in the filter-wheel (most likely u-band in place of y-band). This will affect the quasar counts based on the quality of the data, where quality is directly related to the number of re-visits for the source.

Modelling 'real' light curves
LSST is a photometric project but the 6-channel photometry (see Figure 8) can be effectively used for the purpose of reverberation mapping and estimation of time delays. We present some preliminary results from our software in development which allows to produce mock light curves and recover the time delays. The code takes into consideration several key parameters to produce these light curves, namely -(1) the campaign duration of the instrument (10 years); (2) number of visits per photometric band; (3) the photometric accuracy (0.01-0.1 mag) 5 ; (4) black hole mass distribution 6 ; (5) luminosity distribution 7 ; (6) redshift distribution 8 ; and (7) line equivalent widths (EWs) consistent with SDSS quasar catalogue (Shen et al., 2011). We create continuum stochastic lightcurve for a quasar of an assumed magnitude and redshift from AGN power spectrum with Timmer-Koenig algorithm (Timmer and Koenig, 1995). The code takes as an input a first estimate for the time delay measurement. We utilize the standard R Hβ − L 5100 relation (Bentz et al., 2013) to estimate this value. In the current version of the code, the results for the photometric reverberation method are estimated by adopting only 2 photometric channels at a time and the time delay is estimated using the χ 2 method. We account for the contamination in the emission line (Hβ, MgII, CIV) as well as the in the continuum. The code also incorporates the FeII pseudo-continuum and contamination from starlight i.e. stellar contribution.   The lightcurves with luminosity 10 45 erg s −1 is shifted by factor of 4.5 with respect to the lightcurves with luminosity 10 46 erg s −1 to highlight the differences between the two cases. The corresponding time delays are reported on top of each plot for these two cases of luminosity. The simulations are currently performed using two channels.
In Martínez-Aldama et al. (2019b), we show preliminary results from our code. In that paper, we show the variation in the simulated lightcurves for MgII and Hβ as a function of the redshift and luminosity and their corresponding time-delay distributions. The power spectral distribution for these lightcurves was assumed to have the low-frequency break corresponding to 2000 days. In Figure 10, we show the results from subsequent analysis performed using two-channel setup. We incorporate a higher density of the continuum and the the line contribution owing to the fact that in the DDFs the probing will be denser compared to the rest of the All-Sky survey. The computations are performed for a representative case of black hole mass, M BH = 3 × 10 8 M . The photometric accuracy is kept at 0.01 and the low-frequency break corresponds to between 200-300 days. The EWs are increased by a factor of 3 with respect to the average values (MgII: 47Å and Hβ: 87Å) which were used in the previous results. The corresponding time-delays (with dispersion) recovered from the analyses are reported for the respective emission line (MgII and Hβ) and two representative cases of the luminosity, L 3000 = 10 45 erg s −1 and 10 46 erg s −1 . The dispersions for the reported time-delays are found to be within ∼1-2% limit.

DISCUSSIONS
The MgII and Hβ both are low-Ionization broad-emission lines, however, MgII shows a narrower profile than Hβ profile (Marziani et al., 2013), suggesting that the emission line is emitted in an outer region. InŚniegowska et al. (2018), we have tested that their emissivity profiles using CLOUDY reveals similar emitting zones in the BLR. When compared to their corresponding FeII emission, the maximum MgII and Hβ emission peaks are ∼ 10 7 cm deeper from the face of the BLR cloud. Yet, these emitting zones are extended, spanning across ∼6 orders in cloud depth. There have been recent studies that predict that the most efficient MgII emitting clouds are always near the outer physical boundary of the BLR, while the Balmer line gas is inside this outer BLR boundary (Guo et al., 2019). This explains the low variability shown by MgII in comparison to Hβ. Since MgII usually does not display a "breathing mode", it suggests that there is not a Radius-Luminosity relation for this line. However, the possibility of a global Radius-Luminosity is viable, as long as the outer radius of the BLR scales with the black hole mass (Guo et al., 2019). The Radius-Luminosity relation for MgII ( Figure 5) shows a slope of α = 0.59 ± 0.06, which is in agreement with the idea of a general R-L relation.
As has been shown in this paper, the Eddington ratio plays an important role in the R Hβ − L 5100 . This is straightforward as we see that if the line width is fixed but we recover a lower time delay than expected and subsequently the black hole mass estimated is lower. Hence, the Eddington ratio that we estimate will be higher. This trend of decrease of the time delay with rising Eddington ratio seems systematic Martínez-Aldama et al., 2019a). Being able to constrain their physical parameter space (in terms of the ionisation parameter, cloud density and accretion rate) and link this to the observables, such as emission line ratios and shape of the emission line profiles, would provide answers to this anomaly (Panda et al. in prep).
The next steps in the code development for modeling the light curves, include testing with different methodologies to generate lightcurves e.g. Damped-Random Walk. This also applies to the time-delay estimation where we are testing the consistency of the predicted values with other available methods e.g. Interpolated cross-correlation function (Gaskell and Peterson, 1987) and JAVELIN (Zu et al., 2011(Zu et al., , 2013(Zu et al., , 2016. In the near future, we plan to extend the analysis taking into account contribution from all the 6 channels, including the lag between the images when going from one filter to the next and including the filter replacement time as explained in Sec. 5.2.1. Finally, we plan to create lightcurves incorporating a distribution of the redshift, luminosity, M BH and EWs. A special attention should be paid to realistic estimates of the uncertainties of different time-lag methods (ICCF, DCF, zDCF, JAVELIN, χ 2 , regularity measures-von Neumann estimator, see Zajaček et al., 2019, for a general overview and the comparison) since this is one of the main sources of the scatter of the radius-luminosity relation. The effect of systematic errors on the time-lag uncertainty was analyzed by Yu et al. (2019), who compared JAVELIN and ICCF methods using simulated light curves, with JAVELIN performing better in terms of uncertainty determination. A similar result is reached by Li et al. (2019) who compare ICCF, JAVELIN, and zDCF methods using generated light curves that mock multi-object spectroscopic reverberation mapping (MOS-RM) surveys. They conclude that JAVELIN and ICCF outperform zDCF, with JAVELIN generally introducing the lowest bias into the radius-luminosity relation. However, more comparison and analysis is needed using both the observed reverberation-mapped quasars (e.g. bootstrap algorithm) as well as simulated light curves (generated by both Timmer-Koenig algorithm and damped random walk for comparison) for the traditional (ICCF, DCF, zDCF) and other, novel methods of the time-lag determination (von Neumann, χ 2 methods), with the special focus on how time-delay uncertainties propagate into constraining cosmological parameters (Czerny et al., 2013).

CONCLUSIONS
Reverberation mapping studies have been quite successful to constrain the R BLR − L λ relation for various broad emission lines (Hβ: Grier et al. 2017 and references therein;MgII: Metzroth et al. 2006;Czerny et al. 2013;Lira et al. 2018;Czerny et al. 2019;CIV: Grier et al. 2019 and references therein). In this contribution, we present the monitoring of the intermediate-redshift quasar CTS C30.10 by the SALT telescope over the course of six years and detected the time-delay of MgII-line response of τ M gII = 551 +49 −82 days in the source frame. In combination with other sources where MgII time-delay was detected, we constructed the radius-luminosity relation considering a sample of 117 sources taken from the literature, and fitted it with the general power-law relation log(R BLR /1lt-day) = K + α log (λL λ /10 44 erg s −1 ), with the best-fit coefficients of K = 1.47 ± 0.06 and α = 0.59 ± 0.06, which are within uncertainties consistent with the values of R Hβ − L 5100 relation by Bentz et al. (2013). It is thus possible to use the radius-luminosity relation also towards higher redshifts using MgII and potentially CIV line. Under the assumption that the time-delay can be measured within 5% uncertainty, about 625 sources are needed in future surveys to constrain the cosmological parameters with 2% uncertainty.
Along R Hβ − L 5100 relation there is an effect of the accretion rate, which induces a departure from the expected value, i.,e., super-Eddington sources show a time delay shorter than the expected. This effect can be corrected, recovering the classical low scatter in the relation. With the estimation of the luminosity distance, we estimated the cosmological parameters. The results are in agreement with the ΛCDM model within 2σ confidence level, which is still not suitable for cosmological results. Unfortunately, there are still some systematic errors that have to be corrected. The current sample size for such studies has now reached 100 and future reverberation campaigns promise to increase this sample by manifolds and get accurate cosmological results. Also, an improvement in the method for the determination of the time delay and a better understanding of the AGN properties (e.g. inclination angle) are required to decrease the errors and get better results (Panda et al., 2019).
LSST is one of the forerunners in such future campaigns which will cover a wide range of wavelength (3050Å λ 11000Å) and this instrument alone will provide ∼5 orders increase in the number of reverberation-mapped quasars. With LSST we will have no dearth in the number of quasars, yet if we are to find answers to some of the long-standing questions in astronomy with the help of quasars, we would require to improve the quality of the observed data from the instrument. In this paper, we provide the foundation for a filtering algorithm that will be useful for the community to account for the "good" quasars especially for the purpose of variability studies and photometric reverberation mapping. The suitability of these good quasars is based on three fundamental criteria -the emission line intensity, the band overlap and the line width. We also show the first results from our code designed to create synthetic, stochastic lightcurves incorporating fundamental quasar properties and the instrument's performance as per Ivezić et al. (2019). This will be useful -first for preparing a mock catalog of quasar lightcurves, and second, to test real data that will be obtained from the LSST in the near future.