Cosmological tests with the FSRQ gamma-ray luminosity function

Abstract

The extensive catalogue of gamma-ray selected flat-spectrum radio quasars (FSRQs) produced by Fermi during a four-year survey has generated considerable interest in determining their gamma-ray luminosity function (GLF) and its evolution with cosmic time. In this paper, we introduce the novel idea of using this extensive database to test the differential volume expansion rate predicted by two specific models, the concordance Λ cold dark matter (ΛCDM) and R_h = ct cosmologies. For this purpose, we use two well-studied formulations of the GLF, one based on pure luminosity evolution (PLE) and the other on a luminosity-dependent density evolution (LDDE). Using a Kolmogorov–Smirnov test on one-parameter cumulative distributions (in luminosity, redshift, photon index and source count), we confirm the results of earlier works showing that these data somewhat favour LDDE over PLE; we show that this is the case for both ΛCDM and R_h = ct. Regardless of which GLF one chooses, however, we also show that model selection tools very strongly favour R_h = ct over ΛCDM. We suggest that such population studies, though featuring a strong evolution in redshift, may none the less be used as a valuable independent check of other model comparisons based solely on geometric considerations.

methods: statistical, quasars: general, cosmology: theory, large-scale structure of Universe, gamma-rays: general

1 INTRODUCTION

The discovery of quasars at redshifts ≳ 6 (Fan et al. 2003; Jiang et al. 2007; Willott et al. 2007; Jiang et al. 2008; Willott et al. 2010a; Mortlock et al. 2011; Venemans et al. 2013; Banados et al. 2014; Wu et al. 2015) suggests that ≳ 10^9 − 10 M_⊙ supermassive black holes emerged only ∼900 Myr after the big bang, and only ∼500 Myr beyond the formation of Population II and Population III stars (Melia 2013a). Such large aggregates of matter constitute an enduring mystery in astronomy because these quasars could not have formed so quickly in ΛCDM without an anomalously high accretion rate (Volonteri & Rees 2006) and/or the creation of unusually massive seeds (Yoo & Miralda-Escudé 2004); neither of these has actually ever been observed. For example, Willott et al. (2010b) have recently demonstrated that no known high-z quasar accretes at more than ∼1 − 2 times the Eddington rate (see fig. 5 in their paper; see also Melia 2014).

This paper will feature two specific cosmologies – the aforementioned ΛCDM (the ‘standard’ or concordance model) and another Friedmann–Robertson–Walker solution known as the R_h = ct Universe (Melia 2007; Melia & Shevchuk 2012; Melia 2016). Our focus will be to explain the luminosity function (LF) of these quasars, particularly as they evolve towards lower redshifts. Part of the motivation for this comparative study is that, unlike ΛCDM, the R_h = ct model does not suffer from the time compression problem alluded to above (Melia 2013a). In this cosmology, cosmic reionization (starting with the creation of Population III stars) lasted from t ∼ 883 Myr to ∼2 Gyr (6 ≲ z ≲ 15), so ∼5 − 20 M_⊙ black hole seeds formed (presumably during supernova explosions) shortly after reionization had begun, would have evolved into ∼10¹⁰ M_⊙ quasars by z ∼ 6 − 7 simply via the standard Eddington-limited accretion rate. The R_h = ct Universe has thus far passed all such tests based on a broad range of cosmological observations, but already, this consistency with the age–redshift relationship implied by the early evolution of supermassive black holes suggests that an optimization of the quasar LF might serve as an additional powerful discriminator between these two competing expansion scenarios. The class of flat-spectrum radio quasars (FSRQs) is ideally suited for this purpose.

FSRQs are bright active galactic nuclei (AGNs) that belong to a subcategory of blazars. These represent the most extreme class of AGNs, whose radiation towards Earth is dominated by the emission in a relativistic jet closely aligned with our line of sight. The discovery of gamma-ray emission from these sources was an important confirmation of the prediction by Melia & Konigl (1989) that the particle dynamics in these jets ought to be associated with significant high-energy emission along small viewing angles with respect to the jet axis. It is still an open question exactly what powers the jet activity, but it is thought that the incipient energy is probably extracted from the black hole's spin, and is perhaps also related to the accretion luminosity. Major mergers might have enhanced the black hole growth rate and activity, which would have occurred more frequently in the early Universe. In this context, the blazar evolution may be connected with the cosmic evolution of the black hole spin distribution, jet activity and major merger events themselves, all of which may be studied via the LF and its evolution with redshift.

Recently, the Fermi Gamma-ray Space Telescope has detected hundreds of blazars from low redshifts out to z = 3.1, thanks to its high sensitivity (Abdo et al. 2010a). Based on the previous analysis of the FSRQ gamma-ray luminosity function (GLF), it is already clear that the GLF evolution is positive up to a redshift cut-off that depends on the luminosity (see, e.g. Padovani et al. 2007; Ajello et al. 2009; Ajello et al. 2012). But all previous work with this sample ignored a very important ingredient to this discussion – the impact on the GLF evolution with redshift from the assumed cosmological expansion itself. Our main goal in this paper is to carry out a comparative analysis of the standard ΛCDM and R_h = ct models using the most up-to-date sample of 408 FSRQs detected by the Fermi Large Area Telescope (LAT) over its four-year survey. We wish to examine the influence on the results due to the assumed background cosmology and, more importantly, we wish to demonstrate that the current sample of gamma-ray emitting FSRQs is already large enough for us to carry out meaningful cosmological testing. Throughout this paper, we will be directly comparing the flat ΛCDM cosmology with Ω_m = 0.315 and H₀ = 67.3 km s⁻¹ Mpc⁻¹, based on the latest Planck results (Planck Collaboration XXXI 2014), and the R_h = ct Universe, whose sole parameter – the Hubble constant – will for simplicity be assumed to have the same value as that in ΛCDM. We will demonstrate that these data already emphatically favour R_h = ct over ΛCDM, even without an optimization of H₀ for R_h = ct.

