Constraints on redshifts of blazars from extragalactic background light attenuation using Fermi-LAT data

The extragalactic high-energy $\gamma$-ray sky is dominated by blazars, which are active galactic nuclei with their jets pointing towards us. Distance measurements are of fundamental importance yet for some of these sources are challenging because any spectral signature from the host galaxy may be outshone by the non-thermal emission from the jet. In this paper, we present a method to constrain redshifts for these sources that relies only on data from the Large Area Telescope on board the Fermi Gamma-ray Space Telescope. This method takes advantage of the signatures that the pair-production interaction between photons with energies larger than approximately 10 GeV and the extragalactic background light leaves on $\gamma$-ray spectra. We find upper limits for the distances of 303 $\gamma$-ray blazars, classified as 157 BL Lacertae objects, 145 of uncertain class, and 1 flat-spectrum-radio quasar, whose redshifts are otherwise unknown. These derivations can be useful for planning observations with imaging atmospheric Cherenkov telescopes and also for testing theories of supermassive black hole evolution. Our results are applied to estimate the detectability of these blazars with the future Cherenkov Telescope Array, finding that at least 21 of them could be studied in a reasonable exposure of 20 h.


INTRODUCTION
Active galactic nuclei that have jets pointing in the direction of our line of sight are known as blazars.These sources dominate the extragalactic sky at the energies detected by the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope, typically around 1 GeV (Ackermann et al. 2016;Abdollahi et al. 2020).
In this paper, we propose a methodology based solely on LAT data and EBL attenuation to set upper limits on the redshifts of blazars.This work is organized as follows.Section 2 includes a brief review on the theoretical aspects about -ray attenuation.In Section 3, we describe the -ray sample, data reduction, and methodology.This is followed by Section 4, where our results are shown and discussed.Finally, we summarize and conclude in Section 5.

THEORETICAL BACKGROUND ON 𝛾-RAY ATTENUATION
We refer the reader to other detailed papers on the theoretical background on -ray attenuation such as Dwek & Krennrich (2013) and Finke et al. (2022), and briefly summarize the main aspects here.
There is a pair production interaction between  rays and EBL photons that produces an optical depth (, ) that is analytically given by where  and  are respectively the observed energy and redshift of the -ray source,   is the photon-photon pair production cross section and ( ′ ,  ′ ) is the proper number density of EBL photons with rest-frame energy  ′ at redshift  ′ given by (2) where   is the EBL spectral intensity.The cross section   depends on the relative rest-frame energies of the -ray photon ( ′ ), the rest mass energy of the electron    2 , and the EBL photon energy ( ′ ) as where   ℎ is the energy threshold for photon-photon pair production, which is given by and variable  = (1 − cos ) connects the energy threshold to the angle of interaction .
While -ray absorption due to the EBL does not have a fixed onset and strongly depends on redshift, a pragmatic choice of 10 GeV as a threshold can be made when considering blazars in the  = 2 − 3 range.The observational effect is that this EBL attenuation leaves a characteristic signature on the -ray spectra of blazars (e.g., Ackermann et al. 2012;Abdollahi et al. 2018).

