Constraints on Extragalactic Background Light using TeV Observations of BL Lacertae objects

The extragalactic background light (EBL) is the diffuse photon field formed mainly by thermal processes resulting from structure formation. It exists today at wavelengths from the cosmic optical background (COB) at $\lambda_{\rm EBL}=0.1-10\mu m$ to the cosmic infrared background (CIB) at $\lambda_{\rm EBL}=10-1000\mu m$. Thus, EBL’s strength and spectral shape hold fundamental information about the history of star formation and cosmological evolution (Inoue et al., 2013; Saldana-Lopez et al., 2021). In addition, the EBL photons can attenuate the distant VHE (E $\geq$ 100GeV) $\gamma$-rays via the $\gamma\gamma\rightarrow e^{+}e^{-}$ process that forms an insurmountable obstacle to the high-energy cosmic astronomy (Partridge & Peebles, 1967; Peebles, 1993; Hauser et al., 1998; Dwek & Krennrich, 2013). Therefore, the knowledge of EBL intensity and spectra is necessary and relatively important but remains controversial.

Directly measuring the EBL faces several challenges, and the most difficult task is removing the pollution of the foreground photons like zodiacal light from the Milky Way and the diffuse galactic light (Hauser et al., 1998; Hauser & Dwek, 2001; Mattila, 2006). These attempts are engrossed in the optical (Bernstein, 2007; Lauer et al., 2021) and the IR(Bernstein, 2007), while the reliability of data sets needs to be in general agreement. In addition, calibration issues also exist in practice. A strict lower limit on EBL intensity can be acquired by collecting the light emitted from resolved galaxies (Madau & Pozzetti, 2000; Hauser & Dwek, 2001). This approach, available in many bands, including UV, optical/NIR far IR, and submillimeter, favours the local EBL determination. The acquired EBL intensity in optical and near-IR are generally below direct photometry estimates. Basic descriptions of EBL and its measurement are elaborated in the context of recent studies by Saldana-Lopez et al. (2021).

An indirect way to constrain the EBL is to obtain the EBL optical depth as the functions of wavelength and redshift if both the intrinsic and observational $\gamma$-rays are available from distant extragalactic sources. Blazars are the most common VHE gamma-ray sources, classified into flat-spectrum quasars (FSRQ) and BL Lacs based on the width of the emission line. Besides, blazars exist at various redshifts, and their VHE spectrum is modified by EBL absorption, becoming the natural candidate for studying EBL. This approach has been tested by several authors (Stecker & de Jager, 1993; Dwek & Slavin, 1994; Mazin & Raue, 2007; Finke & Razzaque, 2009; Yang, & Wang, 2010; Domínguez et al., 2013, 2019; Fermi-LAT Collaboration et al., 2018; Desai et al., 2019; Longhua et al., 2020; Singh & Meintjes, 2020; Singh et al., 2021), and has been remarkably successful in practice. However, this method still faces the uncertainty of the intrinsic spectral shape and optical depth model in the fitting process. In fact, the VHE $\gamma$-rays of blazars can be generated through one-zone leptonic synchrotron self-Compton model (SSC) (Zhang et al., 2012), hadronic process (Zheng et al., 2016) or electromagnetic cascades (Chen et al., 2015). The intrinsic/unattenuated energy spectra of blazars are usually assumed to constrain the EBL (Domínguez et al., 2019; Fermi-LAT Collaboration et al., 2018; Desai et al., 2019).

Unlike previous studies focusing on each distant object (e.g., Acciari et al. (2019)), we only consider three main types of TeV BL Lacs. To get lower statistics and the FIR limits due to a pile-up at energies, we use two near and well-measured BL Lacs, Mkn421 and Mkn501. Besides, we take one distant source, 1ES0229-200, and two sources with the hard spectra, 1ES1101-232 and 1ES1218+304, because their TeV emissions are more subject to EBL absorption, and their hard spectra are natural candidates to study the lower limit on EBL. We obtain the VHE data compiled by (Biteau & Williams, 2015) and derive the updated data in an online database ¹¹1https://gamma-cat.readthedocs.io/. Note that we exclude the $Fermi$-LAT data in this work since the origin of photons in the GeV bands may not be the same as in the TeV band, and fitting intrinsic spectra also requires the (quasi)-simultaneous observations.

