High-$z$ Universe probed via Lensing by QSOs (HULQ) I. Number Estimates of QSO-QSO and QSO-Galaxy Lenses

It is unclear how galaxies and their central supermassive black holes (SMBHs) co-evolve across cosmic time, especially for the non-local universe ($z \gtrsim 0.5$). The High-$z$ Universe probed via Lensing by QSOs (HULQ) project proposes to utilize quasi-stellar object (QSO) host galaxies acting as gravitational lenses (QSO lenses) to investigate this topic. This paper focuses on the feasibility of this project, i.e., whether sufficiently large numbers of QSO lenses are expected to be found in various concurrent and future imaging surveys. We find that $\sim 440$ QSO lenses will reside in the Hyper Suprime-Cam Wide survey (HSC/Wide), which is expected to be the most prolific concurrent survey, with this number being boosted by one to two orders of magnitude (to $\sim 10000$) with upcoming surveys such as that conducted with the Large Synoptic Survey Telescope (LSST). We discuss several methods of how to study the redshift evolution of the $M_{\rm BH} - \sigma_*$ relation, which is a stand-out illustration of the co-evolution. In addition, we demonstrate how the intimacy of lensed images to the bright deflector QSO for most systems will affect the \textit{detectability} of QSO lenses. We estimate that only $\sim 82$ and $900$ will be detectable for HSC/Wide and LSST, respectively; the decrease is significant yet still yields an acceptable sample for the main objective. This decrease will be less of a problem for space-based imaging surveys, for their small point spread function FWHMs will allow detections of lensed images lying relatively close to the deflector QSO, and thus unveil the less massive yet more numerous QSO hosts.

