The First VHE Activity of OJ 287 and the Extragalactic Background Light

Abstract

The BL Lacertae (BL Lac) object OJ 287 underwent an intense X-ray activity phase, exhibiting its brightest recorded X-ray flare in 2016–2017, characterized by much softer X-ray spectra and, concurrently, its first-ever recorded very-high-energy (VHE) emission (100–560 GeV), reported by the VERITAS observatory. Broadband spectral energy distribution reveals a new jet emission component similar to high-synchrotron-peaked BL Lac objects, thereby implying the soft X-ray spectrum for the synchrotron emission. Using the advantage of simultaneous X-ray and VHE spectral information, as well as the source being a low-synchrotron-peaked BL Lac object, we systematically explored the extragalactic background light (EBL) spectrum by demanding that the VHE spectrum cannot be harder than the X-ray spectrum. We used three different phenomenological forms of the EBL spectral shape (power-law, parabola, and polynomial) motivated by current constraints on the EBL with the Bayesian Monte Carlo approach to infer the credible EBL range. Our study favors an almost flat power-law spectral shape and is consistent with previous studies. The other spectral forms capable of capturing curvature though result in a better statistics value; the improvement is statistically insignificant given the additional parameters.

1. Introduction

The Universe as apparent today has gone through many different phases of evolution since beginning from a very hot and dense phase—the Big Bang. Each of these different evolutionary phases is expected to be associated with a distinct radiation signature, e.g., the cosmic microwave background (CMB) is the relic associated with the decoupling of radiation–matter in the matter-dominated phase ($\mathrm{z}$∼1100). One of the milestones in this evolutionary sequence is gravity locally becoming dominant, marking the start of structure formation, i.e., the first stars with copious optical and ultraviolet emission. A significant part of this radiation is expected to be re-processed to infra-red energies through interactions with matter/dust [1]. The extragalactic background light (EBL) is this cumulative radiation in the wavelength range of 0.1–1000 $\mathsf{\mu }$m [2] produced through this cosmic history, from both resolved and unresolved sources. The inferred spectral shape of the EBL is bimodal with a peak at around 1 $\mathsf{\mu }$m, attributed to be from direct starlight emission (but see [3]) and the other peak at around 100 $\mathsf{\mu }$m argued to be from re-emission of absorbed starlight by the dust that encodes the formation and evolution history of galaxies and the structure formation of the Universe.

Direct measurement of the EBL is a very challenging task, primarily due to contamination by the Zodiacal and galactic light [3,4]. However, lower limits on the EBL intensity have been provided by the integrated light from resolved galaxies [5,6,7,8], and upper limits have been derived by direct observations of the EBL with space- and ground-based instruments sensitive in specific bands [9,10,11]. Most of the EBL is produced by discrete galactic and primordial stellar sources. Studying the fluctuations in their number and clustering properties has provided the spatial fluctuations in the EBL intensity that provide a limit on its intensity through fluctuation analysis [12]. These limits indicate a bimodal structure of the EBL between the UV and far-IR with a peak at ∼1 $\mathsf{\mu }$m and the other at ∼200 $\mathsf{\mu }$m, but the absolute spectral energy distribution (SED) of the EBL remains uncertain and highly debatable. Several theoretical models predicting the SED of the EBL at $z=0$ have been developed (see Figure 1; Franceschini et al. [13], Gilmore et al. [14], Finke et al. [15], Kneiske and Dole [16]), but the lack of observations between the two peaks and the observational challenges have impeded its measurement. Backward-evolution models start with the determination of the luminosity density, $\mathcal{L}(\lambda ,z=0)$, in the local Universe and then evolve it with redshift at different wavelengths [13,17]. Forward-evolution models use the cosmic star formation rate (CSFR) to calculate $\mathcal{L}(\lambda ,z=0)$ determining the distribution of the energy over wavelengths [15]. Semi-analytic models calculate $\mathcal{L}(\lambda ,z=0)$ by taking physical processes into account for the formation and evolution of the structure of the Universe [14].

A complimentary and powerful probe to put limits on the EBL is by analyzing the signatures of the EBL on the very-high-energy (VHE; E > 100 GeV) gamma-ray spectrum of distant blazars [18,19,20], provided the intrinsic VHE spectrum of the source is known. The VHE gamma-rays while traversing the intergalactic medium interact with the EBL photons and get absorbed by the process of electron–positron pair production [21]. The absorption of VHE photons by the EBL steepens the VHE spectra, hence providing the imprints of the EBL (e.g., [22,23]). Dis-entangling the EBL-induced attenuation signature in the observed spectrum and the intrinsic spectrum offers a powerful tracer of EBL density. Thanks to the tremendous overall advancement in observational gamma-ray astronomy, as well as blazar’s SED modeling, several studies have attempted to constrain the EBL based on the attenuation of gamma-rays (e.g., [4,24,25,26]). For example, using the H 1426 + 428 spectrum, along with the spectra of already discovered TeV blazars, limits on the EBL were derived [27,28]. Similarly, HESS observations of 1ES 1101-232 and H 2356-309 were used to estimate the EBL in optical–NIR wavelengths by assuming the intrinsic spectrum cannot be harder than ${\mathrm{\Gamma }}_{int}$ = 1.5 [29]. However, subsequent EBL studies argued for a much harder particle spectrum of distant blazars [30], and now, many sources have shown a very hard photon spectrum, implying a harder particle spectrum (${\mathrm{\Gamma }}_{int}$∼1), and magnetic re-connection has been argued to be the most prominent driver. Contrary to this, studies have also demonstrated that, in certain scenarios, the hardness can be well reproduced, e.g., internal gamma-ray absorption through interaction with narrow-band strong photon fields [31,32] without requiring a hard particle spectrum.