The outline of this paper is as follows. In Section 2, we will summarize the observational data, specifically the 3FGL catalogue (Acero et al. 2015), and describe how the gamma-ray luminosity is determined for each specific model. Section 3 will provide an account of the critical differences between these two cosmologies that directly impact the calculation of the GLF, and we discuss the currently preferred ansatz for this LF based on the most recent analysis of these data in Section 4. We present and discuss our results in Section 5, and conclude in Section 6.

2 OBSERVATIONAL DATA AND SOURCE SAMPLE

The third Fermi-LAT source catalogue (3FGL) provided by Acero et al. (2015) lists 3303 sources detected by Fermi-LAT during its four years of operation. These data include the source location and its spectral properties. A subset of these is the third LAT AGN catalogue (3LAC; Ackermann et al. 2015), containing 1591 AGNs of various types located at high Galactic latiude, i.e. ∣b∣ ≥ 10^o. Most of the detected AGNs are blazars, which consist of 467 FSRQs, 632 BL Lacs, 460 blazar candidates of uncertain type (BCUs), and 32 non-blazar AGNs. Removing the entries in 3LAC for which the corresponding gamma-ray sources were not associated with AGNs, had more than one counterpart or were flagged for other reasons in the analysis, Ackermann et al. (2015) reduced the AGN catalogue to a ‘clean’ sample of 1444 sources, including 414 FSRQs, 604 BL Lacs, 402 BCUs and 24 non-blazar AGNs. The energy flux distribution of all the Fermi sources may be seen in fig. 18 of Acero et al. (2015). The flux threshold in 3FGL is ≃ 3 × 10⁻¹² erg cm⁻² s⁻¹, lower than the value (≃ 5 × 10⁻¹² erg cm⁻² s⁻¹) in 2FGL and (≃ 8 × 10⁻¹² erg cm⁻² s⁻¹) in 1FGL. The sample above the 3FGL flux threshold is essentially complete (see fig. 18 of Acero et al. 2015). Note also that all of the FSRQs catalogued in 3FGL have measured redshifts. In this paper, we have chosen to use only the FSRQs from 3LAC, and not the BL Lacs, because of their greater redshift coverage and better sample completeness, both of which strengthen our statistical analysis.

Our sample is therefore comprised of the 414 FSRQs detected by Fermi with a test statistic (TS) ≥25 and latitude |b| ≥ 10^o. Their gamma-ray fluxes S_γ and photon indices Γ in the energy range 0.1–100 GeV are obtained from 3FGL¹ and their redshifts are from table 7 of Ackermann et al. (2015).² We calculate the gamam-ray luminosity L_γ of an FSRQ using the expression:  

\begin{equation} L_{\gamma }=\frac{4\pi D^2_{\rm L}(z)S_{\gamma }}{1+z}K, \end{equation}

(1)

where D_L(z) is the luminosity distance at redshift z and K is the K-correction for the observed fluxes using S_γ = S_{γ, obs}(1 + z)^{Γ − 1}. The corresponding photon flux F_γ (in units of photons cm⁻² s⁻¹) in the 0.1–100 GeV energy band used in this paper is readily obtained from the following expressions (see, e.g. Ghisellini, Maraschi & Tavecchio 2009; Singal, Ko & Petrosian 2014):  

\begin{eqnarray} F_{\gamma }= \left\lbrace \begin{array}{lcl}S_{\gamma } \frac{2-\Gamma }{1-\Gamma } \frac{1}{E_{0.1}}\frac{1-10^{3(1-\Gamma )}}{1-10^{3(2-\Gamma )}} & & {\rm (if\, \Gamma \ne 2.0)} \\\\ S_{\gamma } \frac{1}{1-\Gamma } \frac{E_{0.1}^{1-\Gamma }}{\ln (10^3)} (10^{3(1-\Gamma )}-1) & & {\rm (if\, \Gamma = 2.0)}\;, \end{array} \right. \end{eqnarray}

(2)

where E_0.1 = 1.602 × 10⁻⁴ erg (corresponding to 0.1 GeV) is the lower energy limit.

Clearly, a determination of L_γ requires the assumption of a particular cosmological model. A detailed description of the differences between the two models we are considering here, the concordance ΛCDM cosmology and the R_h = ct Universe, may be found in Melia (2012a,b, 2013a,b). The quantity most relevant to the analysis in this paper is the luminosity distance D_L(z), which in ΛCDM is given as  

\begin{eqnarray} D_{{\rm L}}^{\Lambda {\rm CDM}}(z)&=&\frac{c}{H_0} \frac{(1+z)}{\sqrt{|\Omega _k|}} {\rm sinn}\Bigg \lbrace \mid \Omega _{k}\mid ^{1/2} \nonumber \\ &&\left.\times \int _0^{z} \frac{{\rm d}z}{\sqrt{(1+z)^2(1+\Omega _m z)- z(2+z) \Omega _{\Lambda }}} \right\rbrace , \end{eqnarray}

(3)

where c is the speed of light and H₀ is the Hubble constant at the present time. In this equation, Ω_m ≡ ρ_m/ρ_c and Ω_Λ ≡ ρ_Λ/ρ_c are, respectively, the energy density of matter and dark energy written in terms of today's critical density (⁠|$\rho _{\rm c}\equiv 3c^2H_0^2/8\pi G$|⁠), and Ω_k is the spatial curvature of the Universe, appearing as a term proportional to the spatial curvature constant k in the Friedmann equation. Also, sinn is sinh  when Ω_k > 0 and sin  when Ω_k < 0. For a flat Universe with Ω_k = 0, equation (3) simplifies to the form (1 + z)c/H₀ times the integral. In the R_h = ct Universe, the luminosity distance is given by the much simpler expression:  

\begin{equation} D_{{\rm L}}^{R_{\rm h}=ct}(z)= R_{\rm h}(t_0)(1+z)\ln (1+z)\,, \end{equation}

(4)

where R_h(t₀) = c/H₀ is the gravitational horizon (equal to the Hubble radius) at the present time.

The luminosities in these two models are simply related according to the expression  