Masses of a considerable number of these SMBHs (M BH ) have been measured via multiple approaches; gas and stellar kinematics around the SMBH are used mainly for nearby, resolved galaxies (e.g., Kormendy & Richstone 1995;Ferrarese & Merritt 2000;Gebhardt et al. 2000a), and reverberation mapping is the primary method for distant QSOs (Blandford & McKee 1982;Peterson et al. 2004;Bentz et al. 2009;Shen et al. 2015a). Based on the broad line region size − luminosity (R − L) relation (Wandel et al. 1999;Kaspi et al. 2000Kaspi et al. , 2005 from reverberation mapping studies, and assuming virial equilibrium, M BH can be estimated with a single spectrum; this technique is commonly known as the single epoch method (Woltjer 1959;Vestergaard & Peterson 2006;Vestergaard & Osmer 2009;Kim et al. 2010Kim et al. , 2015bJun et al. 2015).
There are well-established relations between M BH and properties of the bulge component of galaxies, among which the most famous is the M BH − σ * relation (Ferrarese & Merritt 2000;Gebhardt et al. 2000a;Kormendy & Ho 2013), where σ * denotes the stellar velocity dispersion of the bulge. The tightness of such correlations imply that SMBHs play an important role in galaxy evolution (Woo et al. 2010), but details remain veiled (Peng 2007;Jahnke & Macciò 2011 propose a different scenario).
To explore how this so-called "co-evolution" actually took place, it is crucial to examine the evolution of the forementioned correlations, if at all. Unfortunately, this is especially difficult at high redshifts, starting from z ∼ 1, where decomposition of galaxies from their central QSOs becomes a delicate task due to their small angular sizes (Peng et al. 2006;Decarli et al. 2010;Merloni et al. 2010). This turns out to be more complicated due to the utterly strong luminosities of QSOs causing them to outshine their host galaxies even for closer targets (Gebhardt et al. 2000b), which makes detecting and studying the host galaxies a challenging task.
In the case that these QSO host galaxies act as strong gravitational lenses (hereafter deflectors, to avoid confusion with lens systems) and distort the shapes of sources behind, the masses of these QSO hosts can be estimated using gravitational lensing (GL). The galaxies are much heavier than their hosted QSOs, albeit being less luminous, and any sort of GL features around QSOs will be predominantly due to the mass of the galaxies. Consequently, careful modeling of these GL features will allow us to measure the masses of the QSO host galaxies to high accuracy.
As mentioned above, QSOs are discovered at low to high redshifts, and their M BH are relatively easy to determine with a single spectrum. Therefore, QSO host galaxies that act as GL deflectors (hereafter QSO lenses) yield a unique method of studying galaxies and their central SMBHs. Furthermore, a sufficiently large sample of QSO lenses at various redshifts will allow us to examine the co-evolution of SMBHs and their host galaxies.
In addition, galaxy-scale GL systems have been used to constrain the density profiles of galaxies, in joint analyses with stellar kinematics or weak lensing (Koopmans et al. 2006;Gavazzi et al. 2007). QSO lenses will provide knowledge about the mass profiles of QSO host galaxies, and allow us to inspect whether QSO hosts are different from normal galaxies with quiescent SMBHs.
There has only been one group studying QSO lenses so far; the authors use Sloan Digital Sky Survey (SDSS) spectra of QSOs to look for emission lines originating from sources at redshifts further than those of the QSOs themselves. The first QSO lens, found from SDSS Data Release 7 (DR7) QSO spectra (Schneider et al. 2010), was presented by Courbin et al. (2010) (hereafter C10). be incomparable to ground-based surveys. These surveys have the following characteristics: deep limiting magnitudes ( 23 mag) and high image qualities (pixel scale & seeing 1 ′′ ) to enable the detection of the faint lensed images separated from the deflector QSO, and large areas ( 100 deg 2 ) that are required to compensate for both the low number density of QSOs and the intrinsically low probability of galaxies acting as deflectors. Therefore, a survey for QSO lenses is likely conceivable nowadays.
The objective of the High-z Universe probed via Lensing by QSOs (HULQ) project is to search for QSO lenses at various redshifts to explore the evolution of SMBH-galaxy correlations. To check the feasibility of this project, this paper aims to estimate the number of QSO lenses that are detectable in various recent and upcoming surveys. Future papers will discuss new findings of QSO lenses, and how these discoveries will affect our understanding of the co-evolution of SMBHs and their host galaxies.
This paper is organized as follows. Section 2 demonstrates how the number of QSO lenses is calculated, and the data used in these calculations are described in Section 3. Results are shown in Section 4, and further issues, such as factors that may affect our estimates, are discussed in Section 5. Section 6 summarizes our work. The AB magnitude system is used throughout this paper (Oke & Gunn 1983). We use a standard ΛCDM cosmology model with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7, which is supported by recent observational studies (Im et al. 1997;Planck Collaboration et al. 2018).

METHODOLOGY
The number of QSO lenses can be calculated following the subsequent arguments (Dobler et al. 2008). We define the Einstein volume (V Ein ) of a deflector to be the region within which a source must be located to be lensed by the deflector. The probability of a specific object to act as a deflector is equivalent to the number of sources that are lensed to be brighter than the depth of the image and located within V Ein for that deflector, which in turn is equal to the probability that a sufficiently bright source is within V Ein of the deflector, summed for all possible sources.
Expressing this in equation form, for an object at redshift z d with velocity dispersion σ, the probability that this object is a deflector to a background source, f QL , is where z s is the redshift of the source, L lim is the limiting luminosity of the imaging data that is being searched, µ 2 is the magnification of the second brightest image, n s (z s , L lim /µ 2 ) is the number density of sources at z s detectable as lensed images with the second brightest image brighter than L lim , u is the angular position vector from the deflector in the source plane, and mult indicates that this integral is calculated for the multiply-imaged region only, i.e., within V Ein . The first component depends on the cosmology, and following Hogg (1999), can be expressed as where D os is the angular diameter distance from the observer to the source redshift, so the entire component is simply a function of z s . To compute the multiply-imaged region, we need to know the parameters z d and σ of the plausible deflectors, and to obtain n s and z s , it is necessary to determine which samples can be used as likely sources. Therefore, to calculate the above equation, we need to constrain the deflector and source populations.
Since the main subject of this study is QSO lenses, it is natural to use QSO host galaxies as the deflector sample. Thus, we use the velocity dispersion function (VDF) of QSO host galaxies as a function of redshift (z d ), which is explained in Section 3.2. For the source population, the source luminosity functions (LFs) with respect to redshift are required. We assume that QSOs and galaxies are the predominant sources; these are the only samples with sufficient luminosities and number densities. These two populations are described in detail in Section 3.1.
To calculate n s , we integrate the LFs over luminosities brighter than the flux limit at a specific z s . As shown in Eq. 1, the flux limit is altered by µ 2 , the magnification of the second brightest image; this is because multiple images (more than two) are required to confirm the object as a lens system. It is possible for a source to simply be distorted into only one image and still be confirmed as a lens system, but we assume that all sources are points, and ignore this possibility. Assuming singular isothermal spherical mass distributions, the multiply-imaged region is within the Einstein radius (θ Ein ), expressed as where D ds is the angular diameter distance from the deflector redshift to the source redshift, so θ Ein is a function of z d , z s and σ. Two images with magnification µ 1,2 = 1 ± 1 (u/θ Ein ) are created, where u is the size of u. Thus, Eq. 1 can be expressed as which is a simple double integral over z s and u, of a function of z d , σ, and L lim .
To summarize, f QL can be calculated for a given deflector with z d and σ, when searching in an imaging survey with limiting magnitude m lim . The only prior information required is the redshift dependence of the source LFs, which is used to calculate the number density of sources brighter than L lim /|µ 2 | at angular displacement u from the deflector, which is integrated over u and z s .
From f QL , it is possible to calculate the surface number density of QSO lenses, n QL , i.e., the number of QSO lenses in a given area in the sky. We simply add up f QL for all QSO host galaxies within the sky area to obtain n QL . Because f QL is a function of z d and σ, and assuming that all QSOs are homogeneously distributed, VDFs of QSO hosts can be used to simplify the task. Expressing this in equation form, where Φ σ (z d ; σ ′ ) is the VDF at z d with velocity dispersion σ ′ and dV d is the differential comoving volume for a unit area in the sky at z d . Thus, n QL can be calculated for a survey depth m lim as long as the source and deflector populations are well defined.

SOURCE AND DEFLECTOR POPULATIONS
As discussed in the previous section, the redshift evolution of the VDF of the deflector population, and that of the source LFs are required to calculate n QL . In this section, we illustrate which deflector and source populations are used for our calculations.
3.1. Source Population 3.1.1. QSOs as Sources We use the SDSS DR9 QSO LFs from Ross et al. (2013), who have performed several fits for the redshift evolution of the four parameters (Φ * , M * i , α, and β) of QSO LFs in the observed i-filter, following a double power-law function of the form where the redshift dependence is from the four parameters. Among these, we select the parameterizations that fit the data well and are more or less continuous at the redshift 2.2 break; the pure luminosity evolution (PLE) model for z < 2.2, and the luminosity evolution and density evolution (LEDE) model for z > 2.2 (orange and black lines in Figure 15 of Ross et al. 2013). The two fits are slightly offset at the z = 2.2 break, so we make the fits continuous by fixing the z = 2.2 parameters to be those from the z < 2.2 fit, since this is the fit that is more trustworty, as can be seen from the prediscussed figure of Ross et al. (2013). This functional form is shown in Table 1. Although the z > 2.2 fit is for up to redshift 3.5 only, which was the limit for the data available at that time, we extend this up to redshift 7, which is where the most distant quasars are being found, while also being roughly the limit of Lyman-break sources that are detectable in the i-filter. The actual QSO LFs are shown in the top panels of Figure 1. The magnitudes given here are for the observed i-filter for QSOs at z = 2, which is the standard introduced in Richards et al. (2006). These magnitudes need to be K-corrected to the observed i-filter at the source redshift, which can be done assuming a power-law SED for QSOs, using the prescription described in Richards et al. (2006) as follows: where M i (z = z ′ ) is the absolute magnitude observed in the i-filter for QSOs at z = z ′ , α ν is the spectral slope for f ν (i.e., f ν ∝ ν αν ), and K cont (z s ) and K em (z s ) are the K-corrections due to the power-law continuum and QSO emission lines, respectively, as defined in Richards et al. (2006). The total K-correction (K cont + K em ) is given in Table 4 of Richards et al. (2006).

Galaxies as Sources
The QSO LFs given in Ross et al. (2013) are optimal for our predictions, in that the LFs are given for a certain observed wavelength range (namely the observed i-filter for objects at z = 2); we wish to calculate the expected number of QSO lenses observable in a certain filter. On the contrary, galaxy LFs in the literature are usually given with respect to their restframe wavelengths. Thus, we compiled multiple galaxy LFs for various wavelengths at various redshifts, that correspond to the observed i-filter wavelengths. For instance, the V -filter LF is used for z s ∼ 0.4 QSOs, ∼ 1500Å for z s ∼ 4, and so on. The list of compiled galaxy LFs is given in Table 2.
We compute the K-corrected absolute magnitudes of the LF break, M * i (z = z s ), as follows: where M F (z = z ′ ) is the absolute magnitude observed in filter F emitted from sources at redshift z ′ , and λ F,z ′ is the central wavelength of filter F redshifted from z ′ so that λ F,z ′ = λ F,0 /(1 + z ′ ). The final approximation holds because the filter wavelengths are selected to be redshifted to the restframe i-filter wavelength, i.e., λ i,zs ≈ λ F,0 .
Since a continuous functional form (with respect to redshift) is desired for computational purposes, we fit the redshift evolution of the three parameters (log Φ * , M * i (z = z s ), and α) of the compiled galaxy LFs, assuming the evolution to be linear, which is the simplest case possible. The linear fit does not suit the log Φ * data well, so we employ a two-piece linear fit for it. The redshift evolution and parameterizations are shown in Figure 2 and Table 1. We can see in Figure 2 that our parameterized fits the data reasonably well, so these choices are justified. The galaxy LFs are shown in the bottom panels of Figure 1.

Deflector Population
The deflector population in question is constrained to QSO host galaxies. However, although the VDF for normal galaxies in the local universe is relatively well-known (e.g., Choi et al. 2007;Sohn et al. 2017), that for active galaxies, and moreover its cosmic evolution, is poorly studied. As an alternative, it is possible to infer the VDF from the active black hole mass function (BHMF) using the M BH − σ * relation. Yet, again our understanding of the redshift evolution of the BHMF is unsatisfactory.
Thus, two approaches are used to construct the deflector VDFs. Our main approach is to infer the VDF from the QSO LFs using empirical relations, and the second approach is to derive it simply from the observed QSO population; we will call these VDF 1 and VDF 2, respectively.
For the first approach, we contemplate that the black hole mass of a QSO can be estimated from its luminosity. Theoretically this is true to some extent; as the single epoch method suggests, more massive SMBHs accrete more mass, and thus they are more luminous. We check this using SDSS DR7 QSOs; Figure 3(a) shows the black hole masses from Shen et al. (2011), plotted against the z = 2 absolute magnitude for the i-filter (M i (z = 2), hereafter simply M i ). We can see that a negative correlation does exist, albeit a weak one.
To translate the LF into a BHMF, we conduct the following steps. First, the number of QSOs in a specific magnitude bin is obtained from the LF. Then, we fit the M BH distribution of these QSOs using a Gaussian function, as can be seen in Figure 3 (b). From this, the fraction of QSOs within a specific M BH bin among the QSOs in the specified magnitude bin can be estimated by integrating the normalized fitted Gaussian function over the M BH bin. Finally, we multiply this fraction with the number of QSOs in the magnitude bin, and integrate this over all magnitudes to calculate the number of QSOs in the prespecified M BH bin. Expressing this in equation form, where µ = log (M BH /M sun ), Φ M BH ∆µ is the BHMF between µ and µ + ∆µ, and G(µ, σ µ ; M i ) is the normalized Gaussian function fit for the magnitude bin at M i , with a mean value of µ and standard deviation of σ µ , both of which we can expect to differ for each M i bin.
We can also postulate that µ depends on redshift. However, Figure 4(a) demonstrates that the variation in µ for different redshift bins is relatively small compared to that for different magnitude bins, so a single value of µ is used for all redshifts in a given M i bin; the result is shown in Figure  4(b). In addition, σ µ is roughly similar for all magnitude bins when the QSO sample in the bin is sufficiently large, so this is also fixed to 0.3 dex for future calculations.
Subsequent to DR7, later generations of SDSS (BOSS/eBOSS; Dawson et al. 2013Dawson et al. , 2016 have delved into the more distant and less luminous populations of QSOs. Since QSO luminosities are directly linked to the SMBH masses, as shown above, our work, which employs the redshifts and masses of QSO host galaxies, may be significantly affected by these recent changes. Thus, we use the most recent SDSS DR14 QSO catalog (Pâris et al. 2018) to check whether these recent updates affect our results. The black hole masses of these QSOs are measured via the single epoch method, and will be published in a subsequent paper (Taak et al. in preparation). Fortunately, Figure 4(b) shows that the two QSO samples give similar results in all magnitude bins. Since the DR7 QSOs are the more complete sample, we choose to use the DR7 fit.
The same method is used to translate the BHMF to the VDF of the QSO host galaxies, using the M BH − σ * relation from Kormendy & Ho (2013), in the form of where ω = log(σ/km s −1 ). ω, the mean of the normalized Gaussian function, is obtained from a one-to-one transformation of the M BH bin using the M BH − σ * relation. As for σ ω , the standard deviation of the Gaussian, the scatter of the relation is used: 0.28 dex in the M BH direction, or 0.063 dex in the σ * direction. Figure 5 shows the BHMF derived from the QSO LF of Ross et al. (2013), and the VDF in turn from the BHMF, for various redshifts. For the transformation from the QSO LF to the BHMF, since the scatter of the correlation (σ µ ) may be uncertain, two different scatters are used. As can be expected, the use of a larger scatter yields a larger BHMF in the high-mass bins, but the low-mass bins are not affected by much.
Since VDF 1 is derived from a QSO sample corrected for completeness, we compare these results with VDF 2, using the individual SDSS DR7 and DR14 QSO black hole masses translated to velocity dispersions, with the scatter of the relation applied, in the form of Eq. 10. No completeness corrections were made for VDF 2. This is shown in Figure 6, which demonstrates that most of the QSOs residing in high-σ * hosts are observed by SDSS, although the more numerous QSOs in low-σ * hosts are not.
Thus, following the above steps, the VDF of QSO host galaxies as a function of redshift is estimated sequentially from the QSO LF. This is used, along with the source population parameters, to calculate f QL .

RESULTS
In this section, we present the statistics of the probabilities of a QSO host galaxy to be a QSO lens (f QL ), calculated using the method and data shown in Sections 2 and 3. How n QL , the surface number density of QSO lenses, depends on m lim of the surveys in question is shown. We also derive the probability distribution functions (PDFs) of f QL depending on parameters of the deflector and source populations, such as z d , z s , and M BH . Finally, the number of expected QSO lenses for several surveys are computed.
4.1. f QL , n QL and Probability Distribution Functions Figure 7 shows n QL versus limiting magnitude m lim . First of all, as expected, n QL increases for deeper limiting magnitudes. A straightforward explanation is that deeper surveys are able to detect fainter images, which originate from fainter sources, so the number of sources magnified to above the limiting flux increases, and so does the number of lens systems. Second, the two source populations cross at about m lim ≈ 15.5 mag. This simply means that below this m lim , QSO sources are more dominant, whereas at the faint end, there are more galaxies than QSOs that will be lensed. Finally, for a limiting magnitude of 26.4 mag, which is the depth of the first data release of the HSC/Wide survey in the i-filter, there is a ∼ 2.6 dex difference between galaxy and QSO sources. This offset will come up again in the following plots. It should be noted that the lower limits of the limiting magnitudes shown here (m lim ∼ 15-20) are shallower than that for the SDSS QSO sample (m lim ≈ 21 mag), so for these limiting magnitudes, the QSO deflector sample itself will not be fully detected, thus decreasing n QL further. Yet, concurrent surveys have surpassed this limiting magnitude range, and this will be relatively unimportant for the surveys that we discuss here. Figure 8 shows f QL with respect to parameters of the deflector QSO host galaxies, for m lim = 26.4 mag. The left plot is f QL versus M BH ; we can see that there is a positive correlation between the two. It is possible to infer what the slope of this correlation should be; first, although the actual relation is complicated by the dependence on u, we can say that f QL is roughly proportional to the "crosssection" for sources to be lensed, or θ 2 Ein . Also, θ Ein ∼ σ 2 * , though this is again not straightforward due to z s . Finally, from the M BH − σ * relation, σ * ∼ M 0.23 BH , so f QL ∼ M 0.90 BH , which is quite similar to the measured slopes of 0.92-0.95 for all four cases shown. In addition, there is a constant ∼ 2.6 dex offset between galaxy and QSO sources, as was previously shown for Figure 7. This demonstrates that M BH does not play a critical role in determining this offset. The right panel shows f QL as a function of z d ; there is a negative correlation. Simply put, for more distant deflectors, the Einstein volume is smaller, and there are fewer sources that can be lensed. There is also the same ∼ 2.6 dex offset in this diagram, so this offset is almost entirely due to the differences in the galaxy and QSO source population densities, not because of deflector properties. The offset is not completely constant; this will be discussed in the next paragraph in detail. Figure 9 shows the probability distribution of QSO lenses in the top plots, as functions of M BH and z d , again for m lim = 26.4 mag, and the cumulative PDF in the bottom plots. From the peaks of the PDFs, we can infer that most of the QSO lenses will be QSOs with log (M BH /M sun ) ≈ 9 and z d ≈ 0.7. It is also interesting to note that about half of all the QSO lenses will have deflector QSOs at z d 1, justifying the main objective of the HULQ project; the furthest QSO lens confirmed so far is at z d ≈ 0.3, and this project will enable us to discover and examine the relatively high-redshift QSO host population. Also, the similarity of each of the two plots in all panels confirm that differences in the source population do not affect the overall PDFs for the deflector parameters. However, if we look at the plots closely, the M BH panels are almost perfectly identical, while those for z d are slightly different. This is due to the intrinsic differences in the redshift distribution of the sources; the galaxy number density falls at a steeper rate for higher redshifts when compared to the QSO number density, so the ratio of QSOs that are lensed by high-redshift (z d ∼ 4 -5) deflectors compared to those lensed by relatively low-redshift (z d ∼ 1 -2) deflectors is larger than that for galaxy sources.

Expected Number of QSO Lenses for Various Surveys
As illustrated in the introduction, with sufficiently wide and deep surveys, the present is a good time to start searching for QSO lenses. It is necessary to find out whether this search will be prolific, and which surveys will yield the most QSO lenses to optimize our methods. The calculation of the number of QSO lenses for each individual survey is quite straightforward: n QL for the limiting magnitude of the survey in question is multiplied by its areal coverage. These numbers for several surveys are given in the fourth column of Table 3, along with some survey statistics. The number of detectable QSO lenses is also shown in the final column this table. As will be discussed later in Section 5.4, these numbers are much smaller than the number of all QSO lenses, due to the image qualities of ground-based surveys.
We can see that among the surveys that have been initiated or completed, HSC/Wide is expected to yield the most QSO lenses with ∼ 440. Although the numbers in Table 3 cannot be simply added up, since some of the surveys overlap with each other, the sample of QSO lenses in currently available imaging data is predicted to be a few dozens, and is expected to reach three digits in the upcoming several years, even with the exception of LSST.
Since QSO lenses are a deflector-selected sample, discoveries inevitably require searching for lensed features around confirmed QSOs. Most of the confirmed QSOs are from the SDSS quasar catalogs, and located in the northern hemisphere, which is also the focus of current and ongoing imaging surveys. Fortunately, several spectroscopic surveys, such as the next phase of SDSS (SDSS-V; Kollmeier et al. 2017) and the 4-metre Multi-Object Spectroscopic Telescope (4MOST; Merloni et al. 2019) surveys plan to confirm large numbers of QSOs in the southern hemisphere. The main LSST survey, which is scheduled to survey mainly the southern hemisphere, will therefore enable the unveiling of a large number of previously unknown QSO lenses. Space missions, however, are not expected to contain new QSO lenses, since their limiting magnitudes are shallower than ground-based surveys. Yet, their importance in terms of the detectability of QSO lenses will be emphasized in Section 5.4.
As expected, small surveys below ∼ 100 deg 2 provide only a small number of QSO lenses, despite their relatively deeper imaging. The predicted number depends on the depth and area of the survey; deeper limiting magnitudes lead to larger n QL , and the number of SDSS QSOs in a survey is proportionate to its area. In fact, it is possible to determine the relative importance of depth and area, from the slope of the total n QL from Figure 11. Consider two surveys, one survey having one tenth of the imaging area of the other but reaching deeper, with identical total exposure times. The smaller survey will have a limiting magnitude 1.25 mag deeper, so assuming that the number of QSOs is proportionate to the area of the survey, n QL must increase by a factor of 10 for a 1.25 mag deeper survey. This results in a slope of (log 10)/1.25 = 0.8 on the log n QL versus m lim plot, which is achieved only for m lim < 16 mag. Thus, for a given total exposure time, it is more advantageous to increase the imaging area than to probe deeper to maximize the number of QSO lenses. For m lim = 25 mag, which is a typical limiting magnitude for concurrent surveys, the slope becomes 0.46, indicating that for a survey 10 times smaller, the depth must be (log 10)/0.46 = 2.2 mag deeper to preserve the number of QSO lenses. For example, the difference in the depths of HSC/Wide DR1 and KiDS is 2.2 mag, and the area coverage of KiDS is 15 times that of HSC/Wide DR1. Thus, the yield from KiDS is 1.5 times that of HSC/Wide DR1. This is the reason for the small surveys having such few QSO lenses; most small surveys that are a set along with larger surveys have areas less than one tenth of their counterparts, but do not reach 2.2 mag deeper limiting magnitudes. Consequently, to search for QSO lenses, it is usually more advantageous to examine larger surveys, and depth is the secondary option.
To summarize, the number of QSO lenses is expected to reach several dozens in currently released data, and will be increased by ongoing surveys within the next few years; LSST will boost that number significantly during its execution. Therefore, the number of QSO lenses is expected be increased from a few tens with currently available data, by several orders of magnitude within a decade.

Comparison with Literature: Galaxy Deflectors
We compare our results with previous studies in the literature discussing galaxy deflectors. Although there are no numerical studies that focus on QSO host galaxies acting as deflectors, Oguri & Marshall (2010) (hereafter OM10) have presented the anticipated number of QSOs lensed by galaxies, while Collett (2015) (hereafter C15) showed the expected number of galaxy-galaxy lens systems, both for various imaging surveys. Here an A-B lens system indicates that the source B is lensed by the deflector A. We show the deflector and source redshift distributions for potential QSOgalaxy and QSO-QSO lens systems in Figure 12, and compare these with the analogous distributions from these two papers.
The left panel of Figure 12 shows the deflector redshift distributions for galaxy and QSO sources, and the right panel displays the source redshift distributions. We can see that the source population does not affect the potential deflector population much, inferred from the similarity of the distributions. The deflectors have a peak slightly below z = 1, while the sources are peaked at z ≈ 2. These characteristics are similar to those of the distributions shown in Figures 5 and 6 of OM10, and Figure  1 of C15, which are overplotted with dashed lines. Overall, the distributions are similar to those found in the literature, barring a few notable differences; the PDF of the QSO host deflector redshifts peaks at a further redshift compared to those of normal galaxy deflectors from the literature, and the PDF of the QSO source redshift has a pointy peak in contrast to the other source redshift distributions. The former can be explained in part by the redshift distributions of the deflector populations; as can be seen in Figure 14(b), which shows the relative number of QSOs to those of galaxies at a certain redshift, the ratio increases monotonically from z = 0 to z ≈ 0.7. This demonstrates that normal galaxies are much more abundant than QSO host galaxies in the local universe when compared with z ≈ 0.7, so the redshift distribution of QSO host deflectors should shift to a higher redshift than that for normal galaxy deflectors. The latter is simply due to the evolution of the QSO source population. According to the QSO LFs from Ross et al. (2013), as can be seen in Figure 1, the QSO number density increases monotonically up to z = 2.2, then decreases for higher redshifts, resulting in an unnatural pointed peak at z = 2.2, which is propagated directly to the z s PDF.
In addition, we plot the Einstein radius distribution of the potential lens systems in Figure 13. The θ Ein distribution from C15, overplotted here, is largely similar to our results; their PDF peaks at θ Ein ≈ 0. ′′ 4, and falls by a factor of ∼ 10 by θ Ein ≈ 1. ′′ 5. Our results give a slightly smaller θ Ein peak because the deflector population is different; their galaxy VDF drops off for smaller σ * , since it is observationally obtained, while our QSO host galaxy VDF (VDF 1) continuously increases. Notice that there is an overall ∼ 2.5-dex difference between the two distributions; this will be explained in the following paragraphs.
The lensed rate f lensed , which is the ratio of sources lensed by the deflectors to all sources, is discussed in several papers. Pindor et al. (2003) suggest f lensed is around 10 −3.4 , and OM10 proclaim that it is independent of m lim , at 10 −3.5 . We plot f lensed for our calculations, for both galaxy and QSO sources in Figure 14(a), along with the value of 10 −3.5 given by OM10. Overall, f lensed for QSO sources is roughly constant over our full m lim range, which is consistent with OM10, but there is a ∼ 2.5-dex difference.
We conjecture that this 2.5-dex difference is mostly due to the difference in the number of the two deflector samples, QSO hosts versus all galaxies. We compare the numbers of these two samples; to elaborate, the number of galaxies used as deflectors in OM10 is integrated up to the redshift in question, and the ratio of the number of SDSS DR14 QSOs within this redshift to the aforementioned number of galaxies is calculated as a function of redshift. This graph is plotted in Figure 14(c); the logarithm of the ratio fluctuates around −2-−2.5, thus confirming that the ratio of the numbers of the deflector populations is the dominant factor in the offset. This infers that the mean probabilities for galaxies and QSO host galaxies to be deflectors should be roughly similar.
In addition, we can conduct a sanity check of our calculations of the number of QSO lenses (Section 4.2) by comparing with the number of galaxy-galaxy lens systems presented in C15. Details of the source population are different for C15; they use a mock source catalog that is complete up to i ∼ 27.5 mag, in contrast to our use of LFs with no source luminosity bounds, and a cut of i = 26.8 mag for lensed images for the LSST survey. However, considering that magnification will boost the source flux by some factor, C15 should give comparable results to our calculations for LSST. The number of galaxy-galaxy lens systems in C15 is 11 million for the entire sky. Scaling for the area of the LSST survey gives 4.8 million galaxy-galaxy lenses, and decreasing this by 2.5 dex results in ∼ 15000 QSO lenses for the LSST coverage. This is similar to ∼ 9700 resulting from our calculations. OM10 also present the number of galaxy-QSO lens systems for several surveys, many of which overlap with surveys in Table 3. For instance, they predict ∼ 8100 galaxy-QSO lens systems, which is ∼ 2.7 dex greater than our results of 15 QSO-QSO lens systems. Numbers for other surveys also agree within a factor of 2-3 (with the exception of PS1/MDS, which reached a shallower depth than what was expected 10 years ago, and is different by a factor of ∼4). In all, the number of QSO lenses inferred from the number of galaxy-galaxy and galaxy-QSO lens systems are in line with our results, supporting the robustness of our calculations.
Overall, our results agree well with previous studies. Distributions of various parameters of the QSO lens population are in good agreement with the literature. The lensed rate, f lensed , and the expected number of QSO lenses are also consistent with other papers, assuming that QSO hosts and normal galaxies have similar lens statistics. The main reason behind this is that for the low-z d (z d 1.5) and high-σ * (σ * 200 km s −1 , which corresponds to log (M BH /M sun ) 8.5) regimes, where most of the QSO lenses lie (see Figure 9), the VDFs for QSO host galaxies are surprisingly similar to those of normal galaxies, apart from a normalization offset, as can be seen in Figure 6.

Comparison with Literature: Imaging versus Spectroscopy
We also compare our results with n QL observed with spectroscopy from the literature. As mentioned in the introduction, the sole comparison can be made with the results of C10 and M19. To reiterate, these studies have looked for candidate QSO lenses using SDSS QSO spectra, searching for emission lines that originate from redshifts further than that of the QSO, in the QSO-subtracted residual spectra. The differences between our work and these studies is that i) we are using imaging data, as opposed to spectral data, so the depths of the data will have distinct meanings; ii) imaging allows us to look for lensed images at large angular distances from the deflector QSO, while spectra can only detect those within the fiber aperture; and iii) imaging requires at least two separate lensed images to confirm as a lens system, but spectra accumulates flux from all images lying within the aperture. Therefore, the number of QSO lenses estimated in this work must be adjusted for the fiber size, the depth of the spectra, and the total magnification to compare with the number of spectroscopically observed QSO lenses.
First, we recalculate f QL for the SDSS fiber size. As an initial calculation, we simply alter Eq. 4 to consider only lensed images within the radius of the fiber as follows: where µ tot is the sum of magnifications of all images inside R SDSS , which is the SDSS fiber radius, which we simply assume to be 1 ′′ , the radius of the BOSS fiber, with which most QSO spectra were taken. We have assumed that the fibers are centered on the targeted deflector QSOs, and the lensed images have PSFs in the shape of two-dimensional δ-functions. For comparison, we also calculate the above for fiducial fiber sizes of 0. ′′ 5 and 0. ′′ 1 radii. The results are illustrated in Figure 15; we can see that larger fibers are able to gather light from lensed images that are further away from the deflector QSO, so they have larger f QL , and consequently yield more QSO lenses. It is worth noting that the 0. ′′ 5 fiber and imaging produces almost identical results. For an SIE lens, as we have used throughout this paper, the dimmer image is inside the Einstein radius, while the brighter image is outside it. Since the peak of the θ Ein distribution is at ∼ 0. ′′ 5 (Figure 13), for the 0. ′′ 5 fiber, usually only the dimmer image is inside it. Thus µ tot becomes just µ 2 for these systems, and Eq. 11 becomes identical to Eq. 4. The actual results differ slightly since not all lens systems have θ Ein = 0. ′′ 5, but the overall results should be quite similar, as is shown here.
The predictions made above are correct if we were to observe with a very large telescope in an atmosphere-free environment, effectively giving very small PSF FWHMs; in reality, atmospheric seeing and telescope optics cause the images to be blurred according to a point-spread function (PSF), and light that belongs to lensed images within the fiber may be pushed outside of the fiber aperture, or vice versa. We take this effect into account, assuming a seeing of 1. ′′ 5 (the mean seeing for SDSS spectroscopy time; Gunn et al. 2006), and use only the light that remains in the fiber after a Gaussian-seeing convolution. The results are shown as dashed lines in Figure 15; even for the largest fibers, a non-negligible portion of the image flux is removed, so the number densities drop inevitably by design. Also, the lensed images are not located at the center of the fiber, so this further reduces the light fraction within it. The decrease is more significant when the apertures are small, since for smaller fibers, more light will be expelled from the fiber.
To compare with the number of QSO lenses from C10 and M19, we need to constrain the deflector sample further. C10 used the DR7 sample, with an additional constraint of z d < 0.7, while M19 used the DR12 sample, also with the z d < 0.7 constraint; we will call these Samples 1 and 2. We plot the expected n QL as functions of m lim for these two deflector samples in Figure 16(a). We also plot the results for the DR14 QSO catalog with the same redshift cut of z d < 0.7 applied, which we label as Sample 3. Sample 4 is simply the full DR14 sample, shown for comparison purposes, and these results are plotted as well. We can see that n QL for Sample 4 is the greatest, and by a significant margin over the others; this can be mostly attributed to the ratio of the number densities of the deflector QSO samples. We divide n QL by the number density of QSOs for each sample to obtain the average f QL ( f QL ), to simply compare the properties of the deflectors; these results are shown in Figure 16(b).
f QL decreases in the order of Samples 1, 3, 2 and 4, and this is due to differences in both their M BH and z d distributions, which are seen in Figure 10, and illustrated in detail in the following paragraphs. Sample 4 yields the smallest f QL , thus highlighting the importance of the deflector redshifts; as discussed for Figure 8, more distant deflectors have smaller Einstein radii, and have a smaller range for the source redshift (the source must be located behind the deflector), so their Einstein volumes are smaller, and so are their f QL . A large portion of Sample 4 are located beyond the redshift cut for the other three samples, which leads to the large deviation of f QL for this sample from those for the rest.
The sample selection for pre-BOSS (SDSS-I/II) QSOs are different from those for BOSS (SDSS-III) and eBOSS (SDSS-IV) QSOs, and the latter surveys have a deeper limiting magnitude, allowing them to probe the low-M BH regime. Since f QL ∼ M 0.90 BH , we also show the average of M 0.90 BH for each sample, to use as a proxy for f QL , in Figure 10(a). We can see that Sample 1 has a larger M 0.90 BH than the other two, explaining why f QL is larger for Sample 1 compared to Samples 2 and 3.
f QL for Samples 2 and 3 are somewhat trickier. They have similar M 0.90 BH , and their redshift ranges are identical. However, a closer look at Figure 10(b) reveals that when compared to Sample 2, Sample 3 has a steeper slope at the low-redshift end (z d 0.3), and a more gradual slope at the high-redshift end (z d 0.5) of the sample. This means that the fraction of QSOs at the low-redshift end is greater for Sample 3, and so should their Einstein radii. Consequently, Sample 3 should have a larger f QL than Sample 2.
Second, the limiting magnitude of the lensed images for the SDSS QSO spectra must be determined. Among the three confirmed lens systems given in Courbin et al. (2012), we select SDSS J1005+4016, since it is the only system with a lensed image that is deblended relatively easily from the deflector QSO and its host galaxy in HST F814W-filter imaging. We first remove the bulge component of the host using the ellipse task of the Image Reduction and Analysis Facility (IRAF). Since the bar of the host passes through the image, we measure the fluxes of two components: the upper bar with the lensed image, and the lower bar. We then take the difference of the two fluxes to obtain the flux of the lensed image, assuming that the bar is symmetric, and then convert this flux to magnitude, which is determined to be 21.1 mag in the F814W-filter.
The SDSS spectrum of this system tells us that the spectral lines used to determine this as a QSO lens candidate, [O II], Hβ, and the [O III] doublet, are close to the detection threshold, with S/N < 6. Thus, we can infer that the "limiting magnitude" of SDSS spectra is ∼ 21 mag in the F814W-filter, and almost identical for the i-filter also. Another sanity check can be done with the i-magnitude distribution of SDSS QSOs shown in Figure 17; it peaks slightly above ∼ 20 mag and falls off steeply after ∼ 21 mag, so again this is roughly where the detection threshold is. Although the detection limits for imaging filters and that for discrete emission lines must be different, this value will give a rough estimate of the detection limit of the SDSS spectra.
It is worth noting the case where a source is not multiply imaged, but the single lensed image is within the fiber. These cases will be classified as QSO lens candidates when using spectroscopy, but will not be true QSO lenses. Figure 15 shows these cases, or "false detections", for the BOSS fiber in violet. For the SDSS spectral "limiting magnitude" of m lim ≈ 20 -22 mag, the number density of false detections are similar to those of actual QSO lenses. This implies that roughly half of the QSO lens candidates discovered via spectroscopy are not bona fide QSO lenses.
Finally, we take the value of f QL from Figure 16(b) that corresponds to m lim ≈ 20 -22 mag. Sample 1 gives us −4.4 < log f QL < −3.4, and Sample 2 gives −4.8 < log f QL < −3.7. The values given from C10 and M19 are −3.2 and −3.5, but these are calculated assuming that all of their candidates (14 and 9) are actual QSO lenses. As discussed above, it is likely that about half of these candidates are non-lenses. In addition, in the case of C10, considering that one of the four observed candidates was confirmed as a non-lens, the number of actual QSO lenses ranges from three to 6.5, which gives −3.9 < log f QL < −3.5 and is more in line with our results for Sample 1. Since no QSO lenses were confirmed in M19, their log f QL ≈ −3.5 is simply an upper limit, which is likely decreased by half to log f QL −3.8, so this is also consistent with the Sample 2 results. Thus, our calculations agree with the f QL values from the spectroscopic searches.