Source selection
For our analysis, we select all blazars in the 4FGL-DR2 catalogue (Ballet et al. 2020).There are several definitions in use to distinguish BL Lacertae objects (BL Lacs) and flat spectrum radio quasars (FS-RQs), but the most common one is as follows.Blazars are classified as BL Lacs if the equivalent width of their emission lines is < 5 Å and FSRQs otherwise (Marcha et al. 1996;Healey et al. 2008).BL Lacs are expected to be the best targets for EBL studies since these are the sources exhibiting a hard -ray spectrum and unlikely to have an external photon field, e.g., broad line region (BLR), that may lead to additional -ray attenuation (e.g., Domínguez & Ajello 2015).Although, Abdollahi et al. (2018) show that internal absorption is unlikely to be significant for FSRQs, a result confirmed by Costamante et al. (2018).In the LAT catalogues there are also sources classified as blazar candidates of uncertain type (BCUs), which are objects that have a broad-band spectral energy distribution typical of blazars but do not have a spectral classification (Ajello et al. 2020).
Emissions from blazars can present variability producing changes of flux and modifications in the spectrum (e.g., Peñil et al. 2020Peñil et al. , 2022;;Otero-Santos et al. 2022).These effects may lead to inconclusive results for information extracted from individual spectra using -ray attenuation (see for instance the section on Variability in the Supplementary Material by Abdollahi et al. 2018), therefore in our approach we will distinguish between variable and non-variable blazars.This distinction is made using the Variability_Index column in 4FGL-DR2 and using the condition lower or higher than 18.48 as explained by Abdollahi et al. (2020).A higher value corresponds to a 99% probability of the source being variable.
Before making any calculations, we expected that our results would be significant only for BL Lacs and BCUs that are probably BL Lacs, since these sources typically have the strongest EBL signature.However we also ran the method on FSRQs as a sanity check.The initial number of sources in our sample was (non-variable/variable blazars) 621/510 BL Lacs, 866/446 BCUs, and 129/565 FSRQs.