\begin{equation} L_{\gamma }^{R_{\rm h}=ct}(z) = L_{\gamma }^{\Lambda {\rm CDM}}(z) \times \left(\frac{D_{{\rm L}}^{R_{\rm h}=ct}(z)}{D_{{\rm L}}^{\Lambda {\rm CDM}}(z)}\right)^2\;. \end{equation}

(5)

Fig. 1 shows the resulting luminosity-redshift distribution for both ΛCDM (left-hand panel) and R_h = ct (right-hand panel). In this figure, the solid curves are calculated using equation (1) for the threshold flux S_{γ, limit} = 3.0 × 10^{− 12} erg cm⁻² s⁻¹, and a fixed photon index Γ = 2.44, which is the mean of the Γ distribution from 3LAC (Ackermann et al. 2015). For a given observed flux, S_γ, the inferred luminosity depends on the assumed cosmology through the model-dependent luminosity distance, D_L, as indicated in equation (5). The data points in Fig. 1 are therefore slightly different for the two models being compared here. However, since we are assuming the concordance model with Planck parameters (see below), it turns out that the ratio |$D_{{\rm L}}^{\Lambda {\rm CDM}}/D_{{\rm L}}^{R_{\rm h}=ct}$| is very close to 1 over the redshift range 0 ≲ z ≲ 3. One may see this in fig. 3 of Melia (2015), which plots this ratio for several values of the matter density Ω_m. For the Planck matter density Ω_m = 0.308 ± 0.012, |$D_{\rm L}^{\Lambda {\rm CDM}}$| is just a few per cent bigger than |$D_{{\rm L}}^{R_{\rm h}=ct}$| over the entire redshift range considered in this paper.

Figure 1.

The 3LAC FSRQ luminosity-redshift distribution for both ΛCDM (left-hand panel) and R_h = ct (right-hand panel). The solid curves are calculated using equation (1) for the threshold flux S_{γ, limit} = 3.0 × 10^{− 12} erg cm⁻² s⁻¹, and a fixed photon index Γ = 2.44, which is the mean of the Γ distribution from 3LAC (Ackermann et al. 2015).

3 MODEL COMPARISON

Given a fixed background cosmology, we constrain the model parameters of the FSRQs GLF (see Section 4) using the method of maximum likelihood evaluation (see, e.g. Chiang & Mukherjee 1998; Narumoto & Totani 2006; Ajello et al. 2009; Abdo et al. 2010b; Zeng, Yan & Zhang 2014). The likelihood function |$\mathcal {L}$| is defined by the expression  

\begin{equation} \mathcal{L} = \textrm{exp}(-N_{\rm{exp}})\; \prod\nolimits^{N_{\rm obs}}_{i=1} \frac{{\rm d}^3 N(z_i,L_{\gamma,i},\Gamma_i)}{{\rm d}z\,{\rm d}L_\gamma\,{\rm d}\Gamma}, \end{equation}

(6)

where d³N/dz dL_γ dΓ is the space density of FSRQs, which generally depends on the LF, ρ(z, L_γ); the intrinsic distribution of photon indices with a Gaussian dependence, |${\rm d}N/{\rm d}\Gamma \propto \rm {exp}(-\frac{(\Gamma -\mu )^2}{2\sigma ^2})$|⁠, where μ and σ are the Gaussian mean and dispersion, respectively; and the comoving volume element per unit redshift and unit solid angle, dV_com/dz. This space density may be expressed as  

\begin{eqnarray} \frac{{\rm d}^3 N}{{\rm d}z\,{\rm d}L_{\gamma }\,{\rm d}\Gamma } &=& \frac{{\rm d}^2 N}{{\rm d}L_\gamma \,{\rm d}V_{\rm com}} \times \frac{{\rm d}N}{{\rm d}\Gamma } \times \frac{{\rm d}V_{{\rm com}}}{{\rm d}z} \nonumber \\ &=& \rho (z,L_\gamma ) \times \frac{{\rm d}N}{{\rm d}\Gamma } \times \frac{{\rm d}V_{\rm com}}{{\rm d}z}\;. \end{eqnarray}

(7)

In equation (6), the quantity N_exp is the expected number of FSRQ gamma-ray detections,  

\begin{eqnarray} N_{\rm exp} &=& \int ^{\Gamma _{\rm max}}_{\Gamma _{\rm min}} \int ^{z_{\rm max}}_{z_{\rm min}} \int ^{L_{\gamma ,\rm max}}_{L_{\gamma ,\rm min}} \frac{{\rm d}^3 N}{{\rm d}z\,{\rm d}L_{\gamma }\,{\rm d}\Gamma } \nonumber \\ &&\quad \times\, \omega (F_{\gamma }, \Gamma ) {\rm d}z\,{\rm d}\Gamma {\rm d}L_{\gamma }\;. \end{eqnarray}

(8)

Based on the properties of the gamma-ray emitting FSRQ sample, we take Γ_min = 1.0, Γ_max = 3.0, z_min = 0.0, z_max = 6.0, L_{γ, min} = 1.0 × 10⁴³ erg s⁻¹, and L_{γ, max} = 1.0 × 10⁵² erg s⁻¹. The quantity ω(F_γ, Γ) is the detection efficiency, and represents the probability of detecting a FSRQ with the photon flux F_γ and photon index Γ (e.g. Atwood et al. 2009; Abdo et al. 2010c; Ajello et al. 2012; Zeng, Yan & Zhang 2013; Di Mauro et al. 2014a,b), where F_γ is strongly dependent on the photon index Γ (obtained from equations 1 and 2), and is also a function of the luminosity L_γ and redshift z. To estimate ω(F_γ, Γ), we use the method provided by Di Mauro et al. (2014a). The efficiency ω(F_γ, Γ) at a photon flux |$F_{\gamma }^k \in (F_{\gamma }^{k,{\rm min}}, F_{\gamma }^{k,{\rm max}})(k={1 \cdot \cdot \cdot \cdot N_1})$| and photon index Γ^l ∈ (Γ^{l, min}, Γ^{l, max})(l = 1 · · · ·N₂) may be estimated as  