Evolution of the M BH − σ * Relation
There is much controversy about whether the correlations between the properties of SMBHs and their hosts evolve with redshift. It is of general belief that for galaxies with similar masses (or velocity dispersions), the central SMBHs at high redshifts (z ∼ 1) are heavier than their counterparts in the local universe (Peng et al. 2006;Woo et al. 2006Woo et al. , 2008Treu et al. 2007;Ding et al. 2017). However, whether this trend is sufficiently noteworthy to call it so much as an "evolution" of such relations is of significant debate (Shen et al. 2015b;Sexton et al. 2019). Since the hosts' velocity dispersions, not masses, are thought to be the main driver for co-evolution (e.g., Silk & Rees 1998), and the main interest of GL studies is towards velocity dispersions, due to its use when calculating Einstein radii, this section focuses on the possible evolution of the M BH − σ * relation only.
We adopt the trend from Woo et al. (2008), whose proposed evolution is the most extreme among various studies, to check how much such a change at higher redshifts affects our results. The redshift dependence of the offset of log M BH relative to the local relation is as follows: ∆ log M BH = (3.1 ± 1.5) log(1 + z) + 0.05 ± 0.21.
For our calculations, we ignore the y-intercept term for agreement with the z = 0 population. We use this offset from the local relation, whilst ignoring the scatter, and translate it into the negative offset of the velocity dispersion, as: We can see that this does not bring a large difference to σ * ; at z = 1, where most of our deflector QSOs are located (Fig. 9b), σ * is decreased by 0.2 dex, and even at z = 5, the decrease is about 0.5 dex. Considering that f QL is roughly proportional to the square of θ Ein , which is in turn proportional to the square of σ * , a large portion of the deflector QSO hosts will have θ Ein smaller by ∼ 0.4 dex when evolution is taken into account, and thus, we can expect f QL to be decreased roughly by ∼ 0.8 dex.
We use this newly-calibrated σ * to calculate a revised set of θ Ein , and in turn, corrected f QL , which is plotted as a function of m lim in Figure 18. For m lim = 26.4 mag, f QL is corrected by a factor of ∼ 0.7 dex, which is similar to the above estimate. Considering that the evolution of the M BH − σ * relation that we adopted is the most extreme case, it is reasonable to conclude that the actual correction is much smaller for the surveys we discussed, giving only a minor effect (decrease of a factor of ∼ 2 -3) on the expected number of QSO lenses, and indicating that the feasibility of the HULQ project is secure. Now, then, will it be possible for us to infer the evolution of the M BH − σ * relation from any confirmed QSO lenses from HULQ? If there are a sufficient number of QSO lenses for some redshift range, the M BH − σ * relation at that redshift can be drawn, with M BH measurements from singleepoch spectra and σ * from lens analysis. If this is possible for various redshifts, the overall evolution of the relation can be determined.
In the case that the QSO lens sample is plentiful but are distributed rather haphazardly in redshift space, statistics of these lenses can be used to infer the redshift evolution of the relation. First, as was shown in Figure 18, the overall number density of QSO lenses will change depending on the direction and amplitude of the evolution. Also, Figures 12 and 13 show the z d and θ Ein distributions for when the evolution is taken into account, and compare them with those for the no-evolution scenario. For the z d distribution, the peak is at a redshift that is significantly lower than that for when the relation is fixed. Similarly, the θ Ein distribution peaks at θ Ein ≈ 0. ′′ 15, a factor of ∼ 2 smaller than that for a constant relation, and decreases faster for larger θ Ein . In addition, the number density of QSO lenses for large values of θ Ein (i.e., θ Ein 1 ′′ ) decreases most significantly (by a factor of ∼ 20 for this most extreme evolution scenario), so if we can constrain the QSO lens number density accurately, it will be possible to figure out how the M BH − σ * relation evolves. Finally, if the θ Ein distribution of the QSO lenses is sufficiently wide, the slope of this distribution offers us another method of checking the evolution of the relation.

