Source Count Distribution of Fermi LAT Gamma-Ray Blazars Using Novel Nonparametric Methods

We utilized a sample from the Fermi-LAT 14-year Source Catalog by adjusting the flux detection threshold, enabling us to derive the intrinsic source count distribution $dN/dF_{25}$ of extragalactic blazars using nonparametric, unbinned methods developed by Efron and Petrosian and Lynden-Bell. Subsequently, we evaluated the contribution of blazars to the extragalactic gamma-ray background. Our findings are summarized as follows: (1) There is no significant correlation between flux and spectral index values among blazars and their subclasses FSRQs and BL Lacs. (2) The intrinsic differential distributions of flux values exhibit a broken-power-law form, with parameters that closely match previous findings. The intrinsic photon index distributions are well described by a Gaussian form for FSRQs and BL Lacs individually, while a dual-Gaussian model provides a more appropriate fit for blazars as a whole. (3) Blazars contribute 34.5\% to the extragalactic gamma-ray background and 16.8\% to the extragalactic diffuse gamma-ray background. When examined separately, FSRQs and BL Lacs contribute 19.6\% and 13\% to the extragalactic gamma-ray background, respectively.


Introduction
The extragalactic gamma-ray background (EGB) refers to diffuse gamma-ray radiation originating beyond our Milky Way galaxy.This background consists of contributions from various astrophysical sources, such as active galactic nuclei (AGNs), star-forming galaxies, and potentially dark matter annihilation or decay.Understanding the composition and origin of the EGB is crucial for advancing our knowledge of high-energy astrophysics and cosmology.A prominent category of astrophysical objects, blazars, which encompass both flat-spectrum radio quasars (FSRQs) and BL Lacertae objects, represent a dominant class of radio-loud AGNs detected by the Large Area Telescope (LAT) of the Fermi Gamma-Ray Space Telescope (see, e.g., [1][2][3]).These blazars are highly variable and luminous, emitting intense gamma-ray radiation due to relativistic jets nearly aligned with Earth's line of sight (see, e.g., [4][5][6][7][8]).
Generally, the contribution of blazers to the EGB is typically estimated using the luminosity function and its evolution with various parameters, as seen in, e.g., [9][10][11][12][13][14][15].However, determining these distributions often requires variables such as redshift (z), which are not readily available for most blazars, especially BL Lacs objects.Conversely, blazars' contribution to the EGB can be directly estimated if the source count distributions of gamma-ray emitters (logS-logN) are provided, as seen in, e.g., [16][17][18][19].Examining the population statistics of gamma-ray sources can provide valuable insights into the cumulative emissions from these objects and their overall impact on the EGB.
Determining the distributions of flux and photon spectral index values is complex due to the limitations of Fermi-LAT in relation to its threshold flux values.For instance, blazars with harder spectra are more easily detected at lower flux levels.Additionally, each detected source has a different threshold flux, determined by a detection test statistic (TS), because the background flux varies with the position on the sky [20].Observational biases have been analyzed and mitigated using Monte Carlo simulations (see, e.g., [16,21]).Moreover, Refs.[17,18] employed a method developed by Efron and Petrosian (EP) [22] and Lynden-Bell [23] to directly determine the distributions from observational data, accounting for complex observational selection biases, initially identifying the intrinsic correlations between variables and subsequently establishing the univariate distribution of each variable.The advantage of these methods lies in their avoidance of ad hoc assumptions and their direct extraction of the mono-variate distributions from data in a nonparametric manner, without resorting to binning.Nonparametric methods have successfully contributed to the study of AGNs (e.g., [17,18,[24][25][26][27]) and gamma-ray bursts (e.g., [28][29][30][31]).
In this paper, we apply nonparametric methods to the truncated data detected by Fermi-LAT in order to ascertain the intrinsic correlations between photon flux and the photon index, followed by analyzing their intrinsic distributions.In contrast to Refs.[17,18], we have employed a novel sample selection method to procure suitable samples and assess the detection efficiency of blazars.We note that the study in [18] focuses on variable luminosity and redshift, analyzing the evolution of the gamma-ray and optical luminosity functions of blazars.In contrast, here, we do not deal with cosmological evolution, and instead focus solely on a detailed analysis of the intrinsic distribution of flux and the photon index.This paper is organized as follows: Section 2 details our selection of an available sample of data from the Fermi-LAT extragalactic catalog.Section 3 focuses on analyzing these data to obtain the intrinsic distribution of flux and photon spectral index values and subsequently calculating the contribution of blazars to the EGB.Section 4 provides a brief summary and discussion.

Sample
The fourth Fermi-LAT AGN catalog, 4LAC [32], comprises 2863 AGNs of various types detected by Fermi-LAT operating between 4 August 2008 and 2 August 2016, situated at galactic latitude |b| ≥ 10 o .After excluding entries in 4LAC not associated with known AGNs, those with multiple counterparts, or flagged for other reasons in the 4FGL FITS file 1 , a clean sample of 2614 sources emerged, including 591 FSRQs, 1027 BL Lacs, 941 blazars of unknown type (BCUs), and 55 non-blazar AGNs.Analysis of the energy flux distributions for high-latitude sources in 4LAC (Figure 4 of [32]) has suggested that the detection threshold approximates a truncation line where the LAT energy flux threshold S 25,lim rises with increasing spectral index Γ.This line was thus chosen for use in this study because the updated sample from the Fermi group does not provide the expected detection threshold.
An updated version (4FGL−DR4, for Data Release 4) of the fourth full catalog of LAT sources, based on 14 years of survey data in the 50 MeV-1 TeV energy range, is now available online 2 .Initially, a sample was selected from the LAT 14-year Source Catalog, available as a FITS file 3 , under the conditions of a source significance over 4.3 σ (TS > 25, fitted with a power-law function) and analysis flags equal to 0, resulting in 2796 blazars, including 619 FSRQs, 1239 BL Lacs, and 938 BCUs.For our analysis, we required the photon flux ranging from 100 MeV to 100 GeV and the photon index.Although the photon flux is not provided in the sample, it can be derived using the following expression: where E 1 = 1.602 × 10 −4 erg (corresponding to 100 MeV) represents the lower energy limit.F 25 and S 25 4 , respectively, denote the photon flux and the energy flux from 100 MeV to 100 GeV, and Γ represents the absolute value of the photon index.The scatter diagram in the left panel of Figure 1 illustrates the photon flux F 25 versus the power-law index Γ, with the dashed line representing the expected detection threshold with 8-year Fermi-LAT data (F 25,lim vs. Γ lim ).This curve underscores the detection threshold's substantial dependence on the power-law index (from Figure 16 of [2]).Due to the significant number of BCUs, they may be classified as either FSRQs or BL Lacs, leading to substantial classification incompleteness in the samples of both types; i.e., causing the classification fraction of source , where i represents FSRQs or BL Lacs, to be low.To mitigate this incompleteness of the FSRQ and BL Lac sample, we introduced a simple translation of the detection threshold by constructing a shift factor α, such that F25,lim = αF 25,lim for a fixed Γ lim .All sources with F 25 > F25,lim ( or S 25 > Ś25,lim ) were selected until the classification fraction reached 0.9 by varying the value of α.Consequently, we obtained two samples: 234 FSRQs corresponding to α = 5.0 and 466 BL Lacs corresponding to α = 4.0.The FSRQ sample is designed to include 90% of sources identified as FSRQs by excluding known BL Lacs, while the BL Lac sample is designed to encompass 90% of the sources identified as BL Lacs by removing known FSRQs.For the blazar sample, we used a more conservative limit with α = 2.87 and a fraction N fsrq +N bll N fsrq +N bll +N Bcu = 0.9.This resulted in the selection of 1003 blazars, ensuring that 90% of all sources could be identified as either FSRQs or BL Lacs.Therefore, we utilize three different samples in the following analysis.The right panel of Figure 1 shows the relationship between the fraction of sources and the translation factor.Γ lim ), which reflects a larger dependence of the detection threshold on the power-law index (from Figure 16 of [2]).Right: The fraction of source N bll +N Bcu vs. the shift factor α. The black, red, and blue lines show the fraction of blazar, FSRQ, and BL Lac objects with the shift factor, respectively.When the shift factor α = 2.87, the fraction of blazars reaches 0.9, while the shift factors of FSRQs and BL Lacs need to be equal to 5.0 and 4.0.

Analysis and Results
The aim of this study is to determine the intrinsic source count distribution of blazars and evaluate their contribution to the EGB.The previous section provides us with a description of the sample, which includes bi-variate truncated data.This requires investigating the independence between the photon flux and photon index.If these two variables are independent, their distributions, represented as G(F 25 , Γ), can be separated into two distinct distributions: ψ(F 25 ) and h(Γ).However, in scenarios like gamma-ray bursts, there may be a persistent correlation between luminosity and the photon index due to the intrinsic relationship between flux and photon index values, even after mitigating cosmic influences [33,34].Therefore, to isolate the intrinsic distributions within our sample, it is crucial to eliminate the correlation between flux and the photon index.

Correlations
Here, we employ the Efron-Petrosian method to assess the correlation between two variables.This method presents a variation of the Kendall Tau statistic test tailored for truncated data, utilizing a test statistic to examine the independence of two variables within a dataset, as follows: where R i is the normalized ranks of the sources within their associated sets, defined as sources with Γ j < Γ i and F j > F i,lim , for ranking in F. Further details on this test can be found in [22,24,27,29].Under the null hypothesis of independence, the expected ranks and variances are E i = 1/2 and V i = 1/12.Our findings indicate that for both subsamples, τ(k = 0) > 5, rejecting the null hypothesis at a significance level exceeding 5σ.To investigate the correlation between photon flux and photon index, we follow the methodology delineated in [17].In this approach, a new parameter, termed the correlation-reduced photon index (Γ cr ), is introduced, defined as follows: This transformation is essentially a simple coordinate rotation.The parameter β is empirically determined such that F 25 and Γ cr are independent (τ = 0), allowing the combined distribution to be expressed as G(F 25 , Γ) = ψ(F 25 ) × ĥ(Γ cr ).Once the photon flux distributions are determined, the intrinsic distribution of the photon index (h(Γ)) is obtained through integration over flux:

Distributions
Upon acquiring the truncated data with independence in the F 25 − Γ cr plane, we apply the Lynden-Bell method, as developed by Lynden-Bell [23], to derive the cumulative mono-variate distributions Ψ(F 25 ) and H(Γ cr ).This method provides a nonparametric and point-by-point description of the cumulative distributions, conventionally sorted in descending order by photon flux and ascending order by photon index (see also, e.g., Zeng et al. [18], Petrosian et al. [29]).Then, their distributions are represented as follows: where N j represents the number of objects in the associated set of object i (with variables F 25,i , Γ cr,i ) consisting of sources with F 25,j > F 25,i and Γ cr,j < Γ cr,max,i (or F min,j < F 25,i ); and M j denotes the number of objects in the associated set of object i comprising sources with Γ cr,j < Γ cr,i and F 25,j > F min,i (or Γ cr,max,j > Γ cr,i ).The cumulative distribution of flux and photon index values can also be expressed as integrals of the independent distributions ψ(F 25 ) and ĥ(Γ cr ), as follows: Figure 3 illustrates the cumulative distributions of photon flux F 25 and index Γ cr corrected for the photon flux limit, alongside the corresponding raw cumulative distributions of F 25 and Γ cr , denoted as N(> F 25 ) and N(< Γ cr ), respectively.The derivatives of the cumulative distributions yield the following differential distributions: , and ĥ(Γ cr ) = dH(Γ cr ) dΓ cr .
Their distributions can be obtained directly by differentiation of the numerical histograms, as described in Appendix B of [29].However, a point-by-point derivative may result in noise, which can be mitigated by smoothing.For flux, employing interpolation smoothing provides the differential distributions shown by the points with error bars in the left panel of Figure 4. To facilitate comparison with the results of previous studies (e.g., [16,17]), we assume that the intrinsic distribution of flux follows a broken power-law model: where A is the normalization and F b is the flux break.β 1 and β 2 are the power-law slopes above and below the break, respectively.Parameters obtained via Markov Chain Monte Carlo (MCMC) fitting of the cumulative distribution of photon flux are listed in Table 1, computed using Equation ( 6).This method necessitates providing the likelihood function L, where we define the χ 2 = −2lnL and Here, Y i represents N(> 3 × 10 −8 )/10.38 sr otherwise.The truncation of data below F 25 < 3 × 10 −8 photons cm −2 s −1 is significant, warranting error scaling relative to Ψ(F 25 = 3 × 10 −8 ).This error treatment method was initially introduced in [17], providing comprehensive details.Notably, the total sky coverage considered is 10.38 sr, corresponding to |b| > 10 o .The left panels of Figures 3 and 4 illustrate that these distributions are well characterized by a broken power-law form for both blazars and their subsets.
Regarding the photon index, Equation ( 4) calculates the intrinsic distribution of the photon index, h(Γ), using the parameterized distribution ψ(F 25 ) and numerical differentiation, resulting in ĥ(Γ cr ).The right panel of Figure 4 demonstrates the intrinsic distributions of the photon index and the raw observed distributions for blazars, FSRQs and BL Lacs.It should be noted that the best-fit value of the correlation parameter β is utilized here.A Gaussian model was applied to fit the intrinsic distribution of the photon index.Interestingly, the distribution for blazars deviates from a Gaussian distribution, while those for FSRQs and BL Lacs are well described.Notably, the intrinsic photon index values' distribution for blazars seems better represented by two Gaussian distributions, suggesting a non-Gaussian nature.
The best-fit model of the intrinsic source count distribution is a broken power-law model for blazars, FSRQs, and BL Lacs, with the break points found at F b = 8.39 +1. 16 −1.07×  4, represented by the thin and the thick solid lines.The best-fit parameters are shown in the Table 1, indicating general agreement with the findings of Refs.[16,17].
We note that β 1 = 2.31 +0.05 −0.06 for BL Lacs deviates from the expected value of 2.5 for a Euclidean distribution in the local universe.Additionally, the break flux F b in the source count distribution for BL Lacs is lower than that reported in [16,17].Abdo et al. [16] proposed that this break is of cosmological origin, coinciding with changes in population evolution sign.BL Lac objects demonstrate positive cosmological evolution at high luminosities and negative evolution at low luminosities [10], potentially leading to a source count distribution following a three-component power law.In the left panels of Figures 3 and 4, it can be seen that there appears to be a break flux at F b ∼ 6.0 × 10 −8 ph cm −2 s −1 .

Detection Efficiencies
Once the intrinsic differential photon flux distribution has been accurately determined, the detection efficiency of the sample can be calculated as follows: Here, we utilize the theoretical distribution (dN/dF 25 ) th rather than the smoothed differential distribution; specifically, we employ the best-fit distribution based on a broken power-law model.The observed distribution, (dN/dF 25 ) obs , encompasses the entire sample from 14 years of Fermi-LAT operations.Due to significant incompleteness in the entire sample for FSRQs and BL Lac objects when analyzed separately, we focus solely on calculating the detection efficiency for all 2796 blazars.Both (dN/dF 25 ) th and (dN/dF 25 ) obs are illustrated in the left panel of Figure 4.The results are depicted in Figure 5 as a black pentagram, along with the black curve described by ω(F) = (1 + (F/F 0 ) −w 1 ) −w 2 , where F 0 = 2.88 × 10 −9 ph cm −2 s −1 , w 1 = 2.25, and w 2 = 1.04, fitted using the least-squares method [15].Generally, detection efficiency is computed by providing lists of the simulated and detected sources for each efficiency simulation (see, e.g., [16,21,35]).Our method is also frequently employed for estimating detection efficiency, as seen in [11,15].It is noteworthy that our result is similar to that of [21], since the ratio √ 14/ √ 10 = 1.18 closely approximates to unity.9)) between the distribution of detected sources and the theoretical distribution.The red squares denote data obtained from [16].The purple dots indicate the detection efficiency of the LAT high-latitude survey (|b| > 20 o ) with the simulated skies from [21], which is compared with our results.

