Most Lensed Quasars at z > 6 are Missed by Current Surveys

The discovery of the first strongly lensed quasar at z > 6 (Fan et al. 2018) marks a breakthrough nearly two decades after the first theoretical prediction of their existence (Wyithe & Loeb 2002). Fan et al. (2018) report the detection of J0439+1634, a z = 6.51 quasar whose best-fit model predicts three images and a total magnification μ = 51.3 ± 1.4. This lensed quasar is the first example of this class detected during the reionization epoch (Barkana & Loeb 2001; Fan et al. 2006).

Wyithe & Loeb (2002) suggested the possibility of a high fraction (up to ∼1/3) of lensed sources among high-z quasars, due to the large lensing optical depths reached at z ≳ 6. This fraction of the then-known population of high-z quasars would have had their observed flux magnified by factors μ ≳ 10. Thus far, none of the ∼150 quasars observed at z ≳ 6 (Bañados et al. 2018) have shown evidence for lensing.

As suggested by Fan et al. (2018), this lack of lensed sources could be accounted for by a selection bias. The selection of z ≳ 6 quasars requires a non-detection (or at least a strong drop) in the drop-out band at shorter wavelengths than the Lyman break (912 Å). In their detection paper, Fan et al. (2018) point out that the presence of a lens galaxy along the line of sight of the lensed quasar injects flux into the drop-out bands. This effect is very relevant for partly resolved images, a likely occurrence for high-z sources. This can thus constitute a relevant issue in the selection of z ≳ 6 lensed quasars. As the drop-out bands also contain transmission windows of the intergalactic medium (IGM), Fan et al. (2018) caution that great care is needed in separating the continuum of the lens galaxy from the transmission spikes of the IGM.

This new discovery has several key implications. First, it suggests the likely existence of a population of more modestly lensed quasars that are thus far undetected (Fan et al. 2018). Second, if this population of undetected quasars does exist, it could potentially impact in a very significant way the luminosity and mass functions of the earliest populations of supermassive black holes (BHs). In this regard, Wyithe & Loeb (2002) pointed out that the inclusion of a magnification bias in the high-z quasar surveys could have a major effect on the abundance of massive halos (M_h ≳ 10¹³ ${M}_{\odot }$), for which the mass function is very steep.

In this Letter we derive the theoretical consequences of the discovery of the first strongly lensed quasar in the epoch of reionization (Fan et al. 2018). Throughout the Letter we use the latest values of the cosmological parameters from the Planck Collaboration et al. (2018).

We employ the formalism developed in Pei (1995) to analytically compute the magnification probability distribution, P(μ), due to cosmologically distributed galaxies. We define μ as the total magnification of a source at a redshift z_s, due to lenses at redshift z'. The probability distribution of magnifications larger than μ is

$$\begin{eqnarray}&&P(\gt \mu )={\int }_{\mu }^{\infty }P(\mu ^{\prime} )d\mu ^{\prime} .\end{eqnarray} \tag{ 1 }$$

Following Pei (1993, 1995), we distinguish between the magnification μ, relative to a filled beam in a smooth universe, and the magnification A, relative to an empty beam in a clumpy universe. Calling $\bar{A}$ the mean magnification, the two variables are related by $\mu \equiv A/\bar{A}$. We define the moment function $Z(s| {z}_{s})$ as

$$\begin{eqnarray}\begin{array}{rcl}Z(s| {z}_{s}) & = & {\displaystyle \int }_{0}^{{z}_{s}}{dz}^{\prime} {\displaystyle \int }_{1}^{\infty }{dA}\rho (A,z^{\prime} | {z}_{s})\\ & & \times \,({A}^{s}-1-0.4\mathrm{ln}(10)A+s).\end{array}\end{eqnarray} \tag{ 2 }$$

