Cosmic reionization after Planck II: contribution from quasars
MMNRAS , 000–000 (0000) Preprint 26 September 2017 Compiled using MNRAS L A TEX style file v3.0
Cosmic reionization after Planck II: contribution from quasars
Sourav Mitra (cid:63) , T. Roy Choudhury † and Andrea Ferrara ‡ University of the Western Cape, Bellville, Cape Town 7535, South Africa National Centre for Radio Astrophysics, TIFR, Post Bag 3, Ganeshkhind, Pune 411007, India Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa, Italy
26 September 2017
ABSTRACT
In the light of the recent
Planck downward revision of the electron scattering opticaldepth, and of the discovery of a faint AGN population at z > , we reassess the actual contri-bution of quasars to cosmic reionization. To this aim, we extend our previous MCMC-baseddata-constrained semi-analytic reionization model and study the role of quasars on globalreionization history. We find that, the quasars can alone reionize the Universe only for modelswith very high AGN emissivities at high redshift. These models are still allowed by the recentCMB data and most of the observations related to H I reionization. However, they predict anextended and early He II reionization ending at z (cid:38) and a much slower evolution in themean He II Ly- α forest opacity than what the actual observation suggests. Thus when we fur-ther constrain our model against the He II Ly- α forest data, this AGN-dominated scenario isfound to be clearly ruled out at 2- σ limits. The data seems to favour a standard two-component picture where quasar contributions become negligible at z (cid:38) and a non-zero escape frac-tion of ∼ is needed from early-epoch galaxies. For such models, mean neutral hydrogenfraction decreases to ∼ − at z = 6 . from ∼ . at z = 10 . and helium becomes doublyionized at much later time, z ∼ . We find that, these models are as well in good agree-ment with the observed thermal evolution of IGM as opposed to models with very high AGNemissivities. Key words: dark ages, reionization, first stars – intergalactic medium – quasars: general –cosmology: theory – large-scale structure of Universe.
Over the past few decades, various observational studies have beenperformed in order to constrain the epoch of hydrogen reionization,among them, observations from high-redshift QSO spectra (Beckeret al. 2001; White et al. 2003; Fan et al. 2006) and CMB by Wilkin-son Microwave Anisotropy Probe (WMAP) (Komatsu et al. 2011;Bennett et al. 2013) or Planck (Planck Collaboration et al. 2014,2016a,b) have been proven to be most useful. These observationsalong with theoretical models and simulations aim to help under-stand the nature of the reionization process and the sources respon-sible for it. Although it is commonly believed that, the IGM becameionized by the UV radiation from star forming high redshift galax-ies and the bright active galactic nuclei (AGN) populations, theirrelative contributions remain poorly understood (Giallongo et al.2012, 2015; Madau & Haardt 2015; Sharma et al. 2016; Finlatoret al. 2016). The rapidly decreasing number density of QSOs andAGNs at z > (Hopkins et al. 2007) usually leads to the assump- (cid:63) E-mail: [email protected] † E-mail: [email protected] ‡ E-mail: [email protected] tion that the high-redshift galaxies might have dominated hydrogenreionization at z (cid:38) (Barkana & Loeb 2001; Loeb & Barkana2001; Loeb 2006; Choudhury & Ferrara 2006a; Choudhury 2009).Based on this idea, several semi-analytic models have been put for-ward using a combination of different observations (Choudhury &Ferrara 2005; Wyithe & Loeb 2005; Gallerani et al. 2006; Dijkstraet al. 2007; Samui et al. 2007; Iliev et al. 2008; Mitra et al. 2011;Kulkarni & Choudhury 2011; Mitra et al. 2012, 2013; Cai et al.2014).However, recent measurements from Planck brought up animportant question regarding the contribution from high-redshiftsources. Unlike the WMAP results, they found relatively smallvalue for integrated reionization optical depth, τ el = 0 . ± . from Planck 2015; (Planck Collaboration et al. 2016a) or . ± . from Planck 2016; (Planck Collaboration et al.2016b), which corresponds to a sudden reionization occurring atmuch lower redshift, mean z reion ≈ . − . (Planck Collabora-tion et al. 2016c) and thus reduces the need for high- z galaxies asreionization sources (Robertson et al. 2015; Bouwens et al. 2015b;Mitra et al. 2015). This also seems to explain the rapid drop ob-served in the space density of Ly- α emitters at z ∼ (Mesingeret al. 2015; Choudhury et al. 2014; Schenker et al. 2014). Using c (cid:13) a r X i v : . [ a s t r o - ph . C O ] S e p Mitra, Choudhury & Ferrara the Planck constraints along with high- z galaxy luminosity func-tion (LF) data (Bouwens & Illingworth 2006; Oesch et al. 2012;Bradley et al. 2012; Oesch et al. 2014; Bowler et al. 2014; McLeodet al. 2014; Bouwens et al. 2015a), one can obtain relatively lowerand almost constant escape fraction ( f esc ) of ionizing photons fromgalaxies (but also see Price et al. 2016), which was somewhat in-evitably higher and evolving towards higher redshift in order tomatch the WMAP data (Mitra et al. 2013). In particular, recentlyMitra et al. (2015) found that a constant ∼ f esc at z (cid:62) issufficient to ionize the Universe and can as well explain most of thecurrent observational limits related to H I reionization. A similarconstraint on f esc has been reported by Bouwens et al. (2015b) andKhaire et al. (2016). Likewise, using an improved semi-numericalcode, Hassan et al. (2016) require f esc ≈ − , independentof halo mass and redshift, to simultaneously match current keyreionization observables. Ma et al. (2015) too found a non-evolvingand considerably lower ( < ) escape fraction using their high-resolution cosmological zoom-in simulation for galaxy formation.Such small leakages from star-forming galaxies again leads to thefact to reconsider their contributions for reionization.On the other hand, the conventional scenario of rapidly de-creasing ionizing QSO background beyond z (cid:38) is also chang-ing, owing to various recent claims on finding significantly higher(than previously expected) number of faint AGNs at higher red-shifts ( z > ) which are able to contribute a rather steep luminos-ity function (Glikman et al. 2011; Civano et al. 2011; Giallongoet al. 2012, 2015). Interestingly, Khaire et al. (2016) found that f esc needs to increase by a factor of ∼ from z ∼ . to 5.5in order to match the measured photoionization rates in case thesefainter AGNs are not taken into account. As a result, it has beenproposed recently that the reionization may have been driven en-tirely by the QSOs (Madau & Haardt 2015; Khaire et al. 2016).These models have no difficulty in matching the τ el constraintsfrom Planck. In addition, they lead to two significant features ofreionization which might resolve some puzzling observations ofIGM: a late hydrogen reionization completing around z ∼ , whichexplains the remarkable flatness observed in the H I photoioniza-tion rate (Becker & Bolton 2013) and a distinguishing feature ofearly He II reionization ending by z ∼ , which seems to be con-sistent with the measured high-transmission regions of He II Ly- α forest out to z ∼ . (Worseck et al. 2016). However, several recentclaims make this AGN-only reionization picture somewhat contro-versial. Kim et al. (2015) claimed that the number of spectroscopi-cally confirmed faint AGNs at z (cid:38) might not be high enough toreionize the Universe alone (see Finlator et al. 2016 as well). Also,D’Aloisio et al. (2017) pointed out a possible tension between theseAGN-dominated models and observed thermal history of the IGM.More recently, based on the SDSS high-redshift quasar sample at z ∼ , Jiang et al. (2016) find that the observed quasar populationis not enough to support the AGN-dominated reionization scenario.Thus, the possibility of a larger QSO emissivity at higher redshiftscombined with the lower τ el data from Planck demands more rigor-ous and careful investigation on the actual contribution from star-forming galaxies and QSOs to reionization.In this paper, we extend our previous Markov Chain Monte-Carlo (MCMC)-constrained reionization model with Planck (Mitraet al. 2015) to explore the possibility of quasar-only reionizationscenario by comparing different AGN emissivities. The paper isstructured as follows. In the next section, we review the main fea-tures of our semi-analytical reionization model. We consider fourdifferent cases to compare their relative quasar contributions. InSection 3, we present the MCMC constraints on reionization sce- nario obtained from each model. We compare our predictions withthe observation related to the thermal evolution of IGM in Sec-tion 4. Finally, we summarize and conclude our main findings inSection 5. We assume a flat Universe with Planck cosmologicalparameters: Ω m = 0 . , Ω Λ = 1 − Ω m , Ω b h = 0 . , h = 0 . , σ = 0 . , n s = 0 . and Y P = 0 . (Planck Collaboration et al. 2016a). In the following sections, we develop a robust statistical approachin order to constrain a simple reionization model against selecteddatasets. We perform a detailed MCMC analysis over our modelparameters related to the mean free path of ionizing photons andstellar and quasar emissivities. The method is built on a series of ourearlier papers (Mitra et al. 2011, 2012, 2013, 2015) with the addi-tion of different quasar contribution models. We study their effectson hydrogen and helium reionization separately and try to see howmuch information can be gained on reionization sources by addinghelium reionization data in our analysis. Later we also check theconsistency of such models with observed IGM thermal evolution.Throughout this work, we use a data-constrained semi-analytical reionization model, which is based on Choudhury & Fer-rara (2005) and Choudhury & Ferrara (2006b). Below we summa-rize its key features. • The model adopts a lognormal a probability distribution func-tion P (∆) ( ∆ being the overdensity of corresponding regions) atlow densities, changing to a power-law at high densities (Choud-hury & Ferrara 2005). The ionization and thermal histories of allhydrogen and helium regions are evaluated self-consistently bytreating the IGM as a multi-phase medium. It accounts for the IGMinhomogeneities in a very similar manner proposed by Miralda-Escud´e et al. (2000), where reionization is said to be complete onceall the low-density regions with overdensities ∆ < ∆ crit are ion-ized, where the critical density ∆ crit is determined from the meanseparation of the ionizing sources (Miralda-Escud´e et al. 2000;Wyithe & Loeb 2003). According to this prescription, the meanfree path of photons is computed from the distribution of high den-sity regions: λ mfp ( z ) = λ [1 − F V ( z )] / (1)where F V is the volume fraction of ionized regions, i.e., F V (∆ i ) = (cid:90) ∆ i d∆ P (∆) (2)and λ is a normalization parameter which we take as a free param-eter. We assume that, ∆ i does not evolve significantly with time inthe pre-overlap stage; it is equal to a critical value ∆ crit . In thiswork, it has been taken as a free parameter. Once ∆ crit is known,we follow the evolution of the ionized volume filling factor of thecorresponding regions, Q i , until it becomes unity and followingthat, we enter the post-overlap stage. • Stellar contribution : Reionization is assumed to be drivenby PopII stellar sources which have sub-solar metallicities and aSalpeter IMF in the mass range − M (cid:12) . We neglect the contri-butions from the metal-free PopIII sources as the impact of those onreionization is likely to be insignificant due to low electron scatter-ing optical depth data from Planck (Mitra et al. 2015). The modelalso calculates the suppression of star formation in low-mass haloesi.e. radiative feedback self-consistently from the thermal evolution MNRAS000
Over the past few decades, various observational studies have beenperformed in order to constrain the epoch of hydrogen reionization,among them, observations from high-redshift QSO spectra (Beckeret al. 2001; White et al. 2003; Fan et al. 2006) and CMB by Wilkin-son Microwave Anisotropy Probe (WMAP) (Komatsu et al. 2011;Bennett et al. 2013) or Planck (Planck Collaboration et al. 2014,2016a,b) have been proven to be most useful. These observationsalong with theoretical models and simulations aim to help under-stand the nature of the reionization process and the sources respon-sible for it. Although it is commonly believed that, the IGM becameionized by the UV radiation from star forming high redshift galax-ies and the bright active galactic nuclei (AGN) populations, theirrelative contributions remain poorly understood (Giallongo et al.2012, 2015; Madau & Haardt 2015; Sharma et al. 2016; Finlatoret al. 2016). The rapidly decreasing number density of QSOs andAGNs at z > (Hopkins et al. 2007) usually leads to the assump- (cid:63) E-mail: [email protected] † E-mail: [email protected] ‡ E-mail: [email protected] tion that the high-redshift galaxies might have dominated hydrogenreionization at z (cid:38) (Barkana & Loeb 2001; Loeb & Barkana2001; Loeb 2006; Choudhury & Ferrara 2006a; Choudhury 2009).Based on this idea, several semi-analytic models have been put for-ward using a combination of different observations (Choudhury &Ferrara 2005; Wyithe & Loeb 2005; Gallerani et al. 2006; Dijkstraet al. 2007; Samui et al. 2007; Iliev et al. 2008; Mitra et al. 2011;Kulkarni & Choudhury 2011; Mitra et al. 2012, 2013; Cai et al.2014).However, recent measurements from Planck brought up animportant question regarding the contribution from high-redshiftsources. Unlike the WMAP results, they found relatively smallvalue for integrated reionization optical depth, τ el = 0 . ± . from Planck 2015; (Planck Collaboration et al. 2016a) or . ± . from Planck 2016; (Planck Collaboration et al.2016b), which corresponds to a sudden reionization occurring atmuch lower redshift, mean z reion ≈ . − . (Planck Collabora-tion et al. 2016c) and thus reduces the need for high- z galaxies asreionization sources (Robertson et al. 2015; Bouwens et al. 2015b;Mitra et al. 2015). This also seems to explain the rapid drop ob-served in the space density of Ly- α emitters at z ∼ (Mesingeret al. 2015; Choudhury et al. 2014; Schenker et al. 2014). Using c (cid:13) a r X i v : . [ a s t r o - ph . C O ] S e p Mitra, Choudhury & Ferrara the Planck constraints along with high- z galaxy luminosity func-tion (LF) data (Bouwens & Illingworth 2006; Oesch et al. 2012;Bradley et al. 2012; Oesch et al. 2014; Bowler et al. 2014; McLeodet al. 2014; Bouwens et al. 2015a), one can obtain relatively lowerand almost constant escape fraction ( f esc ) of ionizing photons fromgalaxies (but also see Price et al. 2016), which was somewhat in-evitably higher and evolving towards higher redshift in order tomatch the WMAP data (Mitra et al. 2013). In particular, recentlyMitra et al. (2015) found that a constant ∼ f esc at z (cid:62) issufficient to ionize the Universe and can as well explain most of thecurrent observational limits related to H I reionization. A similarconstraint on f esc has been reported by Bouwens et al. (2015b) andKhaire et al. (2016). Likewise, using an improved semi-numericalcode, Hassan et al. (2016) require f esc ≈ − , independentof halo mass and redshift, to simultaneously match current keyreionization observables. Ma et al. (2015) too found a non-evolvingand considerably lower ( < ) escape fraction using their high-resolution cosmological zoom-in simulation for galaxy formation.Such small leakages from star-forming galaxies again leads to thefact to reconsider their contributions for reionization.On the other hand, the conventional scenario of rapidly de-creasing ionizing QSO background beyond z (cid:38) is also chang-ing, owing to various recent claims on finding significantly higher(than previously expected) number of faint AGNs at higher red-shifts ( z > ) which are able to contribute a rather steep luminos-ity function (Glikman et al. 2011; Civano et al. 2011; Giallongoet al. 2012, 2015). Interestingly, Khaire et al. (2016) found that f esc needs to increase by a factor of ∼ from z ∼ . to 5.5in order to match the measured photoionization rates in case thesefainter AGNs are not taken into account. As a result, it has beenproposed recently that the reionization may have been driven en-tirely by the QSOs (Madau & Haardt 2015; Khaire et al. 2016).These models have no difficulty in matching the τ el constraintsfrom Planck. In addition, they lead to two significant features ofreionization which might resolve some puzzling observations ofIGM: a late hydrogen reionization completing around z ∼ , whichexplains the remarkable flatness observed in the H I photoioniza-tion rate (Becker & Bolton 2013) and a distinguishing feature ofearly He II reionization ending by z ∼ , which seems to be con-sistent with the measured high-transmission regions of He II Ly- α forest out to z ∼ . (Worseck et al. 2016). However, several recentclaims make this AGN-only reionization picture somewhat contro-versial. Kim et al. (2015) claimed that the number of spectroscopi-cally confirmed faint AGNs at z (cid:38) might not be high enough toreionize the Universe alone (see Finlator et al. 2016 as well). Also,D’Aloisio et al. (2017) pointed out a possible tension between theseAGN-dominated models and observed thermal history of the IGM.More recently, based on the SDSS high-redshift quasar sample at z ∼ , Jiang et al. (2016) find that the observed quasar populationis not enough to support the AGN-dominated reionization scenario.Thus, the possibility of a larger QSO emissivity at higher redshiftscombined with the lower τ el data from Planck demands more rigor-ous and careful investigation on the actual contribution from star-forming galaxies and QSOs to reionization.In this paper, we extend our previous Markov Chain Monte-Carlo (MCMC)-constrained reionization model with Planck (Mitraet al. 2015) to explore the possibility of quasar-only reionizationscenario by comparing different AGN emissivities. The paper isstructured as follows. In the next section, we review the main fea-tures of our semi-analytical reionization model. We consider fourdifferent cases to compare their relative quasar contributions. InSection 3, we present the MCMC constraints on reionization sce- nario obtained from each model. We compare our predictions withthe observation related to the thermal evolution of IGM in Sec-tion 4. Finally, we summarize and conclude our main findings inSection 5. We assume a flat Universe with Planck cosmologicalparameters: Ω m = 0 . , Ω Λ = 1 − Ω m , Ω b h = 0 . , h = 0 . , σ = 0 . , n s = 0 . and Y P = 0 . (Planck Collaboration et al. 2016a). In the following sections, we develop a robust statistical approachin order to constrain a simple reionization model against selecteddatasets. We perform a detailed MCMC analysis over our modelparameters related to the mean free path of ionizing photons andstellar and quasar emissivities. The method is built on a series of ourearlier papers (Mitra et al. 2011, 2012, 2013, 2015) with the addi-tion of different quasar contribution models. We study their effectson hydrogen and helium reionization separately and try to see howmuch information can be gained on reionization sources by addinghelium reionization data in our analysis. Later we also check theconsistency of such models with observed IGM thermal evolution.Throughout this work, we use a data-constrained semi-analytical reionization model, which is based on Choudhury & Fer-rara (2005) and Choudhury & Ferrara (2006b). Below we summa-rize its key features. • The model adopts a lognormal a probability distribution func-tion P (∆) ( ∆ being the overdensity of corresponding regions) atlow densities, changing to a power-law at high densities (Choud-hury & Ferrara 2005). The ionization and thermal histories of allhydrogen and helium regions are evaluated self-consistently bytreating the IGM as a multi-phase medium. It accounts for the IGMinhomogeneities in a very similar manner proposed by Miralda-Escud´e et al. (2000), where reionization is said to be complete onceall the low-density regions with overdensities ∆ < ∆ crit are ion-ized, where the critical density ∆ crit is determined from the meanseparation of the ionizing sources (Miralda-Escud´e et al. 2000;Wyithe & Loeb 2003). According to this prescription, the meanfree path of photons is computed from the distribution of high den-sity regions: λ mfp ( z ) = λ [1 − F V ( z )] / (1)where F V is the volume fraction of ionized regions, i.e., F V (∆ i ) = (cid:90) ∆ i d∆ P (∆) (2)and λ is a normalization parameter which we take as a free param-eter. We assume that, ∆ i does not evolve significantly with time inthe pre-overlap stage; it is equal to a critical value ∆ crit . In thiswork, it has been taken as a free parameter. Once ∆ crit is known,we follow the evolution of the ionized volume filling factor of thecorresponding regions, Q i , until it becomes unity and followingthat, we enter the post-overlap stage. • Stellar contribution : Reionization is assumed to be drivenby PopII stellar sources which have sub-solar metallicities and aSalpeter IMF in the mass range − M (cid:12) . We neglect the contri-butions from the metal-free PopIII sources as the impact of those onreionization is likely to be insignificant due to low electron scatter-ing optical depth data from Planck (Mitra et al. 2015). The modelalso calculates the suppression of star formation in low-mass haloesi.e. radiative feedback self-consistently from the thermal evolution MNRAS000 , 000–000 (0000) osmic reionization after Planck II of IGM (Choudhury & Ferrara 2005; also see Okamoto et al. 2008;Sobacchi & Mesinger 2013). The production rate of ionizing pho-tons in the IGM is then computed as ˙ n ph ( z ) = n b N ion d f coll d t ≡ n b (cid:15) II N γ d f coll d t (3)where, f coll is the collapsed fraction of dark matter haloes, n b is thetotal baryonic number density in the IGM, N ion is the number ofphotons entering the IGM per baryon in stars and (cid:15) II is the productof the star formation efficiency (cid:15) ∗ and escape fraction f esc of theionizing photons from stars. We treat (cid:15) II as a free parameter in thisanalysis. The model also incorporates the contribution of quasars by com-puting their ionizing emissivities. In our earlier works (Mitra et al.2011, 2012, 2015), we adopted a fixed model for the AGN contri-bution and calculated the number of ionizing photons from QSOsby integrating their observed luminosity functions at z < (Hop-kins et al. 2007; hereafter, H07), where we assumed that they havenegligible effects on IGM at higher redshifts (Choudhury & Fer-rara 2005). The comoving QSO emissivity (cid:15) ν at ( (cid:15) H07912 inunits of erg s − Hz − Mpc − ) is then obtained by integrating theQSO luminosity function at each redshift assuming its efficiencyto be unity (but also see Faucher-Gigu`ere et al. 2009), which isidentical to integrating the B-band luminosity functions with theconversion L ν (912˚A) = 10 . erg s − Hz − ( L B /L (cid:12) ) (Schir-ber & Bullock 2003; Choudhury & Ferrara 2005; Richards et al.2006; Hopkins et al. 2007). However, this standard picture of neg-ligible QSO contributions at z > has been challenged by the re-cent findings of numerous faint AGN candidates at higher redshifts(Giallongo et al. 2015). Based on this, Madau & Haardt (2015)proposed an AGN comoving emissivity fit which is significantlyhigher at high- z . The faint AGN population alone can dominate thereionization process in their model. We shall refer this as MH15model.However, it has been argued that the observations of faint endQSO LFs at high- z ( z > ) are very uncertain (Yoshiura et al.2016). So, rather than taking a fixed AGN model, in this workwe aim to get the constraints on high-redshift LF slope backwardby varying it as an additional free parameter. This will help us togain more insight on the relative contributions of the stellar andAGN components allowed by the current data. For that, we take amodel with comoving AGN emissivity same as that provided byH07 model for z (cid:54) , while its redshift evolution at z > is ac-companied by an exponential fall depending on the free parameter β , so that this β will now determine the high- z slope: (cid:15) ( z ) = (cid:15) H07912 ( z ) , for z (cid:54) (cid:15) H07912 ( z = 2) × e β ( z − , for z > (4)This method is quite similar to what described in the recent studyby Yoshiura et al. (2016), where they try to constrain the contri-bution of high- z galaxies and AGNs to reionization by varying f esc and β (although the definition of β is different in their case)while simultaneously satisfying the Planck electron scattering op-tical depth data and the redshift for completion of hydrogen orhelium reionization. However, note that our analysis, in principle,should give a more robust constraint on these parameters since it isalso constrained by many other observables. We assume the quasarshave a frequency spectrum (cid:15) ν ∝ ν − . for frequencies above thehydrogen ionization edge. The above parametrization of the quasar emissivity can fairlyreproduce both the H07 ( β = − . ) and MH15 ( β = − . )cases (see Fig. 1). Note that, at lower redshifts ( z (cid:46) ) the co-moving AGN emissivity for our β = − . model is lower thanthe original MH15 fit. In that sense, we somewhat underestimatethe AGN contributions in compare to their model at those redshifts.But we shall see that our overall conclusions will remain unaffectedby this reconciliation. Our main goal is to investigate whether afaint AGN population alone can dominate the reionization process,or equivalently find out what are the acceptable ranges of β allowedby current observations related to reionization. We hence treat β asa free parameter along with our other model parameters. Although, the above model can be constrained by a variety of ob-servational datasets, we check that most of the relevant constraintscome from a subset of those data (Mitra et al. 2011). So, to keepthis analysis simple we obtain the constraints on reionization using mainly three datasets:(i) observed HI Ly α effective optical depth τ eff , HI at . (cid:54) z (cid:54) , as obtained from the QSO absorption spectra. The values adoptedfor τ eff , HI in this paper are based on the data tabulated in Fan et al.(2006) and Becker et al. (2013). We bin the data into redshift binsof width ∆ z = 0 . and compute the mean in each bin. The cor-responding uncertainties are estimated using the extreme values of τ eff , HI along different lines of sight ;(ii) redshift evolution of Lyman-limit systems (LLS), d N LL / d z over a wide redshift range . < z < from the combined dataof Songaila & Cowie (2010) and Prochaska et al. (2010) with theerrors calculated using the quadrature method;(iii) Thomson scattering optical depth τ el data ( . ± . )from Planck 2016 (Planck Collaboration et al. 2016b).We also impose somewhat model-independent upper limits on theneutral hydrogen fraction x HI at z ∼ − from McGreer et al.(2015) as a prior to our model.As mentioned earlier, the free parameters of this model are (cid:15) II , λ , ∆ crit and β ; all the cosmological parameters are fixed at theirbest-fit Planck value. In principle, (cid:15) II can have a dependence onredshift z and halo mass M , but for simplicity, we assume (cid:15) II to beindependent of z or M throughout this work. This is also motivatedby the results from our earlier work (Mitra et al. 2015), where wefind that both f esc and (cid:15) ∗ are almost non-evolving with redshift.Unlike our previous works (Mitra et al. 2011, 2012, 2013), we al-low ∆ crit to be a free parameter. In our models ∆ crit sets the meanfree path of ionizing photons and its value depends on the typi-cal separation between the ionizing sources (Miralda-Escud´e et al.2000). Since QSOs are relatively rarer sources, the implied ∆ crit is expected to be relatively higher for QSO-dominated models thanfor the stellar-dominated ones. For our simplified mean free pathprescription, this value usually turns out to be ∼ − depend-ing on the density profile of the halo (Choudhury & Ferrara 2005; From here on, H07 and MH15 refer to the models with β = − . and β = − . respectively. In our earlier works (Mitra et al. 2011, 2012, 2013) we have used thehydrogen photoionization rate measurements (Bolton & Haehnelt 2007;Becker & Bolton 2013) instead of the τ eff , HI used in this work. We havechecked and found that the main conclusions of the paper remain unchangedirrespective of which of the two data sets is used in the analysis.MNRAS , 000–000 (0000) Mitra, Choudhury & Ferrara
Choudhury 2009). Keeping this in mind, we allow it to vary only inthe range [20 , . We also impose a prior β (cid:54) to ensure that thequasar emissivities do not diverge at high redshifts. As we will see in the following section(s), the constraints on theQSO emissivity models are relatively weaker when we consideronly those datasets related to the hydrogen reionization. The AGNscan produce significantly more hard ionizing photons compared tostellar sources and should ionize He II more efficiently. Thus thehelium reionization history should have a more decisive depen-dence on β than the H I reionization (Yoshiura et al. 2016). Ob-servations of the He II Ly α forest have already been used to studythe He II reionization (Gleser et al. 2005; Dixon & Furlanetto 2009;McQuinn et al. 2009; Khaire & Srianand 2013).Given this, we include one more observable related to heliumreionization, namely the effective optical depth of He II Ly α ab-sorption τ eff , HeII data at . < z < . from Worseck et al.(2016). This observation manifests surprisingly low He II absorp-tion at z ∼ . , which might indicate that the bulk of the heliumwas ionized much earlier (at z > . ). This is somewhat in tensionwith the general findings from current numerical radiative transfersimulations (McQuinn et al. 2009; Compostella et al. 2013, 2014),where the He II reionization ends at relatively later epoch, z ≈ (McQuinn 2016; D’Aloisio et al. 2017). Thus, it would be inter-esting to see what sources could support such an early reionizationwithin our semi-analytic formalism.The observed values of τ eff , HeII along different sightlinesshow a large scatter particularly at z (cid:38) (Worseck et al. 2016).In order to adapt them into our likelihood analysis, we have binnedthe data points within redshift intervals of z = 0 . and calculatedthe mean. The errors are calculated using the extreme values of τ eff , HeII along different lines of sight.The inclusion of the τ eff , HeII would imply additional calcula-tions in our theoretical model. We briefly summarize the additionalsteps which have been included for this purpose, the details can befound in Choudhury & Ferrara (2005). • As in the case of hydrogen reionization, we assume the He II in the low-density ∆ < ∆ crit , HeII to be ionized first, where ∆ crit , HeII is the critical overdensity similar to the ∆ crit used forhydrogen reionization. It is essentially determined by the mean sep-aration between helium ionizing sources and can, in principle, bedifferent from ∆ crit . However, we find that our results are relativelyinsensitive to the exact chosen value of ∆ crit , HeII , hence we avoidintroduction of one more free parameter in our model and simplyuse ∆ crit , HeII = ∆ crit . • We use relations analogous to equations (1) and (2) to esti-mate the mean free path λ mfp , HeII for He II ionizing photons. Giventhe emissivity (cid:15) ν and λ mfp , HeII , it is straightforward to calculatethe photoionization rate for He II (assuming λ mfp , HeII to be muchsmaller than the horizon size, which holds true for z (cid:38) ). • We evolve the average fraction of different ionization statesof helium, along with that of hydrogen and the temperature, whichcan then be used for calculating the Ly α optical depth τ HeII forhelium. The effective optical depth is simply obtained from the av-erage value of the corresponding transmitted flux and is given by τ eff , HeII = − ln (cid:10) e − τ HeII (cid:11) . (5) z † ( e r g s − H z − M p c − ) β = − . β = − . without HeII datawith HeII dataH07 β = − . MH15 β = − . Figure 1.
MCMC constraints on ionizing comoving AGN emissivity fortwo different cases: (i) without taking He II data - white dashed lines forbest-fit result with orange shaded region around it for 2- σ C.L., (ii) withHe II data - solid blue curve for the best-fit with shaded cyan region for2- σ limits. The H07 (or β = − . ) and MH15 (or β = − . ) mod-els are also shown here by dotted and dot-dashed curves respectively.Theblack points with errorbars are the observed data inferred from QLFs byGiallongo et al. (2015). We perform an MCMC analysis over all the parameter space { (cid:15) II , λ , ∆ crit , β } using the above mentioned datasets. We employa code based on the publicly available COSMOMC (Lewis & Bri-dle 2002) code and run a number of separate chains until the usualGelman and Rubin convergence criterion is satisfied. This methodbased on the MCMC analysis has already been developed in ourprevious works (Mitra et al. 2011, 2012, 2015).While presenting the results, we clearly distinguish betweentwo cases: (i) “without He II data” where we obtain constraints us-ing only hydrogen reionization data (Section 2.2), and (ii) “withHe II data” where we include the τ eff , HeII data as well (Section2.3). This is done to emphasize the significance of the He II reion-ization data in constraining the QSO emissivity models.The results from our current MCMC analysis for these twomodels are summarized in Table 1, where the first four rows are forthe free parameters of the model and the last one ( τ el ) correspondsto the derived parameter. For comparison, we also show here H07( β = − . ) and MH15 ( β = − . ) cases. For that, we fixthe ∆ crit = 60 (following our earlier works) and choose somecombination of (cid:15) II and λ (without running the MCMC) so thatthey can reasonably match all the observed datasets for HI Ly α effective optical depth, redshift evolution of LLSs and also produce τ el which is allowed by Planck 2016. We have plotted the MCMCconstraints on the comoving AGN emissivities in Fig. 1 and theother quantities of interest in Fig. 2. The constraints on the fractionof hydrogen photoionization rates contributed by QSOs (comparedto its total, i.e., QSOs + galaxies, value) are shown in Fig. 3. Forcomparison, we also show the H07 and MH15 models in these threefigures. MNRAS000
MCMC constraints on ionizing comoving AGN emissivity fortwo different cases: (i) without taking He II data - white dashed lines forbest-fit result with orange shaded region around it for 2- σ C.L., (ii) withHe II data - solid blue curve for the best-fit with shaded cyan region for2- σ limits. The H07 (or β = − . ) and MH15 (or β = − . ) mod-els are also shown here by dotted and dot-dashed curves respectively.Theblack points with errorbars are the observed data inferred from QLFs byGiallongo et al. (2015). We perform an MCMC analysis over all the parameter space { (cid:15) II , λ , ∆ crit , β } using the above mentioned datasets. We employa code based on the publicly available COSMOMC (Lewis & Bri-dle 2002) code and run a number of separate chains until the usualGelman and Rubin convergence criterion is satisfied. This methodbased on the MCMC analysis has already been developed in ourprevious works (Mitra et al. 2011, 2012, 2015).While presenting the results, we clearly distinguish betweentwo cases: (i) “without He II data” where we obtain constraints us-ing only hydrogen reionization data (Section 2.2), and (ii) “withHe II data” where we include the τ eff , HeII data as well (Section2.3). This is done to emphasize the significance of the He II reion-ization data in constraining the QSO emissivity models.The results from our current MCMC analysis for these twomodels are summarized in Table 1, where the first four rows are forthe free parameters of the model and the last one ( τ el ) correspondsto the derived parameter. For comparison, we also show here H07( β = − . ) and MH15 ( β = − . ) cases. For that, we fixthe ∆ crit = 60 (following our earlier works) and choose somecombination of (cid:15) II and λ (without running the MCMC) so thatthey can reasonably match all the observed datasets for HI Ly α effective optical depth, redshift evolution of LLSs and also produce τ el which is allowed by Planck 2016. We have plotted the MCMCconstraints on the comoving AGN emissivities in Fig. 1 and theother quantities of interest in Fig. 2. The constraints on the fractionof hydrogen photoionization rates contributed by QSOs (comparedto its total, i.e., QSOs + galaxies, value) are shown in Fig. 3. Forcomparison, we also show the H07 and MH15 models in these threefigures. MNRAS000 , 000–000 (0000) osmic reionization after Planck II Parameters best-fit with 95% C.L. fixed β -model (No MCMC)without He II data with He II data H07 ( β = − . ) MH15 ( β = − . ) β − .
06 [ − . , . − .
43 [ − . , − . − . − . (cid:15) II × .
51 [0 . , .
94] 3 .
52 [2 . , .
56] 3 .
30 1 . λ .
13 [1 . , .
71] 3 .
35 [1 . , .
18] 3 .
85 4 . crit .
14 [25 . , .
97] 52 .
63 [25 . , .
96] 60 . . τ el .
058 [0 . , . .
062 [0 . , . .
062 0 . Table 1.
Best-fit value and 95% C.L. errors of the model parameters (above four) and derived parameter (bottom row) obtained from the current MCMCanalysis for two cases: with and without He II data. We also show the H07 and MH15 models with β = − . and − . respectively, where we fix ∆ crit = 60 and choose some (cid:15) II and λ so that they can fairly match all the datasets considered here for hydrogen reionization. z τ e l without HeII datawith HeII dataH07 ( β = − . )MH15 ( β = − . ) (a) Planck15 d N LL / d z (b) 5 10024681012 τ e ff , H I (c) z τ e ff , H e II (d)5 10 z Q H II (e) 5 10 z -4 -3 -2 x H I z Q H e III (g) 2 3 4 5 6 z -4 -3 -2 x H e II Figure 2.
The MCMC constraints on various quantities related to reionization. Different panels indicate: (a) electron scattering optical depth τ el , (b) redshiftevolution of Lyman-limit systems d N LL / d z , (c) effective H I Ly α optical depth τ eff , HI , (d) effective He II Ly α optical depth τ eff , HeII , (e) volume fillingfactor of ionized hydrogen Q HII , (f) neutral hydrogen fraction x HI ( z ) , (g) He III ionized volume filling factor Q HeIII and (h) global He II fraction x HeII . Thesolid blue and dashed white lines correspond to the mean evolution obtained from the models with He II data and without He II data respectively. The shadedregions refer to the 2- σ confidence limits around the mean value of the corresponding models. The black points with errorbars are the current observationallimits on reionization, see the main text for references. For comparison, we also plot the H07 (i.e. β = − . ) and MH15 ( β = − . ) models. II data” The 2- σ or 95% confidence limits on our model parameters withoutHe II data are given in the second column of Table 1 and shownby the shaded orange regions in Figs 1–3. The best-fit model isshown by the dashed white curves in these figures. We can see fromthe table that the data allows a wide range of β values spanningfrom ∼ − (signifying that the quasar contributions are almostnegligible at z > ) to ∼ (signifying that the quasars dominateeven at higher redshifts and a very little contribution is coming fromgalaxies). This is also evident from Fig. 1 where we see that theshaded orange region is remarkably wide. Both the H07 (dottedblack curve) and MH15 (dot-dashed black curve) cases are wellinside this allowed region at z (cid:38) .From the MCMC constraints on various quantities related toreionization shown in Fig. 2, we find that the allowed 2- σ confi-dence limits are relatively narrower for low- z regime and increaseat z (cid:38) . This is expected as most of the datasets used in this workexist only at low redshifts z (cid:54) , whereas the higher- z epoch is poorly constrained (Mitra et al. 2011, 2015). The reionization op-tical depth for the best-fit model is in a good agreement with thePlanck 2016 data (black point with errorbar in panel a). We alsoshow the earlier Planck 2015 limit by dotted red lines. In panel b,we plot the evolution of Lyman-limit systems which again matcheshe combined observational data points from Songaila & Cowie(2010) and Prochaska et al. (2010) at < z < quite reason-ably for the constraints obtained. We plot the evolution of effectiveoptical depth of H I Ly α absorption τ eff , HI in panel c. The MCMCconstraint obtained is in general agreement with the observationalmeasurements from Fan et al. (2006); Becker et al. (2013) in theinterval < z < . From the evolution of volume filling factor Q HII of H II regions (panel e) and the global neutral hydrogen frac-tion x HI (panel f), we find that the hydrogen reionization history isreasonably well constrained, in spite of the quasar emissivity beingallowed to take such wide range of values. This is because any vari-ation in β is appropriately compensated by a similar change in (cid:15) II making sure that the models agree with the observations. In fact, MNRAS , 000–000 (0000)
Mitra, Choudhury & Ferrara z Γ Q S O P I / Γ T o t a l P I without HeII datawith HeII dataH07 ( β = − . )MH15 ( β = − . ) Figure 3.
Same as Fig. 2, but showing the constraints on the ratio of pho-toionization rates coming from quasars to its total (galaxies + quasars) valueat different redshifts. the allowed range of (cid:15) II , as can be seen from Table 1, is also quitewide between ∼ × − − ∼ × − . We also find that themodel can match the current observed constraints on x HI quite rea-sonably within their errorbars. Note that the match is quite impres-sive, given the fact that we did not include these datasets, exceptthe McGreer et al. (2015) data at z ∼ − , as constraints in thecurrent analysis. The observational limits (black points) are takenfrom various measurements by Fan et al. (2006) (filled circle), Mc-Greer et al. (2015) (open triangle), Totani et al. (2006); Chornocket al. (2013) (filled triangle), Bolton et al. (2011); Schroeder et al.(2013) (filled diamond), Ota et al. (2008); Ouchi et al. (2010) (opensquare), Schenker et al. (2014) (filled square).Moving on to quantities related to He II reionization, we plotthe effective optical depth τ eff , HeII in panel d. The points with er-ror bars are from Worseck et al. (2016), binned appropriately for theMCMC analysis. Note that, these data points are not taken into ac-count for the results presented in this subsection. Clearly, ignoringthese data points leads to the 2- σ allowed regions that are consider-ably wide which is also related to the fact that β is allowed to take awide range of values. We show the evolution of He III ionized vol-ume fraction Q HeIII and the global He II fraction x HeII in panelsg and h, respectively. The case Q HeIII = 1 implies that helium isdoubly reionized, and the x HeII in that case is determined by theresidual He II fraction in the He III regions. The global He II reion-ization history mainly depends on the contribution from QSOs andthus leaves a distinguishing feature for different β models - higherthe value of β , earlier the He II reionization occurs. The best-fitmodel without He II data or the AGN-dominated MH15 model in-dicates that, the helium becomes doubly reionized at z (cid:38) , how-ever, the 2- σ limits allow a wide range of reionization redshifts. Forexample, He II reionization is allowed to be completed as early as z ∼ (implying a high β ∼ ) and as late as z ∼ . (implying alow β ∼ − ).From Fig. 3, we find that the fraction of the photoionizationrate contributed by the quasars can take any value between and ∼ at z > . Interestingly, the data allow the fraction to haveconsistently high value ∼ over the full redshift range 50% ( ∼ of the IGM at z (cid:39) . z (cid:39) . . This shouldbe compared with our predictions of the volume filling factor ofHe III regions for the best-fit model which gives Q HeIII ∼ at z (cid:39) . , remarkably consistent with that given by Worseck et al.(2016). The corresponding value at z (cid:39) . is Q HeIII ∼ forour model which is slightly higher than the Worseck et al. (2016)value, however, the two results are fully consistent if we accountfor the statistical uncertainties in our analysis. The H07 model alsoyields a very similar evolution as the best-fit blue curve. Modelswith β < − β > − . imply a very late (early) reionizationhistory and fail to match observational data. However, an improvedsurvey with much larger sample of z > He II sightlines will beneeded to put a tighter constraint on β in future.The fraction of photoionization rate contributed by thequasars, as shown in Fig 3, is significantly small. It can take valuesat most ∼ ( ∼ ) at z = 6(4) , thus implying that quasardominated reionization models at z > are strongly disfavoured.For the best-fit model with He II data and H07, the AGN contribu-tions become negligible at z (cid:38) , which is in stark contrast againstthe MH15 case where the AGNs dominate ( (cid:38) of total value)even at higher redshifts. MNRAS000 50% ( ∼ of the IGM at z (cid:39) . z (cid:39) . . This shouldbe compared with our predictions of the volume filling factor ofHe III regions for the best-fit model which gives Q HeIII ∼ at z (cid:39) . , remarkably consistent with that given by Worseck et al.(2016). The corresponding value at z (cid:39) . is Q HeIII ∼ forour model which is slightly higher than the Worseck et al. (2016)value, however, the two results are fully consistent if we accountfor the statistical uncertainties in our analysis. The H07 model alsoyields a very similar evolution as the best-fit blue curve. Modelswith β < − β > − . imply a very late (early) reionizationhistory and fail to match observational data. However, an improvedsurvey with much larger sample of z > He II sightlines will beneeded to put a tighter constraint on β in future.The fraction of photoionization rate contributed by thequasars, as shown in Fig 3, is significantly small. It can take valuesat most ∼ ( ∼ ) at z = 6(4) , thus implying that quasardominated reionization models at z > are strongly disfavoured.For the best-fit model with He II data and H07, the AGN contribu-tions become negligible at z (cid:38) , which is in stark contrast againstthe MH15 case where the AGNs dominate ( (cid:38) of total value)even at higher redshifts. MNRAS000 , 000–000 (0000) osmic reionization after Planck II As far as the other quantities are concerned, we find that thebest-fit value of ∆ crit is slightly larger for the model without He II data, as it allows higher contribution from the QSOs than the modelwith He II data. However, both the models support a much broaderrange for this parameter; ∆ crit ∼ − within the 2- σ limits.Thus, our results do not vary considerably for the choice of ∆ crit aslong as it is within this limit. We also find that, the electron scatter-ing optical depths τ el for all the best-fit models are in good agree-ment with the Planck 2016 value. Interestingly, when we includeHe II data, the best-fit model predicts slightly higher τ el than theother models, as it allows the highest stellar contributions amongthem.Interestingly, we find that for the best-fit He II model which isconsistent with all the datasets considered in this paper, the stellarcomponent of the comoving ionizing emissivity (in units of s − Mpc − ) can be fitted quite well by fitting functions of the form log ˙ n stellarion ( z ) = 59 . 786 e − . z − . 844 e − . z (6)which, when combined with the QSO emissivities given in Sec-tion 2.1, should be useful in computing the reionization history. A different way of constraining the quasar emissivity and the asso-ciated He II reionization is by studying the thermal history of theIGM. The AGN-dominated models can lead to an additional heat-ing of the IGM as they produce significantly more hard ionizingphotons with energies in excess of 4 Ry compared to stellar sources.Thus the observational measurements on the IGM temperature canalso be used to discriminate between different source models andput additional constraints on reionization (Zaldarriaga et al. 2001;Trac et al. 2008; Furlanetto & Oh 2009; Cen et al. 2009; Raskuttiet al. 2012; Padmanabhan et al. 2014).Note that a careful calculation of the thermal evolution of theIGM is computationally much more expensive than that for thesemi-analytic models described in the earlier sections and henceis not suitable for including in the MCMC analysis. Instead, wecalculate the thermal histories for our best-fit models (along withthe H07 and MH15 cases) and compare our predictions with theobservations from Becker et al. (2011).Our model computes the IGM temperature ( T ) evolution self-consistently and separately for each of the three regions, i.e., (i)completely neutral, (ii) regions where hydrogen is ionized and he-lium is singly ionized and (iii) where both species are fully ion-ized. For neutral regions, the temperature is usually assumed to de-crease adiabatically while in the ionized regions, it is calculated as(Choudhury & Ferrara 2005; Upton Sanderbeck et al. 2016): d T d t = − H ( t ) T + 23 T ∆ d∆d t − T (cid:80) i X i dd t (cid:88) i X i + 23 k B n b (1 + z ) d E d t (7)where we have defined X i ≡ n i m p ρ b (∆) (8)for three independent species X HI , X HeI , X HeIII and ρ b is theproper mass density of baryons. The first term on the right-handside accounts for the adiabatic cooling of gas due to cosmic expan-sion. The second term calculates the adiabatic evolution due to col-lapsing overdensities. The third term describes the rate of changein the number of particles of the system and in the last term on the right-hand side, d E/ d t represents the net heating rate per baryonwhich encodes all other possible heating and cooling processes.For most purposes, it is sufficient to take into account only the pho-toheating, recombination cooling and Compton cooling off CMBphotons (Choudhury & Ferrara 2005).In order to calculate the temperature at a given density ∆ andredshift z , we track the thermal evolution of a large number of den-sity elements ( ∼ , ), generated according to the IGM proba-bility distribution. It is this step which increases the computationalcost of the analysis. During hydrogen reionization, our model as-sumes that all regions with ∆ > ∆ crit remain neutral, while afraction Q HII of ∆ < ∆ crit regions are ionized (Miralda-Escud´eet al. 2000). In the post-reionization Q HII = 1 era, the effect ofionizing radiation is to increase the value of the density thresholdabove which the IGM is neutral. At every time step, we calculatethe increase in the ionized fraction and randomly assign the corre-sponding number of density elements as newly ionized. FollowingUpton Sanderbeck et al. (2016) and D’Aloisio et al. (2017) (andthe references therein), each of these density elements is instanta-neously heated, on top of the uniform photoheating background, toa temperature of , K. We follow an identical procedure forHe II reionization as well and assume the density elements to beadditionally heated by , K when they reach their He II reion-ization redshifts.This method allows us to calculate, at any given redshift, thedependence of the gas temperature T (∆) on the overdensity ∆ . Wecan fit the T (∆) relation using a power-law of the form T (∆) = T ∆ γ − , (9)where T is the temperature of the mean density ( ∆ = 1 ) gas and γ is the slope of the temperature-density relation. We can estimatethe values of T and γ at every redshift.In the left panel of Fig. 4, we show the redshift evolution ofIGM temperature T (∆) at different overdensities ∆ at redshiftswhere the data points from Becker et al. (2011) exist. Differentcurves are for different quasar emissivity models, while the pointswith error bars are the observed data (Becker et al. 2011). The firstpoint to note is that all the models, irrespective of their quasar emis-sivity, underpredict the temperature at z < . This could be becauseof assumptions in our semi-analytic model (e.g., the effective spec-tral index of the radiation from the quasars could be harder thanwhat we assumed, which would lead to additional heating), or be-cause of additional heating sources like the blazars (Puchwein et al.2012), cosmic rays or intergalactic dust absorption (Upton Sander-beck et al. 2016). Given the fact that our calculations of the tem-perature do not match the observations at z < , any conclusionsdrawn in this section should be interpreted with caution.We find that the H07 model and our best-fit model with He II data are in excellent agreement with the observations at z (cid:38) .On the other hand, models with a more quasar contribution, i.e.,the MH15 and the best-fit model without He II data, predict con-siderably higher temperature than what is observed at these red-shifts. These models have considerably earlier He II reionizationand hence the IGM is significantly photoheated at (cid:46) z (cid:46) .The same trends can be seen from the evolution of tempera-ture at the mean density T ( z ) ( right panel). The two peaks seenin all the models correspond to the hydrogen and He II reionizationrespectively. Also shown are the data points from the same obser-vation by Becker et al. (2011), but adjusted accordingly using thetemperature-density relation of the best-fit model with He II data(solid blue lines). Models with higher β heat up the IGM at muchearlier epoch than the what the observations show. Our main find- MNRAS , 000–000 (0000) Mitra, Choudhury & Ferrara z T ( ∆ ) / K ∆ = without HeII data (best-fit)with HeII data (best-fit)H07 ( β = − . )MH15 ( β = − . ) z T / K Figure 4. Thermal history for four different reionization models with different AGN emissivities; model with β = − . (orange dashed), β = − . (blue solid), β = − . (black dotted) and β = − . (black dot-dashed line). Left panel : evolution of IGM temperature at different overdensities printedat the bottom of the plot. The data points shown here are from the Ly α forest temperature measurements by Becker et al. (2011). Right panel : temperatureof gas at mean IGM density ( ∆ = 1 ). Data points are inferred from the same temperature measurements, but corrected for the appropriate γ (slope of thetemperature-density relation) of our best-fit model with He II data (solid blue lines). ings are in excellent agreement with the recent study by D’Aloisioet al. (2017), where they conclude that the AGN-dominated sce-nario struggles to match the temperature measurements. The new Planck results along with the recent discovery of signifi-cantly high number of faint AGNs from multiwavelength deep sur-veys at higher redshifts (Fiore et al. 2012; Giallongo et al. 2015)demand careful investigation of the actual contribution from star-forming galaxies and quasars as reionization sources. In addition,various current studies suggest that the globally-averaged escapefraction of ionizing photons is very small ∼ a few percent (Madau& Dickinson 2014; Ma et al. 2015; Mitra et al. 2015; Hassan et al.2016). Based on this idea, a QSO dominated reionization scenariohas been explored in many contemporary studies (Madau & Haardt2015; Khaire et al. 2016; Yoshiura et al. 2016). On the other hand,Kim et al. (2015) and Jiang et al. (2016) claimed that the number ofspectroscopically identified faint AGNs may not be high enough tofully account for the reionization of the IGM. Recently, D’Aloisioet al. (2017) also indicates a possible drawback of such high AGNemissivity models. This makes the claim of quasar-dominated sce-nario debatable. Thus, it would be interesting to see how our data-constrained semi-analytical reionization model can distinguish be-tween the models with different comoving AGN emissivities. Inthis work, we have expanded our previous model (Mitra et al. 2015)and evaluated the impact of QSOs on reionization for differentcases. As the observations on AGN emissivities at z > are stillvery uncertain, we model the AGN contribution in a way that itsevolution at high redshifts ( z > ) will solely be determined by asingle parameter β . We check that, β = − . and β = − . canreasonably restore our old model with the observed quasar luminos-ity functions from Hopkins et al. (2007) and a quasar-dominated model with best-fit AGN emissivity from Madau & Haardt (2015)respectively. We denote the former case as H07 and the later one asMH15. We try to constrain the value of β against different obser-vations related to hydrogen and helium reionization: measurementsof effective H I Ly α optical depth τ eff , HI , redshift evolution of LLS d N LL / d z , reionization optical depth from Planck 2016 and He II Ly α forest data. • We find that if we do not include He II data in our analysis,a wide range of β -model can be allowed by the current observa-tions. Both the quasar-dominated ( β ∼ ) and stellar-dominated( β ∼ − ; negligible QSOs at z > ) scenarios are possible withinthe 2- σ confidence limits. Any variation in the quasar emissiv-ity is appropriately compensated by the radiation from the stellarsources, thus ensuring that the constraints arising from hydrogenreionization data are not violated. It is thus not possible to putany meaningful constraints on the quasar emissivity using hydro-gen reionization alone. • In contrast, using He II data in the analysis can put a tighterconstraint on β and rules out the possibility of quasar dominatedmodels of hydrogen reionization. This model supports the decade-old two-component view of cosmic reionization; quasars contributesignificantly at lower redshifts and later the star-forming earlygalaxies take over. An escape fraction of to from thosegalaxies at z (cid:38) are needed to provide the appropriate ionizingflux. The fraction of photoionization rate contributed by the quasarscannot be larger that ∼ ∼ at z = 6(4) . • The models where AGNs throughout dominate the reioniza-tion process (e.g. MH15 model) favour an early helium reionization(completed at z (cid:38) ) history. They produce a much more slowlyevolving He II effective optical depth than observation (Worsecket al. 2016) suggests. Also such models are inconsistent with theobservational measurements of IGM temperature (Becker et al. MNRAS000 Thermal history for four different reionization models with different AGN emissivities; model with β = − . (orange dashed), β = − . (blue solid), β = − . (black dotted) and β = − . (black dot-dashed line). Left panel : evolution of IGM temperature at different overdensities printedat the bottom of the plot. The data points shown here are from the Ly α forest temperature measurements by Becker et al. (2011). Right panel : temperatureof gas at mean IGM density ( ∆ = 1 ). Data points are inferred from the same temperature measurements, but corrected for the appropriate γ (slope of thetemperature-density relation) of our best-fit model with He II data (solid blue lines). ings are in excellent agreement with the recent study by D’Aloisioet al. (2017), where they conclude that the AGN-dominated sce-nario struggles to match the temperature measurements. The new Planck results along with the recent discovery of signifi-cantly high number of faint AGNs from multiwavelength deep sur-veys at higher redshifts (Fiore et al. 2012; Giallongo et al. 2015)demand careful investigation of the actual contribution from star-forming galaxies and quasars as reionization sources. In addition,various current studies suggest that the globally-averaged escapefraction of ionizing photons is very small ∼ a few percent (Madau& Dickinson 2014; Ma et al. 2015; Mitra et al. 2015; Hassan et al.2016). Based on this idea, a QSO dominated reionization scenariohas been explored in many contemporary studies (Madau & Haardt2015; Khaire et al. 2016; Yoshiura et al. 2016). On the other hand,Kim et al. (2015) and Jiang et al. (2016) claimed that the number ofspectroscopically identified faint AGNs may not be high enough tofully account for the reionization of the IGM. Recently, D’Aloisioet al. (2017) also indicates a possible drawback of such high AGNemissivity models. This makes the claim of quasar-dominated sce-nario debatable. Thus, it would be interesting to see how our data-constrained semi-analytical reionization model can distinguish be-tween the models with different comoving AGN emissivities. Inthis work, we have expanded our previous model (Mitra et al. 2015)and evaluated the impact of QSOs on reionization for differentcases. As the observations on AGN emissivities at z > are stillvery uncertain, we model the AGN contribution in a way that itsevolution at high redshifts ( z > ) will solely be determined by asingle parameter β . We check that, β = − . and β = − . canreasonably restore our old model with the observed quasar luminos-ity functions from Hopkins et al. (2007) and a quasar-dominated model with best-fit AGN emissivity from Madau & Haardt (2015)respectively. We denote the former case as H07 and the later one asMH15. We try to constrain the value of β against different obser-vations related to hydrogen and helium reionization: measurementsof effective H I Ly α optical depth τ eff , HI , redshift evolution of LLS d N LL / d z , reionization optical depth from Planck 2016 and He II Ly α forest data. • We find that if we do not include He II data in our analysis,a wide range of β -model can be allowed by the current observa-tions. Both the quasar-dominated ( β ∼ ) and stellar-dominated( β ∼ − ; negligible QSOs at z > ) scenarios are possible withinthe 2- σ confidence limits. Any variation in the quasar emissiv-ity is appropriately compensated by the radiation from the stellarsources, thus ensuring that the constraints arising from hydrogenreionization data are not violated. It is thus not possible to putany meaningful constraints on the quasar emissivity using hydro-gen reionization alone. • In contrast, using He II data in the analysis can put a tighterconstraint on β and rules out the possibility of quasar dominatedmodels of hydrogen reionization. This model supports the decade-old two-component view of cosmic reionization; quasars contributesignificantly at lower redshifts and later the star-forming earlygalaxies take over. An escape fraction of to from thosegalaxies at z (cid:38) are needed to provide the appropriate ionizingflux. The fraction of photoionization rate contributed by the quasarscannot be larger that ∼ ∼ at z = 6(4) . • The models where AGNs throughout dominate the reioniza-tion process (e.g. MH15 model) favour an early helium reionization(completed at z (cid:38) ) history. They produce a much more slowlyevolving He II effective optical depth than observation (Worsecket al. 2016) suggests. Also such models are inconsistent with theobservational measurements of IGM temperature (Becker et al. MNRAS000 , 000–000 (0000) osmic reionization after Planck II • On the other hand, the model constrained by He II data has nodifficulties in matching both the τ eff , HeII and temperature data. Thebest-fit model predicts a relatively late helium reionization endingat z ≈ . τ eff , HeII increases rapidly from z = 2 . to z = 3 . bya factor of ∼ . Our old H07 case is well within the 2- σ limits ofsuch models.As a final note, we should mention some caveats of ourmodel. The star-formation efficiency and escape fraction param-eters should depend on the halo mass and/or redshifts (howeverweakly). For simplicity, we take them as a constant throughout thiswork. A detailed probe of parameter space using the principal com-ponent analysis, at the expense of larger computational time, mightbe useful in this case. We notice that, our simplified λ mfp prescrip-tion starts to break down when it becomes comparable to the Hub-ble radius at low-redshift regime ( z (cid:46) ). Perhaps a more accuratemean free path calculation will be needful. In our model, followingMiralda-Escud´e et al. (2000), we assume that reionization is saidto be complete once all the low-density regions with ∆ < ∆ crit are ionized, while the higher density gas is still neutral. This postu-late might not be fully correct, as in a real scenario, the ionizationfactor can depend on the local intensity of ionizing photons and/orthe degree of self-shielding (Miralda-Escud´e et al. 2000; Choud-hury 2009). Furthermore, a proper account of the quasar proximityeffect (where the ionization rate around a quasar exceeds the back-ground flux) in our model might provide additional insights (Bajt-lik et al. 1988; Miralda-Escud´e & Rees 1994; Padmanabhan et al.2014; D’Aloisio et al. 2017). Finally, recent observations (Beckeret al. 2015) indicate a rapid dispersion in Ly- α forest effective op-tical depth at z > . D’Aloisio et al. (2016) argued that match-ing such large scatter requires a much shorter (factor of (cid:38) ) spa-tially averaged mean free path than the actual measured value fromWorseck et al. (2014) at high redshifts, assuming galaxies to be thedominant sources of ionizing photons (Davies & Furlanetto 2016).However, they also claimed that these opacity fluctuations mightbe driven by the residual inhomogeneities in the IGM temperaturearising from a patchy reionization process rather than by the ion-ization background itself. A more rigorous modeling of IGM alongwith a careful estimation of mean free path is needed in our anal-ysis to fully include these effects. Also, our model for temperatureevolution is probably not very accurate at low redshifts ( z < ).Additional sources of heating like blazars might also be required,which we have not considered in this work.Nonetheless, future discovery of faint quasars at high redshiftsmight be worthwhile to clearly distinguish between the models withdifferent AGN contribution. A refined reionization model alongwith the improved measurements on intergalactic He II Ly α ab-sorption spectra and high-redshift temperature measurements willbe needed to put a more decisive constraint on the helium reioniza-tion scenario. ACKNOWLEDGEMENTS We thank the anonymous referee for constructive suggestions thatimproved this paper. We would also like to thank Vikram Khairefor useful discussions. REFERENCES Bajtlik S., Duncan R. C., Ostriker J. P., 1988, ApJ, 327, 570Barkana R., Loeb A., 2001, Phys. Rep., 349, 125Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023Becker R. H., et al., 2001, Astron.J., 122, 2850Becker G. D., Bolton J. S., Haehnelt M. G., Sargent W. L. W., 2011, MN-RAS, 410, 1096Becker G. D., Hewett P. C., Worseck G., Prochaska J. X., 2013, MNRAS,430, 2067Becker G. D., Bolton J. S., Madau P., Pettini M., Ryan-Weber E. V., Vene-mans B. P., 2015, MNRAS, 447, 3402Bennett C. L., et al., 2013, ApJS, 208, 20Bolton J. S., Haehnelt M. G., 2007, MNRAS, 382, 325Bolton J. S., et al., 2011, MNRAS, 416, L70Bouwens R., Illingworth G., 2006, New Astronomy Reviews, 50, 152Bouwens R. J., et al., 2015a, ApJ, 803, 34Bouwens R. J., Illingworth G. D., Oesch P. A., Caruana J., Holwerda B.,Smit R., Wilkins S., 2015b, ApJ, 811, 140Bowler R. A. A., et al., 2014, arXiv:1411.2976,Bradley L. D., et al., 2012, ApJ, 760, 108Cai Z.-Y., Lapi A., Bressan A., De Zotti G., Negrello M., Danese L., 2014,ApJ, 785, 65Cen R., McDonald P., Trac H., Loeb A., 2009, ApJL, 706, L164Chornock R., et al., 2013, ApJ, 774, 26Choudhury T. R., 2009, Current Science, 97, 841Choudhury T. R., Ferrara A., 2005, MNRAS, 361, 577Choudhury T. R., Ferrara A., 2006a, preprint,( arXiv:astro-ph/0603149 )Choudhury T. R., Ferrara A., 2006b, MNRAS, 371, L55Choudhury T. R., Puchwein E., Haehnelt M. G., Bolton J. S., 2014,arXiv:1412.4790,Civano F., et al., 2011, ApJ, 741, 91Compostella M., Cantalupo S., Porciani C., 2013, MNRAS, 435, 3169Compostella M., Cantalupo S., Porciani C., 2014, MNRAS, 445, 4186D’Aloisio A., McQuinn M., Davies F. B., Furlanetto S. R., 2016, preprint,( arXiv:1611.02711 )D’Aloisio A., Upton Sanderbeck P. R., McQuinn M., Trac H., Shapiro P. R.,2017, MNRAS, 468, 4691Davies F. B., Furlanetto S. R., 2016, MNRAS, 460, 1328Dijkstra M., Wyithe J. S. B., Haiman Z., 2007, MNRAS, 379, 253Dixon K. L., Furlanetto S. R., 2009, ApJ, 706, 970Fan X., et al., 2006, AJ, 132, 117Faucher-Gigu`ere C.-A., Lidz A., Zaldarriaga M., Hernquist L., 2009, ApJ,703, 1416Finlator K., Oppenheimer B. D., Dav´e R., Zackrisson E., Thompson R.,Huang S., 2016, MNRAS, 459, 2299Fiore F., et al., 2012, A&A, 537, A16Furlanetto S. R., Oh S. P., 2009, ApJ, 701, 94Gallerani S., Choudhury T. R., Ferrara A., 2006, MNRAS, 370, 1401Giallongo E., Menci N., Fiore F., Castellano M., Fontana A., Grazian A.,Pentericci L., 2012, ApJ, 755, 124Giallongo E., et al., 2015, A&A, 578, A83Gleser L., Nusser A., Benson A. J., Ohno H., Sugiyama N., 2005, MNRAS,361, 1399Glikman E., Djorgovski S. G., Stern D., Dey A., Jannuzi B. T., Lee K.-S.,2011, ApJL, 728, L26Hassan S., Dav´e R., Finlator K., Santos M. G., 2016, MNRAS, 457, 1550Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731Iliev I. T., Shapiro P. R., McDonald P., Mellema G., Pen U.-L., 2008, MN-RAS, 391, 63Jiang L., et al., 2016, ApJ, 833, 222Khaire V., Srianand R., 2013, MNRAS, 431, L53Khaire V., Srianand R., Choudhury T. R., Gaikwad P., 2016, MNRAS, 457,4051Kim Y., et al., 2015, ApJL, 813, L35Komatsu E., et al., 2011, Astrophys.J.Suppl., 192, 18Kulkarni G., Choudhury T. R., 2011, MNRAS, 412, 2781MNRAS , 000–000 (0000) Mitra, Choudhury & Ferrara Lewis A., Bridle S., 2002, Phys. Rev. D, 66, 103511Loeb A., 2006, preprint, ( arXiv:astro-ph/0603360 )Loeb A., Barkana R., 2001, ARA&A, 39, 19Ma X., Kasen D., Hopkins P. F., Faucher-Gigu`ere C.-A., Quataert E., KereˇsD., Murray N., 2015, MNRAS, 453, 960Madau P., Dickinson M., 2014, ARA&A, 52, 415Madau P., Haardt F., 2015, ApJL, 813, L8McGreer I. D., Mesinger A., D’Odorico V., 2015, MNRAS, 447, 499McLeod D. J., McLure R. J., Dunlop J. S., Robertson B. E., Ellis R. S.,Targett T. T., 2014, arXiv:1412.1472,McQuinn M., 2016, ARA&A, 54, 313McQuinn M., Lidz A., Zaldarriaga M., Hernquist L., Hopkins P. F., DuttaS., Faucher-Gigu`ere C.-A., 2009, ApJ, 694, 842Mesinger A., Aykutalp A., Vanzella E., Pentericci L., Ferrara A., DijkstraM., 2015, MNRAS, 446, 566Miralda-Escud´e J., Rees M. J., 1994, MNRAS, 266, 343Miralda-Escud´e J., Haehnelt M., Rees M. J., 2000, ApJ, 530, 1Mitra S., Choudhury T. R., Ferrara A., 2011, MNRAS, 413, 1569Mitra S., Choudhury T. R., Ferrara A., 2012, MNRAS, 419, 1480Mitra S., Ferrara A., Choudhury T. R., 2013, MNRAS, 428, L1Mitra S., Choudhury T. R., Ferrara A., 2015, MNRAS, 454, L76Oesch P. A., et al., 2012, ApJ, 745, 110Oesch P. A., et al., 2014, ApJ, 786, 108Okamoto T., Gao L., Theuns T., 2008, MNRAS, 390, 920Ota K., et al., 2008, ApJ, 677, 12Ouchi M., et al., 2010, ApJ, 723, 869Padmanabhan H., Choudhury T. R., Srianand R., 2014, MNRAS, 443, 3761Planck Collaboration et al., 2014, A&A, 571, A16Planck Collaboration et al., 2016a, A&A, 594, A13Planck Collaboration et al., 2016b, A&A, 596, A107Planck Collaboration et al., 2016c, A&A, 596, A108Price L. C., Trac H., Cen R., 2016, preprint, ( arXiv:1605.03970 )Prochaska J. X., O’Meara J. M., Worseck G., 2010, ApJ, 718, 392Puchwein E., Pfrommer C., Springel V., Broderick A. E., Chang P., 2012,MNRAS, 423, 149Raskutti S., Bolton J. S., Wyithe J. S. B., Becker G. D., 2012, MNRAS,421, 1969Richards G. T., et al., 2006, ApJS, 166, 470Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015, ApJL,802, L19Samui S., Srianand R., Subramanian K., 2007, MNRAS, 377, 285Schenker M. A., Ellis R. S., Konidaris N. P., Stark D. P., 2014, ApJ, 795, 20Schirber M., Bullock J. S., 2003, ApJ, 584, 110Schroeder J., Mesinger A., Haiman Z., 2013, MNRAS, 428, 3058Sharma M., Theuns T., Frenk C., Bower R., Crain R., Schaller M., SchayeJ., 2016, MNRAS, 458, L94Sobacchi E., Mesinger A., 2013, MNRAS, 432, 3340Songaila A., Cowie L. L., 2010, ApJ, 721, 1448Totani T., et al., 2006, Pub. Astron. Soc. Japan, 58, 485Trac H., Cen R., Loeb A., 2008, ApJL, 689, L81Upton Sanderbeck P. R., D’Aloisio A., McQuinn M. J., 2016, MNRAS,460, 1885White R. L., Becker R. H., Fan X., Strauss M. A., 2003, AJ, 126, 1Worseck G., et al., 2014, MNRAS, 445, 1745Worseck G., Prochaska J. X., Hennawi J. F., McQuinn M., 2016, ApJ, 825,144Wyithe J. S. B., Loeb A., 2003, ApJ, 586, 693Wyithe J. S. B., Loeb A., 2005, ApJ, 625, 1Yoshiura S., Hasegawa K., Ichiki K., Tashiro H., Shimabukuro H., Taka-hashi K., 2016, preprint, ( arXiv:1602.04407 )Zaldarriaga M., Hui L., Tegmark M., 2001, ApJ, 557, 519 MNRAS000