\begin{equation} \omega (F_{\gamma }, \Gamma ) = (1+\eta ) \frac{N_{\rm sources}^{k,l}}{\Delta \Omega \int _{F_{\gamma }^{k,{\rm min}}}^{F_{\gamma }^{k,{\rm max}}} \frac{{\rm d}N}{{\rm d}F_{\gamma }} {\rm d}F_{\gamma } \int _{\Gamma ^{l,{\rm min}}}^{\Gamma ^{l,{\rm max}}} \frac{{\rm d}N}{{\rm d}\Gamma } {\rm d}\Gamma }\;, \end{equation}

(9)

where ΔΩ is the solid angle associated with |b| > 10^o, η is the incompleteness of the sample representing the ratio of unassociated sources to the total number of sources, and |$N_{\rm sources}^{k,l}$| is the number of selected sources. The integrated values of dN/dF_γ and dN/dΓ are, respectively, the intrinsic flux and index distributions of the sources. Note that the effect of the photon index on the detection efficiency is fully accounted for in our method, the details of which are discussed in the appendix of Di Mauro et al. (2014a). Here we assume that the intrinsic flux distributions in the low-flux band of FSRQs and blazars have the same power-law index, because the distribution logN–logS of FSRQs is flatter than that of blazars at low fluxes (see fig. 14 of Abdo et al. 2010c). Fig. 2 shows the efficiency for this sample of FSRQs evaluated using equation 9, in the region Γ ∈ [1.5, 3.2], compared with the 1FGL (Abdo et al. 2010c) and 2FGL (Di Mauro et al. 2014a) samples. This comparison shows that this method of evaluation is a simplified and effective way to find the efficiency, though the correct method would include a simulation and an estimate of the number of detected sources versus the number of simulated objects in different flux bins (Abdo et al. 2010c; Ackermann et al. 2016).

While maximization of the likelihood function |$\mathcal {L}$| is an appropriate and reliable method for optimizing the parameters of a given model, to determine statistically which of the models is actually preferred by the data it is now common in cosmology to use several model selection tools (see, e.g. Melia & Maier 2013, and references cited therein). These include the Akaike Information Criterion, |${\rm AIC}\equiv -2\ln \mathcal {L}+2n$|⁠, where n is the number of free parameters (Liddle 2007); the Kullback Information Criterion, |${\rm KIC}=-2\ln \mathcal {L}+3n$| (Cavanaugh 2004); and the Bayes Information Criterion, |${\rm BIC}=-2\ln \mathcal {L}+(\ln N_{\rm obs})n$|⁠, where N_obs is the number of data points (Schwarz 1978). A more quantitative ranking of models can be computed as follows. When using the AIC, with AIC_α characterizing model |$\mathcal {M}_\alpha$|⁠, the unnormalized confidence in |$\mathcal {M}_1$| is given by the Akaike weight exp( − AIC₁/2). The relative probability that |$\mathcal {M}_1$| is statistically preferred is  

\begin{equation} P(\mathcal {M}_1) = \frac{\rm {exp}(-AIC_1/2)}{\rm {exp}(-AIC_1/2)+\rm {exp}(-AIC_1/2)}\;. \end{equation}

(10)

The difference Δ_AIC ≡ AIC₂ − AIC₁ determines the extent to which |$\mathcal {M}_1$| is favoured over |$\mathcal {M}_2$|⁠. For Kullback and Bayes, the likelihoods are defined analogously. In using these model selection tools, the outcome Δ_AIC (and analogously for KIC and BIC) is judged to represent ‘positive’ evidence that model 1 is to be preferred over model 2 if Δ_AIC > 2.0. If 2 < Δ_AIC < 6, the evidence favouring model 1 is moderate, and it is very strong when Δ_AIC > 10.

In this paper, we have chosen to compare two specific models: the concordance ΛCDM cosmology, with Planck-measured prior values for all its parameters, and the R_h = ct Universe, whose sole parameter – the Hubble constant H₀ – is, for simplicity, assumed to have the same value as that in ΛCDM. The total number of free parameters in this study is therefore limited in both cases to the formulation of the GLF (see Section 4), which is common to both models. In other words, the model selection statistic we will be using, i.e. Δ (for AIC, KIC, or BIC, as the case may be), depends solely on the quantity |$W\equiv -2\ln \mathcal {L}$| and in fact, for this reason, all three of these information criteria share the same values of Δ. As it turns out, W represents the χ² distribution. Transforming to the standard expression, we may therefore write  

\begin{equation} W=-2 \,\displaystyle\sum\limits_i^{N_{\rm obs}}\;\ln \frac{{\rm d}^3 N(z_i,L_{\gamma,i},\Gamma_i)}{{\rm d}z\,{\rm d}L_{\gamma}\,{\rm d}\Gamma} + 2 N_{\rm exp}. \end{equation}

(11)

For each model, we optimize the GLF parameters using the Markov Chain Monte Carlo (MCMC) technology, which is widely applied to give multidimensional parameter constraints from observational data. In practice, this means we will find the parameter values that minimize W, which yields the best-fitting parameters and their associated 1σ errors. We have adapted the MCMC code from cosraymc (Liu et al. 2007), which itself was adapted from the cosmomc package (Lewis & Bridle 2002). Additional details about the MCMC method may be found in Gamerman (1997) and Mackay (2003). We remark here that this way of posing the maximum likelihood problem is different from that chosen by Ajello et al. (2009, 2012, 2014) but, as pointed out in Ajello et al. (2009), who tested these various approaches, one gets exactly the same results using these different formulations, so there is no preference for one over the other, except in terms of convenience.

One of the principal differences between the two models affecting the value of W is the space density of FSRQs due to its dependence on the luminosity distance D_L(z). The comoving differential volume is given as follows:  

\begin{equation} \frac{{\rm d}V_{\rm com}}{{\rm d}z} = D_{\rm com}^2 \frac{{\rm d}D_{\rm com}}{{\rm d}z}\;, \end{equation}

(12)

where D_com ≡ D_L/(1 + z) is the comoving distance. For flat ΛCDM and R_h = ct we have, respectively,  