The quantity $\rho (A,z^{\prime} | {z}_{s})$ is the mean number of lenses in the range ($z^{\prime}$, $z^{\prime} +{dz}^{\prime}$) with magnifications in the range (A, A + dA) for a given source at a redshift z_s. The probability distribution function for μ reads

$$\begin{eqnarray}&&P(\mu )={\mu }^{-1}{\int }_{-\infty }^{+\infty }{ds}\exp [-2\pi {is}\mathrm{ln}\mu +Z(2\pi {is}| {z}_{s})].\end{eqnarray} \tag{ 3 }$$

With some analytical adjustments, this integral can be efficiently computed via a fast Fourier transform algorithm.

Assuming that the population of lenses is composed of galaxies with flat rotation curves modeled as a truncated singular isothermal sphere, their surface mass density profile can be written as

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)=\displaystyle \frac{{\sigma }_{v}^{2}}{2{GR}}{\left(1+\displaystyle \frac{R}{{R}_{T}}\right)}^{-2},\end{eqnarray} \tag{ 4 }$$

where σ_v is the line-of-sight velocity dispersion, R is the projected radius, R_T is the radius containing half of the projected mass, and G is the gravitational constant. The total mass is finite and equal to ${M}_{T}=\pi {\sigma }_{v}^{2}{R}_{T}/G$. For such a population of lenses, calling δ(x) the Dirac delta function

$$\begin{eqnarray}&&\rho (A,z^{\prime} | {z}_{s})=2\tau (z^{\prime} | {z}_{s}){\int }_{0}^{\infty }{dff}\delta [A-A(f)].\end{eqnarray} \tag{ 5 }$$

The coefficient $\tau (z^{\prime} | {z}_{s})$ is the strong lensing optical depth between source and lens galaxy at $z^{\prime}$ (see Pei 1995). The factor f is defined as $f=l/{({a}_{\mathrm{cr}}{R}_{T})}^{0.5}$, where l is the unperturbed impact parameter, and a_cr is the critical impact parameter for double images. In Equation (5), $A(f)={A}_{+}(f)$ for $f\gt 1/r$ and $A(f)={A}_{+}(f)+{A}_{-}(f)$ for f ≤ 1/r, with

$$\begin{eqnarray}&&{A}_{\pm }(f)=\displaystyle \frac{1}{2}\left(\displaystyle \frac{{f}^{2}\pm {fr}+2}{f{\left[{\left(f\pm r\right)}^{2}+4\right]}^{1/2}}\pm 1\right),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&r={\left[\displaystyle \frac{{R}_{T}{{cH}}_{0}}{4\pi {\sigma }_{v}^{2}D({z}_{s},z^{\prime} )}\right]}^{1/2}.\end{eqnarray} \tag{ 7 }$$

Here, $D({z}_{s},z^{\prime} )$ is the angular diameter distance between source and lens at $z^{\prime}$, and c is the speed of light. The size parameter r has dimensions of length^−1/2 and is usually expressed as a function of the dimensionless parameter $F=3{{\rm{\Omega }}}_{G}/[2{r}^{2}D({z}_{s},z^{\prime} )]$, with Ω_G the density parameter of galaxies. Following Pei (1995), we assume F ∼ 0.05 and we check that our results are unchanged throughout the range 0.01 < F < 0.1, which covers the full domain of interest.

2.1. Distribution of Lens Galaxies

The theoretical framework described thus far to make predictions for the lensing probability P(>μ) has to be supplemented with a realistic cosmological distribution of foreground (z ≲ 6) galaxies, assumed to be early type (E/S0) and distributed uniformly in space. The ultraviolet (UV) luminosity function for galaxies is modeled as a simple Schechter function (Schechter 1976):

$$\begin{eqnarray}&&{\rm{\Phi }}(L)=\displaystyle \frac{{{\rm{\Phi }}}_{\star }}{{L}_{\star }}{\left(\displaystyle \frac{L}{{L}_{\star }}\right)}^{{\alpha }_{g}}\exp \left(-\displaystyle \frac{L}{{L}_{\star }}\right),\end{eqnarray} \tag{ 8 }$$