Data reduction and methodology
For each source in our sample, we develop a LAT analysis as follows: (i) We perform an analysis in the 1 GeV to 2 TeV range without an EBL model, in order to determine the parameters that define the spectral function such as index and flux, also known as spectral parameters, of the sources around the source of interest.The spectral functions are either PowerLaw or LogParabola models depending on the best fit provided by 4FGL (see below for details).
(ii) We perform an analysis in the 1-10 GeV range without an EBL model, using the source spectral parameters (other than the source of interest) frozen to the values from step (i).The upper limit of 10 GeV is adopted considering the fact that the EBL absorption is negligible below this energy (Ackermann et al. 2012;Abdollahi et al. 2018).
(iii) We extrapolate the spectrum from (ii) to 2 TeV, and include absorption with the EBL model, assuming a redshift of  = 0.01, and record the log-likelihood value.We repeat this step for all redshifts between  = 0.01 and  = 3.00 in 100 steps and save the loglikelihood value at each redshift.
(iv) We compare the log-likelihood values from step (iii) with the one from step (i) to create a TS profile.The peak of this profile gives the most likely value for the redshift.We take this redshift to be an upper limit.Because there can be a component in the spectral curvature of the sources that is not produced by EBL attenuation but which is intrinsic to their emission.
A detailed description of this analysis is as follows.We perform a standard likelihood analysis on the Fermi-LAT data covering a given time period (first 14 years of telescope time), energy range (from 1 GeV to 2 TeV), and a selected region of interest (ROI) using fermiPy (Wood et al. 2017).The 4FGL catalogue of Fermi-LAT detected sources Data Release 2 (Abdollahi et al. 2020;Ballet et al. 2020) is used to model the -ray sky including the latest interstellar emission model (gll_iem_v07.fits)and the standard template for the isotropic emission (iso_P8R3_SOURCE_V2_v1.txt).We have in the likelihood fitting all 4FGL-DR2 sources lying within the size of the circular region of interest (ROI) plus a radius .The ROI size as well as  depend on the minimum energy ( min ) of the analysis.Since the point spread function of the Fermi-LAT improves considerably at 1 GeV compared to that at the usually adopted  min = 100 MeV, we considered an ROI with radius of 7 • and  = 3 • .Spectral parameters associated with the sources lying outside the ROI are kept fixed, whereas those within the ROI are allowed to vary during the likelihood fitting.For each source, we use the same spectral model as adopted in the 4FGL-DR2 catalog.After the first round of the optimization, the ROI is scanned to search for unmodeled -ray sources by generating test statistic (TS) maps.The maximum likelihood TS is defined as TS = 2log( 1 - 0 ), where  0 and  1 denote the likelihood values of the null hypothesis (i.e., no source at the position of interest) and that of the alternative hypothesis (i.e., the existence of a point source), respectively (Mattox et al. 1996).When an excess emission (TS> 25) is found, it is added to the sky model following a power-law spectral model and, iteratively, a second set of TS maps is generated.A final likelihood fit is performed to optimize all the free parameters in the ROI once all excesses above the background are identified and inserted to the model to optimize the sky model as best as possible.This procedure is necessary since our data are taken in longer exposure than those from which the 4FGL-DR2 sky maps were constructed.
Next, we freeze the spectral parameters of all sources to their bestfitted values except for the source of interest.Then, we apply the EBL attenuation to the spectral model of the source of interest.For this purpose, we update the spectral model of the blazar either to EBLAtten::PowerLaw2 or EBLAtten::LogParabola depending on whether the 4FGL-DR2 spectrum of the blazar reveals a significant curvature or not, which is given by the Spectrum_Type column (see Abdollahi et al. 2020;Ballet et al. 2020, for details).In particular, the EBL attenuation is inserted in the spectral PowerLaw2 or LogParabola models 1 as exp[−(, )], where  is the EBL optical depth as a function of the energy of the -ray photon and the redshift, .In the first step, we use  = 0, i.e., no EBL, and perform the likelihood fitting in the energy range 1−10 GeV.This exercise enables us to determine the intrinsic shape of the -ray spectrum, including any possible curvature.Then, we extrapolate the fitted spectrum to 2 TeV and allow the EBL optical depth to vary in the likelihood fit by changing the redshift value from  = 0.01 to  = 3 in 100 linear steps.In every step, the log-likelihood value is calculated thus permitting us to generate a likelihood profile as a function of redshift for each blazar in the sample.
Later, we create a different TS profile by comparing this likelihood as a function of redshift with the likelihood from the null hypothesis of no EBL.Therefore, this TS is related to the EBL strength.The peak of the TS profile (TS  ) corresponds to the maximum-likelihood value of the redshift of the target blazar, which, due to previous considerations on the possible intrinsic curvature of the spectra, should be interpreted as an upper limit.We remove from our results (1) all sources with TS  < 1 since these correspond to noisy profiles and 2) those sources whose likelihood profiles do not show a peak in the 0.01−3 redshift range.We determine confidence limits from the 1 https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html .An example of one of the TS profiles.This is the case of 4FGL J0802.3-0942 with TS  = 38.3.Note that negative TS, but not negative TS  , are possible at the highest redshifts and they mean that the alternate hypothesis of having a point source with -ray spectrum attenuated by EBL at a given redshift is disfavoured w.r.t. to that of having no EBL.This may be an indication that the EBL model overestimates the attenuation at these redshifts.
TS profile.Figure 1 shows an example of one likelihood profile, in particular, for 4FGL J0802.3-0942.
For estimating the impact on the redshifts due to the uncertainties of the intrinsic spectrum, we implement a Markov Chain Monte Carlo simulation for a subset of test blazars.For each blazar, we generate 100 random values from a Gaussian distribution centered on the bestfit intrinsic slope with a standard deviation equal to its uncertainty.We then rerun the LAT spectral analysis with these values fixed as the intrinsic spectrum.The resultant redshift distributions are consistent with our original results within the statistical uncertainties.For example, for J0540.5+5823, the redshift of TS  obtained is =1.0±0.3, while from our simulation, we find a redshift distribution whose mean is 1.1 and standard deviation 0.4.Similar consistency is observed for the other blazars we tested, indicating that the impact of the statistical uncertainties on the redshift limits is minimal.These results hold at a 68% confidence level, reinforcing the robustness of our original findings.