The first systematic exploration of the EBL signature in blazar spectra at VHEs was presented in H.E.S.S. Collaboration et al. [20], reporting the EBL signature at the $8.8\sigma$ level with respect to the model of Franceschini et al. [13] and constraining it in $0.3\mathsf{\mu }\mathrm{m}<\lambda <17\mathsf{\mu }\mathrm{m}$. This approach was similar to the one introduced in the Fermi large area telescope (Fermi-LAT) collaboration EBL work [19]. The Fermi-LAT team subsequently expanded this work, exploiting improved exposure data of 739 blazars and found attenuation due to the EBL in the entire explored redshift range $0.03<z<3.1$ for all twelve models that were tested [33]. Using this, they also reconstructed EBL evolution, as well as the cosmic star formation history with the latter being consistent with other independent studies. It should, however, be noted that, for sources exhibiting strong spectral changes, using long-exposure data may introduce false EBL signatures.

Blazars currently have been divided into three primary spectral classes: low- synchrotron-peaked (LSP), intermediate-synchrotron-peaked (ISP), and high-synchrotron-peaked (HSP), based on the band/frequency in which the low-energy component of their bi-model SED attains the maximum [34]. This scheme is an extension of the previous existing classification designating BL Lacertae (BL Lac) objects as low-peaked (LBLs), intermediate-peaked (IBLs), and high-peaked (HBLs) with LSP now representing FSRQs and LBLs. The majority of the blazars that have exhibited frequent VHE activities belong to the HBL/HSP class and are strong X-ray emitters. The X-ray (or low-energy component of the SED) is understood to be synchrotron, while high-energy is primarily inverse Comptonization (IC) (e.g., [35,36]). A simultaneous X-ray and VHE spectral measurement can be used to study the EBL, except that the HBL VHE can be significantly steepened due to the Klein–Nishina effect, thereby overestimating the EBL [37]. Such systematic exploration, however, is expected to not affect the outcome when using a sample of blazars at different redshifts. LBL’s gamma-rays, on the other hand, are primarily due to the IC of infra-red photons and are relatively less prone to the IC effect, except that only a few counted ones are known1, and thus, at best, this can provide a constraint on the limited range. Nonetheless, this can be used to confirm and cross- check the current EBL limits.

OJ 287 is a BL Lac-type object located at z = 0.306 [38] and is one of the three BL Lac sources detected at VHEs [39]. BL Lacs are jetted active galactic nuclei (AGN) with a very weak or a complete absence of emission line features in their optical spectrum [equivalent width $<5$A] and are generally accepted to lack the standard accretion disk and infra-red torus. They emit an entirely jet-dominated highly variable continuum covering the entire accessible electromagnetic spectrum from radio to GeV–TeV gamma-rays. The broadband emission exhibits a characteristic broad dual-humped spectral energy distribution [34]. In the broadband SED-based categorization, OJ 287 is classified as a low-energy-peaked BL lac (LBL/LSP) object, and its gamma-ray emission is primarily attributed to the inverse Comptonization of the IR photon field [40]. The inverse Compton interpretation of MeV–GeV gamma-rays is also strongly favored by the SED modeling of candidate neutrino blazars, indicating a sub-dominant hadronic contribution (e.g., [36]).

OJ 287 exhibited its first-ever very-high-energy (VHE; E > 100 GeV) $\gamma$-ray activity in 2016–2017, reported by the Very Energetic Radiation Imaging Telescope Array System (VERITAS), which continued for a very long time [39]. The broadband study of the VHE period along with previous and later activity episodes revealed an additional HBL-like emission component, establishing the soft X-ray emission to synchrotron [41], as is the case in HBLs. This VHE emission can be explained via the HBL/HSP scenario or the IC of an IR photon field. Given that gamma-ray emission of LBLs/LSPs is due to the IC-IR [40,42], we believe that the VHE is also likely of the same origin.

Figure 1. EBL SED from direct and indirect measurement. Data points are from direct measurements [26], and the curves are the theoretical models2 predicting limits on the EBL.

Figure 1. EBL SED from direct and indirect measurement. Data points are from direct measurements [26], and the curves are the theoretical models2 predicting limits on the EBL.