where L_⋆ is the break luminosity, Φ_⋆ is the number density of galaxies of luminosity L_⋆, and α_g is the slope of the faint end. Several surveys (e.g., Schmidt et al. 2014; Coe et al. 2015) suggest that the Schechter function correctly models the distribution of galaxies up to z ∼ 6, and possibly up to z ∼ 8 (Bouwens et al. 2014). We employ previous results on the UV luminosity function for galaxies, which are customarily divided into two large redshift ranges. For z < 1 we follow Beifiori et al. (2014), while for z ≳ 1 we follow Bernardi et al. (2010) and Mason et al. (2015). As shown by previous studies (e.g., Wyithe et al. 2011; Mason et al. 2015) most of the lensing optical depth for sources at z ≳ 6 originates from lens galaxies at z ≲ 1.5. The dependence of the velocity dispersion σ_v on the magnitude of the galaxy is modeled assuming the Faber–Jackson relation (Faber & Jackson 1976), i.e., $L\propto {\sigma }_{v}^{4}$. As the peak of star formation activity is reached at z ∼ 2 (Madau & Dickinson 2014), the stellar velocity dispersion σ_v is expected to increase with redshift, at least for z ≲ 3. This redshift range includes the population of lens galaxies that contribute the most to the lensing optical depth. We model the redshift evolution of the stellar velocity dispersion as ${\sigma }_{v}{(z)\propto (1+z)}^{\gamma }$. Following Beifiori et al. (2014), we use the value γ = 0.18 ± 0.06, indicating a mild evolution. This value is consistent with other studies (e.g., Mason et al. 2015), but smaller than what found by van de Sande et al. (2013). As discussed in Mason et al. (2015), measurements of velocity dispersions at z ≳ 0.5 are very difficult, and large uncertainties in lens models result from a lack of knowledge of how the dark matter contained in galaxies (traced by σ_v) evolves with redshift.

We are now in a position to make theoretical predictions for the population of lensed z > 6 quasars. Our model for the population of sources assumes a double power-law shape for the z > 6, quasar luminosity function, with a faint-end slope α = 1.3 (Manti et al. 2017) and a variable bright-end slope, as discussed in the following.

3.1. The Population of Lensed z > 6 Quasars

The resulting probability distribution function P(>μ) for a source at z = 6.51 is shown in Figure 1. Defining $P(\mu \gt 2,{z}_{s}\gt 6)\equiv {P}_{0}$ as the probability of strong lensing (μ > 2, leading to multiple images, e.g., Comerford et al. 2002), the plot shows the result for P₀ = 0.20, valid if the slope of the bright end of the quasar luminosity function β, Φ(L) ∝ L^−β, is β = 3.6 (Yang et al. 2016), and for P₀ = 0.04, valid if β = 2.8 (Jiang et al. 2016). Note that for β ≳ 4.5 we obtain P₀ ≳ 0.92. Very recently Kulkarni et al. (2018) predicted the bright-end slope of the z ∼ 6 quasar luminosity function to be as high as $\beta \approx {5.05}_{-0.76}^{+1.18}$.

Zoom In Zoom Out Reset image size

In our case, the observed frequency of lensed z > 6 quasars is P_obs ∼ 1/150 ≈ 7 × 10⁻³ (Bañados et al. 2018), with one source at z_s = 6.51 having a magnification factor μ ≈ 51. If we assume that this source is drawn randomly from a population smoothly distributed in μ, we require that P(μ > 50) ≥ P_obs. Depending on the value of the slope β, the predictions for P(μ > 50) vary, and we obtain P(μ > 50) = P_obs for β ≈ 3.2 (see Figure 1).