Contribution to the EGB and IGRB
Having derived the source count distribution and the detection efficiency, we now estimate the contribution of blazars to the EGB and the isotropic diffuse gamma-ray background (IGRB).The composition of the IGRB is dominated by undetected Fermi sources.We adopt a computational approach similar to that of Singal et al. [17] and Zeng et al. [18], which avoids the necessity for redshift evolution data or reliance on the luminosity function.The differential background photon intensity from blazars can be determined as follows: where f Γ (E) is the spectrum of sources.In the range of 100 MeV to 100 GeV, most blazars show a power-law spectrum, and the observed spectrum consists of a power law with exponential cutoff due to the optical depth τ(E, z), f z) , where z is the redshift.Ignoring the effect of the optical depth τ(E, z) in the E 1 = 100 MeV to 100 GeV range, we can write ), where Substituting this into Equation (10), we obtain with where P(E) is extracted from Equation (10) and encompasses the integration of the spectral index, which characterizes this particular spectral shape property.We can obtain an analytic expression for P(E) assuming a Gaussian distribution for Γ with mean value of Γ 0 and dispersion σ, which gives According to the analysis in [18], this expression can be written as where ) and b = 9.8σ, and the second derivative Therefore, the integrated intensity photon energy above 100 MeV can be written as Fermi-LAT has previously measured the intensities of the EGB and IGRB above 100 MeV as 1.13 ± 0.07 × 10 −5 ph cm −2 s −1 sr −1 and 7.2 ± 0.6 × 10 −6 ph cm −2 s −1 sr −1 , respectively, using 50 months of LAT data from 5 August 2008 and 6 October 2012 [36].According to Figure 8 in [37], the minimum photon flux detected from blazars by LAT over 48 months is approximately 5.0 × 10 −10 ph cm −2 s −1 with a spectral index Γ ∼ 1.5.Thus, employing Equation (17) with F min = 5.0 × 10 −10 ph cm −2 s −1 , F max = 1.0 × 10 −5 ph cm −2 s −1 , and the best-fit parameters in Table 1, we can estimate the contribution of blazars to the EGB and IGRB.The calculated contributions are 3.90 × 10 −6 ph cm −2 s −1 sr −1 for the EGB and 1.21 × 10 −6 ph cm −2 s −1 for the IGRB, representing 34.5 and 16.8 percent of the intensities observed by Fermi-LAT, respectively.The predicated EGB and IGRB spectra of blazars, using Equation (12) with the different intrinsic distributions of the indexes, are shown in the left panel of Figure 6.Extrapolating to zero flux, we find that blazars in total can produce a photon output of 8.27 × 10 −6 ph cm −2 s −1 sr −1 to the EGB and 5.72 × 10 −6 ph cm −2 s −1 to the IGRB.
The detection efficiencies for FSRQs and BL Lacs are not obtained in this work; we solely estimated their contribution to the EGB.Our results indicate the contribution from FSRQs to be 2.22 × 10 −6 ph cm −2 s −1 sr −1 , which accounts for 19.6 percent of Fermi's EGB intensity, and the contribution from BL Lacs to be 1.47 × 10 −6 ph cm −2 s −1 sr −1 , representing 13 percent of the Fermi EGB intensity.Consequently, their combined contribution amounts to 32.6 percent of the Fermi EGB intensity.The right panel of Figure 6 illustrates the predicated EGB spectra of FSRQs, with z = 1.0, and BL Lacs, with z = 0.3.When using Equation (17) with F min = 0 to assess their maximum contributions to the EGB, FSRQs contribute 2.36 × 10 −6 ph cm −2 s −1 sr −1 and BL Lacs contribute 1.77 × 10 −6 ph cm −2 s −1 sr −1 .When the jointly fitted spectral index distribution is utilized, the intensity of the EGB increases slightly due to the harder spectrum.As shown in right panel of Figure 6, the contribution from FSRQs exceeds that of BL Lacs at energies below approximately 600 MeV, whereas the contribution from BL Lacs dominates the EGB at energies around 600 MeV.Notably, this boundary energy is higher than the value reported in [16] (300 MeV).According to Equation (17), with E min = 50 GeV and F min = 0, the total integrated flux from BL Lacs is found to be 1.98 × 10 −9 ph cm −2 s −1 sr −1 , which accounts for about 78% of the EGB above 50 GeV.This estimate is lower than that reported in [15], at about 100%, which was obtained through employing the gamma-ray luminosity function of BL Lacs, and lower than that reported in [19], at 84 +16 −14 %, which was obtained after the authors analyzed a sample from the second catalog of hard Fermi-LAT sources (2FHL) consisting of 90% BL Lacs.Furthermore, the contribution of FSRQs to the EGB is approximately one-tenth that of BL Lacs, about 1.5 × 10 −10 ph cm −2 s −1 sr −1 ; this disparity can be attributed to the larger spectral index and significant absorption effects at higher redshifts.
It is noteworthy that here, we utilize the EBL model of [38] to calculate τ(E, z), which requires a redshift parameter.However, we do not consider redshift evolution of sources but aim to select characteristic redshifts to represent the average evolutionary effects or assume no evolution with a constant redshift.Base on the redshift distribution of 3LAC [37] and 4LAC [32], FSRQs exhibit a broad peak around z = 1, while the predominant peak occurs at z = 0.3 for BL Lacs.Therefore, we employ specific redshift values of 0.5, 1.0, and 0.3 to calculate τ(E, z) for the entire dataset, for FSRQs, and for BL Lacs, respectively.As depicted in Figure 6, the predicted EGB spectrum exhibits an earlier cutoff compared to the observed data, suggesting that the characteristic redshift of blazars and BL Lacs may be lower than those aforementioned values.The predicated EGB (solid line) and IGRB (dashed line) spectra for blazars, obtained by using the optical depth model of [38] with z = 0.5.The black and blue lines represent the intrinsic distribution of indexes used in the calculations as single Gaussian and double Gaussian, respectively.Right: The predicated EGB spectra for FSRQs with z = 1.0 (dashed lines) and BL Lacs with z = 0.3 ( dotted lines).The black solid line is the total EGB spectrum of the subset of blazars.The blue lines represent the results of the EGB spectrum calculated using the jointly fitted spectral index distribution.The observed EGB and IGRB data are taken from [36].