In this work, we have carried out a systematic investigation of the reported VHE spectrum of OJ 287 to explore the EBL imprints it carries and, thereby, infer the EBL. In Section 2, we revisit the basics of the photon–photon pair production mechanism and the relation between the intrinsic and observed spectrum. In Section 3, we describe the constraints from the observations on the VHE spectral shape imposed by the general physical principles of the blazar’s radiative processes. The methodology/formalism to constrain the free parameters and the selection criteria of the models is explained in Section 4. The results are discussed in Section 5. A flat $\mathrm{\Lambda }$ cold dark matter ($\mathrm{\Lambda }$CDM) cosmology is used in this work with Hubble constant $H_0=70\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$; matter density parameter ${\mathrm{\Omega }}_m$ = 0.27; and curvature density parameter ${\mathrm{\Omega }}_{\mathrm{\Lambda }}$ = 0.73.

2. Attenuation of Gamma-Ray Photons by Pair Production

The VHE $\gamma$-rays traversing the Universe are attenuated via pair production through interaction with low-energy EBL photons, ${\gamma }_{VHE}+{\gamma }_{EBL}\to e^{-}+e^{+}$, and the following condition must be satisfied by the energies of two photons for the creation of the electron–positron pair:

$$\begin{array}{c}\hfill E_{\gamma }{ϵ}_{th}(1-\mu )=2m_e^2{c}^4\end{array}$$

(1)

where $E_{\gamma }$ and ${ϵ}_{th}$ are the energies of VHE $\gamma$-rays and the threshold energy for EBL photons, respectively, and $\mu =cos\varphi$, such that $\varphi$ is the angle between the momenta of two photons and $m_e$ is the rest mass of an electron.

The cross-section for the $\gamma -\gamma$ interaction [21] is given by

$$\begin{array}{c}\hfill {\sigma }_{\gamma \gamma }(E_{\gamma },ϵ,\mu ,z)=\frac{3{\sigma }_T}{16}(1-{\beta }^2)\left[2\beta ({\beta }^2-2)+(3-{\beta }^4)ln\left(\frac{1+\beta }{1-\beta }\right)\right]\end{array}$$

(2)

where ${\sigma }_T$ is the Thompson cross-section and $\beta$ is given by

$$\begin{array}{c}\hfill \beta ={\left(1-\frac{2m_e^2{c}^4}{E_{\gamma }ϵ(1-\mu )}\right)}^{1/2}=\sqrt{1-\frac{{ϵ}_{th}}{ϵ}}\end{array}$$

(3)

The optical depth for pair creation of a $\gamma$-ray photon emitted by a source at a redshift z via interacting with the EBL field [2] is given by

$$\begin{array}{c}\hfill {\tau }_{\gamma \gamma }(E_{\gamma },z)={\int }_0^{z}d{z}^{\prime }\frac{d\mathit{l}}{d{z}^{\prime }}{\int }_{-1}^{1}d\mu \frac{1-\mu }{2}{\int }_{{ϵ}_{th}^{\prime }}^{\infty }dϵ{\mathit{n}}_{ϵ}(ϵ,z^{\prime }){(1+z^{\prime })}^3{\sigma }_{\gamma \gamma }({\beta }^{\prime },z^{\prime })\end{array}$$

(4)

where l is the proper distance such that

$$\begin{array}{c}\hfill \frac{d\mathit{l}}{dz}=c\frac{d\mathit{t}}{dz}=\frac{c}{H_0(1+z){[{\mathrm{\Omega }}_m{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}]}^{1/2}}\end{array}$$

(5)

assuming a flat $({\mathrm{\Omega }}_K=0)$ and matter-dominated $({\mathrm{\Omega }}_R<<1)$$\mathrm{\Lambda }$CDM Universe. ${\mathit{n}}_{ϵ}(ϵ,z)=\frac{d\mathit{n}(ϵ,z)}{dϵ}$ is the specific comoving number density of the EBL. ${(1+z)}^3$ is the conversion term from the specific to the proper number density. In this context, the threshold condition leading to pair production is modified to

$$\begin{array}{c}\hfill {ϵ}_{th}^{\prime }=\frac{2{\left(m_e{c}^2\right)}^2}{E_{\gamma }ϵ(1-\mu )(1+z)}\end{array}$$

(6)

with $\beta$ becoming ${\beta }^{\prime }=\sqrt{1-\frac{{ϵ}_{th}^{\prime }}{ϵ}}$.

Calculating the pair opacity of $\gamma$-rays requires the knowledge of the EBL comoving specific photon number density ${\mathit{n}}_{ϵ}(ϵ,z)$ as a function of redshift. The conversion of specific comoving intensity to comoving specific photon number density is given by

$$\begin{array}{c}\hfill {ϵ}^2{\mathit{n}}_{ϵ}(ϵ,z)=\frac{4\pi }{c}\nu I_{\nu }(\nu ,z)=2.62\times {10}^{-4}\nu I_{\nu }(\nu ,z)\end{array}$$

(7)

where $ϵ$ is the energy of EBL photons and the coefficient in the second term is calculated for $ϵ$ in eV, ${\mathit{n}}_{ϵ}$ in ${\mathrm{cm}}^{-3}{\mathrm{eV}}^{-1}$, and $\nu I_{\nu }$ in $\mathrm{nW}{\mathrm{m}}^{-2}{\mathrm{sr}}^{-1}$.