\begin{equation} \frac{{\rm d}D_{\rm com}^{\Lambda {\rm CDM}}}{{\rm d}z} = \frac{c}{H_0} \frac{1}{\sqrt{(1+z)^2(1+\Omega _{\rm m} z)- z(2+z) \Omega _{\Lambda }}}\;, \end{equation}

(13)

and  

\begin{equation} \frac{{\rm d}D_{\rm com}^{R_{\rm h}=ct}}{{\rm d}z} = \frac{c}{H_0} \frac{1}{1+z}.\; \end{equation}

(14)

For illustration, we show in Fig. 3 the space density (equation 7) as a function of redshift z, using the two formulations of the GLF described below, i.e. for pure luminosity evolution (PLE) and for luminosity-dependent density evolution (LDDE). The choice of parameters in this example is based on the discussion in Section 5.

Figure 3.

Space density of gamma-ray emitting FSRQs (equation 7) as a function of z for the concordance ΛCDM and R_h = ct cosmologies, using two formulations of the GLF (described in Section 4): pure luminosity evolution and luminosity-dependent density evolution. This illustration corresponds to two fixed parameters: L_γ = 1.0 × 10⁴⁸ erg s⁻¹ and Γ = 2.44.

4 THE GAMMA-RAY LUMINOSITY FUNCTION

The so-called PLE formulation for the GLF is motivated by the observed space density of radio-quiet AGNs, peaking at intermediate redshifts that correlate with source luminosity (see, e.g. Ueda et al. 2003; Hasinger, Miyaji & Schmidt 2005). This peak may correspond to the combined effect of black hole growth and the falloff in fuelling activity. The conventional formulation for the space density of this GLF (see, e.g. Ajello et al. 2012) has the form  

\begin{eqnarray} \rho (L_{\gamma },z)&=&\rho (L_{\gamma }/e[z]) \nonumber \\ &=&\frac{A\,e(z)}{\ln 10\;L_{\gamma }}\left[ \left(\frac{L_{\gamma }/e(z)}{L_{\ast }}\right)^{\gamma _1}\hspace{-7.22743pt}+\left(\frac{L_{\gamma }/e(z)}{L_{\ast }}\right)^{\gamma _2}\right],\quad \end{eqnarray}

(15)

where e(z) = (1 + z)^κe^z/ξ is the evolution factor correlated with source luminosity, A is a normalization factor, L_* is the evolving break luminosity, γ₁ is the faint-end slope index, γ₂ is the bright-end slope index, and κ and ξ represent the redshift evolution. Including the additional parameters μ and σ characterizing the (Gaussian) photon-index distribution (see discussion following equation 6) then results in a total of eight parameters that need to be optimized in our PLE analysis.

However, while the PLE GLF generally provides a good fit to the observed redshift and luminosity distributions, it is a very poor representation of the observed log N–log S (Ajello et al. 2012). Closer scrutiny of the values of κ and ξ in different redshift bins suggests that there is a significant shift in the redshift peak, with the low- and high-luminosity samples peaking at ∼1.15 and ∼1.77, respectively.

Since the simple PLE GLF may not be a completely adequate fit to the Fermi data, and since the redshift peak apparently evolves with luminosity, it is also beneficial to consider a GLF with LDDE (see Ueda et al. 2003; Ajello et al. 2012). In this formulation, the GLF evolution is decided by a redshift cut-off that depends on luminosity. The space density for this GLF is given by the following expression  

\begin{eqnarray} \rho (L_{\gamma },z)&=&\frac{A}{\ln 10\;L_{\gamma }} \left[\left(\frac{L_{\gamma }}{L_{\ast }}\right)^{\gamma _1}+ \left(\frac{L_{\gamma }}{L_{\ast }}\right)^{\gamma _2}\right]^{-1} \nonumber \\ &&\times \left[\left(\frac{1+z}{1+z_c(L_{\gamma })}\right)^{p_1}+ \left(\frac{1+z}{1+z_c(L_{\gamma })}\right)^{p_2}\right]^{-1},\qquad \end{eqnarray}

(16)

with  

\begin{equation} z_c(L_{\gamma })\equiv z_{c}^{\ast }(L_{\gamma }/10^{48})^{\alpha }\;, \end{equation}

(17)

where A is a normalization factor, L_* is the evolving break luminosity, γ₁ and p₁ are the faint-end slope indices, γ₂ and p₂ are the bright-end slope indices, |$z_{c}^{\ast }$| is the redshift peak with luminosity 10⁴⁸ ergs s⁻¹, and α is the power-law index of the redshift-peak evolution.

Previous studies based solely on the concordance ΛCDM model (see, e.g. Ajello et al. 2012) have shown that the LDDE provides a good fit to the LAT data and can reproduce the observed distribution quite well. The log-likelihood ratio test strongly favours it over PLE. In this paper, we will use both formulations of the GLF, just to be sure that we are not biasing our results prematurely with an ansatz for ρ(L_γ, z) that is too specific. As it turns out, both the PLE and LDDE formulations give completely consistent results when it comes to model selection. We will therefore conclude that the results of our model comparison using the gamma-ray emitting FSRQs is not at all dependent on assumptions concerning the form of the GLF.

5 RESULTS AND DISCUSSION

5.1 Pure luminosity evolution

We begin with the PLE assumption and optimize the GLF parameters by maximizing the likelihood function |$\mathcal {L}$| using the MCMC method. The best-fitting parameters and one-dimensional (1D) probability distributions are shown in Fig. 4 for ΛCDM and Fig. 5 for R_h = ct. The mean-fit parameter values and their 1σ confidence levels are listed in Table 1.

Figure 4.

The 1D probability distributions of the GLF parameters and their 1σ statistical uncertainties for the concordance ΛCDM cosmology, assuming PLE.

Figure 5.

The 1D probability distributions of the GLF parameters and their 1σ statistical uncertainties for the R_h = ct cosmology, assuming PLE.

Figure 6.