For the BL Lacs, the intrinsic spectra in the TeV regime are characterized by $dN/dE\propto E^{-\Gamma_{i}}$. One can expect the same photon index $\Gamma_{i}=1.5$ as the synchrotron emission spectrum if the inverse-Compton(IC) occurs in the Thomson regime and $\Gamma_{i}\geq 1.5$ if the scattering happens in Klein-Nishima regime. Furthermore, if the VHE comes from the hadronic process, then $\Gamma_{i}=2.0$ (Mazin & Raue, 2007). However, Katarzyński et al. (2006) argue that if synchrotron emission and IC do not occur in the region of electron acceleration, the electron spectrum may become truncated due to propagation, such as $\Gamma_{i}<1.0$. Besides, at least some distant TeV BL Lacs, e.g., 1ES 1101-232 and 1ES 0229+200, exhibit a much harder spectral index. Krennrich et al. (2008) point out that their spectral index $\Gamma_{i}$ is 1.28 or harder even under the weakest EBL model. The high density of EBL photons will lead to more absorption, requiring the harder intrinsic spectrum to get the measured spectrum that is not very soft. Considering the above mentioned, we can set $\Gamma_{i}$ as 1.5 and 1.0 in the fitting procedure to determine the EBL density’s lower limits.

In this paper, a flat $\Lambda$CDM cosmology with $\Omega_{\Lambda}=0.7$, $\Omega_{M}=0.3$, $H_{0}=70\rm{km\ s^{-1}\ Mpc^{-1}}$ is assumed. Besides, the unprimed and primed quantities represent the quantities in the observer and the co-moving frames respectively.

The distant $\gamma$-ray photon of energy $E_{\gamma}$ at $z$ will interact with a background photon of energy $\epsilon$ by the cross section $\sigma_{\gamma\gamma}\left(E_{\gamma},\epsilon,\mu\right)$:

where $\sigma_{\rm T}$ is the Thompson cross section, $\epsilon_{th}(E_{\gamma},\ \mu)$ is the threshold energy dictated by kinematics, $\mu$ are the cosine of the angle for the incident photons and the mass of electron, respectively. The corresponding optical depth is given by

where $dl/dz^{\prime}$ is the differential distant of redshift written as

$\epsilon_{\mathrm{th}}^{\prime}=\epsilon_{\mathrm{th}}(E_{\mathrm{\gamma}}^{\prime},\mu)$, $E_{\mathrm{\gamma}}^{\prime}=E_{\mathrm{\gamma}}(1+z^{\prime})$, and $n_{\rm{EBL}}$ is the EBL density in the observer frame given by

Following Dwek & Krennrich (2005), we set the EBL intensity in $\rm{nWm^{-2}sr^{-1}}$ as the 9 degrees of polynomial:

The form of the intrinsic spectrum $(\frac{dN}{dE})_{\rm int}$ used in this paper is taken from very general considerations about the emission mechanisms. The SEDs of blazars are commonly thought to be originated from the leptonic model, e.g., the SSC (Band & Grindlay, 1985). In addition, the VHE spectrum may result from hadronic scenarios like emission from leptonic models demonstrating a smooth and concave spectrum. Furthermore, the spectrum may have an addition of cut-off components caused by the Klein-Nishina effect or a cut-off in the underlying electron distribution. Considering the above reasons, the intrinsic spectrum, following H. E. S. S. Collaboration et al. (2013), can be expressed as a power law (PL), $N_{0}(E/E_{0})^{-\Gamma}$, a power law with the exponential cut-off (PLE), $N_{0}(E/E_{0})^{-\Gamma}{\rm exp}(-E/E_{cut})$, a log parabola (LP), $N_{0}(E/E_{0})^{{\rm-a-bln}(E/E_{0})}$, and a log parabola with the exponential cut-off (LPE), $N_{0}(E/E_{0})^{{\rm-a-bln}(E/E_{0})}{\rm exp}(-E/E_{cut})$. In the above formulas, $\Gamma$ and $N_{0}$ are the photon index and flux normalization, $a$ is the photon index at $E_{0}$, and $b\geq 0$ is the curvature parameter. Here, the reference energy $E_{0}$ is set to the decorrelation energy of the spectrum. In this work, following Biteau & Williams (2015), $E_{0}$ is fixed to the central value of the energy range of each spectrum and obtained by $\sqrt{E_{\rm{min,VHE}}E_{\rm{max,VHE}}}$.

The distant VHE $\gamma$-rays arriving at Earth can be attenuated by the EBL photons, and the observed spectrum is given by

where $\alpha$ is the scale factor determined by spectral fitting.