The pair opacity leads to the steepening of the observed VHE spectrum from the intrinsic VHE $\gamma$-ray spectrum ${\left(\frac{d{N}_{\gamma }}{dE}\right)}_{int}$ emitted from the source due to EBL absorption, and the observed spectrum ${\left(\frac{d{N}_{\gamma }}{dE}\right)}_{obs}$ is related to intrinsic one as:

$$\begin{array}{c}\hfill {\left(\frac{d{N}_{\gamma }}{dE}\right)}_{obs}={\left(\frac{d{N}_{\gamma }}{dE}\right)}_{int}{e}^{-\tau (E_{\gamma },z)}\end{array}$$

(8)

Thus, in this approach, knowledge of the intrinsic spectrum is required to estimate the EBL. It should be noted that the kinematics of the process and the energy dependency of the pair cross-section (${\sigma }_{\gamma \gamma }$; Equation (2)) can change the observed VHE spectrum from a simple power-law to more complex with steepening followed by a hardening depending on $E_{\gamma }$, and thus, a correct understanding of the EBL is essential to explore the VHE and TeV emission mechanisms (e.g., [43]).

3. Observational Constraints

VERITAS detected an almost persistent, marginally variable VHE $\gamma$-ray emission from OJ 287 over a relatively extended duration from December 2016–March 2017 at ∼$10\sigma$ above the background in the energy range 100 GeV–560 GeV [39]. The total time-averaged VHE spectrum is consistent with a power-law photon spectral index (${\mathrm{\Gamma }}_{VHE}$) of 3.49 ± 0.28. The total time-averaged flux above 150 GeV for the entire dataset is $(4.61\pm 0.62)\times {10}^{-12}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$, which corresponds to ∼$(1.3\pm 0.2)\%$ of the flux of the Crab Nebula above the same threshold [44].

The reported VHE activity followed in the declining phase, after the brightest-ever-recorded X-ray flare of the source [41,45]. Detailed studies and investigations of the spectral and temporal behavior establish it as a jet emission component as well, but of a different spectral class (HBL-like) [41,46] than the traditional spectral class of the source [34,40]. Such an appearance of a new spectral class is very rare in blazars.

As stated, for OJ 287, an external photon field at a temperature of $250K$ as soft target photons up-scattered by the highly relativistic electrons present in the jet by the IC mechanism has been argued to explain the high energy part of the spectrum [40] in the LBL spectral state of the source while the 2015 gamma-ray spectral change was explained by the IC-BLR [47,48]. The photons emitted from the broad-line region (BLR) and the infrared (IR) emission from the dusty torus have typical photon energies ${ϵ}_{BLR}\simeq 10\mathrm{eV}$ and ${ϵ}_{IR}\simeq 0.3\mathrm{eV}$ (∼$1200\mathrm{K}$), respectively. For the IC scattering of soft photons ${ϵ}_{soft}$, the Lorentz factor of ultra-relativistic electrons that can produce the VHE emission is ${\gamma }_{e,IC}\simeq {[(1+z){ϵ}_{VHE,obs}/\left({\delta }_{source}{ϵ}_{soft}\right)]}^{1/2}$, where ${\delta }_{source}$ is the Doppler factor and ${ϵ}_{VHE,obs}$ is the typical VHE photon energy taken to be 100 GeV for our purposes. The scattering occurs in the Thomson regime if ${\gamma }_{e,IC}{ϵ}_{soft}/\left(m_e{c}^2\right)<1$, which modifies to

$$\begin{array}{c}\hfill {ϵ}_{soft}<\frac{{\delta }_{source}{\left(m_e{c}^2\right)}^2}{(1+z){ϵ}_{VHE,obs}}.\end{array}$$

(9)

If the emission region is located within the external field extent, the energy density of the external field in the comoving frame is boosted by a factor of ${\mathrm{\Gamma }}_{source}^2$ and the photon energy increases by a factor of ${\mathrm{\Gamma }}_{source}$, where ${\mathrm{\Gamma }}_{source}$ is the bulk Lorentz factor. Hence, the Thomson limit is ${ϵ}_{ER}<1.9\mathrm{eV}\times ({\delta }_{source}/{\mathrm{\Gamma }}_{source})$. For our study, we have taken bulk Lorentz factor ${\mathrm{\Gamma }}_{source}=16.5\pm 4.0$ and Doppler factor ${\delta }_{source}=18.9\pm 6.4$ computed by Jorstad et al. [49] from 17 epochs at a 7 mm wavelength with the Very Long Baseline Array (VLBA) from 1998–2001. As a result, only IR photons are scattered in the Thomson regime.

The external field up-scattered to gamma-rays via the IC can, in turn, also interact with the emitted VHE photons. The VHE gamma-ray photons upon interaction with the ER field are destroyed in the photon–photon pair production process when the condition in Equation (1) is satisfied. The BLR radiation can absorb VHE photons of observed energy:

$$\begin{array}{c}\hfill {ϵ}_{VHE,BLR,obs}\gtrsim \frac{{\left(m_e{c}^2\right)}^2}{(1+z){ϵ}_{BLR}}\simeq 20{\left(\frac{{ϵ}_{BLR}}{10.2eV}\right)}^{-1}\mathrm{GeV}\end{array}$$