Top: the predicted cumulative distributions (in luminosity, redshift and photon index) from the luminosity functions with best-fitting parameters, versus the observed distributions of the corresponding quantities, assuming a PLE GLF, for both ΛCDM and the R_h = ct Universe. Bottom: to bring out the model differences more clearly, we also plot in the bottom three panels the ratio of predicted to observed distributions for the luminosity, redshift and photon index. The line definitions correspond to those in the upper panels. In these lower three panels, the top is for luminosity, the middle is for photon index, and the lowest is for redshift.

Table 1.

Optimized parameters (with 1σ errors) of the PLE GLF for ΛCDM and R_h = ct.

PLE	H0a	Ωm	log10Ab	γ1	log10L*	γ2	κ	ξ	μ	σ	W =
											|$-2\ln \mathcal {L}$|
ΛCDM	67.3	0.315	|$-10.24^{+0.36}_{-0.34}$|	|$0.57^{+0.06}_{-0.06}$|	|$47.75^{+0.31}_{-0.31}$|	|$1.90^{+0.53}_{-0.47}$|	|$4.76^{+0.63}_{-0.62}$|	|$-0.53^{+0.07}_{-0.07}$|	|$2.44^{+0.01}_{-0.01}$|	|$0.18^{+0.01}_{-0.01}$|	89 106
Rh = ct	67.3	–	|$-10.11^{+0.40}_{-0.39}$|	|$0.56^{+0.08}_{-0.08}$|	|$47.64^{+0.36}_{-0.37}$|	|$1.81^{+0.52}_{-0.46}$|	|$5.08^{+0.64}_{-0.65}$|	|$-0.48^{+0.06}_{-0.06}$|	|$2.44^{+0.01}_{-0.01}$|	|$0.18^{+0.01}_{-0.01}$|	88 962
KS Test	PD(Lγ)	PD(z)	PD(Γ)
ΛCDM	98.4 per cent	67.5 per cent	72.9 per cent
Rh = ct	98.6 per cent	73.4 per cent	74.9 per cent

^aIn units of km s⁻¹ Mpc⁻¹.

^bIn units of Mpc⁻³ erg⁻¹ s.

Ajello et al. (2012) examined whether the PLE GLF could adequately account for the Fermi data, though strictly only for the ΛCDM cosmology, and concluded that it was not a good representation of the log N–log S distribution. To see whether this is still true for our sample, and also to test whether this deficiency is also present for the R_h = ct Universe, we apply the Kolmogorov–Smirnov (KS) test for the predicted one-parameter cumulative distributions using the measured populations as individual functions of redshift, luminosity and photon index. The theoretical one-parameter distributions are calculated as follows:  

\begin{eqnarray} N(<z) &=& \int ^{\Gamma _{\rm max}}_{\Gamma _{\rm min}} \int ^{z}_{z_{\rm min}} \int ^{L_{\gamma , \rm max}}_{L_{\gamma , \rm min}} \frac{{\rm d}^3 N}{{\rm d}z\,{\rm d}L_{\gamma }\,{\rm d}\Gamma } \nonumber \\ &&\qquad \times\, \omega (F_{\gamma }\,\Gamma )\, {\rm d}z\,{\rm d}\Gamma \, {\rm d}L_{\gamma }\;, \end{eqnarray}

(18)

 

\begin{eqnarray} N(<L_{\gamma }) &=& \int ^{\Gamma _{\rm max}}_{\Gamma _{\rm min}} \int ^{z_{\rm max}}_{z_{\rm min}} \int ^{L_{\gamma }}_{L_{\gamma , \rm min}} \frac{{\rm d}^3 N}{{\rm d}z\,{\rm d}L_{\gamma }\,{\rm d}\Gamma } \nonumber \\ &&\qquad \times\, \omega (F_{\gamma }\,\Gamma )\, {\rm d}z\,{\rm d}\Gamma \, {\rm d}L_{\gamma }\;, \end{eqnarray}

(19)

 

\begin{eqnarray} N(<\Gamma ) &=& \int ^{\Gamma }_{\Gamma _{\rm min}} \int ^{z_{\rm max}}_{z_{\rm min}} \int ^{L_{\gamma , \rm max}}_{L_{\gamma , \rm min}} \frac{{\rm d}^3 N}{{\rm d}z\,{\rm d}L_{\gamma }\,{\rm d}\Gamma } \nonumber \\ &&\qquad \times\, \omega (F_{\gamma }\,\Gamma )\, {\rm d}z\,{\rm d}\Gamma \, {\rm d}L_{\gamma }\;. \end{eqnarray}

(20)

In addition, the source-count distribution is given by the expression  

\begin{eqnarray} N(> F_{\gamma }) &=& \int ^{\Gamma _{\rm max}}_{\Gamma _{\rm min}} \int ^{z_{\rm max}}_{z_{\rm min}} \int ^{L_{\gamma , \rm max}}_{\rm {Max}(L_{\gamma }[F_{\gamma }, z, \Gamma ], L_{\gamma , \rm min})} \frac{{\rm d}^3 N}{{\rm d}z\,{\rm d}L_{\gamma }\,{\rm d}\Gamma }\nonumber \\ && \qquad \times\, \omega (F_{\gamma }\,\Gamma )\, {\rm d}z\,{\rm d}\Gamma \, {\rm d}L_{\gamma }\;, \end{eqnarray}

(21)

where L_γ(F_γ, z, Γ) is the luminosity of a source at redshift z with photon index Γ, having a flux F_γ. The corresponding curves, together with the binned data with which they are compared, are shown in Figs 6and 7 for both ΛCDM and R_h = ct.

At least visually, the predicted cumulative distributions appear to match the data quite well, aside from the source-count distribution. Indeed, the PLE GLF passes the KS test in three cumulative distributions (luminosity, redshift, photon index) for both cosmologies, though only at a modest level of confidence in the case of Γ. The R_h = ct Universe does better than the concordance ΛCDM model for all the distributions. As we shall see shortly, the results of our KS comparison between the cumulative distributions for PLE and LDDE are somewhat mixed. Certainly in the case of z, the LDDE GLF passes the KS test with a significantly higher level of confidence, where it reaches ∼98 per cent in the case of R_h = ct, compared to only ∼74 per cent for PLE. However, the KS test results for the cumulative distributions in L_γ are very similar between PLE and LDDE. Note that the left-hand panels in Fig. 7 show that the predicted distributions are a very poor representation of the observed logN–logS. Based solely on the cumulative redshift distribution, together with the relatively poor source-count distribution, we do confirm the result in Ajello et al. (2012), that the LDDE GLF appears to be a better representation of the Fermi data than the PLE GLF.