Detectability of Lensed Images
The calculations executed up to this point assume that all lensed images brighter than the limiting magnitude are detectable. In practice, the lensed images will be close (∼ 1 ′′ ) to a bright QSO and usually be much dimmer, so they are likely to be buried within the shot noise of the QSO. In order to take this into account, we consider the actual detectability (f D ), which depends on a number of parameters, e.g., image quality, brightness of the deflector QSO and lensed images, etc. In general, f D is a function of five parameters: the deflector QSO magnitude (m 1 ), the lensed image magnitude (m 2 ), the angular distance from the deflector QSO to the lensed image (θ 12 ), the PSF FWHM of the survey, and the limiting magnitude of the survey. For a given survey, the latter two are fixed, so f D depends on the former three parameters: m 1 , m 2 , and θ 12 .
First, we need to understand how the three parameters affect f D . It is evident that the lensed image is more likely to be detected when the deflector QSO is dim, the lensed image is bright, and the angular distance between them is large. Thus we can expect f D to be monotonically increasing for m 1 and θ 12 , and monotonically decreasing for m 2 .
For a survey with a certain depth and a parameter set of m 1 , m 2 , and θ 12 , f D can be calculated with the following arguments. We construct an image created with the three parameters, and add random noise corresponding to the image quality of the survey, i.e., sky noise, along with Poisson noise from the deflector QSO and lensed images. We fit the image using GALFIT (Peng et al. 2002) with a PSF of the survey. Then, we run SExtractor (Bertin & Arnouts 1996) on the residual images to check whether the lensed images can be detected. This process repeated for a statistically meaningful number of trials will give us the expected f D for the parameter set, for a given survey.
The ideal method is to calculate f D for all possible parameter sets. However, since this is impossible to achieve time-wise, we create an array of the three parameters for which f D will be estimated, for a single survey; HSC/Wide DR1 in this occasion. Regarding m 1 , most of the SDSS DR14 QSOs have magnitudes between 18.5 and 21.5 mag (as can be seen in Figure 17), so we employ three values of m 1 : 19, 20 and 21 mag. Since the limiting magnitude of HSC/Wide DR1 is 26.4 mag, the steps used for m 2 is from m 1 to 26.5 mag, in 0.5 mag steps. Also, θ 12 is varied from 0 to 2 ′′ , with step sizes of one-sixth of a pixel, or 0. ′′ 028. Each configuration is created 500 times, and f D is calculated to be the rate of detection among the 500 PSF-subtracted residual images. The lensed image is defined to be "detected" when the detected lensed image is less than 0.5 mag fainter than the input lensed image magnitude. A PSF obtained from the HSC PSF Picker is used for the PSF subtraction.
This result is shown in Figure 19. The top panel demonstrates that the effect of m 1 is relatively insignificant; f D drops abruptly to 0 for θ 12 0. ′′ 7, regardless of m 1 , for m 2 = 22 mag. For dimmer lensed images, the decline is more gradual, and the cutoff moves to slightly larger θ 12 . Since the dependence on m 1 is weak, it is reasonable to assume that the m 1 = 20 mag plot, shown on the bottom panel, applies for all QSOs.
Using the values of f D calculated as discussed above, we alter Eq. 4 to get the modified f QL as where L lim is the luminosity corresponding to the limiting magnitude of HSC/Wide DR1, and f D is f D for the nearest m 2 and θ 12 . For instance, for a deflector QSO of 20.2 mag, and for deflector-source configurations that give a dimmer image magnitude of 24.3 mag and deflector-image angular distance of 1. ′′ 5, f D corresponding to m 2 = 24.5 mag, and θ 12 = 1. ′′ 512 is used. We then incorporate this into Eq. 5 to obtain the number density of detectable QSO lenses, and multiply it by the area of HSC/Wide DR1 to obtain the number of detectable QSO lenses in HSC/Wide DR1 to be 8.0; accounting for f D has decreased the number of QSO lenses to roughly one-fifth for the image quality of HSC/Wide DR1. It is reasonable to presume that this factor of ∼ 20% should be related to the θ Ein distribution shown in Figure 13 and the seeing of the survey; larger seeing values should naturally lead to lower f D at a given θ 12 , which is linked to θ Ein . We speculate that the fraction of detectable QSO lenses is equivalent to the fraction of QSO lenses with θ Ein greater than the seeing multiplied by some factor, and since the fraction of the θ Ein distribution with θ Ein > 0. ′′ 96, or 1.7 times the seeing of HSC/Wide DR1, is about 20%, we assume that the fraction of detectable QSO lenses is the fraction of QSO lenses with θ Ein greater than 1.7 times the seeing value.
Based on this hypothesis, the number of detectable QSO lenses is calculated for each survey and given in the last column of Table 3. The decrease is significant; around ten are predicted to be detectable in currently available imaging data, and before LSST, HSC/Wide is expected to be most fruitful with 82. Still, LSST is complete, the numbers will hopefully reach one thousand, and space missions are forecast to discover a few thousands more. It becomes clear why no QSO lenses have been discovered with ground-based imaging data so far; the most concurrent large survey is PS1/3π, but its seeing was quite poor, so only ∼ 3% of the hundred or so QSO lens sample in the survey would have been found.
These results also emphasize the importance of space-based surveys for QSO lens discoveries. First, the majority of QSO lenses have θ Ein < 0. ′′ 5, meaning that surveys with seeings less than 0. ′′ 3 are desired to increase the QSO lens sample by factors of several. For instance, CSS-OS is expected to recover ∼ 80% of all QSO lenses with its excellent R EE80 ≈ 0. ′′ 15, so compared to the ∼ 20% for HSC/Wide DR1, which has the best seeing for ground-based surveys, a much larger fraction of QSO lenses can be discovered by space missions. Second, these small-θ Ein systems are necessary if we wish to study the M BH − σ * relation in detail. As explained in Section 5.3, if the M BH − σ * relation is to evolve with redshift as proposed, most of the QSO lenses will have θ Ein 0. ′′ 2. The position of the peak of the PDF and its slope beyond the peak can be used to constrain the magnitude of the evolution of the M BH − σ * relation, and the relative number of QSO lenses with small θ Ein (i.e., the slope of the θ Ein distribution) will be critical in determining this. Finally, this small-θ Ein subsample allows us to probe the low-mass end of the QSO host galaxy population, thus enabling us to establish the M BH − σ * relation over a wider range of both M BH and σ * . To summarize, the upcoming space-based imaging surveys will be a decisive factor in the future of the HULQ project.