Redshift upper limits
There are 344/331 non-variable/variable sources that satisfy the criteria that are given at the end of Section 3.2.From this sample, there are 5/40 non-variable/variable FSRQs, all of them, with redshifts in the 4LAC catalogue (Ajello et al. 2020), except 4FGL J1532.7-1319.However, we note that this source is classified as BCU in the SIMBAD database2 .
Then, we have 338/293 non-variable/variable sources (BL Lacs and BCUs), from which 166/162 have a redshift in the 4LAC catalogue and 172/130 do not have a redshift.As a sanity check we plot in Figure 2, the redshift obtained from our EBL-attenuation methodology versus the one in 4LAC for the 166 sources with redshift and TS  ≥ 1, and also for the 91 sources with redshift and TS  ≥ 4. Figure 3 shows similar information for the variable blazars, 162 with TS  ≥ 1 and 111 with TS  ≥ 4. We can see empirically for both sub-samples that the redshift obtained from our methodology can be considered an upper limit for the redshift of the sources in both cases.This is expected since there can be a component in the spectral curvature of the sources that is not produced by EBL attenuation but which is intrinsic to their emission.We tested whether our results change when they are extracted from different periods for the variable sources, concluding that the redshift upper limit generally varies significantly.However in the vast majority of cases, from the comparison with known redshifts, the result is still an upper limit.Furthermore, the empirical test shown in Figure 3 also supports the validity of our upper limits.
For the non-variable blazars, and the case with TS  ≥ 1 (and also TS  ≥ 4), there is only one source that is away from the one to one line by more than 1, which are also within the statistical expectations.This source is Ton 116 ( 4  = 1.065 vs.    < 0.632 at 84% C.L., TS  = 5.69).The redshift of Ton 116 given by SIMBAD comes from Abazajian et al. (2009).We could not see any feature leading to that redshift and we think that it could have been automatically assigned by the line search pipeline, probably with low significance.Yet there is a lower limit by Paiano et al. (2017c) at  > 0.48 that is compatible with our upper limit of    < 0.632 (84% C.L.).
In the case of the variable blazars there are eight blazars that are away from the one to one line more than 1 (although all within 2), also in agreement with statistical expectations.These sources are: • 3C 66A ( 4  = 0.444 vs.    < 0.150 at 84% C.L., TS  = 1.66).The  4  = 0.444 value comes from very weak lines reported in Miller et al. (1978); Kinney et al. (1991), and Lanzetta et al. (1993).Based on the marginal detection of the host galaxy by Wurtz et al. (1996), Finke et al. (2008) estimated the redshift to be  = 0.321.However, Stadnik & Romani (2014) were not able to detect the host galaxy and, based on this non-detection, gave a lower limit of  > 0.42.In SIMBAD, the reported value for the redshift of this source is      = 0.340, using the HI column density (Neeleman et al. 2016).Based on a comparison of the LAT and MAGIC spectrum, and assuming EBL attenuation, Aleksić et al. (2011) estimated the redshift  < 0.68.Recently, Paiano et al. (2017a) observed the optical spectra of 3C 66A in a spectroscopic campaign carried out at the 10 m Gran Telescopio Canarias, leading to inconclusive results.To this date, the redshift of this source remains a subject of debate.
• 1ES 0502+675 ( 4  = 0.416 vs.    < 0.331 at 84% C.L., TS  = 54.39).As for the previous case, the redshift reported in SIMBAD is different from the value that appears in the 4LAC.In 4LAC, the redshift value  4  = 0.416 is taken from Landt et al. (2002) , while in SIMBAD the value      = 0.314 is selected from Sbarufatti et al. (2005), that is compatible with the upper limit here derived.However, these studies rely on imaging methods, since there is no optical spectrum available for this source.
• TXS 0628−240 ( 4  = 1.238 vs.    < 1.124 at 84% C.L., TS  = 45.70).The redshift comes from Landt (2012), which actually gives a lower limit based on a clear Mg II absorption line.This  > 1.238 limit is similar to our upper limit of    < 1.124 (84% C.L.), likely indicating that the Mg II line actually comes from the BL Lac host galaxy instead of an intercepting gas cloud.
• NVSS J090226+205045 ( 4  = 2.055 vs.    < 1.274 at 84% C.L., TS  = 5.44).The 4LAC redshift is taken from the automatic line finding pipeline given by the SDSS DR6 (Richards et al. 2009), which can be biased for faint lines as are typical in BL Lacs.However, there are no emission/absorption lines in the spectra that allow its confirmation.A more recent upper limit for the redshift of this object is available in the literature:  < 1.21 given by Rau et al. (2012) using UV-to-NIR photometry.
• PKS 0903-57 ( 4  = 0.695 vs.    < 0.451 at 84% C.L., TS  = 7.82).This value is selected from optical spectroscopic observations performed by Thompson et al. (1990).However, these authors select the position of the source by its radio emission, and they mention that the optical counterpart corresponding to the radio emission is a star, so the optical spectra they took corresponds to a source 4 arcsec west of the position of the radio emission, a Seyfert I. Therefore, it would be necessary to re-observe the optical spectra for this object, to see whether the optical counterpart observed to determine the value  = 0.695 is correct or is a different source.
• 87GB 105148.6+222705( 4  = 2.055 vs.    < 1.274 at 84% C.L., TS  = 6.71).As in the previous case, the 4LAC redshift value is taken from the SDSS DR6 (Richards et al. 2009), while this value is assigned automatically with poor significance.Rau et al. (2012) gave an upper limit of  < 1.36, and there is a more recent value given automatically by the SDSS DR133 for this object:  = 0.63, both compatible with the upper limit presented here.
• MG3 J225517+2409 ( 4  = 1.370 vs.    < 0.752 at 84% C.L., TS  = 2.77).The 4LAC redshift value is from Albareti et al. ( 2017), but we did not identify any feature leading to that redshift and we think that it could have been automatically assigned, probably with small significance.
Since the upper limit interpretation is robust for non-variable as well as for variable sources, we do not distinguish between them in the following discussions.Table 1 lists the 303 blazars that do not have a redshift in 4LAC and whose likelihood profile peaks at  < 3. We stress that the 125 sources with 1 ≤ TS  < 4 correspond to a significance lower than 2, which according to Figure 2 seem to be compatible with upper limits as well.At any rate, although we give these low significance results (TS  < 4) in Table 1, we suggest the user be prudent about using them.
Figure 4 shows the 303 sources listed in Table 1 in comparison with our initial source selection.It is clear that these blazars lie in the lower region of the initial sample cloud, where, as expected, they tend to have a lower photon index, i.e., hard spectrum.
In Figure 5, we show the 303 blazars using the redshift upper limit at 84% C.L. along the cosmic -ray horizon (CGRH), this is the energy as a function of redshift for which  = 1.This implies CGRH divides the universe in two different regions for which the surviving flux is larger or smaller than exp(−1) ∼37%.The upper limits that lie around the most opaque regions, such as at  ∼ 2.5 and energies of approximately 200 GeV are poorly constrained.
In order to estimate the uncertainties introduced by the EBL, we also run our pipeline using the lower uncertainties of the optical depths given by Domínguez et al. (2023) based on the model by Saldana-Lopez et al. (2021).These lower uncertainties will translate into larger values of the redshift where the TS profile peaks, which goes from smaller shifts in TS  for the lower redshifts and larger for the higher redshifts.This effect is a consequence of having better constrained the EBL at the lower than at the higher redshifts.The increase in the redshifts where TS  is, in redshift bins of 0.5 from  = 0 to  = 3, is on average approximately 13%, 16%, 21%, 28%, 26% and 30%, respectively.We stress that the redshifts upper limits obtained for other models such as Finke et al. (2010), Domínguez et al. (2011), Franceschini & Rodighiero (2017) and Finke et al. (2022) are within these uncertainties.Any larger intensity EBL will not affect the redshifts upper limits because a more opaque universe leads to lower redshifts for the TS  .Furthermore, we tested that TS  tend to be larger for Domínguez et al. (2023) than for Domínguez et al. (2011) indicating that the former model produces better fits, as expected since it is based on deeper and extended multiwavelength data.