A discrepancy between the frequentist approximation of the probability and the theoretically computed value can be explained by a magnification bias (Turner 1980). The corresponding factor ${ \mathcal B }$ is employed in surveys to account for the fact that lensed sources are brighter than the unlensed population from which they are drawn. When a source with an observed magnitude m is lensed, the probability of detecting it is $\sim { \mathcal B }(m)$ times higher than the probability of detection without the lensing effect. Following Kochanek et al. (2006), we express the magnification bias ${ \mathcal B }(m)$ as

$$\begin{eqnarray}&&{ \mathcal B }{(m)={ \mathcal N }(m)}^{-1}\int d\mu \displaystyle \frac{{dP}}{d\mu }{ \mathcal N }(m+2.5\mathrm{log}\mu ),\end{eqnarray} \tag{ 9 }$$

where ${ \mathcal N }(m)$ is the number count of sources with magnitude m. Assuming a typical power-law form for the number counts in flux ${ \mathcal N }(F)=d{ \mathcal N }/{dF}\propto {F}^{-\alpha }$, the corresponding expression in magnitude is ${ \mathcal N }(m)\propto {10}^{0.4(\alpha -1)m}$. Assuming α = 1.3 as the faint-end slope of the z ∼ 6 quasar luminosity function, the magnification bias ${ \mathcal B }(m)$ varies between a value of ∼25 at m_AB = 19 and a value of ∼2.5 at m_AB = 23. For J0439+1634 (m_AB = 21.04 ± 0.01) we thus expect a bias factor of ∼10. While the detection of J0439+1634 was serendipitous and not necessarily representative of the population from which it was drawn, a comprehensive probabilistic analysis with more examples can inform us about the value of the slope β, once the magnification bias factor is taken into account. For any value of β < 3.6, our predictions for the strong lensing probability P(μ > 2) are lower than the value estimated by Wyithe & Loeb (2002), i.e., P(μ > 2) ∼ 0.3. In this regard, Mason et al. (2015) pointed out that early works (e.g., Wyithe & Loeb 2002; Wyithe et al. 2011) might have overestimated the strong lensing optical depth. The values suggested by Mason et al. (2015), P(μ > 2) ∼ 3%–15%, are consistent with our estimates.

3.2. An Undetected Population of z > 6 Quasars

Figure 1 implies that, independently on the value β of the bright end of the quasar luminosity function, the probability of observing quasars with a magnification μ ≲ 10 is higher than the observed frequency of lensed z > 6 quasars. This indicates the theoretical possibility that some of the observed z > 6 quasars are magnified with factors μ ≲ 10. High-resolution imaging spectroscopy could potentially reveal whether or not a quasar is lensed. Furthermore, diagnostics such as the quasar proximity zone (e.g., Eilers et al. 2017) allow an estimate of the intrinsic luminosity of the quasar. However, it is crucial to note that the maximum image separation for J0439+1634 is ∼02 (Fan et al. 2018), close to the highest-resolution limit obtainable with the Hubble Space Telescope (HST; see Figure 2). Therefore, it is a clear possibility that many of the detected z > 6 quasars are actually lensed, with image separations below the resolution threshold. Following Barkana & Loeb (2000), in Figure 2 we show the distribution in angular separations of lensing images for z = 6.51, ${\rm{\Delta }}\theta =8\pi {({\sigma }_{v}/c)}^{2}D({z}_{s},z^{\prime} )/D({z}_{s},0)$, where $D({z}_{s},0)$ is the angular diameter distance between source and observer. The distribution indicates that a fraction of the lensed sources at z > 6 have image separations ≲01. Thus far, no multiple images have been detected with the HST in samples of ∼10 z > 6 quasars (e.g., Kim et al. 2009; McGreer et al. 2013). This is equivalent to stating that the fraction of μ > 2 magnified quasars at z > 6 is ≲10% in current samples, leading to an upper limit of β ≲ 3.1 (see Figure 1). It is worth noting that, while the diffraction limit for the HST is ∼01, its point-spread function (PSF) is very well characterized and stable. This would, in principle, allow to discern a point-source quasar from a multiply imaged one well below the diffraction limit, possibly at a level ∼0001. For example, Libralato et al. (2018) measured proper motions of stars by identifying the centroid of the PSF with a precision ∼00003. This remark could foster a re-examination of the HST images of z > 6 quasars thus far detected.

