HELIUM REIONIZATION SIMULATIONS. I. MODELING QUASARS AS RADIATION SOURCES

When modeling helium reionization, we employ the RadHydro code, which includes N-body, hydrodynamics, and radiative transfer calculations simultaneously. The code includes a particle mesh (PM) solver for gravity calculations, a fixed-grid Eulerian code for solving hydrodynamics, and radiative transfer solved by performing ray-tracing. For more details on the hydrodynamics portion of the simulation code, see Trac & Pen (2004). For more details regarding the RadHydro code and its application to hydrogen reionization, see Trac & Cen (2007) or Trac et al. (2008).

The simulation strategy we employ for our exploration of helium reionization consists of two steps. First, a high-resolution N-body simulation is run for a given set of initial conditions. Halos are found on-the-fly using the friend-of-friends algorithm, and a corresponding catalog of spherical overdensity halos are saved at even steps in cosmological time (Trac et al. 2015). Then, using the same initial conditions, a medium-resolution simulation using the RadHydro code is run. In order to provide accurate sources of ionizing photons for the radiative transfer calculations, it is necessary to convert the halo catalogs produced from the first simulation into quasar catalogs for the second simulation. Since the resolution of the RadHydro simulations is comparatively low (typically a hydro grid unit is 10–100 ${h}^{-1}$ kpc), the simulations are not able to accurately capture the subgrid, galaxy-level physics to include quasars directly. Thus, either a halo-level scaling relation or observational constraint is needed in order to create a physically reliable sample. Rather than having to rely on scaling relations that require several steps to convert between halo mass and quasar luminosity, we use abundance matching to calculate luminosity as a function of mass, and then use observations to create a population with the proper characteristics.

In order to calibrate the proper quasar properties to use, a suite of 10 N-body P³M simulations with $L\,=\,1$ ${h}^{-1}$ Gpc and 2048³ dark matter particles were run, which corresponds to a particle mass of ${m}_{p}=8.72\times {10}^{9}\ {h}^{-1}\,{M}_{\odot }$. The total volume is thus 10 (${h}^{-1}\,$ Gpc)³; the BOSS measurement of the two-point correlation function in White et al. (2012) has an effective volume of 9.8 (${h}^{-1}\,$ Gpc)³, so the volumes are comparable. Then halo finding was performed which produced the associated halo catalog snapshot every 20 Myr between $2\leqslant z\leqslant 10$. Since only comparatively massive halos serve as hosts for the bright quasars of interest, the simulations have a sufficient resolution to capture the required number of halos.

The first step in our model construction is to define the properties of individual quasars. The two most important of these are the light curve (i.e., L(t)) and the quasar lifetime. The most common model found in the literature for the light curve of quasars is the so-called lightbulb model, in which a quasar emits radiation at a constant luminosity for a lifetime t_q before turning off. Though largely unphysical, this model has the convenience of being simple to implement in calculations. A further simplification is typically made in which it is assumed that t_q is independent of luminosity, so that this quantity becomes a universal property.

A more realistic model of the light curve is to assume an exponential form. This type of model can be motivated physically by noting that it corresponds to Eddington accretion onto the central SMBH. Several variations on this version include an exponential ramp-up to some peak luminosity followed by abrupt turn-off, a symmetric exponential about some peak luminosity, or an exponential ramp-up with a power-law fall-off in luminosity (Hopkins & Hernquist 2009; McQuinn et al. 2009). While these models are more physically motivated, they are slightly more complicated. The approach we outline below is able to reproduce a given luminosity function for quasar light curves of this form.

Specifically, we consider here two classes of quasar light curves: the "lightbulb" model and "exponential" model, defined as:

$$\begin{eqnarray}&&{L}_{\mathrm{lb}}(t)={L}_{\mathrm{peak}}{\rm{\Theta }}(t+{t}_{q}/2-{t}_{0}){\rm{\Theta }}({t}_{q}/2-t+{t}_{0}),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{L}_{\exp }(t)={L}_{\mathrm{peak}}\exp (-| {t}_{0}-t| /\tau ),\end{eqnarray} \tag{ 2 }$$

where ${\rm{\Theta }}(t)$ is the Heaviside theta function and t₀ is the time when the quasar reaches its peak luminosity L_peak. In the exponential case, the parameter τ can be treated as a free parameter in a manner analogous to t_q in the lightbulb case. Nevertheless, we relate τ to t_q, which we will describe in more detail in Section 2.4.

Another consideration is the quasar lifetime itself, which in general need not be a universal property of all quasars. We have parameterized quasar lifetime as a function of luminosity using a power-law form:

$$\begin{eqnarray}&&{t}_{q}(L)={t}_{0}{\left(\displaystyle \frac{L}{{10}^{10}{L}_{\odot }}\right)}^{\gamma },\end{eqnarray} \tag{ 3 }$$

where we vary the values of t₀ and γ. We explore models in which ${10}^{7}\leqslant {t}_{0}\leqslant {10}^{9}\,{\rm{year}}$, and $-0.25\leqslant \gamma \leqslant 0.10$. Positive values of γ imply that brighter quasars have longer lifetimes compared to dimmer ones, and $\gamma =0$ is the case of a universal lifetime for all quasars.

We have discussed considerations for the individual quasars (i.e., light curves and lifetimes), and we wish to connect them to the universal quasar population (i.e., the QLF). In order to do so, we use the concept of a triggering rate $\dot{n}({L}_{\mathrm{peak}},z)$, which dictates the differential number density of quasars that reach their peak luminosity L_peak as a function of luminosity and redshift per unit logarithmic luminosity. Using the formalism outlined in Hopkins et al. (2006), we distinguish between the peak luminosity of a quasar L_peak and the instantaneous luminosity at which it is measured for the construction of the QLF L, and relate the two with the triggering rate $\dot{n}$. Essentially, the triggering rate must be convolved with the light curve of the quasars, since the measured luminosity function reflects a given quasar's current luminosity L rather than its intrinsic peak luminosity L_peak. The result of this convolution is the observed QLF from the intrinsic triggering rate:

$$\begin{eqnarray}&&\phi (L,z)=\int \displaystyle \frac{{dt}(L,{L}_{\mathrm{peak}})}{d\,\mathrm{log}\,L}\dot{n}({L}_{\mathrm{peak}},z)d\,\mathrm{log}\,{L}_{\mathrm{peak}}.\end{eqnarray} \tag{ 4 }$$