5.2 Luminosity-dependent density evolution

We next optimize the 10 parameters of the LDDE GLF, given in equation (16), using the MCMC method to maximize the likelihood function for the same sample of 414 FSRQs that we used for PLE. The 1D probability distributions of these parameters are shown in Fig. 8 for ΛCDM and Fig. 9 for R_h = ct. Their mean-fit values and 1σ confidence levels are listed in Table 2. As before, and specifically to examine which GLF is a better match to the data, we carry out the same KS test as for PLE, and compare the predicted one-parameter cumulative distributions with the data in Fig. 10. The favourable visual impression one gets is confirmed by the confidence levels of the matches, which are quoted in Table 2.

Figure 8.

The 1D probability distributions of the GLF parameters and their 1σ statistical uncertainties for the concordance ΛCDM cosmology, assuming LDDE.

Figure 9.

The 1D probability distributions of the GLF parameters and their 1σ statistical uncertainties for the R_h = ct cosmology, assuming LDDE.

Figure 10.

Top: the predicted cumulative distributions (in luminosity, redshift and photon index) from the luminosity functions with best-fitting parameters, versus the observed distributions of the corresponding quantities for the LDDE GLF, for both ΛCDM and the R_h = ct Universe. Bottom: as in Fig. 6, we also show here the ratio of predicted to observed distributions for luminosity, redshift and photon index, corresponding the distributions in the upper panels. The respective confidence levels are listed in Table 2.

Table 2.

Optimized parameters (with 1σ errors) of the LDDE GLF for ΛCDM and R_h = ct.

LDDE	H0a	Ωm	log10Ab	γ1	log10L*	γ2	|$z_{c}^{\ast }$|	α	p1	p2	μ	σ	W =
													|$-2\ln \mathcal {L}$|
ΛCDM	0.673	0.315	|$-8.78^{+0.14}_{-0.14}$|	|$0.32^{+0.05}_{-0.05}$|	|$47.93^{+0.17}_{-0.17}$|	|$1.71^{+0.23}_{-0.23}$|	|$2.06^{+0.19}_{-0.19}$|	|$0.20^{+0.02}_{-0.02}$|	|$9.71^{+2.53}_{-2.36}$|	|$-4.15^{+0.99}_{-1.00}$|	|$2.44^{+0.01}_{-0.01}$|	|$0.19^{+0.01}_{-0.01}$|	89 047
Rh = ct	0.673		|$-8.61^{+0.14}_{-0.14}$|	|$0.31^{+0.06}_{-0.06}$|	|$47.78^{+0.17}_{-0.17}$|	|$1.70^{+0.21}_{-0.21}$|	|$2.06^{+0.19}_{-0.18}$|	|$0.19^{+0.02}_{-0.02}$|	|$9.64^{+2.31}_{-2.18}$|	|$-4.23^{+0.99}_{-1.00}$|	|$2.44^{+0.01}_{-0.01}$|	|$0.19^{+0.01}_{-0.01}$|	88 903
KS Test	PD(Lγ)	PD(z)	PD(Γ)
ΛCDM	99.1 per cent	78.2 per cent	98.0 per cent
Rh = ct	99.5 per cent	81.2 per cent	98.3 per cent

^aIn units of km s⁻¹ Mpc⁻¹.

^bIn units of Mpc⁻³ erg⁻¹ s.

The R_h = ct cosmology does at least as well as ΛCDM, and usually better, in all the KS tests using the various one-parameter cumulative FSRQ distributions. In both models, the predicted source-count distribution is a better match to the data for LDDE than for PLE (the right-hand panels of Fig. 7), supporting the conclusion drawn earlier by Ajello et al. (2012) that LDDE is favoured over PLE by the measured log N–log S_γ relation.

Figure 7.

The source-count distribution of FSRQs for both cosmologies, assuming PLE (left-hand panels) and LDDE (right-hand panels). The curves are the best-fitting models reported in the text for the ΛCDM (solid) and R_h = ct (dashed) cosmologies. (a) The cumulative distributions in photon flux. (b) Solid circles represent the observed cumulative distribution in photon flux for the 186 FSRQs with TS ≥ 50 and |b| ≥ 15^o reported in Ajello et al. (2012); open circles represent the observed cumulative distribution for the 414 FSRQs in our sample. For comparison, the (additional) thick curves in the lower panels are the intrinsic cumulative distributions assuming ω(F_γ Γ) = 1.0.

To complete our discussion, we also summarize here a comparison of our results with others reported in the literature. Since LDDE appears to be strongly favoured by the data over PLE, we will focus our attention on this particular GLF. Fig. 11 compares the differential local (z=0) and z=1 GLFs with those reported by Ajello et al. (2012) and Singal et al. (2014). Ajello et al. (2012) analysed the LF by using the sample comprised of 186 FSRQs detected by Fermi with TS ≥ 50, |b| ≥ 15^o and F_γ ≥ 10⁻⁸ photons cm⁻² s⁻¹. The results of Singal et al. (2014) were obtained by analysing the sample of 184 FSRQs with TS ≥ 50, |b| ≥ 20^o reported by Shaw et al. (2012). We can see that our distributions have a normalization approximately two times larger than theirs. This is merely a reflection of the fact that our sample (414) is about two times bigger than theirs (186 and 184); other than this obvious difference, our results for the local Universe are virtually identical to theirs. The right-hand plot in Fig. 11 shows some slight differences in the determination of the GLF at z = 1, possibly due to the different redshift distributions of the various samples used for the optimization of the model parameters or the incompleteness of the earlier samples. In this regard, we note from the bottom panels of Fig. 7 that the observed data of our sample are concordant with those of Ajello et al. (2012) at high fluxes, but they clearly differ in the low-flux region. This would confirm the fact that our results should be the same at low redshifts, but differ with Ajello et al. (2012) at high redshifts, as is evident in Fig. 11.