In this paper, four widely accepted EBL templates of Franceschini et al. (2008)(Fr), Finke & Razzaque (2009)(Fi), Domínguez et al. (2011)(Do) and Gilmore et al. (2012)(Gi) are used to fit each source. The first template by Franceschini et al. (2008) is based on data sets available from $Spitzer$ IR observations and ground-based telescopes under an empirical backward evolution model to estimate EBL contribution. The distant observables at near IR bands are obtained from the contribution of the evolution of spheroidal, spiral, and merger galaxies. The rest of the EBL bands are from extrapolation (for shorter wavelengths) or the thermal dust emission from the interstellar galaxy medium (for longer IR wavelengths). The second template proposed by Finke & Razzaque (2009) predicts the intensity of EBL directly from the stellar radiation and re-emission by dust in the stellar medium. In this template, the emission in the shorter wavelength is obtained from the main and off-main sequence stars under the assumption of the star formation rate, initial mass function, and dust extinction. At the larger wavelength, the emission mainly originates from the dust, similar to a combination of three black bodies. The third template by Domínguez et al. (2011) is mainly based on the observations of 6000 galaxies. The EBL’s evolving spectrum comes from the observed rest-frame $K$-band galaxy luminosity. The final template by Gilmore et al. (2012) is based on the semi-analytical approach describing the evolving radiation released by galaxies and quasars starting from initial cosmological conditions.

We perform the spectral fitting with a minimum chi-square given by the LMFIT package (Newville et al., 2016) under MCMC method (Lewis & Bridle, 2002; Mackay, 2003; Foreman-Mackey et al., 2013). LMFIT is a Least-Squares Minimization routine and provides a simple and flexible approach to explore the high-level parameters’ space. We input the $\chi^{2}$, which comes from $\chi^{2}=\sum_{i=1}^{N}(\frac{f_{obs}-f_{model}}{\sigma_{i}})^{2}$, where $f_{obs}$,$f_{model}$ are observed and modelled spectra of TeV BL Lacs, $\sigma_{i}$ is the observed error. We thereby derive the intrinsic spectrum of a chosen source by adopting the function to obtain the lowest value of $\chi^{2}(\alpha)$ when $\alpha=0$. Based on the adopted model, we then fit $\alpha$ that should satisfy $\chi^{2}(1)\geq\chi^{2}(\alpha)$. We derive the final optical depth and uncertainty by calculating the mean of four scaled individual optical depths. Upon those procedures, we employ the MCMC algorithm (Lewis & Bridle, 2002; Mackay, 2003) to do the EBL density fitting based on Equations from Eq.2 to Eq.5.

To get the EBL shapes at the redshift of z=0, we use the TeV data to fit the observed SEDs of 5 BL Lacs based on the assumed intrinsic spectra. At first, following the formalism described in section 2, we derive the best fitting SEDs, the scale factors, and their uncertainties. We present the results in Figure 1 - 4 and the Table 1. Then we determine the nine parameters of the EBL model by fitting average scaled optical depth. We provide the best fitted spectra of the EBL and the 68% confident level (C.L.) results in Figure 5 - 6 and the Table 2. Finally, as the comparison, we also give the upper limits of the EBL density derived by direct measurements (Dwek & Krennrich, 2013), $spitzer$ observations in mid-IR bands (Fazio et al., 2004), galaxy counts (Madau & Pozzetti, 2000), and the lower limits (LW08) at 3.6$\mu m$ band (Levenson & Wright, 2008).

In this paper, we study the EBL property by fitting the observed spectra of five TeV BL Lacs, including close-by and well-studied sources, distant sources, and hard spectral sources. We thoroughly test four shapes of intrinsic gamma-ray spectra to fit the observed SEDs using the commonly used EBL templates. Unlikely the stacking analysis (Fermi-LAT Collaboration et al., 2018; Desai et al., 2019) that requires a large amount of data, we get the average scale factor $\alpha$ of EBL templates and its uncertainty by finding the lowest value of $\chi^{2}(\alpha)$ during the fitting procedures. From Table1, we can see that the scale factors are from 0.8 to 4.0, although their uncertainties are around 50%, which are similar to those obtained by (Sinha et al., 2014). Finally, we get the EBL SED from Eq.2 to Eq.5 based on the average scaled EBL templates. It is stressed that the evolution of the EBL with redshift must be considered because both the universe expansion and the evolution of radiation sources affect the EBL intensity. Therefore, to compare with the observed EBL SED, we convert the EBL SED at the redshift z to one at z=0.

We stress the following points: At first, we carefully select the samples that contain the well-observed gamma-ray data, the moderate redshifts, and the hard $\gamma$ spectra. Commonly, a source with many data can minimize the uncertainty of the derived EBL, and the appropriate distant source can balance EBL’s extinction study and the re-emission of $\gamma$-rays during the propagation process (Zheng et al., 2016). Secondly, the sources with hard gamma-ray spectra could reveal the lower limit of the EBL. Under simple assumptions, we can obtain the lowest limits of the EBL SED and test the radiation mechanism of TeV BL Lacs. Finally, the choice of $\Gamma$ as 1.0 and 1.5 is somewhat arbitrary. However, the previous studies revealed that the hardest index ($\sim 1.28$) is required to fit the gamma-ray spectra for some TeV BL Lacs (Krennrich et al., 2008) and the second $Fermi$-LAT flare catalogue (2FAV) could contain the sources with large uncertainties (Abdollahi et al., 2017).