As explained in Hopkins et al. (2006), $\phi (L)$ is the QLF (i.e., the comoving number density of quasars per logarithmic bin in luminosity), and the quantity ${dt}(L,{L}_{\mathrm{peak}})/d\,\mathrm{log}\,L$ is the amount of time that a quasar spends in a logarithmic luminosity bin. Essentially, the triggering rate can be thought of as analogous to the halo mass function, though with the light curve convolution to account for changes in quasar brightness. In simple cases of the light curve the triggering rate can be solved for analytically: in the case of a lightbulb light curve, ${dt}(L,{L}_{\mathrm{peak}})/d\,\mathrm{log}\,L$ is a delta function at $L={L}_{\mathrm{peak}}$, and so the triggering rate is proportional to the QLF:

$$\begin{eqnarray}&&{\dot{n}}_{\mathrm{lightbulb}}(L,z)=\displaystyle \frac{1}{{t}_{q}}\phi (L,z).\end{eqnarray} \tag{ 5 }$$

In the case of an exponential light curve as defined in Equation (2), we have

$$\begin{eqnarray}&&{\dot{n}}_{\exp }(L,z)=\displaystyle \frac{1}{2\tau }{\left.\displaystyle \frac{d\phi (L,z)}{d\mathrm{log}L}\right|}_{L={L}_{\mathrm{peak}}},\end{eqnarray} \tag{ 6 }$$

where the factor of 2 arises because a quasar will be observed at a luminosity L while its luminosity is increasing and then decreasing. In practice, the QLF is typically reported in magnitude units rather than luminosity. One common convention is to report the quasar's absolute i-band magnitude at z = 2, ${M}_{i}(z=2)$. This quantity is then converted to the specific luminosity at 2500 Å, ${L}_{2500\mathring{\rm A} }$, in cgs units (erg s⁻¹ Hz⁻¹) by using Equation (4) of Richards et al. (2006):

$$\begin{eqnarray}{\mathrm{log}}_{10}\left(\displaystyle \frac{{L}_{2500\mathring{\rm A} }}{4\pi {d}^{2}}\right) & = & -0.4[{M}_{i}(z=2)\\ & & +48.60+2.5\,{\mathrm{log}}_{10}(1+2)],\end{eqnarray} \tag{ 7 }$$

where $d=10\ \mathrm{pc}=3.08\times {10}^{19}$ cm. To find the approximate bolometric luminosity, the relation of Shen et al. (2009) can be used to convert ${M}_{i}(z=2)$ to luminosity in erg s⁻¹:

$$\begin{eqnarray*}&&{M}_{i}(z=2)=90-2.5\,{\mathrm{log}}_{10}(L).\end{eqnarray*}$$

One should note that this relation is approximate, and depends on the assumed spectral energy distribution (SED) of the quasar. Equation (4) is soluble for a few classes of light curves, such as the ones explored here.

The technique of abundance matching has already been applied to populations of galaxies with great success (e.g., Simha et al. 2012; Hearin et al. 2013), and has also been discussed in the context of quasars (e.g., Martini & Weinberg 2001; Porciani et al. 2004; Croton 2009). However, we wish to extend the techniques mentioned above to include different quasar light curves and lifetimes. The methods we outline below are also fairly general, and can be extended to include semi-analytic models as well. We start with the Ansatz for abundance matching of galaxies, namely that the most luminous galaxies are found in the most massive halos. This makes intuitive sense: more massive halos have more dark matter and baryonic matter to eventually convert to stars. Specifically, halo mass is highly correlated to the luminosity in the red bands, which shows the percentage of older stellar mass.

For quasars, we have a similar situation where the most luminous quasars are found in the most massive halos. However, in this case the situation is slightly more complicated because quasars have a lifetime which is much shorter than the period from the halo's formation to the activation of the quasar. Thus, we need to introduce a factor to account for the fact that not all halos host an active quasar. If we assume that the fraction of halos hosting an active quasar is universal (i.e., independent of halo mass or quasar luminosity), we can express abundance matching for quasars, assuming a lightbulb light curve, as:

$$\begin{eqnarray}&&\phi (\gt L)={f}_{\mathrm{on}}{n}_{\mathrm{halo}}(\gt M).\end{eqnarray} \tag{ 8 }$$

Expressed this way, f_on is simply the fraction of halos of a mass M that host an active quasar. Alternatively, we could define this fraction in terms of the quasar lifetime:

$$\begin{eqnarray}&&{f}_{\mathrm{on}}(L,z)=\displaystyle \frac{{t}_{q}(L)}{{t}_{H}(z)},\end{eqnarray} \tag{ 9 }$$

where in some models t_H(z) is formulated as the halo lifetime (Martini & Weinberg 2001), or the Hubble time (Conroy & White 2013). We follow Conroy & White (2013) and use the Hubble time. As we shall see, though, the exact choice for t_H(z) does not strongly affect the results. For the redshifts of interest, for a uniform value of ${t}_{q}=30\,{\rm{Myr}}$, this implies that ${f}_{\mathrm{on}}\sim 0.1 \% \mbox{--}1 \%$.

We can generalize the procedure of abundance matching to different light curves by using the triggering rate. In integral form, we can write abundance matching as equating the cumulative number of quasars above a particular peak luminosity given by the triggering rate with the cumulative number of halos given by the halo mass function. The total number of halos which should host quasars within a time interval ${\rm{\Delta }}t$ is:

$$\begin{eqnarray}&&{\displaystyle \int }_{{\rm{\Delta }}t}{\displaystyle \int }_{L}^{\infty }\dot{n}({L}^{* })d\,\mathrm{log}\,{L}^{* }{dt}\\ &&\quad ={\displaystyle \int }_{{\rm{\Delta }}t}{\displaystyle \int }_{L}^{\infty }\displaystyle \frac{{{dn}}_{\mathrm{halo}}({L}^{* })}{d\,\mathrm{log}\,{M}^{* }}\displaystyle \frac{d\,\mathrm{log}\,{M}^{* }}{d\,\mathrm{log}\,{L}^{* }}\displaystyle \frac{{dP}}{{dt}}d\,\mathrm{log}\,{L}^{* }{dt}\\ &&\quad =\displaystyle \frac{{\rm{\Delta }}t}{{t}_{H}}{\displaystyle \int }_{M}^{\infty }\displaystyle \frac{{{dn}}_{\mathrm{halo}}({M}^{* })}{d\,\mathrm{log}\,{M}^{* }}d\,\mathrm{log}\,{M}^{* }.\end{eqnarray} \tag{ 10 }$$