(10)

and for IR radiation, the observed threshold is

$$\begin{array}{c}\hfill {ϵ}_{VHE,IR,obs}\gtrsim \frac{{\left(m_e{c}^2\right)}^2}{(1+z){ϵ}_{IR}}\simeq 665{\left(\frac{{ϵ}_{IR}}{0.3eV}\right)}^{-1}\mathrm{GeV}\end{array}$$

(11)

which is beyond the VHE energies detected by VERITAS. The observation of the VHE instead indicates that the pair opacity is not high within the source. A plot of the BLR photon number density for different opacity values is shown in Figure 2 (right) for a 100 GeV photon.

The above gamma-ray opacity arguments imply that the VHE emission is solely due to the IC scattering of IR photons from the dusty torus. This is also consistent with the report of the association of the VHE activity with a new radio feature at parsec scales (∼10 pc, [50])—the putative location of the IR torus in the standard radio-loud AGN unification scheme. Detailed and systematic investigations of broadband SEDs during the VHE phase compared to the well-known OJ 287 SED (LBL/LSP; e.g., [46]) show a hardening of the optical–UV spectrum and a much softer X-ray, implying a peak at UV energies (e.g., [39,41,46,48,51]). A similar hardening was observed at MeV–GeV energies, accompanied by the reported softer VHE spectrum, implying an additional peak at GeV energies. A comparison with the well-known blazar spectral sequence [34] establishes that such spectral features are characteristics of HBL/HSP emission as stated above and argued in previous works as well (e.g., [41,51]).

Since the IC-IR occurs in the Thomson regime, in the standard blazar emission paradigm, the intrinsic VHE spectrum should be the same as the corresponding synchrotron spectrum in the case of both originating from the same particle population, i.e., from one emission region, as is the case here with the X-ray and VHE associated with the new HBL-like component. Thus, the VHE spectral index (${\mathrm{\Gamma }}_{VHE}$) cannot be harder than the corresponding synchrotron counterpart, i.e., the X-ray spectral index here (${\mathrm{\Gamma }}_{X-ray}$). A fit of a power-law $\frac{dN}{dE}=N_0{\left(\frac{E}{E_0}\right)}^{-\mathrm{\Gamma }}$ to Swift-XRT data (Figure 2-left) results in a satisfactory fit (${\chi }^2/d.o.f=3.15/5$) with flux normalization $N_0=(3.32\pm 1.02)\times {10}^{-6}{\mathrm{MeV}}^{-1}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$ and spectral index ${\mathrm{\Gamma }}_{X-ray}=2.66\pm 0.05$ at normalization energy $E_0=0.2\mathrm{TeV}$, providing the maximum hardness the intrinsic VHE spectrum can attain, ${\mathrm{\Gamma }}_{VHE}^{max}\simeq {\mathrm{\Gamma }}_{X-ray}=2.66$. It should be noted that the majority of the X-ray spectra during this VHE activity duration and such extremely soft X-ray state have been shown to be log-normal [41,46,52], rather than a simple power-law. However, this log-normal part appears at the high-energy end of the X-ray band and is due to the traditional LBL spectral state of the source [51]. The NuSTAR observation during an intermediate spectral state also supports a simple power-law spectrum up to 30 keV. Beyond this, the spectrum is background-dominated [51]. Thus, within the current context and practical limitations, the extremely soft X-ray spectrum associated with the VHE activity is consistent with a power-law.

The range of TeV detection limits the range over which the EBL can be explored/constrained. This can be seen from the threshold condition (Equation (1) or (6)). For a given VHE energy range, the corresponding wavelength range that can be probed is given by

$$\begin{array}{c}\hfill {\lambda }_{ebl,max}=\frac{h}{2m_{e}c}\times \frac{E_{VHE}}{m_e{c}^2}\times {(1+z)}^2\sim 95\mathsf{\mu }\mathrm{m}\times \left(\frac{E_{VHE}}{20\mathrm{TeV}}\right)\times {(1+z)}^2.\end{array}$$

(12)

For OJ 287, the observed VHE range of 100 GeV to 560 GeV corresponds to an EBL wavelength range of 0.8–2 $\mathsf{\mu }$m.

4. Analysis Method

With the above restriction on the VHE intrinsic spectrum, we performed the maximum likelihood analysis using the Markov chain Monte Carlo (MCMC) algorithm to constrain the EBL flux density. The MCMC code emcee [53], a Python implementation of an affine invariant MCMC ensemble sampler [54], was used to obtain the best fit of the parameters and their posterior probabilities based on the observed VHE data. This approach is based on Bayesian statistics with a set of parameters (hereafter, $\theta$) for the data (hereafter, $\mathcal{D}$). The emcee sampling algorithm generates multiple MCMC chains (‘walkers’) that evolve through the parameter space simultaneously at each step. Assuming that the observed data points are independent of each other and each one follows a Gaussian normal distribution, the likelihood function can be expressed as:

$$\begin{array}{c}\hfill \mathcal{L}\left(\theta \right|\mathcal{D})=\prod _{i=1}^m\frac{exp\left(-{\left({\left(\frac{dN}{dE}\right)}_{i,obs}-{\left(\frac{dN}{dE}\right)}_i\right)}^2/2{\sigma }_{i,obs}^2\right)}{\sqrt{2\pi {\sigma }_{i,obs}^2}}\end{array}$$

(13)

where ${\left(\frac{dN}{dE}\right)}_i$ is the model flux, m is the number of data points, and ${\left(\frac{dN}{dE}\right)}_{i,obs}$ and ${\sigma }_{i,obs}$ are the observed flux and the variance, respectively.

Following Baye’s theorem, the posterior probability is given by:

$$\begin{array}{c}\hfill P\left(\theta \right|\mathcal{D})\propto \mathcal{L}(\mathcal{D}\left|\theta \right)P\left(\theta \right).\end{array}$$

(14)

Due to the simplicity of the models, we used uninformative (flat) priors ($P\left(\theta \right)$) for all parameters with appropriately large limits.

To determine the best EBL shape that describes the near-infrared (NIR) range, three plausible shapes motivated by the current understanding/constraints of EBL flux density have been explored thoroughly and are listed in

Table 1. Functions describing the different EBL forms.

Name	Abbreviation	Function	Parameters to Evaluate
Power Law	PL	${\varphi }_0{\left(E\right)}^{-\mathrm{\Gamma }}$	${\varphi }_0$, $\mathrm{\Gamma }$
Parabola	PB	$-4h{(\lambda -{\lambda }_0)}^2+k$	$h,k$
Polynomial	PM	$a{\lambda }^3+b{\lambda }^2+c\lambda +d$	$a,b,c,d$

. Since the EBL range implied by the OJ 287 VHE spectrum is very small, we adopted the power-law (PL) shape [55]. However, the first EBL peak (Figure 1) resides within the constraint range, and thus, complex models parabola (PB, with axis ${\lambda }_0=1.4\mathsf{\mu }$m) and polynomial (PM) were considered to study the possible curvature information present.

The fitting procedure consists of the following steps: (i) One of the three functional forms was chosen to model the EBL. (ii) We initialized the code with 64 walkers exploring the parameter space in 10,000 steps each, with a burn-in of 2000 steps for each walker. For chains to converge we sampled for a large number of steps with different walkers initialized differently and then ran with different random number seeds to create posterior probability inferences that are substantially the same. (iii) Finally, the posterior probability and the likelihood profile were computed as a function of the model parameters ($\theta$) to estimate the statistical uncertainty associated with the maximum likelihood estimates (MLEs).

To select the EBL model most compatible with the observational data, we used various numerical procedures to obtain the goodness-of-fit and to determine the corresponding best value of the fitting parameters. In general, the least chi-squared (${\chi }^2$) method is sufficient for comparing models, and the model with the least ${\chi }^2$ value implies that it has a more acceptable fit with the data than the others. However, a reduced chi-squared (${\chi }_{red}^2$) is used extensively in testing the goodness-of-fit, so we have considered ${\chi }_{red}^2$ for a better interpretation of the result, which is defined as follows:

$$\begin{array}{c}\hfill {\chi }_{red}^2=\frac{{\chi }_{min}^2}{N-K},\end{array}$$

(15)

where N is the number of data points, K is the number of free parameters for the model, and for Gaussian errors, ${\chi }_{min}^2=-2ln\left({\left(\mathcal{L}\left(\theta \right|\mathcal{D})\right)}_{max}\right)$. When the variance of the measurement error is ${\chi }_{red}^2\approx 1$, this implies that the measurement error is in accord with the error variance. Having the value of ${\chi }_{red}^2\gg 1$ indicates a poor model fit. A ${\chi }_{red}^2\ll 1$ indicates that the model has over-fit the data are or the variance of the measurement error has been overestimated. Since the degree of freedom (N-K) is small, ${\chi }^2$ efficiency reduces greatly when comparing more complex models; to make a fair comparison, we considered two advanced information criteria: the Akaike Information Criterion (AIC) by Akaike [56] and the Bayesian Information Criterion (BIC) by Schwarz [57], defined as:

$$\begin{array}{c}\hfill AIC={\chi }_{min}^2+2K\end{array}$$

(16)

$$\begin{array}{c}\hfill BIC={\chi }_{min}^2+Kln\left(N\right).\end{array}$$

(17)

Now, after considering a reference model (denoted by R), for any model M, the difference $∆$AIC$=|{\mathrm{AIC}}_M-{\mathrm{AIC}}_R|$ will get us the following conclusions. As explained in Burnham and Anderson [58], (i) if $∆$AIC$\le 2$, then there is significant support for the $M^{th}$ model, and $M^{th}$ is probably the proper description of the data, (ii) $4\le ∆\mathrm{AIC}\le 7$ implies less support for the $M^{th}$ model over the reference model, and, (iii) $∆$AIC$\ge 10$ means that the $M^{th}$ model essentially has no support, in principle, this model is of no use. For BIC model selection criteria, following the guidelines of Raftery [59], if the difference in the BICs is $0<∆\mathrm{BIC}\le 2$, this implies weak evidence, for $2\le ∆\mathrm{BIC}\le 6$, this implies positive evidence, and $6\le ∆\mathrm{BIC}\le 10$ is strong evidence, whereas more than 10 is very strong evidence against the model with a higher BIC.