Figure 11.

Differential local (left-hand side) and z=1 (right-hand side) gamma-ray luminosity functions for FSRQs assuming LDDE. The filled and dashed curves represent the results in this paper. Stars are the results of Singal et al. (2014), who restricted their analysis solely to FSRQs with a ≳ 7σ detection threshold in the first-year catalogue of the Fermi LAT. The dotted curves are the FSRQ LFs reported by Ajello et al. (2012).

Figure 2.

Efficiency for the sample of FSRQs evaluated using equation (9), in the region Γ ∈ [1.5, 3.2] for this work (indicated by red stars), and the detection efficiency of the 1FGL (Abdo et al. 2010c) and 2FGL (Di Mauro et al. 2014a) samples.

5.3 Model comparisons

We now turn to the main goal of our analysis, which is to directly compare these two cosmologies, for which we must use the model selection tools discussed in Section 3 above. Starting with the PLE GLF (Section 5.1), the values of W (from which Δ is calculated) are shown in Table 1. Our optimization procedure shows that |$\Delta \equiv W_{\Lambda {\rm CDM}} -W_{R_{\rm h}=ct}$| (for this assumed PLE LF) is ∼134, well into the ‘very strong’ category. At least for the PLE GLF, the size of our gamma-ray emitting FSRQ sample is already large enough for this statistical assessment to overwhelmingly favour R_h = ct over the concordance ΛCDM model.

As we have seen, the LDDE GLF is a significantly better match to the data than PLE. Here too, the model selection tools very strongly favour R_h = ct over ΛCDM. In the case of LDDE, |$\Delta \equiv W_{\Lambda {\rm CDM}}- W_{R_{\rm h}=ct}=144$|⁠, which again is well into the ‘very strong’ category. The choice of GLF does not appear to have much influence in deciding which of these two cosmologies is favoured by the Fermi FSRQ data. The sample is already large enough for the observations to strongly prefer the differential volume dependence on z predicted by R_h = ct over that in ΛCDM.

6 CONCLUSIONS

The extensive, high-quality sample of gamma-ray emitting FSRQs observed by Fermi has generated considerable interest in identifying the GLF and its evolution with cosmic time. The number density of such objects has changed considerably during the expansion of the Universe, growing dramatically up to redshift ∼0.5 − 2.0 and declining thereafter. Aside from the obvious benefits one may derive from better understanding this evolution as it relates to supermassive black hole growth and its connection to the haloes of host galaxies, its strong dependence on redshift all the way out to z ∼ 3 offers the alluring possibility of using it to test different cosmological models.

In this paper, we have introduced this concept by directly comparing two specific expansion scenarios, chiefly to examine the viability of the method. To do so, we have opted to use prior values for the model parameters themselves, and instead focus on the optimization of the parameters characterizing the chosen ansatz for the LF. In doing so, one may question whether the choice of GLF unduly biases the fit for one model or the other. This is a legitimate concern, and considerable work still needs to be carried out to ensure that one is not simply customizing the GLF for each background cosmology.

For this reason, we have opted in this paper to use two different forms of the GLF, one for PLE and the second for a LDDE, even though earlier work had already established a preference by the data for the latter over the former. We have found that selecting either of these GLFs has no influence at all on the outcome of model comparison tools. In both cases, information criteria, such as the AIC, KIC, and BIC, show quite conclusively that the evolution of the GLF for FSRQs very strongly favours R_h = ct over the concordance ΛCDM model.

Cosmic evolution is now studied using a diversity of observational data, including high-z quasars (Melia 2013a, 2014), Gamma-ray bursts (Wei, Wu & Melia 2013), the use of cosmic chronometers (Melia & Maier 2013; Melia & McClintock 2015), Type Ia supernovae (Wei et al. 2015) and, most recently, an application of the Alcock–Paczyński test using model-independent baryon acoustic oscillation (BAO) data (Font-Ribera et al. 2014; Delubac et al. 2015; Melia & López-Corredoira 2015), among others. The BAO measurements are particularly noteworthy because, with their ∼4 per cent accuracy, they now rule out the standard model in favour of R_h = ct at better than the 99.34 per cent C.L.

In this paper, we have provided a compelling confirmation of these other results by demonstrating that population studies, though featuring a strong evolution in redshift, may also be used to independently check the outcome of model comparisons based purely on geometric considerations. We emphasize, however, that much work still needs to be done to properly identify how to best characterize the number density function for this type of analysis. This would be critically important in cases, unlike ΛCDM and R_h = ct, where cosmological models are so different that an appropriate common ansatz may be difficult to find.

We thank the referee, Mattia Di Mauro, for a careful reading of our manuscript, and for thoughtful comments that have led to an improved presentation, including several clarifying descriptions of the results. We acknowledge the use of cosraymc (Liu et al. 2012) adapted from the cosmomc package (Lewis & Bridle 2002). FM is grateful to Amherst College for its support through a John Woodruff Simpson Lectureship, and to Purple Mountain Observatory in Nanjing, China, for its hospitality while part of this work was being carried out. LZ acknowledges partial funding support by the National Natural Science Foundation of China (NSFC) under grant No. 11433004. This work was partially supported by grant 2012T1J0011 from The Chinese Academy of Sciences Visiting Professorships for Senior International Scientists, and grant GDJ20120491013 from the Chinese State Administration of Foreign Experts Affairs. This work is also supported by the Key Laboratory of Particle Astrophysics of Yunnan Province (Grant 2015DG035). This work is also partially supported by the Strategic Priority Research Program, the Emergence of Cosmological Structures, of the Chinese Academy of Sciences, Grant No. XDB09000000, and the NSFC grants 11173064, 11233001, and 11233008.

1

" URL: http://fermi.gsfc.nasa.gov/ssc/data/access/lat/4yr_catalog/

2

" URL: http://www.asdc.asi.it/fermi3lac/