This form of our abundance matching equation becomes the central mechanism by which we are able to equate quasar luminosity with host halo mass. In this construction, we have implicitly used the mass-to-light ratio $d\,\mathrm{log}\,M/d\,\mathrm{log}\,L$ to convert halo mass to quasar luminosity. Additionally, we have introduced the factor ${dP}/{dt}$ to represent the probability that an individual halo will host a quasar. We have set this quantity to be equal to $1/{t}_{H}$. Thus, for the case of a lightbulb light curve and a universal quasar lifetime, this formalism reduces to Equation (8). Formally, this expression is an expansion of $\dot{n}({L}^{* },z)$ about z that is first-order accurate to ${\rm{\Delta }}t/{t}_{H}$ (Hopkins et al. 2006). Thus, so long as the time-steps between determining the triggering rate are small compared to t_H (defined either as the Hubble time or the halo lifetime, both several orders of magnitude longer than the typical quasar lifetime), this expression should reproduce the target QLF.

In the exponential case, we are free to choose the parameter τ in any way that we like, as long as it is constant with respect to L (though it may vary with L_peak). We have chosen τ such that $\dot{n}({L}_{\mathrm{peak}})$ is the same between the lightbulb and exponential cases for all luminosities. We accomplish this by equating Equations (5) and (6), and solving for τ in terms of t_q. The expression involves the ratio of the QLF and its derivative. This means that when we perform abundance matching, the same implicit mass-to-light ratio is used in the two cases. Since the halo mass function is the same between the two cases (due to the same population of halos being used), and the functional form of $\dot{n}$ is the same, we must have the same form of $d\,\mathrm{log}\,M/d\,\mathrm{log}\,L$. This has the advantage of allowing us to apply certain intuition from the lightbulb case to the less straightforward exponential case. The downside to this approach is that when exploring the parameter space of quasar lifetimes t_q in the lightbulb case, it is not immediately obvious how this translates to the exponential time constant τ, since we effectively have different values of t_q for different luminosities. For instance, even in cases where t_q is independent of luminosity, τ still changes as a function of L. However, the benefits of being able to interpret the results of the exponential case using the intuition provided by the lightbulb case outweigh the downsides of not exploring parameters in τ directly.

The general procedure is as follows.

• 1.  
  The halo mass found from the halo catalog at redshift z_cat is read and converted to an expected number density in a particular cosmology using the universal mass function described in Tinker et al. (2008). The fitted form of the mass function is used rather than the empirical one from the catalog in order to decrease the variation in number density at the high-mass end, since these quasars are disproportionately important for the reionization process.
• 2.  
  Using Equation (10), the halo number density is converted to an expected quasar number density using a specified QLF.
• 3.  
  The quasar magnitudes are binned into equal intervals in magnitude ${\rm{\Delta }}M$, such that the expected triggering rate $\dot{n}(M,{z}_{\mathrm{cat}})\ldots \dot{n}(M+{\rm{\Delta }}M,{z}_{\mathrm{cat}})$ is found, which is converted from a number density to a total number $\dot{N}(M)$ using the volume of the simulations.
• 4.  
  Within each magnitude bin, each quasar is assumed to have an equal probability of becoming active. Each quasar candidate is randomly turned on with probability $1/{\dot{N}}_{\mathrm{bin}}(M)$.
• 5.  
  To ensure that the volume self-consistently follows the merging of the underlying host halos, the quasars are propagated forward using a halo merger tree. By design, the halo catalog snapshots are made at times that are shorter than the expected lifetimes of the quasars. This approach allows for halos hosting quasars to be tracked throughout the simulation. In most cases, an active quasar from time step $i-1$ in a progenitor halo passes to the single descendent halo at time step i. Additionally, this halo hosting an active quasar is not eligible to host a new quasar. This approach covers the majority of halos for the majority of time steps. However, there are several special cases related to merger events worth discussing. Specifically, when two progenitor halos merge into a single descendent and one of them is hosting an active quasar, the descendent halo inherits the active quasar. If a single active progenitor halo splits to form two descendent halos, the larger halo retains the quasar. In the case of a merger between two active quasar halos, only the larger quasar survives. These cases represent a comparatively low number of instances of our total population evolution, and do not strongly influence our conclusions.

Throughout this work, we use a series of QLFs as determined at different epochs. For relatively low-redshift ($2\lesssim z\lesssim 3$), we use the QLFs as determined by R13 from the BOSS survey, specifically the high-z stripe 82 sample (S82) form which includes luminosity evolution and density evolution (LEDE). Above a redshift of 3, the QLF has been measured at $z\sim 3.2$ and $z\sim 4$ by M12 using data from COSMOS.¹ At $z\sim 5$, the QLF has been measured by M13 using data from the SDSS.² Although these works use slightly different values for cosmological parameters from the ones assumed here, the impact on the reported quantities is minimal.

In order to span the different epochs over which the luminosity function has been measured, it is necessary to combine the different data sets. All of the data sets fit to a double power law form of the QLF, written as:

$$\begin{eqnarray}&&{\rm{\Phi }}(M)=\displaystyle \frac{{\phi }^{* }}{{10}^{0.4(1+\alpha )(M-M* )}+{10}^{0.4(1+\beta )(M-M* )}},\end{eqnarray} \tag{ 11 }$$

where Φ is the comoving number density of quasars of magnitude M per unit magnitude, ${\phi }^{* }$ is the normalization of the QLF, α is the faint-end slope of the luminosity function, β is the steep-end slope (which is reversed from the parameterizations of M12), and M* is the so-called break magnitude where the luminosity function transitions from the faint-end to the steep-end. In most formulations at high-redshift, redshift evolution is incorporated by a change in ${\phi }^{* }$, M*, or both, that is linear in redshift. For the data from R13, the evolution is given by the equations:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\,{\phi }^{* }(z)={\mathrm{log}}_{10}\,{\phi }_{0}^{* }+{k}_{1}(z-2.2),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{M}_{i}^{* }(z)={M}_{0}^{* }+{k}_{2}(z-2.2).\end{eqnarray} \tag{ 13 }$$

For the data in M13, there is linear evolution in $\,{\mathrm{log}}_{10}\,{\phi }^{* }$ as well, given as:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\,{\phi }^{* }(z)={\mathrm{log}}_{10}\,{\phi }_{0}^{* }+{k}_{1}(z-6).\end{eqnarray} \tag{ 14 }$$