5. Results

Figure 3 shows the best-fit EBL curves for different assumed phenomenological models listed in Table 1 along with the $1\sigma$ credible range, extracted from the posterior. The corresponding best-fit parameters are given in

Table 2. The median and $68\%$ credible intervals of free parameters resulting from MCMC fits along with the statistics value.

Model: (Case)	Parameters	${\mathit{\chi }}^2$	${\mathit{\chi }}_{\mathbf{red}.}^2$	AIC	BIC
PL: (i)	$N_{VHE}=Fixed$ $N_{ebl}=5.1_{-2.4}^{+4.0}\times {10}^{-4}$ ${\mathrm{\Gamma }}_{ebl}=3.1_{-1.6}^{+1.3}$	0.434	$0.434/3$	−249.740	−247.938
PB: (i)	$N_{VHE}=Fixed$ $h=9.4_{-4.8}^{+7.2}\times {10}^{-1}$ $k=1.3_{-0.8}^{+1.8}\times 10$	1.204	$1.204/3$	−248.970	−247.168
PM	$N_{VHE}=Fixed$ $a=-0.1_{-0.8}^{+1.0}\times 10$ $b=1.2_{-3.0}^{+3.5}\times 10$ $c=-0.6_{-3.2}^{+3.9}\times 10$ $d=0.3_{-1.1}^{+1.3}$	0.431	$0.431/1$	−247.546	−243.941
PL: (ii)	$N_{VHE}=2.5_{-1.4}^{+6.8}\times {10}^{-10}$ $N_{ebl}=7.9_{-6.3}^{+6.8}\times {10}^{-4}$ ${\mathrm{\Gamma }}_{ebl}=2.1_{-0.8}^{+3.1}$	0.218	$0.218/2$	−248.858	−246.154
PB: (ii)	$N_{VHE}=3.1_{-1.5}^{+4.1}\times {10}^{-10}$ $h=-2.1_{-7.4}^{+8.7}$ $k=2.2_{-1.4}^{+1.7}\times 10$	0.208	$0.208/2$	−248.867	−246.163

.

The best-fit model parameters are derived following the methodology outlined in the previous section. We have three models: power-law (PL), parabola (PB), and a polynomial (PM). The respective functional form is given in Table 1. Of these, the PM model has the maximum number of parameters (four), while the PL and PB have two each. Given this and noting that we have only five VHE data points, there are possibilities in terms of fixing/freezing a few of the parameters in different models based on our current understanding. The simplest is the normalization factor. The analysis is performed by taking two cases into account: (i) Considering that the $\gamma$-ray absorption at $E_{VHE}$∼100 GeV is negligible, hence fixing the normalization of the intrinsic VHE spectrum to the first VHE observed data at $E_{VHE}=0.1187\mathrm{TeV}$, and (ii) keeping the VHE flux normalization factor as a fitting parameter along with the EBL model parameters.

In Table 2 and the best-fit plot (Figure 3), the $1\sigma$ range is the median, the $68\%$ credible intervals (analogous to confidence intervals in frequentist statistics) of the model parameters resulting from the MCMC fits. The table also summarizes the parameters that were fixed or free during the model fits and the best-fit criteria: ${\chi }^2$, ${\chi }_{red}^2$, AIC, and BIC. The $1\sigma$ confidence bands for all the EBL models are in good agreement with the lower and upper limits from direct measurements. Since the data available are very limited, this causes the confidence intervals to be very wide.

Secondly, to study the best possible shape describing the constraint EBL flux density region, we thoroughly examined three shapes (see Table 1). The PB and PM models have the perk of fitting a wide range of curvature despite that this roughly resembles a power-law-like shape (see Figure 3). Using the AIC and BIC analysis, no significant ($∆\mathrm{AIC},∆\mathrm{BIC}$) difference is observed between the results of the models, and we concluded that all the EBL models fit the constraint region well. For the PL model, the AIC and BIC are smaller than the others, implying that the power-law is the best description for the constrained EBL flux density region. The corresponding corner plot is in the bottom-right panel in Figure 3.

Thirdly, the predicted EBL by PL (i) model fitting suggests the increase in the EBL flux density beyond wavelength $\lambda$∼$1.4\mathsf{\mu }$m, not in agreement with our present knowledge of the EBL. However, it is in between the constraints derived with galaxy counts and the direct measurements (see Figure 3a). For PL (ii) scenario, it predicts roughly flat EBL density in the range $\lambda$∼(0.8–2) $\mathsf{\mu }$m, which is in agreement with the concave-like SED of the EBL at around $\lambda$∼$1.4\mathsf{\mu }$m (see Figure 3d). The $1\sigma$ contour lies in between the constraints derived from the galaxy counts and the direct measurements with a peak amplitude at $1.4\mathsf{\mu }$m of $\nu I_{\nu }=17.700n{\mathrm{Wm}}^{-2}{\mathrm{sr}}^{-1}$. The $1\sigma$ confidence band of the intrinsic VHE spectrum for PL (ii) is in agreement with the constrained intrinsic VHE spectrum as for case (i) (see Figure 4).