We derive the final optical depth and uncertainty by calculating the mean of four scaled individual optical depths. Upon these procedures, we employ MCMC algorithm (Lewis & Bridle, 2002; Mackay, 2003) to fit the EBL SED based on Equations from Eq.2 to Eq.5. We can see that the maximum EBL SED is roughly consistent with the published measurement, around 50 $\rm{nWm^{-2}sr^{-1}}$. We also find that the ”bump-like” spectral feature in our derived EBL SED around 1 eV band is similar to ones obtained by (Korochkin et al., 2020), demonstrating that the excess features in EBL flux could have physical significance. Furthermore, inspired by theory and observation, we use the fixed spectral index to detect the lower limit of the EBL SED that is most likely compatible with direct observations and galaxy count. Interestingly, we know that two spectral windows mostly favour finding a faint extragalactic background; one is near the infrared band at 3.5 $\mu m$. We find that this faint value is also higher than the derived EBL density by 1ES 1101-232 under the case of $\Gamma_{i}=1$, implying this TeV BL Lacs may have much harder spectra. We summarize the results as following:

(1) Four TeV BL Lacs of MKn421, MKn501, 1ES1101-232, and 1ES0229+200 can well constrain the upper limits of EBL SED. However, for 1ES0229+200 with lower $\alpha$, the EBL SED is higher than that obtained by another source with large scale factor, such as 1ES1101-232, because the fitting optical depths have large uncertain (the error to be $\geq$ 1.6) and the EBL templates are only given by Finke & Razzaque (2009);

(2) The upper limits of EBL density is around 50 $\rm{nWm^{-2}sr^{-1}}$. at 1 eV, and the lower density could be at ones between given by $spitzer$ (Fazio et al., 2004) and galaxy counts (Madau & Pozzetti, 2000);

(3) 1ES1218+304 favours the scenario that the EBL is more transparent to VHE $\gamma$-rays. However, 1ES1101-232 with a similar redshift does not favour this collusion;

(4) Using the fixed spectral index to limit the EBL, for the 1ES1101-232, we doubt that its spectral index could be harder than that we considered before and below 1.0.

From above, the behavior of EBL is still complex, and the EBL SEDs constrained by TeV objects vary with redshifts and directions from our results. Mirizzi & Montanino (2009) singled out that axion-like particles (ALPs) could affect the EBL distribution in different sky directions. However, ALPs work at the redshift z$>$0.2 (Mirizzi & Montanino, 2009). The EBL inhomogeneity could exist at only the order of 1%, which is much smaller than the results in this paper. The proper reason is that the evolution of the EBL with redshift must be fully considered. In fact, the evolution of the EBL number density $n_{\rm EBL}$ is scaled redshift as $(1+z)^{3-k}$ rather than the empirical $(1+z)^{3}$, and Singh et al. (2021) point out that $k$ relates to redshift. Another mechanism to explain our results is that the origin of TeV $\gamma$-rays is still unclear, the intrinsic spectrum used here does not fully cover, and the VHE could be much harder than that in previous studies when the spectral index is fixed at 1.0. Besides, the recent VHE observations of distant gamma-ray sources show that the EBL is more transparent to $\gamma$-rays than we thought before.

The authors would like to thank the anonymous referee for the helpful comments, which improved the manuscript significantly. The authors gratefully acknowledge the financial supports from the National Natural Science Foundation of China (grants 12063005, 12063006 and U2031111), the Special Basic Cooperative Research Programs of Yunnan Provincial Undergraduate Universities (grants 2019FH001-012, 2019FH001-076 and 202001BA070001-031), the Program for Innovative Research Team (in Science and Technology) in University of Yunnan Province (IRTSTYN), the Science Research Foundation of Yunnan Education Department of China (grants 2020J0649), the program for Reserve Talents of Young and Middle-aged Academic and Technical Leaders in Yunnan Province (grants 202205AC160087) and the National Natural Science Foundation of Yunnan Province (grants 202101AU070010). The authors (QLH) gratefully acknowledge the financial supports from the Hundred Talents Program of Yuxi (grants 2019).

The data underlying this article are derived from gammacat: https://gamma-cat.readthedocs.io, and will be shared on reasonable request to the corresponding author.