To combine the R13, M12, and M13 data sets into a single set of quantities, we first assume that the results from R13 are accurate for redshifts $z\leqslant 3.5$. This is the nominal limit of the LEDE fits, and though there are small differences between the fit QLF and the binned data, overall the fits are excellent. To incorporate the results at higher redshifts, we cast the four parameters of the QLF (${\phi }^{* }$, M*, α, and β) as quantities that have linear evolution in redshift. We define these parameters as:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\,{\phi }^{* }(z)={\mathrm{log}}_{10}\,{\phi }_{0}^{* }+{c}_{1}(z-3.5),\end{eqnarray} \tag{ 15a }$$

$$\begin{eqnarray}&&{M}^{* }(z)={M}_{0}^{* }+{c}_{2}(z-3.5),\end{eqnarray} \tag{ 15b }$$

$$\begin{eqnarray}&&\alpha (z)={\alpha }_{0}+{c}_{3}(z-3.5),\end{eqnarray} \tag{ 15c }$$

$$\begin{eqnarray}&&\beta (z)={\beta }_{0}+{c}_{4}(z-3.5).\end{eqnarray} \tag{ 15d }$$

These parameterizations are applied to redshifts where $z\gt 3.5$. The constant values are defined to be equal to the values of R13 at z = 3.5, and the values for the slopes (c₁–c₄) are allowed to take on a range of values. The range is generally chosen such that the values for the different parameters brackets the range of best-fit values provided by the highest redshift (M13) data. The fiducial values for the slopes are taken to be ones that reasonably reproduce the high-redshift measurements. Table 1 shows the fiducial values for the slopes, as well as the range of values for the parameters at $z\sim 5$ used in the parameter space exploration in Section 5. For a complete discussion on selecting the parameters for the QLF, see Appendix A.

Table 1.  A List of The QLF Parameters of the Data Sets Incorporated

Data Sets	z	$\,{\mathrm{log}}_{10}({\phi }^{* })$ a	${M}_{0}^{* }$ b	k1c	k2	α	β
R13	2.2–3.5	$-{5.93}_{-0.01}^{+0.02}$	$-{26.57}_{-0.02}^{+0.04}$	$-{0.689}_{0.027}^{0.021}$	$-{0.809}_{0.166}^{0.033}$	$-{1.29}_{-0.03}^{+0.15}$	$-{3.51}_{-0.18}^{+0.09}$
M12	3.2	$-6.{58}_{-0.79}^{+0.26}$	−27.03 ± 0.68	⋯	⋯	−1.73 ± 0.11	−2.98 ± 0.21
M12 d	4	$-7.{12}^{+0.62}$	−27.13 ± 2.99	⋯	⋯	−1.72 ± 0.28	−2.6 ± 0.63
M13 e	5	$-{8.47}_{-0.24}^{+0.20}$	$-{28.70}_{-0.33}^{+0.27}$	−0.47	⋯	$-{2.03}_{-0.14}^{+0.15}$	−4.00
M13	5	$-{7.63}_{-0.25}^{+0.30}$	$-{27.34}_{-0.49}^{+0.60}$	−0.47	⋯	−1.50	$-{3.12}_{-0.41}^{+0.28}$
M13	5	$-{7.93}_{-0.03}^{+0.03}$	−27.88	−0.47	⋯	−1.80	−3.26

Notes.

^a ${\phi }^{* }$ has units of Mpc⁻³ mag⁻¹. ^b ${M}_{0}^{* }={M}_{i}(z=2)={M}_{1450}-1.486$. ^ck₁ and k₂ are defined for models with redshift evolution in Equations (12)–(14). ^dThe authors of M12 provide a value for ${\phi }_{0}^{* }$ at $z\sim 4$ where the reported error is greater than the value itself. Since this value must be positive, the resulting lower-bound is unphysical. We reproduce the value and upper-bound here for completeness, but do not include this value directly when determining the values of the QLF. See Appendix A for further details. ^eIn M13, the authors provide three fits, each with at least one parameter held constant. Values without error ranges indicated correspond to the parameters held fixed for a particular fit.

Download table as:  ASCIITypeset image

Table 2.  A List of the Parameters Used in Equations (15a)–(15d) Based on the Data Listed in Table 1

Parameter	Fiducial Value	Parameter Range
$\,{\mathrm{log}}_{10}\,{\phi }_{0}^{* }$	−6.82	⋯
c1	−0.790	$[-1.10,-0.536]$
${M}_{0}^{* }$	−27.6	⋯
c2	−0.238	$[-0.716,0.170]$
${\alpha }_{0}$	−1.29	⋯
c3	−0.324	$[-0.493,-0.140]$
${\beta }_{0}$	−3.51	⋯
c4	0.0333	$[-0.327,0.260]$

Note. These provide a fit to the luminosity function through redshift, and ensure that the abundance of quasars matches observations as nearly as possible. For additional details on the parameters and the fitting procedure, see Appendix A.

Download table as:  ASCIITypeset image

Table 1 lists the parameters that we include from the measurements of R13, M12, and M13. The values from M12 are not included in the fitting procedure directly, and serve primarily as a consistency check due to their comparatively large error bars. The parameters from M13 are determined at $z\sim 5$, and the ones from M12 are determined at $z\sim 4$ and $z\sim 3.2$. Note that the authors of M13 provide three independent fits to their data, which are all incorporated into the final QLF parameterization. (See Appendix A for more details.) For the measurements from R13, whose fiducial LEDE model includes redshift evolution in ${\phi }^{* }$ and M*, the model is valid over a range of redshift, from $2.2\leqslant z\leqslant 3.5$. For the purposes of generating our quasar catalogs, we are interested in exploring the QLF until z = 2. For the sake of simplicity, we simply extend the LEDE model from R13 to this redshift. Although the LEDE fit is ostensibly not valid below z = 2.2, we expect helium reionization to be largely finished by this redshift, and so the precise form of the QLF at $z\sim 2$ is not of fundamental importance to our study. Also, for the value of M*, it is necessary to convert to a single magnitude system. As explained in Section 2.3, we use ${M}_{i}(z=2)$, the absolute i-band magnitude at z = 2. The QLFs of M12 and M13 use M₁₄₅₀, which is related to ${M}_{i}(z=2)$ by ${M}_{i}(z=2)={M}_{1450}-1.486$ (Richards et al. 2006; Ross et al. 2013, Appendix B). Note that this conversion assumes that the quasar SED follows a power-law with an effective spectral index of $\alpha =0.5$ (using the convention that ${f}_{\nu }(\nu )\propto {\nu }^{-\alpha }$). Modifying the spectral index α changes the magnitude conversion, so care must be taken when converting between magnitude systems. See Appendix A for further discussion.