6. Discussion and Conclusions

The gamma-ray emission in BL Lacs with the LBL SED is primarily due to IC scattering of IR photons (e.g., [40,42,60]) and thus, for such sources, if detected at VHE, the softening due to the KN is expected to be considerably less compared to the HBLs/HSPs that dominate VHE detected sources. This makes these sources an excellent EBL estimator and can provide a good re-check of the existing models and constraints on the EBL. However, only three such BL Lacs have been detected at VHE to date, limiting such tests over a very limited range of EBL energies. One of these is OJ 287, which has exhibited an almost non-variable VHE activity over a rather extended duration [39,61]. The activity was coincident with the observation of a new jet feature in the parsec-scale jet [50] at ∼10 pc, indicating the new HBL-like emission component being associated with it and supporting the IC-IR origin of gamma-rays. The extensive observation at the optical to X-ray spectra provides excellent simultaneous data to explore the EBL. Using this and the VERITAS VHE spectrum data of OJ 287 in the range from 100 to 560 GeV, we extracted the maximum possible and the best functional form explaining the spectral shape of the EBL in the NIR range. Using the AIC and BIC to compare the models, we concluded that the power-law (PL) model remains the best explanation for the constraint EBL range. Meanwhile, the predicted EBL flux density for all models is in accord with the limits derived by direct measurements (Figure 3). The assumption we made on the intrinsic VHE spectrum is the up-scattering of soft photons (${ϵ}_{IR}$), and its spectral shape is related to the corresponding synchrotron spectrum i.e., X-ray here. The derived limit in the NIR region indicates that the intergalactic space is less transparent to $\gamma$-rays than predicted by the integrated light of resolved galaxies, as has been reported in other studies.

The other BL objects with LBL SED detected at VHEs are OT 081 [62,63] and AP Librae [64]. However, the OT 081 broadband SED associated with the VHE activity is unavailable to date [63], while the AP Librae VHE SED is argued to be due to IC scattering of CMB radiation from the extended jet [65], while MeV–GeV (Fermi-LAT) from the IC of the IR field has been recently observed [60]. The VHE emission seems to be the continuity of the gamma-ray spectrum from the Fermi-LAT [60]. However, the corresponding synchrotron component is missing and, thus, the intrinsic shape of the VHE spectrum, which is essential in the current approach. Furthermore, being a very nearby source, for the energy range of VHE detection and the uncertainty on the VHE data, the gamma–gamma opacity is too small to be significantly reflected in the VHE spectrum.

Approaches exploiting the Bayesian methodology have been used to obtain the EBL intensity in the optical band by using 259 VHE spectra from the Spectral TeV Extragalactic Catalog (STeVECat) from 56 extragalactic sources with $z>0.01$. Being a fully Bayesian framework enabled improvement in the typical resolution of gamma-ray measurements of the EBL to around $15\%$ [66]. Another recent work by Genaro et al. [67] used the Bayesian with Hamiltonian Monte Carlo approach on 65 VHE spectra of 36 AGNs observed by the Imaging Atmospheric Cherenkov Telescopes (IACTs), simultaneously predicting the EBL flux density in the mid-IR region and intrinsic spectral parameters. For constraining the IR part of the EBL, the authors have considered three blackbody emission components at 40 K, 70 K, and 450 K, representing different dust contributions. They found that the observation of the ∼20 TeV VHE from the nearby HBL Mkn 501 during a flare (HEGRA/1997) is a crucial input in constraining the far-IR region.

Estimating the EBL using VHE is still evolving, primarily due to large uncertainties in the VHE data due to the counting statistic and detection of mostly low luminous nearby sources (HBLs) providing input over a very limited energy range of the EBL. However, the indirect constraints on the NIR by studying the EBL imprints on TeV $\gamma$-rays will continue to improve as more sources over a greater range of redshifts are observed. Next-generation TeV $\gamma$-ray telescopes, such as the Cherenkov Telescope Array (CTA) [68,69], with an order of magnitude better sensitivity and much-expanded energy range from ∼20 GeV up to ∼300 TeV, will be able to detect fainter sources at higher redshift, thereby significantly expanding the sample of VHE blazars over both the energy and redshift range. The low-energy part of the gamma-ray spectrum will allow an unhindered view of the unaffected spectrum, while the high-energy part will reflect the EBL signature and, especially, the expected peculiar features like flattening, cutoff, etc., (e.g., [43]) serving as a strong test of the existing EBL models, a tighter constraint on the EBL spectrum, as well as a constraint on the contribution from faint galaxies. The coverage will allow a thorough reconstruction of the EBL spectrum over the entire UV-to-IR from the current existing VHE samples that are constrained primarily to the UV, optical, and mid-IR [67], its evolution, and thereby, indirectly an improved and better reconstruction of cosmic star formation history to higher redshifts than the current one (e.g., like in [33]).