Scatters of Scaling Relations Used for VDFs
In Section 3.2, scaling relations were used to translate the QSO LFs to BHMFs, which in turn were transformed to host galaxy VDFs. As can be seen in Figure 5, the size of the scatter of the relations determines the number of QSO host galaxies in the high-σ * end, which are most likely to be QSO lenses. In this section, we demonstrate the importance of these scatters in calculating n QL .
Scatters of 0.3 dex and 0.28 dex in the M BH direction for the M i − M BH and M BH − σ * relations, respectively, were used in Section 3.2. We employ scatters of 0.5 dex for both relations for comparison, for which the results are shown in Figure 21. As expected, a larger scatter pushes more QSO hosts to higher σ * , so more QSO lenses are anticipated. Both increments in the scatters boosts n QL by a similar amount (∼ 0.15 dex at m lim = 26.4 mag). Therefore, even if the scatters were inaccurately estimated, n QL will not change by a large amount.

Upper limits to σ *
It is yet uncertain whether galaxies have a limit on velocity dispersions; some studies hypothesize that a hard upper limit to σ * exists (∼ 400-450 km s −1 ; Bernardi et al. 2008;Salviander et al. 2008), while others have found galaxies beyond such limits (510 km s −1 ; van Dokkum et al. 2009), albeit at high redshifts. Regardless, the number of galaxies with σ * > 400 km s −1 are rare, and it is possible to presume that such an upper limit exists. In this section, we test how such an upper limit affects our results.
We modify the VDF so that all QSO hosts with σ * larger than some upper limit have a velocity dispersion corresponding to that upper limit. Similar to the discussion in 5.5.1, the high-σ * end is altered, although in the opposite direction, so we can expect n QL to show the opposite behavior as was shown above. Figure 21 shows that this is indeed the case; n QL decreases when such an upper limit is introduced. However, even an exaggerated cut of σ * = 300 km s −1 does not result in a significant change in n QL (0.04 dex at m lim = 26.4 mag), and an impractical upper limit of σ * = 200 km s −1 is necessary for the reduction to be conspicuous (0.2 dex at m lim = 26.4 mag). Thus, it is safe to say that upper limits to the QSO host galaxy velocity dispersions do not affect our results at a meaningful level.