Figure 1 shows the combined QLF from R13, M12, and M13 (which at this epoch is essentially that of R13), as well as two different quasar models at $z\sim 2.4$. We can see that there is generally very good agreement between the constructed quasar catalog and the target luminosity function, as should be expected. The differences between the constructed catalogs and target luminosity function are typically on average $\lesssim 5 \%$, which is comparable to or smaller than the uncertainties in the luminosity function itself at these redshifts. At high luminosities (${M}_{i}\lesssim -28$), though, there are some comparatively large differences that can arise between the predicted and empirical luminosity functions. This deviation is largely due to Poisson shot-noise introduced by the rarity of the objects. For objects in this luminosity range, there are typically only a few objects (${ \mathcal O }(10)$) in the entire 1 (${h}^{-1}$ Gpc)³ volume. At the dim end of the QLF, there can be insufficient halos of a particular mass given the mass resolution of our simulation. The minimum halo mass is ${M}_{\mathrm{halo},\min }=4.36\times {10}^{11}\,{h}^{-1}\,{M}_{\odot }$. Since quasars with ${M}_{i}\leqslant -25$ are most important for this study, this does not affect our results significantly.

Zoom In Zoom Out Reset image size

Throughout most of the following analysis, we focus our attention on several models in particular, parameterized in terms of t₀ and γ as in Equation (3). The first four of these models have particularly good agreement with the BOSS measurements. The last two are included to demonstrate how the clustering signal changes as a function of t₀ for a fixed value of γ: models L1, L3, and L4 all have the same γ value. We summarize these models in Table 3.

Table 3.  A List of the Parameters of Some Quasar Models Considered

Model Name	Light Curve	$\,{\mathrm{log}}_{10}({t}_{0}/\mathrm{year})$ a	γ
L1	Lightbulb	7.75	0
L2	Lightbulb	8.25	−0.125
E1	Exponential	7.25	0
E2	Exponential	7.75	−0.15
L3	Lightbulb	7	0
L4	Lightbulb	8.5	0

Note.

^at₀ and γ as defined in Equation (3).

Download table as:  ASCIITypeset image

By construction, our method matches the input QLF at all redshifts, regardless of the individual properties of the underlying quasar population. However, we are not guaranteed to match the observed clustering of quasars. Changing the implicit mass-to-light ratio of Equation (10) through changing the quasar lifetimes will affect how halos are populated with quasars. In general, longer quasar lifetimes lead to quasars of the same luminosity being matched into hosts of larger masses. Since their hosts are more biased, this leads to quasars of the same luminosity showing a larger clustering signal. This is true at all luminosities. We want to match the clustering because it can affect the topology of reionization. There can also be spatial correlations present in the radiation field as a result of reionization, which are important for making measurements of the BAO from the Lyα forest (e.g., White et al. 2010; Slosar et al. 2013).

Here, we explore how to include clustering measurements from the two-point correlation function in our quasar catalog. Recent results from the BOSS survey for the clustering of quasars in the redshift range of interest are presented in White et al. (2012). The above work examines the clustering signal of quasars in both 2D-projected and 3D-redshift-space correlation functions at intermediate scales ($3\lesssim s\lesssim 25$ ${h}^{-1}$ Mpc). The authors also introduce luminosity cuts to make the results more robust. For the purposes of this comparison, we consider their selection for which they imposed luminosity cuts on both the bright and faint ends, so that only objects with $-25\geqslant {M}_{i}\geqslant -27$ were considered across the entire redshift range (Sample 4 as defined by the authors). For a fair comparison, we impose similar cuts on our object selection. We also examine the redshift evolution of the results, and compare against the high-z/low-z samples (Samples 5 and 6) as well. See Appendix B for further discussion of these different redshift samples.

We explore the parameter space of available quasar models by examining the lightbulb and exponential light curves defined in Equations (1) and (2), as well as luminosity-dependent quasar lifetimes defined in Equation (3), parameterized by t₀ and γ. For each combination of parameters, we construct a quasar catalog in the manner described above.³ Then, we extract from this catalog all objects that satisfy the magnitude constraints at the central redshift of the survey z = 2.39. This redshift represents the average redshift of quasars chosen in the BOSS sample; the actual quasar objects span in redshift from $2.2\lt z\lt 2.8$. However, as noted in White et al. (2012), the redshift evolution of the signal is weak. Thus, extracting objects from our quasar catalogs at a single redshift rather than a range should have little effect on our overall conclusions. We measure the monopole of the two-point correlation function using the "natural estimator" ξ:

$$\begin{eqnarray}&&\xi (s)=\displaystyle \frac{\langle {DD}(s)\rangle }{\langle {RR}(s)\rangle }-1,\end{eqnarray} \tag{ 16 }$$

where $\langle {DD}(s)\rangle$ is the average number of quasar pairs from the quasar catalog separated by a real-space distance of $[s-{\rm{\Delta }}s/2,s+{\rm{\Delta }}s/2]$, and $\langle {RR}(s)\rangle$ is the number of pairs of points at the same separation drawn from a distribution with Poisson noise.

In order to quantify the statistical uncertainty in our catalog, we ran a suite of 10 N-body simulations with different initial conditions. We then performed our abundance matching procedure on each of the different simulations, including several realizations for each volume. Since our abundance matching procedure stochastically determines which halos should be hosting active quasars at a given time step, we create several quasar catalogs for each individual halo catalog, using a different initial random seed (three realizations per volume for these results). Additionally, we have augmented the effective number of samples by including redshift space distortions along the different principal axes of the simulation. This strategy gives us a total of 90 samples for which to measure the clustering signal. The best estimate for the correlation function $\xi (s)$ for a given radial bin s_i is given by averaging over all of the individual estimates ${\xi }_{k}$:

$$\begin{eqnarray}&&\bar{\xi }({s}_{i})=\displaystyle \frac{1}{N}\displaystyle \sum _{k=1}^{N}{\xi }_{k}({s}_{i}).\end{eqnarray} \tag{ 17 }$$

We then estimate the covariance between the radial bins by computing the entries of the covariance matrix C_ij. We compute the entries of the covariance matrix as (Zehavi et al. 2005):

$$\begin{eqnarray}&&{C}_{{ij}}=\displaystyle \frac{1}{N}\displaystyle \sum _{k=1}^{N}({\xi }_{k}({s}_{i})-\bar{\xi }({s}_{i}))({\xi }_{k}({s}_{j})-\bar{\xi }({s}_{j})).\end{eqnarray} \tag{ 18 }$$

