The first measurement of the quasar lifetime distribution
Ilya S. Khrykin, Joseph F. Hennawi, Gabor Worseck, Frederick B. Davies
MMNRAS , 1–14 (2021) Preprint 10 February 2021 Compiled using MNRAS L A TEX style file v3.0
The first measurement of the quasar lifetime distribution
Ilya S. Khrykin, ★ Joseph F. Hennawi, , Gábor Worseck, Frederick B. Davies , Kavli Insitute for Physics and Mathematics of the Universe (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan Department of Physics, University of California, Santa Barbara, CA 93106-9530, USA Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA, Leiden, The Netherlands Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Str. 24/25, D-14476 Potsdam, Germany Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany Lawrence Berkeley National Laboratory, 1 Cyclotron Rd, Berkeley, CA 94720, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Understanding the growth of the supermassive black holes (SMBH) powering luminousquasars, their co-evolution with host galaxies, and impact on the surrounding intergalacticmedium (IGM) depends sensitively on the duration of quasar accretion episodes. Unfortu-nately, this time-scale, known as the quasar lifetime, 𝑡 Q , is still uncertain by orders of magnitude( 𝑡 Q (cid:39) .
01 Myr − 𝛼 proximity zones in the absorptionspectra of 𝑧 qso ∼ − ∼
30 Myr. Our recent analysis of 22 archival
Hubble Space Telescope
He ii proximityzone spectra reveals a surprisingly broad range of emission timescales, indicating that somequasars turned on (cid:46) (cid:38)
30 Myr. Determiningthe underlying quasar lifetime distribution (QLD) from proximity zone measurements is achallenging task owing to: 1) the limited sensitivity of individual measurements; 2) randomsampling of the quasar light curves; 3) density fluctuations in the quasar environment; and 4)the inhomogeneous ionization state of He ii in a reionizing IGM. We combine a semi-numericalHe ii reionization model, hydrodynamical simulations post-processed with ionizing radiativetransfer, and a novel statistical framework to infer the QLD from an ensemble of proximity zonemeasurements. Assuming a log-normal QLD, we infer a mean (cid:104) log (cid:0) 𝑡 Q / Myr (cid:1) (cid:105) = . + . − . and standard deviation 𝜎 log 𝑡 Q = . + . − . . Our results allow us to estimate the probabilityof detecting very young quasars with 𝑡 Q ≤ . 𝑝 (≤ . ) = . + . − . , which is broadly consistent with recent determination at 𝑧 ∼ Key words: intergalactic medium – quasars: absorption lines – quasars: general – quasars:supermassive black holes – dark ages, reionization, first stars
Observational efforts of the recent decades resulted in the discoveryof dormant supermassive black holes (with typical masses in therange 𝑀 BH (cid:39) –10 𝑀 (cid:12) ) in the centres of all nearby bulge-dominated galaxies (Soltan 1982; Yu & Tremaine 2002; Kormendy& Ho 2013). It is believed that these are the remnants of the oncepowerful and luminous quasar phase that galaxies went through intheir evolution. As quasars, these black holes played an importantrole in the thermal and ionization evolution of the intergalacticmedium (McQuinn et al. 2009; Compostella et al. 2013; Chardinet al. 2017; D’Aloisio et al. 2017; Khrykin et al. 2017; La Plante et al.2017; Puchwein et al. 2019; Kulkarni et al. 2019; Dayal et al. 2020).In addition, many models of galaxy formation and evolution includevarious quasar feedback prescriptions that are believed to have a ★ E-mail: [email protected] significant impact on their host galaxies, for instance, regulating thestar-formation rate (Hopkins et al. 2006; Choi et al. 2018; Habouzitet al. 2019). However, despite their cosmological importance, manyaspects of quasar evolution remain relatively unconstrained.For example, while it is widely accepted that quasars are pow-ered by accretion on to their supermassive black holes, we do not yetfully understand the physical mechanisms responsible for “poweringup” the accretion. Quasar activity is thought to be triggered by ei-ther major galaxy mergers (e.g. Di Matteo et al. 2005; Springel et al.2005; Di Matteo et al. 2005; Hopkins et al. 2005a,b, 2008; Capeloet al. 2015; Steinborn et al. 2018; Bañados et al. 2019), or by seculardisc instabilities (e.g. Goodman 2003; Hopkins & Quataert 2010,2011; Novak et al. 2011; Bournaud et al. 2011; Gabor & Bournaud2013; Anglés-Alcázar et al. 2013, 2017), or by some combination ofthe two. Another puzzling discovery is the existence of SMBHs withmasses up to 𝑀 SMBH (cid:39) 𝑀 (cid:12) in quasars already at 𝑧 (cid:39) © a r X i v : . [ a s t r o - ph . GA ] F e b Ilya S. Khrykin et al. 2014; Wu et al. 2015; Mazzucchelli et al. 2017; Bañados et al.2018; Shen et al. 2019; Wang et al. 2020; Yang et al. 2020a; Wanget al. 2021). Current theoretical models struggle to explain suchhigh masses existing already less than ∼ 𝑡 Q , with a rangeof uncertain estimates to date covering several orders of magni-tude (0 .
01 Myr (cid:46) 𝑡 Q (cid:46) 𝑡 dc . Measurements at 𝑧 (cid:39) 𝑡 dc (cid:39) − 𝑡 dc (cid:39) −
100 Myr also come from the comparison of the integratedquasar luminosity function to the present day SMBH number den-sity (Soltan 1982; Yu & Tremaine 2002). Most importantly, thesemethods only constrain the integrated quasar lifetime, and not theduration of individual accretion episodes 𝑡 Q , which theoretical mod-els suggest could be much shorter (Ciotti & Ostriker 2001; Novaket al. 2011; Anglés-Alcázar et al. 2017; Angles-Alcazar et al. 2020).The impact of the quasar’s ionizing radiation on the surround-ing IGM provides an independent observational probe of the quasarlifetime 𝑡 Q . For example, the time delay between variations in theradiation field and the environment’s response to these changes wasused to argue for quasar variability on short 𝑡 Q (cid:39) . 𝛼 forest in the spectra of 𝑧 (cid:39) Γ . This equilibration time-scale 𝑡 eq (cid:39) Γ − setsthe upper limit on the duration of episodic quasar activity that canbe probed by their proximity zones (Khrykin et al. 2016; Eilers et al.2017; Davies et al. 2020). Moreover, the extent of the quasar prox-imity zones is actually sensitive only to the quasar on-time , 𝑡 on ≤ 𝑡 Q ,as illustrated in Fig. 1. Imagine that the quasar light is detected byan observer at the time 𝑡 obs , then the quasar on-time is defined suchthat the quasar accretion episode began at a time 𝑡 obs − 𝑡 on in thepast. However, this quasar episode may continue and end at a latertime 𝑡 Q , which can be recorded on Earth only if the observer couldconduct observations of the same quasar in the future. For instance,the H i proximity effect at 𝑧 (cid:39) 𝑡 on ≥ 𝑡 eq (cid:39) .
03 Myr, owing to the high post-reionization H i background radiation ( Γ HI (cid:39) − s − ; Becker& Bolton 2013). However, the situation is qualitatively different t = 0 t = t obs t = t Q Timequasar lifetime t Q quasar on − time t on t on ∈ [0; t Q ] (observe quasar) Figure 1.
Illustration of the difference between the quasar lifetime 𝑡 Q andthe lower limit on 𝑡 Q , the quasar on-time 𝑡 on , inferred from the analysis ofindividual quasar proximity zones. at 𝑧 (cid:39) 𝛼 forest far away fromquasars (Fan et al. 2006; Eilers et al. 2018; Bosman et al. 2018;Yang et al. 2020b) provides a better contrast, and analysis of thequasar H i proximity zones at 𝑧 (cid:39) 𝑡 on (cid:39) . 𝑧 (cid:38) 𝑡 on (cid:39) . + . − . Myr (Moreyet al. prep). This is consistent with finding newly turned on quasarswith 𝑡 on ∼ . ∼ 𝛼 absorption spectra toward 𝑧 ∼ 𝛼 forest at 𝑧 (cid:39) (cid:39)
30 Myr owingto the He ii ionizing background at 𝑧 (cid:39) 𝑡 S (cid:39)
45 Myr; Salpeter 1964), the characteristic time for the SMBHevolution. In Khrykin et al. (2019, hereafter Paper I) we presentedthe first fully Bayesian statistical method to infer the quasar on-timesof individual quasars from their He ii proximity zones at 𝑧 (cid:39)
4, andmeasured short lifetimes of order 1 Myr for several quasars. We re-cently expanded our inference to the full sample of He ii-transparentquasars, discovering a surprisingly broad distribution of individualon-times 𝑡 on , ranging from (cid:46) (cid:38)
30 Myr (Worseck et al.2021, hereafter Paper II). The nature of such wide scatter is yetunknown, but it might indicate that the extent of the He ii proximityzones samples a broad underlying distribution of quasar lifetimes.In this work, we introduce a fully Bayesian statistical algorithmto recover the intrinsic quasar lifetime distribution. We re-analysethe sample of He ii proximity zones from Paper I and Paper IIwith similar numerical modelling but with a modified Bayesianapproach to measure for the first time the underlying quasar lifetimedistribution. In contrast to previous analyses of He ii proximity zones(Syphers & Shull 2014; Zheng et al. 2015) that often ignored theapparent degeneracy between the quasar lifetime and the ionizationstate of the IGM (Khrykin et al. 2016), our algorithm instead takesthis degeneracy into account and allows us to fully marginalizeover the poorly constrained ionization state of He ii in the IGM.Moreover, for individual quasars 𝑡 on is often poorly constrainedowing to large errors in proximity zone sizes resulting from largequasar redshift errors (Paper I; Paper II). Our Bayesian approachtakes these errors into account and optimally combines these weak MNRAS000
30 Myr (Worseck et al.2021, hereafter Paper II). The nature of such wide scatter is yetunknown, but it might indicate that the extent of the He ii proximityzones samples a broad underlying distribution of quasar lifetimes.In this work, we introduce a fully Bayesian statistical algorithmto recover the intrinsic quasar lifetime distribution. We re-analysethe sample of He ii proximity zones from Paper I and Paper IIwith similar numerical modelling but with a modified Bayesianapproach to measure for the first time the underlying quasar lifetimedistribution. In contrast to previous analyses of He ii proximity zones(Syphers & Shull 2014; Zheng et al. 2015) that often ignored theapparent degeneracy between the quasar lifetime and the ionizationstate of the IGM (Khrykin et al. 2016), our algorithm instead takesthis degeneracy into account and allows us to fully marginalizeover the poorly constrained ionization state of He ii in the IGM.Moreover, for individual quasars 𝑡 on is often poorly constrainedowing to large errors in proximity zone sizes resulting from largequasar redshift errors (Paper I; Paper II). Our Bayesian approachtakes these errors into account and optimally combines these weak MNRAS000 , 1–14 (2021) uasar lifetime distribution z qso R p z [ M p c ] z qso -1.0-0.50.00.51.01.52.0 l og ( t o n / M y r) Figure 2.
Left:
The sizes of the He ii proximity zone in the spectra of 22 quasars used in this work as a function of redshift. The error bars correspond to themeasured redshift uncertainty (see Table 1 for details).
Right:
Inferred values of quasar on-times from Paper I and Paper II. The blue dots with error bars showthe measurements, while the red and grey arrows indicate 1 𝜎 lower and upper limits, respectively. constraints from individual quasars to infer the underlying lifetimedistribution.This paper is organized in the following way. In Section 2we outline the main properties of the quasars in our sample. Wesummarize our hydrodynamical and radiative transfer simulationsin Section 3, and discuss their results in Section 4. In Section 5we describe our statistical algorithm for measuring the lifetimedistribution and present the results of our inference from MarkovChain Monte Carlo (MCMC). We discuss our findings in Section 6and conclude in Section 7.Throughout this work, we assume a flat Λ CDM cosmologywith dimensionless Hubble constant ℎ = . Ω 𝑚 = . Ω 𝑏 = . 𝜎 = .
8, and 𝑛 𝑠 = .
96, the helium mass fraction is 𝑌 He = .
24, consistent with the latest
Planck results (Planck Collaborationet al. 2018). All quoted distances are in units of proper Mpc if notstated otherwise. Throughout this work we quote all quasar lifetimesand on-times in the form log ( 𝑡 / Myr ) . Our parent sample of 24 2 . (cid:46) 𝑧 qso (cid:46) . Hubble Space Telescope ( HST ) has been described in Paper Iand Paper II, to which we refer the interested reader. As in our previ-ous work, we exclude two quasars with formally negative proximityzone sizes due to associated absorption (HE 2347 − − 𝑧 (cid:39) 𝑅 pz asthe proper distance from the quasar location where smoothed He iitransmission profile drops below the 10 per cent threshold for thefirst time. We adopt a Gaussian filter with FWHM = − . Among the remaining quasars withno [O iii] detection, 6 have broad Mg ii lines with corresponding un-certainties Δ 𝑣 =
273 km s − ; 3 quasars have neither [O iii] nor Mg iilines covered, but 𝐻𝛽 is present in the 𝐾 -band spectra, with uncer-tainties of Δ 𝑣 =
400 km s − . Finally, the remaining 7 quasars do nothave high-quality near-IR spectra, and their redshifts are determinedfrom the C iv line in the optical band. We take into account the lu-minosity dependent blueshift of the C iv line (the Baldwin effect;Baldwin 1977), and find the high uncertainty of Δ 𝑣 =
656 km s − calibrated against the quasars with known redshift (Paper I; Pa-per II). We provide the identified redshifts and associated redshiftuncertainties in Table 1. We refer the interested reader to a de-tailed description of the redshift measurement procedure providedin Paper II.In Paper I we introduced a Bayesian method for measuringthe on-times of individual quasars based on a statistical comparisonbetween the observed and simulated sizes of the He ii proximityzones. We successfully applied this method to measure on-timesof individual He ii-transparent quasars at 2 . (cid:46) 𝑧 (cid:46) . 𝑡 on , whereas in whatfollows we focus on the intrinsic quasar lifetime distribution. In this study, we utilize the set of simulations previously used inPaper I and Paper II, although with modifications. In what follows,we briefly outline the main parameters of our models, and referthe interested reader to our previous papers for a more detaileddescription of the simulations.The impact of quasar radiation on the IGM is modelled through
MNRAS , 1–14 (2021)
Ilya S. Khrykin
Table 1.
Sample of 22 quasars used in this work. From left to right, thecolumns show: quasar name, quasar redshift with measured uncertainty,base-10 logarithm of the He ii ionizing photon production rate, measuredsize of the He ii proximity zone with corresponding 1 𝜎 uncertainty, and theinferred quasar on-time 𝑡 on (Paper I; Paper II). Quasar 𝑧 qso log 𝑄 𝑅 pz 𝑡 on s − Mpc Myr
HS 1700 + . ± . .
34 7 . ± .
01 0 . + . − . HS 1024 + . ± . .
62 9 . ± .
97 9 . + . − . Q 1602 + . ± . .
80 6 . ± .
97 1 . + . − . PC 0058 + . ± . .
19 7 . ± . > . SDSS J0936 + . ± . .
49 8 . ± .
15 11 . + . − . SDSS J0818 + . ± . .
38 2 . ± . < . HE2QS J2157 + . ± . .
72 17 . ± . > . SDSS J1237 + . ± . .
27 1 . ± . < . HE2QS J1706 + . ± . . − . ± . < . HE2QS J2149 − . ± . . − . ± . < . Q 0302 − . ± . .
89 13 . ± . > . HE2QS J0233 − . ± . .
47 4 . ± . < . HS 0911 + . ± . .
74 4 . ± .
19 1 . + . − . HE2QS J0916 + . ± . .
45 3 . ± . < . SDSS J1253 + . ± . .
48 11 . ± . > . SDSS J2346 − . ± . .
79 2 . ± .
77 0 . + . − . HE2QS J2311 − . ± . .
66 1 . ± . < . SDSS J1137 + . ± . .
19 4 . ± . > . HE2QS J1630 + . ± . .
92 8 . ± .
02 2 . + . − . SDSS J1614 + . ± . .
14 2 . ± . < . SDSS J1711 + . ± . .
19 2 . ± . < . SDSS J1319 + . ± . .
82 3 . ± .
98 0 . + . − . two stages. First, we utilize the outputs of cosmological hydrody-namical simulations performed with the Gadget-3 code (2 × particles and 𝐿 box = ℎ − cMpc; Springel 2005) and extract1000 realistic 1D density, velocity and temperature distributions(to which we refer as skewers ) from 𝑧 sim = [ . , . , . ] snap-shots using periodic boundary conditions . The resulting skewershave a physical length of 160 cMpc, and a spatial resolution ofd 𝑟 = .
012 comoving Mpc (d 𝑣 ≈ . − ).Next, we post-process the extracted skewers with a custom 1Dradiative transfer algorithm based on the C -Ray code (Mellemaet al. 2006), which tracks the evolution of H i, He ii, 𝑒 − , and gastemperature (Khrykin et al. 2016, 2017). For each quasar in oursample (see Table 1) we create a set of radiative transfer modelsdepending on the quasar redshift 𝑧 qso , quasar photon productionrate 𝑄 , assumed quasar on-time 𝑡 on , and initial fraction of He ii 𝑥 HeII , (set by the He ii ionizing background prior to quasar activity).We compute the 𝑄 given the observed 𝑖 -band magnitudes ofthe quasars in our sample (see Paper I and Paper II for details),and together with observed 𝑧 qso fix them for each quasar that wesimulate. In order to better capture the cosmological evolution of thedensity field between the redshift of the simulation output 𝑧 sim andthe redshift of quasars 𝑧 qso in our sample, we rescale the simulateddensities 𝜌 sim in every pixel along the skewers as We used 𝑧 sim = . 𝑧 sim = . 𝑧 qso > . 𝑧 sim = . 𝑧 qso ≤ . Starting from the location of the quasars, we create skewers by castingrays through the simulation volume at random angles, and traversing the boxmultiple times. We use the periodic boundary conditions to wrap a skewerthrough the box along the chosen direction. R simpz [Mpc] p (cid:0) R s i m p z (cid:1) log ( t on / Myr) = 0.000 x HeII , = 1.00 R simpz R simpz + σ R datapz N ( µ = 3 . , σ = 1 . R datapz = 3.9 Mpc Figure 3.
Size distribution of the simulated He ii proximity zones forquasar SDSS J1319 + 𝑥 HeII , = .
00. The grey histogram illustrates the 𝑅 simpz distributionwithout redshift uncertainties, while the black one incorporates them (seediscussion in Section 4). The blue dashed line marks the size of the observedHe ii proximity zone 𝑅 datapz for the quasar in question, while the red curveillustrates the Gaussian fit to the distribution. 𝜌 (cid:0) 𝑧 qso (cid:1) = 𝜌 sim × (cid:0) + 𝑧 qso (cid:1) ( + 𝑧 sim ) (1)In contrast to the previously used simulations (Paper I;Paper II), we extended the upper limit of the logarithmi-cally spaced quasar on-time values from log ( 𝑡 on / Myr ) = .
000 to log ( 𝑡 on / Myr ) = . ( 𝑡 on / Myr ) = [− . , . ] with step Δ log ( 𝑡 on / Myr ) = . 𝑥 HeII , = .
01. This allows us tobetter capture the ionization state of the IGM at redshifts 𝑧 qso (cid:46) . 𝑥 HeII , ∈{ . , . , . , . , . , . , . , . , . , . } .This results in a grid of 410 radiative transfer models per quasar(41 values of quasar on-time and 10 values of initial He ii fraction)with 1000 He ii Ly 𝛼 transmission spectra per model, 9020 modelsin total for 22 quasars. We smooth the simulated He ii Ly 𝛼 transmission spectra and mea-sure the sizes of the proximity zones 𝑅 simpz following the same al-gorithm applied to observational data (see Section 2), which resultsin the distribution of 1000 𝑅 simpz for each model. Note, becausethe adopted smoothing scale is larger than observed spectral res-olution and spans several pixels, the resulting extent of the He iiproximity zones only weakly depends on the observational effectslike the signal-to-noise ratio, or spectral resolution (see discussionin Appendix A of Paper II). Therefore, we do not include theseeffects in our numerical simulations. Fig. 3 shows an example ofthe distribution of He ii proximity zone sizes for one model of MNRAS000
01. This allows us tobetter capture the ionization state of the IGM at redshifts 𝑧 qso (cid:46) . 𝑥 HeII , ∈{ . , . , . , . , . , . , . , . , . , . } .This results in a grid of 410 radiative transfer models per quasar(41 values of quasar on-time and 10 values of initial He ii fraction)with 1000 He ii Ly 𝛼 transmission spectra per model, 9020 modelsin total for 22 quasars. We smooth the simulated He ii Ly 𝛼 transmission spectra and mea-sure the sizes of the proximity zones 𝑅 simpz following the same al-gorithm applied to observational data (see Section 2), which resultsin the distribution of 1000 𝑅 simpz for each model. Note, becausethe adopted smoothing scale is larger than observed spectral res-olution and spans several pixels, the resulting extent of the He iiproximity zones only weakly depends on the observational effectslike the signal-to-noise ratio, or spectral resolution (see discussionin Appendix A of Paper II). Therefore, we do not include theseeffects in our numerical simulations. Fig. 3 shows an example ofthe distribution of He ii proximity zone sizes for one model of MNRAS000 , 1–14 (2021) uasar lifetime distribution l o g ( t o n / M y r ) -2.0-1.5-1.0-0.50.00.51.01.52.02.53.0 x H e II , R s i m p z [ M p c ] -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0log ( t on / Myr)0102030 R s i m p z [ M p c ] x HeII , = 0 . x HeII , R s i m p z [ M p c ] log ( t on / Myr) = 2 . Figure 4.
Left: Evolution of the simulated He ii proximity zone sizes 𝑅 simpz with incorporated redshift uncertainties for quasar SDSS J1319 + 𝑅 simpz distributions, while the shaded area illustrates the1 𝜎 variation. Right: 2D slices through the 3D plot on the left. The upper panel shows the evolution of the 𝑅 simpz distribution parameters as a function ofquasar on-time with fixed 𝑥 HeII , = .
05, while the bottom panel shows the same evolution but as a function of initial He ii fraction with fixed quasar on-timelog ( 𝑡 on / Myr ) = . the quasar SDSS J1319 + 𝑅 simpz dis-tribution random errors drawn from the Gaussian distribution withstandard deviation equal to the observed redshift uncertainty, whichfor the case of SDSS J1319 + 𝜎 𝑧 = . 𝜎 𝑅 datapz = .
98 Mpc ( 𝜎 𝑣 =
400 km s − ; see Table 1).Each 𝑅 simpz distribution is well described by a Gaussian, whichis apparent from the red line in Fig. 3. For each quasar and eachof the 410 models parametrized by { log ( 𝑡 on / Myr ) , 𝑥 HeII , } , wefit the resulting 𝑅 simpz (with incorporated redshift uncertainties) witha Gaussian by determining its mean and standard deviation. Thisallows us to easily track the evolution of the distribution of prox-imity zone sizes as a function of the model parameters, i.e, quasaron-time log ( 𝑡 on / Myr ) and initial He ii fraction 𝑥 HeII , (see Sec-tion 3). We illustrate an example of this evolution for the modelsof the quasar SDSS J1319 + 𝑡 eq (Khrykinet al. 2016), at which point the 𝑅 simpz saturates and sensitivity to theon-time diminishes as 𝑅 simpz stops growing with increasing on-time.On the other hand, as Paper I showed, the size of the He ii prox-imity zone is largely insensitive to the initial He ii fraction. Thisarises from the competition between the decrease in the level ofIGM transmission for large initial 𝑥 HeII , fractions, and the increasein transmission due to the thermal proximity effect for high 𝑥 HeII , fractions (Khrykin et al. 2017). Moreover, the very definition ofthe proximity zone size 𝑅 pz , i.e., the location where transmissiondrops below the 10 per cent threshold, reduces sensitivity to theHe ii background (which sets the 𝑥 HeII , ). This is because at 10per cent transmission levels the quasar radiation still dominates thetotal photoionization rate, and the effect of the He ii ionizing back- ground on the transmission becomes prominent at larger distancescorresponding to the transmission levels much lower than 10 percent.We use bivariate spline interpolation to predict the meanand standard deviation of the 𝑅 simpz distributions for values oflog ( 𝑡 on / Myr ) and 𝑥 HeII , at locations between our simulated pa-rameter grid points. We use the dependence of the 𝑅 simpz distributionson the model parameters in the likelihood calculations in the nextsection. In this section, we describe our algorithm for measuring the param-eters of the quasar lifetime distribution. In what follows, we assumethat all inferred on-times of individual quasars (see Table 1) aredrawn from the underlying statistical distribution. Moreover, we as-sume that the quasar lifetime distribution (to which we further referto as QLD) is described by a log-normal distribution 𝑝 (cid:0) 𝑡 Q (cid:1) = 𝑡 Q log 𝑒𝜎 √ 𝜋 × exp (cid:34) − (cid:0) log (cid:0) 𝑡 Q / Myr (cid:1) − 𝜇 (cid:1) 𝜎 (cid:35) , (2)where 𝜇 = (cid:104) log (cid:0) 𝑡 Q / Myr (cid:1) (cid:105) and 𝜎 = 𝜎 log 𝑡 Q are the mean andstandard deviation of the distribution, respectively. In what followswe show how we can infer these parameters. In order to estimate the parameters of the QLD we need to determinethe Bayesian likelihood function for the observed proximity zonesizes L qso (cid:16) 𝑅 datapz ,𝑖 | 𝜇, 𝜎 (cid:17) for each quasar in our sample. MNRAS , 1–14 (2021)
Ilya S. Khrykin
We begin by constructing a grid of parameter values on whichthe likelihood is estimated. For the mean of the QLD 𝜇 we adoptthe range 𝜇 = (cid:104) log (cid:0) 𝑡 Q / Myr (cid:1) (cid:105) = [− . , . ] with step Δ 𝜇 = . 𝜎 we choose 𝜎 = 𝜎 log 𝑡 Q = [ . , . ] dex, with step Δ 𝜎 = .
1. We obtain thevalues of the individual likelihood L qso (cid:16) 𝑅 datapz ,𝑖 | 𝜇, 𝜎 (cid:17) for a singlequasar in our sample at each point on the parameter grid as follows:(i) Given the values of the mean 𝜇 and standard deviation 𝜎 of theQLD we randomly draw 1000 quasar lifetime values log (cid:0) 𝑡 Q / Myr (cid:1) from the resulting log-normal distribution given by eq. (2).(ii) In order to account for the fact that our observations measurethe quasar on-times, 𝑡 on , which represent a random sampling of thequasars light-bulb light curve, we construct a distribution of 𝑡 on times for each value of log (cid:0) 𝑡 Q / Myr (cid:1) drawn from the QLD. Forthat we assume a uniform distribution of on-times 𝑡 on ∼ U (cid:0) , 𝑡 Q (cid:1) (see Fig. 1). We draw 10 𝑡 on values for each 𝑡 Q value, which resultsin 10000 𝑡 on values per location on the QLD parameter grid ( 𝜇, 𝜎 ).(iii) For each 𝑡 on value we need to find the corresponding distri-bution of the proximity zone sizes in our simulations. The 𝑅 simpz dis-tributions, however, also depend on the value of initial He ii fraction(see Fig. 4), i.e., 𝑅 simpz (cid:0) log ( 𝑡 on / Myr ) , 𝑥 HeII , (cid:1) . In order to takethis into account we use the semi-numerical model of the He ii ion-izing background, Γ HeII , from Davies et al. (2017), which provides1000 realizations of Γ HeII at each quasar redshift in our sample.Assuming a density equal to the mean density of the Universe at theredshift in question and that the IGM is in ionization equilibriumwith the predicted Γ HeII , we calculate the distribution of 1000 𝑥 HeII values at this redshift. For each value of the on-time we then ran-domly draw a value of 𝑥 HeII , from this distribution, which results in10000 { log ( 𝑡 on / Myr ) , 𝑥 HeII , } pairs per evaluation of the QLD.By doing so, we effectively perform a Monte-Carlo sampling of theHe ii fraction and marginalize out the unknown ionization state ofthe IGM. Finally, we use bivariate spline interpolation to find themean and standard deviation of the corresponding Gaussian 𝑅 simpz distribution at each parameter location { log ( 𝑡 on / Myr ) , 𝑥 HeII , } (see Section 4).(iv) Using the Gaussian fits at these 10000 { log ( 𝑡 on / Myr ) , 𝑥 HeII , } locations, we draw 500 𝑅 simpz sam-ples from the Gaussian fit at each parameter value and concatenatethem into one combined distribution of the He ii proximity zonesizes corresponding to all { log ( 𝑡 on / Myr ) , 𝑥 HeII , } values. Thisprocedure has then effectively marginalized over the stochasticrelationship between 𝑡 on and 𝑡 Q , as well as stochastic distributionof 𝑥 HeII , in the IGM at each redshift.(v) Finally, we apply Kernel Density Estimation (KDE) to findthe continuous probability density function (PDF) of this combined 𝑅 simpz distribution. We calculate the value of the likelihood by eval-uating the KDE PDF at the observed value of the He ii proximityzone size 𝑅 datapz for the quasar in question (see Table 1);(vi) we repeat steps (i)-(v) for each combination of the QLDparameters.This procedure results in 990 determinations of the likelihoodfunction on the { 𝜇, 𝜎 } parameter grid for each quasar. Next, we ob-tain a joint likelihood for all quasars in our sample simply by takingthe product of the individual (independent) likelihoods, calculatedfollowing the procedure stated above L joint = 𝑓 (cid:16) 𝑅 datapz (cid:17) = 𝑁 qso (cid:214) 𝑖 = L qso (cid:16) 𝑅 datapz ,𝑖 | 𝜇, 𝜎 (cid:17) , (3) σ log t Q -2.0-1.00.01.02.0 h l og ( t Q / M y r) i h log t Q i = 0 . +0 . − . σ log t Q = 0 . +0 . − . -2.0 -1.0 0.0 1.0 2.0 h log ( t Q / Myr) i Figure 5.
Results of the MCMC inference on the mock quasar sample. Thebottom left panel shows the 2D contours, while to top left and bottom rightpanels illustrate the constrained marginalized posterior probabilities of themean 𝜇 and standard deviation 𝜎 of the mock QLD. The red dot shows thevalues of the input QLD parameters where 𝑁 qso =
22 is the number of the quasars in the sample.Finally, we use a bivariate spline to interpolate the L joint giveneq. (3) for any combination of parameter values between the gridpoints in our parameter space. We now can sample this joint like-lihood with MCMC to estimate the posterior distributions of theQLD mean 𝜇 and standard deviation 𝜎 . However, first, we need todetermine the priors for each parameter. We assume a flat logarithmic prior on the mean of the QLD, i.e., 𝜇 = (cid:104) log (cid:0) 𝑡 Q / Myr (cid:1) (cid:105) = [− . , . ] . As discussed in Paper I,the choice of this range is motivated by the current estimates ofthe quasar lifetime (see Martini 2004, for a review). The lowerlimit is set by the existence of the line-of-sight proximity effect inthe H i Ly 𝛼 forest, because it requires the quasar to shine at leastfor 𝑡 eq (cid:39) .
03 Myr (Bajtlik et al. 1988; Khrykin et al. 2016). Onthe other hand, the upper limit is chosen to lie in the upper rangeof estimates provided by the analysis of the quasar clustering andmodels of the SMBH growth (Shen et al. 2007; White et al. 2008;Eftekharzadeh et al. 2015).For the standard deviation of the QLD we also adopt a flatuniform prior in range 𝜎 = [ . , . ] dex. Although other ob-servations do not provide much guidance on what range should beadopted, it is motivated by the large scatter found in the presentestimates of the quasar lifetime and by the variety of methods.However, we note that combined with the prior on the meanof the log-normal distribution, this range allows quasar lifetimeslonger than 100 Myr, which can be problematic for our modellingstrategy. Namely, the post-processing radiative transfer approachthat we use assumes that cosmic structures do not evolve over thetime-scales of the calculations. We also neglect cooling of the gasdue to cosmological expansion and recombination (Hui & Gnedin MNRAS000
03 Myr (Bajtlik et al. 1988; Khrykin et al. 2016). Onthe other hand, the upper limit is chosen to lie in the upper rangeof estimates provided by the analysis of the quasar clustering andmodels of the SMBH growth (Shen et al. 2007; White et al. 2008;Eftekharzadeh et al. 2015).For the standard deviation of the QLD we also adopt a flatuniform prior in range 𝜎 = [ . , . ] dex. Although other ob-servations do not provide much guidance on what range should beadopted, it is motivated by the large scatter found in the presentestimates of the quasar lifetime and by the variety of methods.However, we note that combined with the prior on the meanof the log-normal distribution, this range allows quasar lifetimeslonger than 100 Myr, which can be problematic for our modellingstrategy. Namely, the post-processing radiative transfer approachthat we use assumes that cosmic structures do not evolve over thetime-scales of the calculations. We also neglect cooling of the gasdue to cosmological expansion and recombination (Hui & Gnedin MNRAS000 , 1–14 (2021) uasar lifetime distribution σ log t Q -2.0-1.00.01.02.0 h l og ( t Q / M y r) i h log ( t Q / Myr) i = 0 . +0 . − . σ log t Q = 0 . +0 . − . -2.0 -1.0 0.0 1.0 2.0 h log ( t Q / Myr) i -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0 log ( t Q / Myr) p ( l og ( t Q / M y r)) Figure 6.
Results of the MCMC inference on the observed sample of quasars in Table 1. Similar to Fig. 5, the left-hand side panels show inferred 2D contoursand marginalized posterior probabilities of the QLD parameters. The right-hand side panel illustrates 100 realizations of the QLD based on random drawsof the mean 𝜇 and standard deviation 𝜎 from the MCMC chains. The solid red curve corresponds to the maximum likelihood values of mean and standarddeviation. 𝑧 (cid:39) (cid:39) . ( 𝑡 on / Myr ) = .
000 (see Section 3) and allow the wideprior range on 𝜎 log t Q . Our reasoning for this is that He ii proximityzones are anyway not sensitive to lifetimes 𝑡 on (cid:38) 𝑡 eq (cid:39)
30 Myr at 𝑧 (cid:39) 𝑅 pz saturates at 𝑡 on ∼ 𝑡 eq and does not significantly increase all the wayup to 𝑡 on ∼ 𝑡 on (cid:38)
100 Myr.
Given the expression for the joint likelihood in eq. (3), and the in-terpolation procedure which allows us to evaluate this likelihood atany point in our parameter space, we can now sample this likelihoodusing MCMC and obtain the posterior probability distribution forour model parameters. In what follows we use the publicly avail-able affine-invariant MCMC sampling algorithm emcee (Foreman-Mackey et al. 2013).
First, to test the accuracy of our inference procedure, we apply itto a mock sample of 𝑁 qso =
22 quasars (the same number as inour observed sample listed in Table 1). In order to create the mocksample of He ii proximity zones we perform the following steps.We begin by choosing the mock QLD parameters. We choosethe mean 𝜇 = .
40 and standard deviation 𝜎 = .
50 to describethe mock QLD, in anticipation of the results in the next section.From this distribution, we randomly draw 𝑁 qso =
22 values of the quasar lifetime log (cid:0) 𝑡 Q / Myr (cid:1) . Using the uniform probability dis-tribution for the on-times given the quasar lifetime (see discussionin Section 5.1), we draw one 𝑡 on value for each log (cid:0) 𝑡 Q / Myr (cid:1) .This results in the sample of 𝑁 qso =
22 on-times, i.e., one for eachmock quasar. Further, we assume that quasars in our mock samplemap to the real quasars in our dataset, and assign them a corre-sponding redshift 𝑧 qso and photon production rate consistent withTable 1. We use these redshifts to assign a value of the initial He iifraction 𝑥 HeII , to each quasar by randomly drawing a value fromthe corresponding 𝑥 HeII distribution at redshift 𝑧 qso similar to ourprocedure in Section 5.1. Next, for each quasar in this mock sam-ple, characterized by the values of on-time and He ii fraction, wefind a corresponding distribution of the He ii proximity zone sizes 𝑅 simpz (cid:0) log ( 𝑡 on / Myr ) , 𝑥 HeII , (cid:1) from our radiative transfer simula-tions. We randomly draw one value of 𝑅 mockpz for each quasar fromthe corresponding 𝑅 simpz distribution, which now corresponds to the“observed" sample of the mock He ii proximity zone sizes. In orderto reproduce the accuracy of the observed He ii proximity zoneswe also add the redshift associated uncertainties to each 𝑅 mockpz , forwhich we adopt the values of 𝜎 R PZ listed in Table 1.We then proceed to calculate the joint likelihood of the mocksample of quasars following the discussion in Section 5.1, and sam-ple it with MCMC. Fig. 5 illustrates the results of this inference.The dark blue (light blue) contours correspond to the 68 per cent(95 per cent) confidence regions, while the marginalized posteriorPDFs of the mean and standard deviation of the QLD are shown bythe histograms. The red dot shows the values of input QLD parame-ters. We quote the 50 th percentiles of the 1D marginalized posteriorprobability distributions as the measured values of the inferred pa-rameters, whereas their uncertainties are derived from the 16 th and84 th percentiles.It is apparent from Fig. 5 that our inference algorithm suc-cessfully recovers the input parameters of the mock QLD with ≈ . − . MNRAS , 1–14 (2021)
Ilya S. Khrykin to apply our method to the observed sample of quasars in the nextSection.
The corner plot in the left panel of Fig. 6 shows the result of theMCMC inference of the QLD parameters { 𝜇, 𝜎 } for the 20 quasarsin our observed sample (see Table 1).We measure the mean of the QLD 𝜇 = . + . − . , and theQLD standard deviation 𝜎 = . + . − . from the correspondingmarginalized posterior distributions shown by the histograms in thecorner plot of Fig. 6. We further illustrate our findings in the rightpanel of Fig. 6, where 200 realizations of the QLD based on randomdraws from MCMC samples of the mean and standard deviation areshown. The solid red curve corresponds to the maximum likelihoodvalues of the mean ( 𝜇 (cid:39) .
25) and standard deviation ( 𝜎 (cid:39) . ≈ . − . 𝜎 (cid:39) . 𝑡 on from the corresponding log-normaldistribution of quasar lifetimes 𝑡 Q . In what follows, we examine how the redshift uncertainties affectthe constraining power of our analysis, and discuss our findingsin the context of previous measurements of the quasar on-timesand lifetimes, as well as the significance of our results for SMBHevolution models.
In Paper I we showed that large uncertainties in the estimated quasarsystemic redshifts, 𝑧 qso , modify the distributions of the simulatedHe ii proximity zone sizes, making them wider (see the comparisonbetween the grey and black histograms in Fig. 3), thus signifi-cantly weakening individual constraints on the quasar on-times 𝑡 on .Therefore, for individual measurements, it would be ideal to obtainsystemic redshifts from sub-millimeter observations of CO and/or[C ii] 158 𝜇 m lines, decreasing the uncertainty to Δ 𝑣 (cid:39)
100 km s − (Eilers et al. 2020). In principle, the redshift uncertainties mightalso degrade the constraints on the QLD parameters. In order toquantify the significance of redshift precision on the results of ourinference we perform the following test.Analogous to the discussion in Section 5.4, we begin by creat-ing a mock sample of 22 quasars with the same photon productionrates and redshifts as our observed catalog (see Table 1). For thevalues of the QLD model that we simulate, we choose the medianbest-fitting values of the QLD parameters inferred in Section 5.5, i.e, 𝜇 = log (cid:0) 𝑡 Q / Myr (cid:1) = .
22 and 𝜎 log 𝑡 Q = .
80. In what followswe assume that the same grid of radiative transfer models used forobserved quasar sample (see Section 3) corresponds to this mockdataset.Similarly to the discussion in Section 4 (see also Paper I and σ log t Q -2.0-1.00.01.02.0 h l og ( t Q / M y r) i σ v = 800 km s − h log t Q i = 0 . +0 . − . σ log t Q ≤ σ v = 100 km s − h log t Q i = 0 . +0 . − . σ log t Q = 0 . +0 . − . -2.0 -1.0 0.0 1.0 2.0 h log ( t Q / Myr) i Figure 7.
Results of the MCMC inference on the mock quasars sample.Similar to Fig. 5, the bottom left panel shows the 2D contours, while to topleft and bottom right panels illustrate the constrained marginalized posteriorprobabilities of the QLD parameters. The filled contours (and blue posteriorhistograms) show the result for the “ideal" case of redshift uncertainties( 𝜎 𝑣 =
100 km s − ), while the black curves correspond to the case with“bad" redshift uncertainties ( 𝜎 𝑣 =
800 km s − ). The red dot indicates theinput values of the QLD parameters. Paper II), the impact of the redshift uncertainties is modelled byadding random Gaussian distributed proximity zone size errors withstandard deviation 𝜎 𝑅 pz to the simulated distributions of the He iiproximity zone sizes of each quasar. In order to test the sensitivityof our inference algorithm to the redshift uncertainty, we considertwo cases: 1) “bad" - the redshift uncertainties for each quasar inthe mock sample are comparable to the biggest uncertainties in ourobserved sample (see discussion in Section 2) and correspond tothe velocity precision 𝜎 𝑣 =
800 km s − , and 2) “ideal", with theredshift uncertainties corresponding to 𝜎 𝑣 =
100 km s − , which isthe precision expected from future sub-millimeter observations (itis also comparable to the values inferred from the analysis of the[O iii] line used to measure systemic redshifts of several quasars inour sample). The value of 𝜎 𝑅 pz for each mock quasar is then foundvia 𝜎 𝑅 pz = 𝜎 𝑣 / 𝐻 (cid:0) 𝑧 qso (cid:1) , where 𝐻 (cid:0) 𝑧 qso (cid:1) is the Hubble parameterat the redshift of a mock quasar.Finally, following the discussion in Section 5.1, we performthe calculations on these noisy distributions and compute the jointlikelihood function given by eq. (3) for these two cases. We sampleeach of these likelihoods with MCMC and compare the resultsin Fig. 7. It is apparent from the constraints in Fig. 7 that thedifference between the two cases is not very significant, implyinga weak dependence on the precision with which we infer the QLDparameters on the level of redshift uncertainty. Although this resultmight be surprising in light of the individual measurements of 𝑡 on (Paper I; Paper II) where redshift errors play a large role, this weakdependence on redshift precision arises from the interplay of severalfactors. Our inference algorithm presented in Section 5 is subject toseveral sources of stochasticity that directly influence the resultingsensitivity to the QLD parameters. We illustrate the effect of each MNRAS000
100 km s − , which isthe precision expected from future sub-millimeter observations (itis also comparable to the values inferred from the analysis of the[O iii] line used to measure systemic redshifts of several quasars inour sample). The value of 𝜎 𝑅 pz for each mock quasar is then foundvia 𝜎 𝑅 pz = 𝜎 𝑣 / 𝐻 (cid:0) 𝑧 qso (cid:1) , where 𝐻 (cid:0) 𝑧 qso (cid:1) is the Hubble parameterat the redshift of a mock quasar.Finally, following the discussion in Section 5.1, we performthe calculations on these noisy distributions and compute the jointlikelihood function given by eq. (3) for these two cases. We sampleeach of these likelihoods with MCMC and compare the resultsin Fig. 7. It is apparent from the constraints in Fig. 7 that thedifference between the two cases is not very significant, implyinga weak dependence on the precision with which we infer the QLDparameters on the level of redshift uncertainty. Although this resultmight be surprising in light of the individual measurements of 𝑡 on (Paper I; Paper II) where redshift errors play a large role, this weakdependence on redshift precision arises from the interplay of severalfactors. Our inference algorithm presented in Section 5 is subject toseveral sources of stochasticity that directly influence the resultingsensitivity to the QLD parameters. We illustrate the effect of each MNRAS000 , 1–14 (2021) uasar lifetime distribution R simpz [Mpc]0.00.20.40.60.81.0 p (cid:0) R s i m p z (cid:1) log ( t on / Myr) = 0.25 x HeII , = 0.10 t on fixed R simpz [Mpc] log ( t Q / Myr) = 0.25 x HeII , = 0.10Effect of t Q − t on sampling R simpz [Mpc] h log ( t Q / Myr) i = 0.25 σ log t Q = 0.80 x HeII , = 0.10 x HeII , fixed σ v = 800 km s − σ v = 100 km s − R simpz [Mpc]0.00.10.2 log ( t Q / Myr) = 0.25 x HeII , = Const t Q fixed R simpz [Mpc] h log ( t Q / Myr) i = 0.25 σ log t Q = 0.80 x HeII , = ConstFinal PDF Figure 8.
Different sources of stochasticity in our inference algorithm. The black and blue histograms illustrate two cases of redshift uncertainties correspondingto velocity precision of 𝜎 𝑣 =
800 km s − (“bad") and 𝜎 𝑣 =
100 km s − (“ideal"), respectively. The panels in the top row illustrate: left: the effect of theredshift uncertainties on the simulated He ii proximity zone sizes distribution in the models of a single quasar HE2QS J0233 − middle: the effect ofsampling the 𝑡 Q − 𝑡 on relation on the joint 𝑅 pz distribution (see Section 5.1 for details); right: the effect of sampling the assumed log-normal distribution ofquasar lifetimes on joint 𝑅 pz distribution for a single realization of the QLD. The value of He ii fraction is fixed at 𝑥 HeII , = .
10. The bottom middle panelshows the combined effects of sampling the 𝑡 Q − 𝑡 on relation and variations in the assumed He ii fraction on the joint 𝑅 pz distribution. Finally, the bottom right panel illustrates the same effect as shown in the top right panel, but combined with variations in the He ii fraction. of these sources in Fig. 8 and elaborate on them separately in whatfollows. We note that for simplicity, the discussion that followsfocuses on a single quasar, HE2QS J0233 − 𝑧 qso = . − ( 𝑡 on / Myr ) = .
25 and 𝑥 HeII , = .
10. It is apparent thatthe 𝑅 simpz distributions are significantly different depending on themagnitude of the redshift uncertainty, i.e, the distribution is muchwider in the “bad" case characterized by the velocity precision 𝜎 𝑣 =
800 km s − . In Paper I we illustrated that incorporating largerredshift uncertainties decreases the constraining power of individ-ual proximity zone measurements. Therefore, individual measure-ments of quasar on-times depend significantly on the accuracy ofthe quasar redshifts (Paper I; Paper II, Eilers et al. 2020).However, the situation is different if one tries to infer the quasarlifetime distribution. As we showed in Section 5.1, the measuredindividual quasar on-times represent random sampling of the quasarlight bulb light curve. Mathematically, this sampling is describedby the uniform distribution 𝑡 on ∼ U (cid:2) , 𝑡 Q (cid:3) . The middle panelin the top row of Fig. 8 labelled "Effect of 𝑡 Q − 𝑡 on sampling"shows the distributions of He ii proximity zone sizes constructed by sampling the single quasar lifetime of log (cid:0) 𝑡 Q / Myr (cid:1) = .
25 forthe "bad" and "ideal" redshift scenarios. Similar to the discussion inSection 5.1, we randomly draw 1000 values of 𝑡 on from U (cid:2) , 𝑡 Q (cid:3) ,and for each 𝑡 on value and fixed value of 𝑥 HeII , = .
10, we find thecorresponding 𝑅 simpz distribution from our radiative transfer modelgrid. We then draw 500 𝑅 simpz values from each Gaussian fit andcombine the resulting proximity zone sizes in one joint distributionshown in the middle panel. It is clear that 𝑡 Q − 𝑡 on sampling broadensthe resulting 𝑅 simpz distribution irrespective of the assumed redshiftuncertainty, in contrast with case of fixed 𝑡 on (illustrated in the topleft panel of Fig. 8) where the differences between "bad" and "ideal"redshifts is larger.The foregoing example only captures the stochasticity due to 𝑡 on sampling from a single quasar lifetime. In reality, there is anunderlying QLD, which, as we assumed in Section 5 is describedby a log-normal distribution (see eq. (2)). The top right panel inFig. 8 demonstrates the combined effect of the QLD sampling andthe previously described 𝑡 Q − 𝑡 on sampling on the resulting joint 𝑅 simpz distributions. For this exercise we choose one QLD realizationwith the mean 𝜇 = .
25 and standard deviation 𝜎 = .
80, and keepthe He ii fraction fixed at 𝑥 HeII , = .
10 for all models. As onecan see, the 𝑅 simpz distributions flatten out and become less peaked, MNRAS , 1–14 (2021) Ilya S. Khrykin further reducing the difference between the two cases of redshiftuncertainties.The situation is further aggravated when variations in thepoorly known He ii fraction 𝑥 HeII , in the surrounding IGM arealso taken into account. In the middle panel of the bottom rowof Fig. 8, we show the same 𝑡 Q − 𝑡 on sampling for a single valueof quasar lifetime as before (top middle panel), but now includ-ing variations in 𝑥 HeII , . As described in Section 5.1, the value of 𝑥 HeII , is randomly drawn from the He ii fraction PDF from oursemi-numerical model at the redshift of quasar in question. It is ap-parent that the 𝑅 simpz distributions are significantly broadened whencompared to the case of fixed He ii fraction in the upper middlepanel. The difference between two cases of redshift uncertainties isalmost completely mitigated.Finally, the right panel in the bottom row of Fig. 8 illustratesthe combined impact of all sources of stochasticity on the resultingjoint 𝑅 simpz distributions for the two cases of redshift error precision.It is evident that the two redshift uncertainty scenarios are basicallyindistinguishable. Consequently, the likelihood function evaluatedfrom these distributions will be practically identical (see Section 5.1for details), which explains the large degree of overlap between theposterior probability distributions of the inferred QLD parametersfor these two cases shown in Fig. 7.According to the results presented in Fig. 8, we conclude thatour inference algorithm only weakly depends on the exact quasarredshift uncertainties. This is because these redshift errors areswamped by a much larger scatter produced by the combinationof several other sources of stochasticity (i.e. 𝑡 Q − 𝑡 on sampling,the QLD itself, and fluctuations in 𝑥 HeII , ). It is, therefore, possi-ble to obtain tight constraints on the QLD parameters (see resultsin Section 5.5) even using quasar samples with poorly determinedredshifts. Note, however, that as one can see from the top left panelof Fig. 8, good precision on measured quasar redshifts is never-theless required to accurately infer individual quasar on-times (seealso Paper I; Paper II). This is especially important for the studyof young quasars with exceptionally small proximity zones (Eilerset al. 2020). In order to directly compare the results of our inference presentedin Section 5.5 to the individual measurements of the on-times 𝑡 on from the quasar He ii proximity zones (Paper I; Paper II), we needto convert the inferred QLD into the distribution of on-times.Assuming that the QLD is described by eq. (2), we derive ananalytical formula for the 𝑝 (cid:0) log ( 𝑡 on / Myr ) (cid:1) – the probability den-sity function for the quasar on-times log ( 𝑡 on / Myr ) as a functionof the mean and standard deviation of the QLD – given by 𝑝 (cid:0) log ( 𝑡 on / Myr ) (cid:1) = − 𝜇 − / − 𝜇 √ (cid:20) 𝜎 ln (cid:21) · Erfc (cid:20) log ( 𝑡 on / Myr ) − 𝜇 + 𝜎 ln10 𝜎 √ (cid:21) · ln10 𝑡 on , (4)where Erfc is the complementary error function, and 𝜇 and 𝜎 are themean and standard deviation of the QLD, respectively. We presenta full derivation of this analytical solution in Appendix A, to whichwe refer the interested reader.The top panel of Fig. 9 illustrates the resulting PDFs com-puted using eq. (4). The grey lines show 100 realizations of the 𝑝 (cid:0) log ( 𝑡 on / Myr ) | 𝜇, 𝜎 (cid:1) for 100 pairs of the { 𝜇, 𝜎 } values ran- p ( l og ( t o n / M y r)) M e a s u r e m e n t s -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0log ( t on / Myr) U pp e r / L o w e r li m i t s Figure 9.
Distribution of quasar on-times 𝑝 ( log ( 𝑡 on / Myr )) . The greycurves in the top panel show 100 realizations of the 𝑝 ( log ( 𝑡 on / Myr )) calculated using eq. (4), based on the values of QLD mean 𝜇 and standarddeviation 𝜎 illustrated in the right panel of Fig. 6. The solid red curve illus-trates the 𝑝 ( log ( 𝑡 on / Myr )) distribution corresponding to the maximumlikelihood values of 𝜇 and 𝜎 . The middle panel shows the measured on-timesof individual quasars, while the bottom panel illustrates the upper/lower lim-its on log ( 𝑡 on / Myr ) (see Table 1). The violins show the correspondingposterior distributions inferred from the analysis of individual quasars (seePaper I and Paper II for details). domly drawn from the MCMC chains (see Section 5.3 for details).The solid red line displays the quasar on-time probability distribu-tion function corresponding to the maximum likelihood values ofthe QLD parameters (see red curve in the right panel of Fig. 6), andcalculated using eq. (4). We compare these results to the on-timeconstraints of individual quasars (see Table 1) in the middle andbottom panels of Fig. 9.It is apparent from the top panel of Fig. 9 that our analyticalsolution for the PDF of the quasar on-times has a tail towards lowvalues of 𝑡 on , which results from the uniform sampling of 𝑡 on given 𝑡 Q , i.e., 𝑡 on ∼ U (cid:2) , 𝑡 Q (cid:3) . The existence of this tail may explain thediscovery of very young quasars at 𝑧 qso (cid:39)
6, assuming that quasarsat higher redshift follow the same distribution. Indeed, Eilers et al.(2020) analysed the H i Ly 𝛼 proximity zones in the spectra of 153 MNRAS000
6, assuming that quasarsat higher redshift follow the same distribution. Indeed, Eilers et al.(2020) analysed the H i Ly 𝛼 proximity zones in the spectra of 153 MNRAS000 , 1–14 (2021) uasar lifetime distribution p ( ≤ . . +0 . − . p ( ≤ .
01 Myr) = 0 . +0 . − . -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0log ( t on / Myr)0.00.51.0 p ( ≤ l og ( t o n / M y r)) Figure 10.
The CDF of the quasar on-times. The grey curves are the CDFrealizations corresponding to 100 realizations of the 𝑝 ( log ( 𝑡 on / Myr )) ,illustrated in Fig. 9, while the red curve is the CDF corresponding to themaximum likelihood values of the QLD parameters. We also quote theestimated fractions of quasars with on-time shorter than 0 . .
01 Myr. quasars at 𝑧 qso (cid:39) 𝑡 on ≤ . 𝑝 (≤ . ) (cid:39) . .
10. We can estimatethe young quasar fraction from the inferred distributions of 𝑡 on timesby integrating 𝑝 (≤ . ) = ∫ − −∞ 𝑝 (cid:0) log ( 𝑡 on / Myr ) (cid:1) d log ( 𝑡 on / Myr ) , (5)where 𝑝 (cid:0) log ( 𝑡 on / Myr ) (cid:1) is given by eq. (4). Fig. 10 shows the re-sulting cumulative probability distributions (CDF) calculated fromthe same realizations of the 𝑡 on PDF illustrated in the top panel ofFig. 9. Applying eq. (5) to ∼ 𝑡 on PDFs (based on total number of { 𝜇, 𝜎 } pairs from MCMC sam-ples) we infer 𝑝 (≤ . ) = . + . − . . This is a factor of 2–3higher than the value estimated by Eilers et al. (2020) who found 𝑝 (≤ . ) = . − .
10. However, we emphasize that the valuequoted by Eilers et al. (2020) is likely a lower limit due to possibleincompleteness of their sample, and thus our estimate is broadlyconsistent. Future analysis of a larger sample of observed He iiproximity zones should help to refine this estimate. In addition, theproperties of the underlying QLD might vary with redshift, whichwill also affect the estimated fraction of young quasars.Finally, we note that the short quasar on-times and quasarlifetimes implied by our inferred QLD (log (cid:0) 𝑡 Q / Myr (cid:1) = . + . − . and 𝜎 log 𝑡 Q = . + . − . ; see Section 5.3) appear to be at odds withmeasurements obtained from the analysis of the quasar clustering.Indeed, Shen et al. (2007) estimated the fraction of the age of theUniverse that a galaxy hosts an active quasar, i.e., the duty cycle 𝑓 dc ,at 𝑧 ≥ . 𝑓 dc (cid:39) −
60% (see also White et al. 2008, Shankaret al. 2010). This measurement translates into the broad range ofintegrated quasar lifetimes 𝑡 dc (cid:39) 𝑓 dc 𝑡 H ( 𝑧 ) (cid:39) −
600 Myr, where 𝑡 H ( 𝑧 ) is the age of the Universe at redshift 𝑧 . These time-scales are ≈ −
200 times longer than the lifetimes suggested by our results(see also He et al. 2018). If quasars are to radiate their energy in onlyone continuous "light-bulb" burst lasting for 𝑡 Q = 𝑡 dc (cid:39) 𝑓 dc 𝑡 H ( 𝑧 ) ,then there is no simple way to resolve the apparent discrepancybetween our results and clustering measurements. However, it needsto be noted that the relationship between the estimated clustering strength and the inferred quasar lifetimes is a subject of significantuncertainty due to our lack of knowledge of how quasars populatethe dark matter haloes (Shen et al. 2007; White et al. 2008; He et al.2018). The short quasar lifetimes that our findings point to pose a significantchallenge to the models of SMBH formation and evolution. Indeed,in the simplest model of constant accretion at the Eddington limit,the SMBH mass follows an exponential growth law given by 𝑀 SMBH = 𝑀 seed · exp (cid:18) 𝑡 − 𝑡 seed 𝑡 S (cid:19) , (6)where 𝑀 seed is the initial mass of the SMBH at time 𝑡 seed , 𝑡 S (cid:39) · 𝜆 − 𝜖 /( − 𝜖 ) Myr is the 𝑒 -folding or Salpeter time-scale, 𝜆 = 𝐿 bol / 𝐿 edd is the Eddington ratio, and 𝜖 is the radiative efficiency ofaccretion which we set to the canonical value from general relativity(Thorne 1974) of 𝜖 = . 𝑀 seed (cid:39) 𝑀 (cid:12) , then, according to eq. (6), SMBH must con-stantly accrete for a total of (cid:39)
16 Salpeter 𝑒 -folding time-scales or (cid:39)
730 Myr in order to grow the 𝑀 SMBH (cid:39) M (cid:12) SMBH massesmeasured for our sample of quasars (Paper II). It is well known thatthis poses a problem for quasars at 𝑧 (cid:39) − < (cid:38) 𝑀 (cid:12) SMBHs (Mortlock et al. 2011; Bañados et al.2018; Wang et al. 2020; Yang et al. 2020a; Wang et al. 2021) fromthe stellar-mass seeds.At fact value, the ∼
730 Myr required to grow a ∼ 𝑀 (cid:12) appears inconsistent with both the measured short on-times of indi-vidual quasars and the inferred QLD, which also favors much shorterquasar lifetimes (see Fig. 6). However, we note that He ii proximityzones are not sensitive to this entire growth history of (cid:39)
16 Salpetertimes given that the SMBH is much smaller and the quasar is muchfainter in the distant past. Only the last e-folding time has an influ-ence on the extent of the He ii proximity zone. As such, the mostfair comparison of our assumed "light-bulb" light curve to this ex-ponential light curve is to approximately commensurate 𝑡 Q with the 𝑒 -folding time-scale 𝑡 S , i.e., by simply equating the total area underthese respective curves, that is ∫ 𝐿 ( 𝑡 ) d 𝑡 = 𝑡 S 𝐿 ≈ 𝑡 Q 𝐿 , where 𝐿 is the luminosity of the quasar at the time (redshift) at which weobserve the quasar. This would seem to imply the Salpeter times of 𝑡 𝑆 (cid:39) 𝜖 (cid:39) . 𝑧 (cid:39) 𝑀 SMBH (cid:39) 𝑀 (cid:12) in less than1 Gyr after the Big Bang.However, application of the Soltan argument (Soltan 1982) thatrelates the total UV emission from quasars to the local SMBH pop-ulation, seems to be consistent with 𝜖 = . 𝑧 ∼ MNRAS , 1–14 (2021) Ilya S. Khrykin big of a problem since 𝑧 (cid:39) 𝑧 (cid:39) −
4, and thus our results appear inconflict with the Soltan argument – the short 𝑒 -folding time-scaleand corresponding low radiative efficiency required by proximityzones would seem to imply that there exists a large mass densityof unseen SMBHs in the local universe. There are two potentialmechanisms to resolve this apparent conflict, namely one can in-voke significant periods of UV-obscured growth or more complexlight curves, which we now discuss in turn.With regards to obscuration, suppose that the ≈ 𝜖 = . 𝑡 S =
45 Myr, the obscured to unobscured ratio at comparable lumi-nosity would then need to be 𝑡 S / 𝑡 Q . This would thus imply about50 obscured quasars to every unobscured one. However, observa-tions at 𝑧 (cid:39) 𝑧 ∼ − (cid:39)
50 times larger could have gone unnoticed, implyingthat predominantly obscured black hole growth does not obviouslyresolve the aforementioned tension with the Soltan argument.On the other hand, SMBH growth might not obey the simplecontinuous Eddington-limited exponential evolution described byeq. (6). Indeed, more complex flickering light curves, characterizedby short time-scale variations in the quasar continuum radiation(Ciotti & Ostriker 2001; Novak et al. 2011), could manage to growthe observed SMBHs and remain consistent with our proximity zoneconstraints. For the sake of illustration, let us assume a simple flick-ering light curve whereby quasars turn on and shine for a constantperiods equal to the light bulb 𝑡 Q drawn from the QLD, but then turnoff for extended periods denoted by 𝑡 off . Even for this simple on-offlight curve the space of possibilities is large. We, however, focus onthe simplest variant, which is that the quasars shine as light bulbsfor 𝑡 Q and then are off for 𝑡 off = 𝑡 eq , and that this behavior continuesfor the entire Hubble time 𝑡 H ( (cid:39) . 𝑧 (cid:39) 𝑡 off allows the proximity zone sizes to decay to zero, which occurson the equilibration time-scale 𝑡 eq as shown by Davies et al. (2020)for analogous H i proximity zones. In this picture the effective dutycycle for UV-luminous phases would be 𝑡 dc = 𝑡 Q / 𝑡 off · 𝑡 H , which for 𝑡 Q = −
10 Myr and 𝑡 off = 𝑡 eq ≈
30 Myr implies 𝑡 dc (cid:39) −
500 Myr.Are duty cycles in this range inconsistent with clustering constraintson the duty cycle and SMBH growth? For the former, we note thatthese duty cycles overlap the range implied by quasar clustering,which allow for duty cycles in the range 30 −
600 Myr (Shen et al.2007; White et al. 2008; Shankar et al. 2010; He et al. 2018). Regard-ing SMBH growth, it appears that this would not be enough time togrowth the SMBHs given the standard model of 𝑡 S =
45 Myr, sinceroughly 16 Salpeter times are required, or about 730 Myr.There thus appears to be a tension between the inferred quasarlifetimes and the SMBH evolution models, although we emphasizethat: 1) we have assumed that this on/off behavior also applies toearly phases of SMBH growth when the quasars are much fainter,which our proximity zone observations do not constrain; 2) we haveconsidered a single value of 𝑡 on / 𝑡 off whereas the space of possiblecombinations is large, and light curves could surely be significantly more complex; 3) we have not considered the effects of obscuredgrowth phases, which current constraints on obscured populationssuggest could modify the light curve math at the factor of a fewlevel (but not significantly more).Further careful work on this question is clearly warranted, andwe conclude by noting that it is quite possible that a similar SMBHgrowth problem that has been touted at 𝑧 (cid:38) 𝑧 ∼ · 𝑡 S ≈
730 Myr is required to grow SMBHs, whichseems straightforward given the ≈ . 𝑧 ∼ We have used the measured He ii proximity zone sizes of 𝑁 qso = . (cid:46) 𝑧 qso (cid:46) . 𝑡 Q for the first time,marginalized over the ionization state of He ii in the surroundingintergalactic medium. The main results of our work are as follows:(i) Assuming a log-normal distribution of quasar lifetimes, weinferred the mean of the distribution (cid:104) log (cid:0) 𝑡 Q / Myr (cid:1) (cid:105) = . + . − . ,and the standard deviation 𝜎 log 𝑡 Q = . + . − . (see Figure 6).(ii) We presented a simple analytical expression for the distribu-tion of the quasar on-times 𝑡 on . We showed that according to thisdistribution, the probability of finding quasars with 𝑡 on ≤ . 𝑝 (≤ . ) = . + . − . , broadly consistent with the valuesfound by Eilers et al. (2020) based on the analysis of the H i prox-imity zones at 𝑧 qso (cid:39) 𝑧 qso (cid:39) . .
2, significantlyimproving the precision of our inference algorithm. Finally, weemphasize that there also exists a large sample ( (cid:39) 𝑧 qso (cid:39) MNRAS000
2, significantlyimproving the precision of our inference algorithm. Finally, weemphasize that there also exists a large sample ( (cid:39) 𝑧 qso (cid:39) MNRAS000 , 1–14 (2021) uasar lifetime distribution such analysis opens up the possibility of constraining the redshiftevolution of the quasar lifetime distribution parameters. ACKNOWLEDGEMENTS
IKS thanks Metin Ata, Valeri Vardanyan, and Anna-Christina Eil-ers for useful discussions. Kavli IPMU is supported by World Pre-mier International Research Center Initiative (WPI), MEXT, Japan.JFH acknowledges support for this work that was provided byNASA through grant numbers HST-AR-15014.003-A and HST-GO-16318.002-A from the Space Telescope Science Institute, which isoperated by AURA, Inc., under NASA contract NAS 5-26555.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author.
REFERENCES
Anderson S. F., Hogan C. J., Williams B. F., Carswell R. F., 1999, AJ, 117,56Anglés-Alcázar D., Özel F., Davé R., 2013, ApJ, 770, 5Anglés-Alcázar D., Faucher-Giguère C.-A., Quataert E., Hopkins P. F., Feld-mann R., Torrey P., Wetzel A., Kereš D., 2017, MNRAS, 472, L109Angles-Alcazar D., et al., 2020, arXiv e-prints, p. arXiv:2008.12303Bañados E., et al., 2018, Nature, 553, 473Bañados E., et al., 2019, ApJ, 881, L23Bajtlik S., Duncan R. C., Ostriker J. P., 1988, ApJ, 327, 570Baldwin J. A., 1977, ApJ, 214, 679Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023Bosman S. E. I., Fan X., Jiang L., Reed S., Matsuoka Y., Becker G., HaehneltM., 2018, MNRAS, 479, 1055Bournaud F., Dekel A., Teyssier R., Cacciato M., Daddi E., Juneau S.,Shankar F., 2011, ApJ, 741, L33Capelo P. R., Volonteri M., Dotti M., Bellovary J. M., Mayer L., GovernatoF., 2015, MNRAS, 447, 2123Carswell R. F., Whelan J. A. J., Smith M. G., Boksenberg A., Tytler D.,1982, MNRAS, 198, 91Cen R., Safarzadeh M., 2015, ApJ, 798, L38Chardin J., Puchwein E., Haehnelt M. G., 2017, MNRAS, 465, 3429Choi E., Somerville R. S., Ostriker J. P., Naab T., Hirschmann M., 2018,ApJ, 866, 91Ciotti L., Ostriker J. P., 2001, ApJ, 551, 131Compostella M., Cantalupo S., Porciani C., 2013, MNRAS, 435, 3169Conroy C., White M., 2013, ApJ, 762, 70D’Aloisio A., Upton Sanderbeck P. R., McQuinn M., Trac H., Shapiro P. R.,2017, MNRAS, 468, 4691Dall’Aglio A., Wisotzki L., Worseck G., 2008, A&A, 491, 465Davies F. B., Furlanetto S. R., Dixon K. L., 2017, MNRAS, 465, 2886Davies F. B., Hennawi J. F., Eilers A.-C., 2019, ApJ, 884, L19Davies F. B., Hennawi J. F., Eilers A.-C., 2020, MNRAS, 493, 1330Dayal P., et al., 2020, MNRAS, 495, 3065De Rosa G., et al., 2014, ApJ, 790, 145Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604Eftekharzadeh S., et al., 2015, MNRAS, 453, 2779Eilers A.-C., Davies F. B., Hennawi J. F., Prochaska J. X., Lukić Z., Maz-zucchelli C., 2017, ApJ, 840, 24Eilers A.-C., Hennawi J. F., Davies F. B., 2018, ApJ, 867, 30Eilers A.-C., et al., 2020, ApJ, 900, 37Fan X., et al., 2006, AJ, 132, 117Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125,306Gabor J. M., Bournaud F., 2013, MNRAS, 434, 606 Goodman J., 2003, MNRAS, 339, 937Green J. C., et al., 2012, ApJ, 744, 60Habouzit M., et al., 2019, MNRAS, 484, 4413Haiman Z., Hui L., 2001, ApJ, 547, 27He W., et al., 2018, PASJ, 70, S33Hennawi J. F., Prochaska J. X., 2007, ApJ, 655, 735Hogan C. J., Anderson S. F., Rugers M. H., 1997, AJ, 113, 1495Hopkins P. F., Hernquist L., 2009, ApJ, 698, 1550Hopkins P. F., Quataert E., 2010, MNRAS, 407, 1529Hopkins P. F., Quataert E., 2011, MNRAS, 415, 1027Hopkins P. F., Hernquist L., Martini P., Cox T. J., Robertson B., Di MatteoT., Springel V., 2005a, ApJ, 625, L71Hopkins P. F., Hernquist L., Cox T. J., Di Matteo T., Martini P., RobertsonB., Springel V., 2005b, ApJ, 630, 705Hopkins P. F., Hernquist L., Cox T. J., Di Matteo T., Robertson B., SpringelV., 2006, ApJS, 163, 1Hopkins P. F., Hernquist L., Cox T. J., Kereš D., 2008, ApJS, 175, 356Hui L., Gnedin N. Y., 1997, MNRAS, 292, 27Inayoshi K., Visbal E., Haiman Z., 2019, arXiv e-prints, p. arXiv:1911.05791Khrykin I. S., Hennawi J. F., McQuinn M., Worseck G., 2016, ApJ, 824,133Khrykin I. S., Hennawi J. F., McQuinn M., 2017, ApJ, 838, 96Khrykin I. S., Hennawi J. F., Worseck G., 2019, MNRAS, 484, 3897Kormendy J., Ho L. C., 2013, ARA&A, 51, 511Kulkarni G., Worseck G., Hennawi J. F., 2019, MNRAS, 488, 1035La Plante P., Trac H., Croft R., Cen R., 2017, ApJ, 841, 87Makan K., Worseck G., Davies F. B., Hennawi J. F., Prochaska J. X., RichterP., 2020, arXiv e-prints, p. arXiv:2012.07876Martini P., 2004, in Ho L. C., ed., Coevolution of Black Holes and Galaxies.p. 169 ( arXiv:astro-ph/0304009 )Martini P., Weinberg D. H., 2001, ApJ, 547, 12Mazzucchelli C., et al., 2017, ApJ, 849, 91McQuinn M., Lidz A., Zaldarriaga M., Hernquist L., Hopkins P. F., DuttaS., Faucher-Giguère C.-A., 2009, ApJ, 694, 842Mellema G., Iliev I. T., Alvarez M. A., Shapiro P. R., 2006, NA, 11, 374Merloni A., et al., 2014, MNRAS, 437, 3550Morey K., Eilers A. C., et al., in prepMortlock D. J., et al., 2011, Nature, 474, 616Novak G. S., Ostriker J. P., Ciotti L., 2011, ApJ, 737, 26Planck Collaboration et al., 2018, arXiv e-prints, p. arXiv:1807.06209Polletta M., Weedman D., Hönig S., Lonsdale C. J., Smith H. E., Houck J.,2008, ApJ, 675, 960Puchwein E., Haardt F., Haehnelt M. G., Madau P., 2019, MNRAS, 485, 47Regan J. A., Downes T. P., Volonteri M., Beckmann R., Lupi A., TrebitschM., Dubois Y., 2019, MNRAS, 486, 3892Salpeter E. E., 1964, ApJ, 140, 796Schawinski K., Koss M., Berney S., Sartori L. F., 2015, MNRAS, 451, 2517Schmidt T. M., Worseck G., Hennawi J. F., Prochaska J. X., Crighton N.H. M., 2017, ApJ, 847, 81Schmidt T. M., Hennawi J. F., Worseck G., Davies F. B., Lukić Z., OñorbeJ., 2018, ApJ, 861, 122Shankar F., Weinberg D. H., Miralda-Escudé J., 2009, ApJ, 690, 20Shankar F., Weinberg D. H., Shen Y., 2010, MNRAS, 406, 1959Shen Y., et al., 2007, AJ, 133, 2222Shen Y., et al., 2009, ApJ, 697, 1656Shen Y., et al., 2019, ApJ, 873, 35Smith A., Bromm V., Loeb A., 2017, Astronomy and Geophysics, 58, 3.22Soltan A., 1982, MNRAS, 200, 115Springel V., 2005, MNRAS, 364, 1105Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776Steinborn L. K., Hirschmann M., Dolag K., Shankar F., Juneau S., KrumpeM., Remus R.-S., Teklu A. F., 2018, MNRAS, 481, 341Syphers D., Shull J. M., 2014, ApJ, 784, 42Thorne K. S., 1974, ApJ, 191, 507Ueda Y., Akiyama M., Hasinger G., Miyaji T., Watson M. G., 2014, ApJ,786, 104Venemans B. P., et al., 2013, ApJ, 779, 24Vito F., et al., 2018, MNRAS, 473, 2378MNRAS , 1–14 (2021) Ilya S. Khrykin
Volonteri M., 2010, A&ARv, 18, 279Volonteri M., 2012, Science, 337, 544Wang F., et al., 2020, ApJ, 896, 23Wang F., et al., 2021, arXiv e-prints, p. arXiv:2101.03179White M., Martini P., Cohn J. D., 2008, MNRAS, 390, 1179White M., et al., 2012, MNRAS, 424, 933Worseck G., Davies F. B., Hennawi J. F., Prochaska J. X., 2019, ApJ, 875,111Worseck G., Khrykin I. S., Hennawi J. F., Prochaska J. X., Farina E. P.,2021, arXiv e-prints, p. arXiv:2101.01196Wu X.-B., et al., 2015, in IAU General Assembly. p. 2251223Yang J., et al., 2020a, arXiv e-prints, p. arXiv:2006.13452Yang J., et al., 2020b, ApJ, 904, 26Yu Q., Tremaine S., 2002, MNRAS, 335, 965Zheng W., Syphers D., Meiksin A., Kriss G. A., Schneider D. P., York D. G.,Anderson S. F., 2015, ApJ, 806, 142
APPENDIX A: ANALYTICAL SOLUTION FORDISTRIBUTION OF ON-TIMES
In what follows we derive the analytical solution for the probabil-ity density function of the logarithm of the quasar on-times 𝑡 on , 𝑝 (cid:0) log 𝑡 on (cid:1) . Given the definition of the on-times, they follow theuniform distribution (see Section 2), and the conditional probabilityof 𝑡 on given the quasar lifetime 𝑡 Q is then can be written as 𝑝 (cid:0) 𝑡 on | 𝑡 Q (cid:1) = (cid:40) 𝑡 Q , ≤ 𝑡 on ≤ 𝑡 Q , 𝑡 on < , 𝑡 on > 𝑡 Q . (A1)On the other hand, we assumed a base-10 log-normal distribu-tion of quasar lifetimes 𝑡 Q , for which the probability density functionis given by eq. (2) 𝑝 (cid:0) 𝑡 Q (cid:1) = 𝑡 Q log e 𝜎 √ 𝜋 exp (cid:34) − (cid:0) log (cid:0) 𝑡 Q / Myr (cid:1) − 𝜇 (cid:1) 𝜎 (cid:35) , (A2)where 𝜇 = (cid:104) log (cid:0) 𝑡 Q / Myr (cid:1) (cid:105) and 𝜎 = 𝜎 log t Q are mean and stan-dard deviation of QLD distribution. The rule of total probabilitydictates that 𝑝 ( 𝑡 on ) = ∫ 𝑝 (cid:0) 𝑡 on | 𝑡 Q (cid:1) 𝑝 (cid:0) 𝑡 Q (cid:1) d 𝑡 Q . (A3)Taking into account eq. (A1) and eq. (A2), eq. (A3) becomes 𝑝 ( 𝑡 on ) = ∫ ∞ 𝑡 on 𝑡 Q · log e 𝑡 Q 𝜎 √ 𝜋 exp (cid:34) − (cid:0) log (cid:0) 𝑡 Q / Myr (cid:1) − 𝜇 (cid:1) 𝜎 (cid:35) d 𝑡 Q . (A4)Assuming 𝑢 = log 𝑡 Q and changing the variables, we canre-write eq. (A4) as follows 𝑝 ( 𝑡 on ) = 𝜎 √ 𝜋 ∫ ∞ log 𝑡 on − 𝑢 · exp (cid:34) − ( 𝑢 − 𝜇 ) 𝜎 (cid:35) d 𝑢. (A5)This integral has analytical solution given by 𝑝 ( 𝑡 on ) = − 𝜇 − / − 𝜇 √ (cid:20) 𝜎 ln (cid:21) · Erfc (cid:20) log ( 𝑡 on / Myr ) − 𝜇 + 𝜎 ln10 𝜎 √ (cid:21) , (A6) where Erfc is the complementary error function. Finally, in orderto find the desired probability 𝑝 ( log ( 𝑡 on / Myr )) , we apply thetransformation of the variables and find that 𝑝 ( log ( 𝑡 on / Myr )) = ln10 𝑡 on · 𝑝 ( 𝑡 on ) . (A7)We illustrate this analytical solution for the 𝑝 ( log ( 𝑡 on / Myr )) given eq. (A6) and eq. (A7) in the toppanel of Fig. 9. This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000