Zoom In Zoom Out Reset image size

Additionally, as originally suggested by Fan et al. (2018), the current selection criteria might have missed an important fraction of the quasar population at z > 6 because a lens galaxy in the same line of sight might have contaminated the drop-out band of the spectrum, leading to a misclassification of the source. It is remarkable to note that, if we de-magnified J0439+1634 by a factor ∼5, reducing its total magnification to μ ∼ 10, its flux would have been comparable to the one from the foreground lens galaxy and, thus, the quasar would have been misclassified. This situation could be occurring systematically for z > 6 quasars with μ ≲ 10. As the typical image separation for strongly lensed z > 6 quasars is ≪1'' and the typical angular dimension of a z ∼ 1 lens galaxy is ∼1'' (see Figure 2), the probability of lensing is equivalent to the probability of having the lens galaxy and the quasar images within the same photometric aperture. The presence of contamination from the light of the lens galaxy is thus inevitable in strong lensing situations. It is worth noting that galaxies that are sufficiently massive (M_⋆ ≳ 10¹⁰ M_⊙) to act as efficient gravitational lenses are most common at z ≲ 1 (Conselice et al. 2016).

Next, we estimate the number of z > 6 lensed quasars that could have been missed by current selection criteria. The following calculations are referred to a flux-limited survey with a flux limit of z_AB = 22.0, in agreement with the Wang et al. (2017) survey from which J0439+1634 was drawn. Let ${{ \mathcal N }}_{\mathrm{obs}}$ be the number of observed quasars, ${{ \mathcal N }}_{\mathrm{obs},{\rm{l}}}$ be the number of observed quasars that are lensed, ${{ \mathcal N }}_{\overline{\mathrm{obs}}}$ be the number of quasars that are not observed, and be ${{ \mathcal N }}_{\overline{\mathrm{obs}},{\rm{l}}}$ the number of quasars that are not observed and are lensed. We make the underlying assumption that the only cause for the lack of detection of this additional population of quasars is the fact that they are magnified, and hence the contamination of their drop-out spectra leads to their misclassification. This results in ${{ \mathcal N }}_{\overline{\mathrm{obs}}}={{ \mathcal N }}_{\overline{\mathrm{obs}},{\rm{l}}}$. We define

$$\begin{eqnarray}&&\displaystyle \frac{{{ \mathcal N }}_{\mathrm{obs},{\rm{l}}}}{{{ \mathcal N }}_{\mathrm{obs}}}={P}_{\mathrm{obs}},\end{eqnarray} \tag{ 10 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{{ \mathcal N }}_{\mathrm{obs},{\rm{l}}}+{{ \mathcal N }}_{\overline{\mathrm{obs}},{\rm{l}}}}{{{ \mathcal N }}_{\mathrm{obs}}+{{ \mathcal N }}_{\overline{\mathrm{obs}}}}={P}_{0}.\end{eqnarray} \tag{ 11 }$$

Based on our previous results, P_obs ≈ 1/150, while P₀ depends on the slope β. Solving for ${{ \mathcal N }}_{\overline{\mathrm{obs}}}$:

$$\begin{eqnarray}&&{{ \mathcal N }}_{\overline{\mathrm{obs}}}={{ \mathcal N }}_{\mathrm{obs}}\displaystyle \frac{{P}_{0}-{P}_{\mathrm{obs}}}{1-{P}_{0}}.\end{eqnarray} \tag{ 12 }$$

We define the lensing boost factor as