The correlation matrix entries for our model L1 is plotted in Figure 2. Notice that the diagonal entries dominate, which means that the bins are mostly independent of each other and dominated by shot-noise (Valageas et al. 2011; White et al. 2012). Implicitly, the samples have been treated as being independent, and this is almost surely not the case. Although the 10 volumes as a whole can be treated as being statistically independent, the different realizations based on the same halo catalog are likely correlated. Further, the projections of peculiar velocities along different axes for the same realization are also likely to produce correlated results. However, producing a sufficient number of independent realizations to decrease the noise in the covariance matrix is computationally infeasible. Further, the variance in the clustering signal among quasar catalog realizations for a given $({t}_{0},\gamma )$ pair is comparable to small displacements in the t₀–γ parameter space, so it is necessary to include this source of uncertainty. Since we are interested only in finding models that are consistent with the BOSS measurements which have their own set of observational uncertainties, we feel that this approach produces sufficiently accurate results.

Zoom In Zoom Out Reset image size

Once the entries of the covariance matrix have been computed, the difference vector $\delta ({s}_{i})\equiv {\xi }_{\mathrm{model}}({s}_{i})-{\xi }_{\mathrm{BOSS}}({s}_{i})$ is calculated. The correlation function ${\xi }_{\mathrm{BOSS}}$ is fit to a power law: ${\xi }_{\mathrm{BOSS}}(s)={(s/{s}_{0})}^{\beta }$, where the authors have fixed the value of $\beta =-2$. In order to investigate the impact this choice has on the conclusions, we performed fits on the correlation function measured from our quasar catalogs using two different parameterizations: one where the best-fit value of s₀ was found when fixing $\beta =-2$, and another where the value of s₀ and β were both fit. In the length scales used for our analysis ($3\leqslant s\leqslant 25$), the deviation of β from the fiducial value of −2 was small, typically less than 5%. Furthermore, the values for s₀ were also largely similar between a fixed slope or a varying one, with deviations typically less than 1%. Thus, the choice to set $\beta =-2$ does not strongly bias the results presented here, or the values reported in ${\xi }_{\mathrm{BOSS}}$.

When comparing one of the quasar models with the BOSS results, the ${\chi }^{2}$ value of the model is then given by:

$$\begin{eqnarray}&&{\chi }^{2}={\delta }^{{\rm{T}}}{C}^{-1}\delta .\end{eqnarray} \tag{ 19 }$$

To define the model that fits the BOSS observations best, we want to minimize the ${\chi }^{2}$ value of the model. A two-dimensional space in t₀ and γ is constructed for both of the light curves, and this space is explored using regular grid points. Following the analysis of White et al. (2012), a ${\chi }^{2}$ distribution with nine degrees of freedom is assumed. Using this distribution, the ${\chi }^{2}$ value for a particular model is converted to a confidence interval. An equivalent $n\sigma$ value is computed based on the confidence interval ($1\sigma$ if the enclosed probability is 0.683, $2\sigma$ if it is 0.955, etc.). This statistic demonstrates how "consistent" a particular model is with the BOSS observations.

Figure 3 shows the clustering measurements for several of our well-fitting models compared to the BOSS measurements. The values of these models are given in Table 3. In general, as t₀ increases at a fixed value of γ, the clustering signal increases as well. Compare specifically the L3, L1, and L4 models, which have the same value of γ but have respectively increasing values of t₀. Mathematically, this behavior can be seen from the form of Equation (10): for the same luminosity and mass functions but a larger value of ${f}_{\mathrm{on}}\propto {t}_{q}$, quasars of the same luminosity will shift to more massive host halos. Since the clustering signal increases with the mass, it follows that increasing t₀ will increase the clustering signal. For similar reasons, increasing values of γ for constant values of t₀ are also associated with a stronger clustering signal, since this also effectively increases the quasar lifetime t_q.

Zoom In Zoom Out Reset image size

Figure 4 shows the ${\chi }^{2}$ values in the two-dimensional parameter space t₀ and γ, as defined by Equation (3), for the different light curves. The region of good agreement between the BOSS measurements and our models takes on a linear relationship between $\,{\mathrm{log}}_{10}({t}_{0})$ and γ. Such a relationship can be parameterized as:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({t}_{0}/\mathrm{year})={\mathrm{log}}_{10}({t}_{\mathrm{eff}}/\mathrm{year})+{L}_{0}\gamma .\end{eqnarray} \tag{ 20 }$$

The parameters t_eff and L₀ can be thought of as a characteristic timescale and a characteristic luminosity, respectively. From the functional form of our power-law for quasar lifetime in Equation (3), L₀ can be interpreted as changing the normalization luminosity. This is the luminosity at which all models have the same lifetime, regardless of the value of γ. In other words, the characteristic luminosity of the power law becomes:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({L}_{\mathrm{eff}}/{L}_{\odot })=10-{L}_{0},\end{eqnarray} \tag{ 21 }$$

where L₀ is defined in Equation (20). The parameter t_eff is the characteristic time because all models have this same lifetime at the luminosity L_eff.

Zoom In Zoom Out Reset image size

For the lightbulb model, the best-fit values are $\,{\mathrm{log}}_{10}({t}_{\mathrm{eff}}/\mathrm{year})=7.76$ and $\,{\mathrm{log}}_{10}({L}_{\mathrm{eff}}/{L}_{\odot })=13.29$. (See Table 4 for evolution of these parameters with redshift.) The characteristic luminosity inferred from this value is ${L}_{\mathrm{eff}}={10}^{13.29}\ {L}_{\odot }$, which has a corresponding magnitude of ${M}_{i}=-27.2$. This value is not surprising, given that quasars were selected for the clustering measurements near this magnitude range. More interesting is the value of $\,{\mathrm{log}}_{10}({t}_{\mathrm{eff}}/\mathrm{year})=7.77$, which gives a characteristic lifetime of ${10}^{7.77}=59\,{\rm{Myr}}$. This is a quasar lifetime that is slightly longer than those typically quoted in the literature (Yu & Tremaine 2002; Porciani et al. 2004; Yu & Lu 2004; Conroy & White 2013), which are closer to the Salpeter e-folding timescale or shorter (∼45 Myr for a quasar accreting at Eddington luminosity and a mass conversion efficiency of $\epsilon =0.1$). Although t_eff is slightly higher than these values, it is within a factor of 2.