Summary and Discussion
In this study, we have employed the rigorous nonparametric methods of Efron-Petrosian and Lynden-Bell (EP-L) to calculate the source count distribution of blazars based on gamma-ray photon fluxes and spectral indexes from the Fermi-LAT 14-year AGN catalog with TS ≥ 25 at galactic latitude |b| > 10 o .To address the significant presence of blazars with uncertain classifications, we mitigated the substantial classification incompleteness in subclass samples by adjusting the Fermi-LAT detection threshold based on flux, thus establishing a sample with 90% classification completeness for comprehensive analysis.We investigated the correlation between photon flux and spectral index values across three blazar samples-blazars as a whole, FSRQs, and BL Lacs-and derived their intrinsic distributions in order to estimate their contributions to the EGB.Below is a concise summary of our procedures and results: 1.
Flux-index relations: Addressing the bias towards hard sources in the photon flux limit (Figure 1), we applied the EP method to correct for this bias, resulting in a correlationcorrected spectral index Γ cr .Our analysis revealed no significant correlation for all blazars (β = 0.006 ± 0.033), with only a weak correlation being observed for FSRQs (β = −0.100+0.031 −0.038 ) and BL Lacs (β = −0.017+0.037 −0.045 ) individually.

2.
Intrinsic distributions: Employing the Lynden-Bell method, we derived cumulative distributions and subsequently obtained differential distributions through numerical derivation.The true intrinsic differential distributions of flux exhibited a broken power law for blazars, FSRQs, and BL Lacs, individually.The intrinsic photon index distributions were described well by Gaussian forms for FSRQs and BL Lacs individually; however, when considering the sample as a whole, double Gaussians provided a better fit.Table 1 summarizes the best-fit parameters for the intrinsic flux and photon index distributions, comparing them with the results of previous studies (i.e., [16,17]).