$$\begin{eqnarray}&&{b}_{L}=\displaystyle \frac{{{ \mathcal N }}_{\overline{\mathrm{obs}}}}{{{ \mathcal N }}_{\mathrm{obs}}}=\displaystyle \frac{{P}_{0}-{P}_{\mathrm{obs}}}{1-{P}_{0}},\end{eqnarray} \tag{ 13 }$$

where b_L is the fractional amount of undetected z > 6 quasars (normalized to the observed ones), which only depends on the overall probability of magnification P₀. Note that while a high magnification bias ${ \mathcal B }$ favors the observation of more sources, a high lensing boost factor b_L indicates that we are not observing a large fraction of the existing sources. If we had observed all of the objects that are predicted to be magnified (P_obs = P₀), then b_L = 0. The lensing boost factor is the number by which the observed BH mass function of z > 6 quasars, ${\varphi }_{\mathrm{obs}}({M}_{\bullet },z)$ needs to be corrected by a factor (1 + b_L)

$$\begin{eqnarray}&&{\varphi }_{\mathrm{real}}({M}_{\bullet },z)=(1+{b}_{L}){\varphi }_{\mathrm{obs}}({M}_{\bullet },z).\end{eqnarray} \tag{ 14 }$$

With the values of P₀ considered in this Letter, the lensing boost factor varies within the range 0.03 < b_L < 12, as shown in Figure 3. Also note that any value of P₀ ≳ 0.5 would lead to a lensing boost factor b_L > 1.

Zoom In Zoom Out Reset image size

With the value of β ∼ 5 predicted by Kulkarni et al. (2018), we might be currently missing the vast majority of the z > 6 quasar population. With ${{ \mathcal N }}_{\mathrm{obs}}\approx 150$, the previous range reads

$$\begin{eqnarray}&&5\lesssim {{ \mathcal N }}_{\overline{\mathrm{obs}}}\lesssim 1800.\end{eqnarray} \tag{ 15 }$$

Considering the space number density of z > 6 quasars from the COSMOS-Legacy survey ${{\rm{\Phi }}}_{\mathrm{obs}}\sim 100\,{\mathrm{Gpc}}^{-3}$ (with an X-ray luminosity ${L}_{{\rm{X}}}\gt {10}^{44.1}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, Marchesi et al. 2016) we might be missing a contribution

$$\begin{eqnarray}&&3\,{\mathrm{Gpc}}^{-3}\lesssim {{\rm{\Phi }}}_{\overline{\mathrm{obs}}}\lesssim 1200\,{\mathrm{Gpc}}^{-3}.\end{eqnarray} \tag{ 16 }$$

We have explored the theoretical consequences of the detection of the first lensed z > 6 quasar (Fan et al. 2018). Our results represent an advance relative to previous calculations (such as Wyithe & Loeb 2002) in several respects: (i) we employ an updated, observationally motivated cosmological distribution of lenses and sources, and (ii) we now have at our disposal an observational point to calibrate our calculations.

The probability of a source at z = 6.51 being lensed with magnification μ > 50 depends on the slope of the bright end of the quasar luminosity function, and we obtain P(μ > 50) = P_obs for β ≈ 3.2. A discrepancy between the frequentist approximation of the probability and the theoretical expectation can be explained by taking into account the magnification bias.

Depending on the slope β, current predictions for the lensing probability range within 0.04 < P₀ < 0.92, in striking contrast with the observed value of P_obs ≈ 7 × 10⁻³. It is thus likely that the observed population of z > 6 quasars contains cases with μ ≲ 10.

Additionally, the current selection criteria might have missed a significant fraction of the z > 6 quasars, due to the contamination of their light by lens galaxies. Assuming that this undetected population coincides with the set of magnified objects at z > 6, we predict that the unknown lensed quasars could account for up to about half of the currently detected population for β ≲ 3.6. Remarkably, for the most recently suggested (Kulkarni et al. 2018) value of β ∼ 5, the vast majority of the z > 6 lensed cases of quasars is still undetected. It is important to note that our formalism predicts the fraction of z > 6 that are undetected only because the magnifying effect of a lens galaxy introduces a contamination in their drop-out band. The presence of additional populations of quasars that are undetected for other reasons (e.g., obscuration, Comastri et al. 2015, or inefficient accretion, Pacucci et al. 2017) is not taken into account here.