Application case: detectability predictions for the Cherenkov Telescope Array
The redshift upper limits are used to compute the expected detection significance of the 303 sources listed in The CTA-North and CTA-South IRFs are used for the sources with positive and negative declination, respectively.We assume observations around culmination time to select the most adequate IRF configuration concerning the zenith angle.We have generated pointlike differential sensitivity curves for a total exposure time of 20 hours by scaling the IRFs corresponding to an exposure of 5 hours utilizing gammapy (Deil et al. 2017;Acero et al. 2022).As products from this calculation we obtained the differential flux, number of excess events, and number of background events that would generate a 5 signal in each energy bin.This is then used to obtain the expected number of excess events for each considered source, which is computed by scaling linearly the number of excess events necessary to get a detection significance of 5 with the ratio of the differential flux.This is the differential flux of the sensitivity curve for each energy bin.The obtained number of excess and background events are used to estimate the statistical detection significance of the blazars using the function WStatCountsStatistic of gammapy, which applies the formula derived by Li & Ma (1983).The energy threshold in each case is selected as the energy of the lower bin of the corresponding differential sensitivity curve, i.e. the energy threshold is considered larger for sources that can only be observed with a zenith angle > 50 deg and also for sources observed with CTA-South.
For the source's spectral shape, we consider the result from the LAT data analysis at the redshift upper limit, which we stress is a PowerLaw or LogParabola including the EBL exponential attenuation.This spectral shape is extrapolated to TeV energies considering the EBL absorption following the optical depths by Domínguez et al. (2023) based on the recent model by Saldana-Lopez et al. (2021), consistent with our LAT data reduction.The significance is estimated using as the redshift of the source, the upper limit values at 84% C.L. from Table 1; this is the redshift of the TS  plus upper uncertainties.For these calculations we use the best zenith angle for the observation either from CTA-North or CTA-South.This results in a total of 21 sources with an expected significance of > 5.We stress that this refers to detections in the average/quiescent state of the source derived from the average flux over 14 years of Fermi-LAT telescope time.Note that considering only the redshift of the TS  , which can be a reasonable approach according to Figure 2, there are 41 sources with a significance of > 5.These additional 20 sources are J0143.7-5846,J0322.0+2335,J0333.9+6537,J0334.2-3725,J0338.5+1302,J0338.9-2848,J0540.5+5823,J0600.3+1244,J0709.2-1527,J0812.0+0237,J0826.4-6404,J1240.4-7148,J1253.2+5301,J1454.4+5124,J1518.0-2731,J1553.5-3118,J1610.7-6648,J1704.5-0527,J2104.3-0212 and J2247.8+4413.We stress that, since we use redshift upper limits in our computation, these CTA detection significances can be considered lower limits.If the sources were closer, the EBL attenuation would be less, so the very-high-energy emission ( > 50 GeV) from the sources would be brighter, and they would be easier for CTA to detect.