3.
Detection efficiency: The source count distribution at |b| > 10 o , shown in the left panel of Figure 4, displayed a broken power law for F 25 > 1.0 × 10 −8 ph cm −2 s −1 .Below this flux threshold, the observed distribution dropped rapidly due to the challenges in detecting fainter sources with Fermi-LAT.We calculated the detection efficiency by comparing observed sources with the theoretical distribution in photon flux.

4.
EGB and IGRB: Through the nonparametric determination of intrinsic distributions of photon flux and the spectral index, along with considerations of the detection efficiency, we directly calculated the contributions of all Fermi-LAT blazars to the EGB's and IGRB's intensity and spectra.Our findings indicate that blazars up to the flux threshold of Fermi-LAT, i.e., F min = 5.0 × 10 −10 ph cm −2 s −1 , could explain 34.5% of observed EGB photons and 32.6% of observed IGRB photons.For FSRQs and BL Lacs, we estimated their contributions to the EGB without considering the detection efficiency, revealing that FSRQs and BL Lacs could account for 19.6% and 13% of the observed EGB, respectively.

Figure 1 .
Figure 1.Left: Scatter diagram of the photon flux F 25 vs. the power-law index Γ.The truncation curve (dashed line) shows the expected detection threshold with 8-year Fermi-LAT data (F 25,lim vs. Γ lim ), which reflects a larger dependence of the detection threshold on the power-law index (from

Here, F 0
= 3.0 × 10 −8 photon cm −2 s −1 serves as a fiducial flux, with its specific value being noncritical.It is essential to transform the truncation curve Γ lim in the F 25 − Γ plane using the same parameter β, Γ cr,lim = Γ lim − β × log 10 ( F 25 F 0 ).The left panel of Figure 2 shows the variation in τ with the correlation parameter β, displaying the best-fit values and 1-sigma ranges of β blazars = 0.006 ± 0.033 for blazars, β fsrq = −0.100+0.031 −0.038 for FSRQs, and β bll = −0.017+0.037 −0.045 for BL Lacs.The results indicate a weak correlation or no correlation for all sources, consistent with the findings of Singal et al. [17].The right panel of Figure 2 shows a scatter diagram of the photon flux F 25 vs. the correlation-reduced photon index Γ cr .

Figure 2 .
Figure 2. Left: The value of the test statistic τ as a function of β for blazars, FSRQs, and BL Lacs, which indicates a correlation between the photon index and flux.The values of β for which τ = 0 and ±1 give the best value and 1-sigma range, respectively.A weak correlation or no correlation are indicated with β blazars = 0.006 ± 0.033 for blazars, β fsrq = −0.100+0.031 −0.038 for FSRQs, and β bll = −0.017+0.037 −0.045 for BL Lacs.Right: Scatter diagram of photon flux F 25 vs the correlation-reduced photon index Γ cr .The black, red, and blue curves show the flux threshold limit αF 25 vs. Γ cr (the truncation curve) for blazars, FSRQs, and BL Lacs, respectively, with the different shift factor.The two boxes enclosed by the dashed (red) and dotted (purple) lines, respectively, represent two associated sets (M j and N j ) of the source marked by the green star, as defined in Equation (6).

Figure 3 .
Figure 3. Left: The cumulative photon flux distribution (solid points), Ψ(F 25 ), corrected for the photon flux limit, and compared with the raw histogram of photon flux (open points), N(> F 25 ).The solid curves display the cumulative distributions of flux achieved by integrating the broken power-law distribution (Equation (8)) using the best-fit parameters listed in Table1.Right: The cumulative photon index distribution H(Γ cr ) vs. Γ cr corrected for selection effects (solid points) and compared with the raw cumulative observed photon index (open points), Γ, N(> Γ).
Notes.EP-L: the results of[17]; MA: the results of[16].a The correlation between photon index Γ and flux F 25 -see Equation (3) and Section 3.1.b The intrinsic flux distribution ψ(F 25 ) at break flux F b , in units of cm 2 s 1 sr −1 .c The power-law slope of the intrinsic flux distribution ψ(F 25 ) at fluxes above the break in the distribution.b The flux at which the power-law break in ψ(F 25 ) occurs, in units of 10 −8 ph cm −2 s −1 (0.1 ≤ E ≤ 100 GeV).e The power-law slope of the intrinsic flux distribution ψ(F 25 ) at fluxes below the break in the distribution.f The mean of the Gaussian fit to the intrinsic photon index distribution h(Γ).g The 1σ width of the Gaussian fit to the intrinsic photon index distribution h(Γ).The first and second rows represent the results of individual and joint fits, respectively.Joint fitting entails constraining the widths of the two Gaussians in the blazar sample to equal those in the respective individual BL Lac and FSRQ samples.

Figure 4 .
Figure 4.The differential intrinsic photon flux distribution (left) and photon index distribution (right) for blazars, FSRQs, and BL Lacs, separately.The dN/dF 25 values derive from the differentiation of the cumulative photon flux distribution for corrected and raw data, and the curves derive from Equation (8) with the best-fit parameters.The intrinsic dN/dΓ = h(Γ) values are shown as solid points, which derive from Equation (4), and the raw observed distributions are shown as open points.The Gaussian distribution form is used to fit h(Γ) and show the best-fit curves.The raw distribution of observations with 2796 blazars is also shown in the left panel (the open stars).

Figure 5 .
Figure 5.The detection efficiency ω(F 25 ) for blazars, calculated by determining the ratio (Equation (9)) between the distribution of detected sources and the theoretical distribution.The red squares denote data obtained from[16].The purple dots indicate the detection efficiency of the LAT high-latitude survey (|b| > 20 o ) with the simulated skies from[21], which is compared with our results.

Figure 6 .
Figure6.Left: The predicated EGB (solid line) and IGRB (dashed line) spectra for blazars, obtained by using the optical depth model of[38] with z = 0.5.The black and blue lines represent the intrinsic distribution of indexes used in the calculations as single Gaussian and double Gaussian, respectively.Right: The predicated EGB spectra for FSRQs with z = 1.0 (dashed lines) and BL Lacs with z = 0.3 ( dotted lines).The black solid line is the total EGB spectrum of the subset of blazars.The blue lines represent the results of the EGB spectrum calculated using the jointly fitted spectral index distribution.The observed EGB and IGRB data are taken from[36].
25,i ) and N i represents N(> F 25,i )/10.38 sr, indicating the number of objects above F 25,i per steradian.The uncertainties associated with N i are typically treated as Poisson errors, represented by σ = √ N i /10.38 sr for F 25,i > 3 × 10 −8 photons cm −2 s −1 , and σ = Based on the fact that each source in the blazar sample is either an FSRQ or a BL Lac, a simultaneous joint fit is performed, constraining the means and widths of the two Gaussian distributions in the blazar sample to equal those in the respective individual BL Lac and FSRQ samples.The mean values of the two Gaussians

Table 1 .
The best-fit parameters for Fermi-LAT blazars of photon flux and index intrinsic distribution.