Table 4.  A List of the Best-fit Parameters for Our Quasar Model as a Function of Redshift

Redshift Selection	zeff	Light Curve	${L}_{\mathrm{eff}}^{* }$ a,b	${t}_{\mathrm{eff}}^{* }$ c
High-zd	2.51	Lightbulb	12.92	7.62
		Exp	12.40	7.14
Fiducial	2.39	Lightbulb	13.29	7.77
		Exp	13.05	7.18
Low-z	2.28	Lightbulb	13.17	7.84
		Exp	13.15	7.29

Notes.

^aL_eff and t_eff as defined in Equations (20)–(21). ^b ${L}_{\mathrm{eff}}^{* }={\mathrm{log}}_{10}({L}_{\mathrm{eff}}/{L}_{\odot })$. ^c ${t}_{\mathrm{eff}}^{* }={\mathrm{log}}_{10}({t}_{\mathrm{eff}}/\mathrm{year})$. ^dThe high-z and low-z samples examine the evolution of these parameters with redshift. See Appendix B for further discussion.

Download table as:  ASCIITypeset image

In the exponential model, the best-fit values for t_eff and L_eff defined in Equation (20) are $\,{\mathrm{log}}_{10}({t}_{\mathrm{eff}}/\mathrm{year})=7.18$ and $\,{\mathrm{log}}_{10}({L}_{\mathrm{eff}}/{L}_{\odot })=13.05$. This luminosity implies a slightly dimmer characteristic luminosity (${M}_{i}=-26.6$). As discussed in Section 2.4, there is not a single τ for all quasars for a given value of t_eff: $L\approx {L}^{* }$ quasars have $\tau \approx {t}_{\mathrm{eff}}$, with brighter quasars having $\tau \gt {t}_{\mathrm{eff}}$. However, the difference between τ and t_eff does not differ by more than a factor of 2 in either direction, and so to a good approximation $\tau \sim {t}_{\mathrm{eff}}$, especially for the luminosity range used to match the clustering measurements. Compared to the lightbulb case, the quasars with an exponential light curve have a shorter characteristic lifetime of 15.1 Myr. The characteristic lifetime is smaller for the exponential than in the lightbulb case because quasars do not shut off entirely after a single lifetime, so the time that a quasar is "bright enough" to be included within the luminosity cuts is longer than its lifetime t_eff. This lifetime is about a third of the Salpeter e-folding timescale, which implies that if quasar light curves are roughly exponential, the combination of the measured QLF and the clustering measurements favors quasars that either radiate at luminosities dimmer than their Eddington ratio ($L/{L}_{\mathrm{edd}}\lt 1$), have a mass-conversion efficiency that is less that the fiducial value ($\epsilon \lt 0.1$), or both. Unfortunately, since our model does not track the underlying physics present, we are not able to distinguish between these two cases.

The reason for the different best-fit values between the two models can be understood as follows. By construction, we have fixed the lifetime of the exponential quasars such that their peak luminosity-to-mass ratio is the same as in the case of the lightbulb for a given choice of t₀ and γ. (See Section 2.4 for more details.) However, the mass-to-light ratios for the two light curves are significantly different. This is due to the fact that the observed luminosity for an exponential quasar can be much smaller than its peak luminosity. A particular luminosity range is selected for the clustering measurements, but the clustering of these quasars is tied to their peak luminosity rather than the observed one. Thus, quasars will tend to have higher clustering at a given luminosity in the exponential case compared to the lightbulb, since they spend comparatively little time at or near their peak luminosity. This luminosity selection includes quasars with a higher peak luminosity than the chosen range (and thus a higher clustering signal), so we must also include quasars that have lower mass hosts to match the average clustering signal. This means that there is a larger spread in host mass compared to the lightbulb case. This behavior explains why the characteristic luminosity is slightly smaller for the exponential model compared to the lightbulb: there is an increased number of low-luminosity quasars occupying high-mass hosts.

Figure 5 shows the range of quasar lifetimes as a function of model parameter γ. The quasar lifetime is broadly similar across different model choices. The exponential model has a lower overall value due to the effect discussed above, i.e., that quasars from a higher peak luminosity will be included in the sample, bringing along a higher clustering signal. Since this is true for nearly all the quasars in the sample, there is an overall decrease in the selected lifetime of quasars. The large difference in the span of quasar lifetimes is due to the way that we have defined the quasar lifetime in the exponential model. As discussed in Section 2.4, the exponential lifetime τ is selected such that the same relationship between host mass and quasar peak luminosity exists in the exponential case as in the lightbulb case. Even in a model where for the lightbulb t_q is independent of L (i.e., when $\gamma =0$), the exponential model parameter τ does have luminosity dependence. In general, quasars with luminosities above L_* will have a lifetime longer than an equivalent luminosity in the lightbulb case for the same choice of t₀ and γ in Equation (3), and those with low luminosities will have a shorter lifetime. This choice for our model leads to the spread in lifetimes of a factor of ∼5, as seen in the case of $\gamma =0$. For models in which $\gamma \gt 0$, there is a widening in the range of values. This is due to the fact that brighter quasars live longer than dimmer ones. Since our choice of quasar lifetime already enforces this difference, these models see an increased effect. Conversely, for $\gamma \lt 0$, there are competing effects between brighter quasars being less long-lived due to the choice of γ, but still simultaneously living longer than their lightbulb counterparts due to the choice of t_q. The latter effect wins out, and these quasars end up having a significantly larger spread than in the lightbulb case. Note that the contours in this figure are smooth compared to Figure 4 because these are results lying along the best-fit line, and the figure shows the range in values rather than a single number (i.e., the ${\chi }^{2}$ value) that fluctuates as a function of position in t₀ and γ.

Zoom In Zoom Out Reset image size

The combination of the QLF and clustering measurements produces an important set of constraints on the space of potential quasar models. Here we investigate the implications of these models. One important implication is the mass of a typical halo for a given quasar luminosity. It is trivial to predict this for the case of a lightbulb model, but less straightforward for the case of the exponential model. Here the peak luminosity L_peak is used to define the mass-to-light ratio, since this is the quantity used in our abundance matching approach. This ratio defines a typical mass for quasars, which can be compared with results of previous analysis (e.g., Martini & Weinberg 2001; Shen et al. 2007; White et al. 2012).