SUMMARY AND CONCLUSIONS
We have calculated redshift upper limits for 303 -ray blazars using LAT data and EBL attenuation.The majority of these blazars are classified as BL Lacs.The redshifts derived with our methodology can be useful for studying the evolution of blazars or planning observations with other telescopes, for instance, with imaging atmospheric Cherenkov telescopes such as MAGIC, VERITAS, H.E.S.S and CTA.We find that 21 blazars, for which we have found a redshift upper limit, can be detected in their average/quiescent state with CTA in 20h.
Figure 1.An example of one of the TS profiles.This is the case of 4FGL J0802.3-0942 with TS  = 38.3.Note that negative TS, but not negative TS  , are possible at the highest redshifts and they mean that the alternate hypothesis of having a point source with -ray spectrum attenuated by EBL at a given redshift is disfavoured w.r.t. to that of having no EBL.This may be an indication that the EBL model overestimates the attenuation at these redshifts.

Figure 2 .Figure 3 .
Figure2.Redshift from our EBL methodology for non-variable blazars (Left panel) for the 166 sources with redshift in 4LAC and TS  ≥ 1 and (Right panel) the 91 sources with redshift in 4LAC and TS  ≥ 4. The uncertainties at lower redshifts,  < 0.5, can be larger than the ones at higher redshifts,  ∼ 1, because at the lower redshifts, the EBL affects energies larger than 200-300 GeV where the LAT's effective area is smaller.

Figure 5 .
Figure 5. Highest energy photons of 4LAC blazars (grey circles) as a function of redshift, the cosmic -ray horizon from Domínguez et al. (2023, red band) based on the EBL model by Saldana-Lopez et al. (2021), and the upper limits at 84% C.L. from this work (blue circles).We find upper limits for the redshifts of 303 blazars.