Assumption of Point-Source Lensed Images
Throughout this paper, we have assumed that all lensed images are point sources, regardless of the morphology of the potential sources. Unfortunately, as Figure 7 shows, most of the QSO lenses will have galaxies as sources, which are intrinsically extended, and more so when they become magnified. Therefore, the lensed images will have lower surface brightnesses, and become more difficult to detect; the number of QSO lenses could decrease significantly. In this section, we verify whether this assumption was correct.
For the assumption to hold, the angular size of galaxies must be similar to or smaller than the PSF seeing. The PDF in Figure 9(b) shows that the most common galaxy source redshift is z s ≈ 2, and the fraction of sources at z s < 1 is meager. The mean galaxy size at z = 2 (z = 1) is ∼ 2 (3) kpc (Ribeiro et al. 2016), which corresponds to an angular size of 0. ′′ 24 (0. ′′ 37). So we can reasonably argue that the average galaxy is smaller than the ground-based seeing at z s > 1. Galaxies that are dimmer, which are more important in this argument since they are closer to the detection threshold, must be even smaller (Im et al. 1995). Since these sizes are smaller than the best seeing achievable from the ground, it is safe to say that most galaxy sources can be assumed to be point sources. Obviously, this assumption becomes incorrect when space telescopes are considered, so the numbers given for Euclid and CSST in Table 3 should be thought as upper bounds.
In addition, gravitational lensing causes the lensed images to become more extended, so this affects the detectability issue mentioned in Section 5.4, where we assumed all images to be point sources. The most definite solution is to follow C15: create mock QSO lenses using our predefined models, and test how many of them can be actually found in simulated observations. This is beyond the scope of this paper and deserves a paper of its own, and will not be discussed here further.