We have introduced the lensing boost factor b_L to model in a compact way the effect of an undetected population of z > 6 quasars on the corresponding BH mass functions: ${\varphi }_{\mathrm{real}}({M}_{\bullet },z)=(1+{b}_{L}){\varphi }_{\mathrm{obs}}({M}_{\bullet },z)$. Future high-z quasar surveys, employing improved selection criteria, will be instrumental to setting constraints on b_L and, consequently, determining the corrected quasar luminosity functions.

Some discussion is due on how these predicted lensed z > 6 quasars are classified. It is reasonable to argue that they are misclassified as low-z galaxies. In fact, a contamination of the drop-out bands of z > 6 quasars (outside the IGM transmission windows) would lead to a failure of the photo-z determination, thus artificially decreasing the redshift and misclassifying the quasar as a galaxy. This is the case as long as the following two conditions are met: (i) the flux from the lens galaxy is comparable to or larger than that of the quasars, and (ii) no X-ray emission is detected from the quasar. The first condition is likely met as the lens is usually a massive elliptical galaxy at z ≲ 1.5 (Mason et al. 2015) while the quasar is at z > 6. The difference in luminosity distance would cause the flux from the quasar to be comparable to the flux of the galaxy, even if the former is a hundred times more intrinsically luminous than the latter. Regarding the second condition, it is worth noting that the presence of the lens galaxy in the same line of sight would artificially increase the total gas column density, thus decreasing the probability of an X-ray detection.

In this Letter we considered two possible solutions to the problem of the lack of detection of lensed z > 6 quasars: (i) a fraction of the currently detected z > 6 quasars is lensed, but their images are unresolved; (ii) the lensed z > 6 quasars are misclassified into low-z galaxies, due to contamination effects in their photometry. These two solutions are complementary and not in conflict with each other. Understanding their feasibility requires advanced modeling of the lens population and, possibly, ad hoc modeling of some specific lensing configurations. In fact, the detectability of the lens galaxy ultimately depends on the ratio between its flux and the magnified flux of the quasar. More advanced models are required to provide a final answer to the question: how many lensed z > 6 quasars might we be missing?

An additional population of undetected z > 6 quasars could have implications to our understanding of the epoch of reionization. For instance, it could change our view of the contribution of quasars to the reionization of the universe (e.g., Madau 2017). Furthermore, there is a long-standing search for the origin of the infrared (e.g., Kashlinsky et al. 2007; Yue et al. 2013) and X-ray backgrounds (e.g., Gilli et al. 2007) and their cross correlation (Cappelluti et al. 2013). Our predicted population of undetected z > 6 quasars, with their spectra extended from the infrared to the X-ray, could make an important contributions to these backgrounds. Moreover, if a fraction of the currently detected population of z > 6 quasars is found to be magnified, the value of their BH mass would be consequently decreased, easing the problem of early supermassive BH growth (e.g., Volonteri & Rees 2005; Pacucci et al. 2017, 2018). For these reasons, employing new selection criteria to better understand the population of high-z quasars will be of fundamental importance in the near future. Our results strongly highlight the importance of the discovery by Fan et al. (2018) and will guide future extended searches for lensed quasars.

We thank the anonymous referee for constructive comments on the manuscript. F.P. acknowledges support from the NASA Chandra award No. AR8-19021A and enlightening discussions with Xiaohui Fan, Massimo Meneghetti, and Andrea Ferrara. This work was supported in part by the black hole Initiative at Harvard University, which is funded by a JTF grant.