Figure 6 shows the luminosity of quasars as a function of the mass of the host halo for the lightbulb and exponential models for all combinations of t₀ and γ. This plot shows the mass-to-light ratio of the entire catalog. The weight assigned to the luminosity as a function of mass ${L}_{\mathrm{peak}}(M)$ for a particular model i is given by a ${\chi }^{2}$ likelihood. We also find the ±1 and 2σ values that enclose 68% and 95% of the likelihood.

Zoom In Zoom Out Reset image size

Figure 7 shows the mass range as a function of the model parameter γ. Note that the range is essentially constant with respect to γ. By averaging the median mass across all values of γ, a characteristic mass for the two models can be defined. This characteristic mass is $2.5\times {10}^{12}$ ${h}^{-1}\,{M}_{\odot }$ for the lightbulb model, and $2.3\times {10}^{12}$ ${h}^{-1}\,{M}_{\odot }$ for the exponential model. These values are broadly consistent with previous studies of quasar clustering measurements (e.g., Porciani et al. 2004; Croom et al. 2005; Lidz et al. 2006; Porciani & Norberg 2006; White et al. 2012). Since in all models the same clustering signal of the quasars is being selected, there is an implicit requirement for the hosts to lie within a certain mass range. Additionally, the mass range in the exponential case is significantly larger than that of the lightbulb model. This is again related to the fact that due to the light curve, quasars with a higher clustering signal are included within the luminosity sample, and so there must also be lower-mass hosts included as well to balance the average clustering strength. There is a significantly larger spread above the median mass than below. The reason for this asymmetry is due to the difference in number density: since the high-mass objects are rarer, a comparatively smaller range in low-mass halos is necessary to make the clustering signal equivalent to the lightbulb case. See Section 3.3 for further discussion.

Zoom In Zoom Out Reset image size

In the case of the lightbulb model, the halo mass that corresponds to the selected luminosity range of quasars is relatively tightly constrained. For the models that agree with the BOSS measurements at $1\sigma$, the average halo mass ranges from 1.35 × 10¹² to 4.93 × 10¹² ${h}^{-1}\,{M}_{\odot }$ for hosts of quasars within the magnitude cutoff. For the exponential model, there is a much larger range in halo mass: for the collection of models that agree at $1\sigma$, the mass ranges between 6.69 × 10¹² and 5.85 × 10¹³ ${h}^{-1}\,{M}_{\odot }$, almost an entire order of magnitude (compared to about half an order of magnitude for the lightbulb model). Also note that the mass range is much larger than in the lightbulb case. This fact can be explained by noting that there is evolution in the mass-to-luminosity ratio during the lifetime of the quasar. Further, the e-folding time for these models is comparatively long, with typical values being $\tau \approx 40\,{\rm{Myr}}$. This means that there are high-mass hosts included in the sample of quasars chosen for the clustering measurements whose quasars are not at their peak luminosity. Since these hosts have a bias larger than the value preferred by the BOSS measurements, this sample must necessarily include hosts which have a smaller clustering signal, so that on average, the total bias agrees with BOSS.

There is a systematic shift upward in the mass of the exponential case compared to the lightbulb. This shift is related to the difference in parameter space discussed in Section 3.1. Due to their exponential change in luminosity, the quasars are not typically found near their peak luminosities. Thus, even though by construction the peak quasar luminosity as a function of halo mass is the same for the two models, the effective luminosity for a given mass is reduced in the exponential case due to the light curve evolution. In other words, quasars of the same luminosity in the two different models are found in more massive hosts in the exponential case. This leads to a systematic shift in the preferred mass range for clustering measurements. To sum up: the increased spread in halo host mass for the exponential model compared to the lightbulb is due to inclusion of highly biased hosts in the measurement being balanced out by lower-mass ones, and the systematic shift toward higher mass is due to the effective increase in the mass-to-luminosity ratio related to evolution of quasar luminosity.

To observe the effect that different points in parameter space have on the halo host properties, the mass function of halos hosting an active quasar has been calculated for the fiducial redshift selection. (For the high- and low-redshift selections, see Appendix B.) Figure 8 shows the total halo mass function as well as the mass function of halos hosting quasars within the luminosity range $-25\geqslant {M}_{i}\geqslant -27$. From this analysis, the duty cycle of halo hosts can be extracted, i.e., the fraction of active halos divided by the total number of halos. As discussed in Section 2.2, in the lightbulb model the duty cycle can be directly related to the quasar lifetime at that luminosity. However, here the duty cycle is defined simply as the active fraction of halos.

Zoom In Zoom Out Reset image size

Figure 8 compares the case of the lightbulb and the exponential light curves. When the halo mass function of active halos is examined in the two different cases, it can be seen that the mass range of hosts spanned by an individual model is quite different. In the lightbulb case, there is a very small range in host mass compared to the exponential case. This difference can be explained in terms of which hosts are included in the clustering measurements. In the lightbulb case, since the luminosity is constant as a function of quasar lifetime, the only evolution in the relationship between mass and luminosity comes from mass accretion, which makes up a small fraction of total halo mass over the timescales for which quasars are active. As such, with an essentially static relationship, there is a very strong correlation of mass to light. For a specific model, there is only about a factor of 2 in halo mass included for the quasars in the selected magnitude range. When looking at the duty cycle of quasar hosts, one can see that the fraction is typically 0.5%–1%, with little evolution with mass within a model.

Conversely, in the exponential model, there is evolution for individual halos in the mass-to-light relationship. Most importantly, this implies that massive halos will be included when selecting quasars at a specific luminosity. Since they are more highly clustered (and more biased), smaller, less biased halos must also be included in order to create an average bias consistent with the BOSS measurements. This has the effect of extending the mass range of halos included in the mass function. Note that within a single model, there is a much larger span in halo mass: in some cases, the span is more than an order of magnitude in halo mass. Additionally, the duty cycle is comparable in magnitude to the lightbulb case, though slightly smaller: the ratio of active halos to total halos ranges from 0.05%–1%. There also seems to be a trend in the evolution of the duty cycle: there is a central "typical mass" for a given model, and the duty cycle decreases in both directions. A similar trend was found by White et al. (2012).

Note that one result of this measurement is the fact that the mass range of host halos is significantly more extended in the exponential case than the lightbulb case. Thus, one way to break the degeneracy between the lightbulb and exponential models would be to measure the mass range of underlying host halos, perhaps through using gravitational lensing to independently find the mass of the dark matter halo (Courbin et al. 2012). If the range of masses for quasar hosts is extended, then there would be observational evidence favoring an exponential model (or a model with evolution in the quasar light curve) as opposed to the lightbulb model.