A More Realistic Model
For our calculations, QSO host galaxies are assumed to have singular isothermal spheroid mass distributions, which is a robust yet oversimplified model. Specifically, the two-image configuration predicted for these deflector mass models is critical for this work, in that the magnification of the dimmer image is used as the criterion for classifying QSO lenses. In the case of quad-like systems, the magnifications become more complicated. Once again, simulating mock QSO lenses will deliver more accurate results.

SUMMARY
In this paper, we introduce the HULQ project, which proposes to use QSO lenses to investigate the co-evolution of SMBHs and their host galaxies. To achieve this objective, an abundant sample of QSO lenses are required at various redshifts. We present the methodology and data to calculate the number of QSO lenses expected for various surveys. The main results are as follows.
1. The surface number density of QSO lenses is calculated as a function of the limiting magnitude of the imaging survey. Currently available surveys, such as PS1/3π and publicly released HSC/Wide data are expected to provide ∼ 300 QSO lenses, and this number will be augmented by at least one order of magnitude within the next decade. This justifies the feasibility of the HULQ project: discovering a statistically sufficient number of QSO lenses at various redshifts, with a significant portion at higher redshifts than currently known samples (z d 0.5).
2. The results above are verified by comparison with several studies regarding gravitational lenses in general. In particular, PDFs for several properties of QSO lenses are largely identical to those for gravitational lenses with normal galaxies as deflectors given in the literature. The expected numbers of QSO lenses for various surveys are also in line with those of galaxy lenses when assuming a ∼ 2.5-dex difference between the two, which corresponds roughly to the number ratio of the two deflector populations at low redshifts (one in ∼ 300 galaxies host QSOs). This implies that on average, the probabilities of normal and QSO host galaxies to be deflectors to background sources should be similar.
• In addition, our calculations were modified to be applicable to spectroscopic data. These results also agree well with those given in previous spectroscopic searches for QSO lenses, thus supporting our methods.
3. The effects of the evolution of the M BH − σ * relation is discussed in detail. The most extreme evolution of the relation decreases the number of QSO lenses by a factor of < 5, so the HULQ project is still feasible. Studying the QSO lens sample in detail will allow us to determine the direction and amplitude of the redshift evolution of the relation in many approaches.
4. We simulated the effects of the Poisson noise of the QSO flux on the detectability of lensed images. Based on these simulations, we propose that QSO lenses with θ Ein smaller than ∼ 1.7 times the seeing of a survey cannot be detected, thus decreasing the number of detectable QSO lenses; for HSC/Wide DR1, this factor is ∼ 20%. This demonstrates the importance of high-resolution imaging surveys, and underlines the significance of space-based surveys, both in discovering more QSO lenses and probing the low-mass QSO host galaxies.

5.
We discuss various factors that may affect our results. Uncertainties in the derivation of the VDF do not change our results significantly ( 0.2 dex). The most critical factor is the oversimplification of the lens system design. The most direct solution is to create mock QSO lenses with more realistic models, and measure their detectability more accurately.
To sum up, the future of the HULQ project is promising. Even when many factors are considered, the QSO lens sample discoverable with imaging data is expected to be abundant, especially with the advent of surveys conducted with space-based telescopes. This will provide us with a new method of studying the co-evolution of galaxies and their central SMBHs at high redshifts. Software: SExtractor (Bertin & Arnouts 1996) Treu, T., Woo, J.-H., Malkan, M. A., et al. 2007, ApJ, 667, 117 van Dokkum, P. G., Kriek, M., & Franx, M. 2009, Nature, 460, 717 Vestergaard, M., & Osmer, P. S. 2009, ApJ, 699, 800 Vestergaard, M., & Peterson, B. M. 2006, ApJ, 641, 689 Wandel, A., Peterson, B. M., & Malkan, M. A. 1999 Wang, F., Fan, X., Yang, J., et al. 2017, ApJ, 839, 27 Woltjer, L. 1959, ApJ, 130, 38 Woo, J.-H., Treu, T., Barth, A. J., et al. 2010, ApJ, 716, 269 Woo, J.-H., Treu, T., Malkan, M. A., et al. 2006, ApJ, 645, 900 Woo, J.-H., Treu, T., Malkan, M. A., et al. 2008 −  . n QL versus m lim . Black and red lines indicate whether VDF 1 or VDF 2 (DR14) was used, respectively, and the dotted and dashed lines represent the two source populations, galaxies and QSOs, respectively. The vertical blue dashed line indicates m lim = 26.4 mag, which is the limiting magnitude for the i-filter of HSC/Wide DR1, and the double-ended arrow shows the difference in n QL for the two source populations for VDF 1.    no seeing 1".5 seeing Figure 15. n QL versus m lim , for spectroscopy. The black squares show n QL for imaging (for VDF 2 using SDSS DR14 QSOs), and is equivalent to the red solid line in Figure 7. The blue, green, and red lines are for when the size of the fiber is taken into account, for the BOSS fiber (radius of 1 ′′ ), and fiducial fibers of 0. ′′ 5 and 0. ′′ 1 radii, respectively. Solid and dashed lines are for when seeing is ignored, and a seeing of 1. ′′ 5 (mean seeing at Apache Point Observatory) is used, respectively. The violet dashed line indicates the "false detection" number density described in Section 5.2, for the BOSS fiber with the seeing applied.    ∼ 2200 1000 0.18 10 720 1 5-σ limiting magnitude for point sources, unless noted otherwise.
2 R EE80 for space missions, with the exception of WFIRST.
6 Euclid is expected to observe in the optical wavelengths with the VIS instrument through a single wide filter, with wavelength coverages of 550-900nm, which is assumed to be equivalent to the i-filter here.
8 WFIRST /WFIHLS is not expected to conduct an optical survey; numbers are given assuming that a survey in the i-filter is undertaken. The details of the survey are for the J-filter portion of the survey.