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

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.


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 (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 (Primack et al. 2005).The extragalactic background light (EBL) is this cumulative radiation in the wavelength range of 0.1-1000 μm (Dwek and Krennrich 2013) 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 μm, attributed to be from direct starlight emission (but see (Lauer et al. 2022)) and the other peak at around 100 μ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 (Hauser and Dwek 2001;Lauer et al. 2022).However, lower limits on the EBL intensity have been provided by the integrated light from resolved galaxies (Madau and Pozzetti 2000;Fazio et al. 2004;Frayer et al. 2006;Dole et al. 2006), and upper limits have been derived by direct observations of the EBL with space-and ground-based instruments sensitive in specific bands (Hauser et al. 1998; Kashlinsky et al. 1996;Gispert et al. 2000).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 (Kashlinsky and Odenwald 2000).These limits indicate a bimodal structure of the EBL between the UV and far-IR with a peak at ∼1 μm and the other at ∼200 μ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  = 0 have been developed (see Figure 1; Franceschini et al. (2008); Gilmore et al. (2009); Finke et al. (2010); Kneiske and Dole (2010)), 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, L (,  = 0), in the local Universe and then evolve it with redshift at different wavelengths (Franceschini et al. 2008;Domínguez et al. 2011).Forward-evolution models use the cosmic star formation rate (CSFR) to calculate L (,  = 0) determining the distribution of the energy over wavelengths (Finke et al. 2010).
Semi-analytic models calculate L (,  = 0) by taking physical processes into account for the formation and evolution of the structure of the Universe (Gilmore et al. 2009).
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 (Georganopoulos et al. 2010;Ackermann et al. 2012; H. E. S. S. Collaboration et al. 2013), 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 (Gould and Schréder 1967).The absorption of VHE photons by the EBL steepens the VHE spectra, hence providing the imprints of the EBL (e.g., (Vassiliev 2000;Mankuzhiyil et al. 2010)).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., (Stecker and de Jager 1993;Hauser and Dwek 2001;Dwek and Krennrich 2005;Mazin and Raue 2007)).For example, using the H 1426 + 428 spectrum, along with the spectra of already discovered TeV blazars, limits on the EBL were derived (Costamante et al. 2004;Kneiske et al. 2004). 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 Γ  = 1.5 (Aharonian et al. 2006).However, subsequent EBL studies argued for a much harder particle spectrum of distant blazars (Krennrich et al. 2008), and now, many sources have shown a very hard photon spectrum, implying a harder particle spectrum (Γ  ∼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 (Aharonian et al. 2008;Zacharopoulou et al. 2011) 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. (2013), reporting the EBL signature at the 8.8 level with respect to the model of Franceschini et al. (2008) and constraining it in 0.3 μm <  < 17 μm.This approach was similar to the one introduced in the Fermi large area telescope (Fermi-LAT) collaboration EBL work (Ackermann et al. 2012).
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 <  < 3.1 for all twelve models that were tested (Fermi- LAT Collaboration et al. 2018).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: lowsynchrotron-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 (Abdo et al. 2010).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., (Böttcher et al. 2013;Gao et al. 2019)).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 (Böttcher et al. 2008).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 (Sitko and Junkkarinen 1985) and is one of the three BL Lac sources detected at VHEs (O'Brien 2017).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 (Abdo et al. 2010).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 (Kushwaha et al. 2013).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., (Gao et al. 2019)).
OJ 287 exhibited its first-ever very-high-energy (VHE; E > 100 GeV) -ray activity in 2016-2017, reported by the Very Energetic Radiation Imaging Telescope Array System (VERITAS), which continued for a very long time (O'Brien 2017).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 (Kushwaha et al. 2018), 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 (Kushwaha et al. 2013;Arsioli and Chang 2018), we believe that the VHE is also likely of the same origin.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 Λ cold dark matter (ΛCDM) cosmology is used in this work with Hubble constant  0 = 70 km s −1 Mpc −1 ; matter density parameter Ω  = 0.27; and curvature density parameter Ω Λ = 0.73.

ATTENUATION OF GAMMA-RAY PHOTONS BY PAIR PRODUCTION
The VHE -rays traversing the Universe are attenuated via pair production through interaction with lowenergy EBL photons,     +    →  − +  + , and the following condition must be satisfied by the energies of two photons for the creation of the electron-positron pair: where   and   ℎ are the energies of VHE -rays and the threshold energy for EBL photons, respectively, and  = cos , such that  is the angle between the momenta of two photons and   is the rest mass of an electron.
The cross-section for the  −  interaction (Gould and Schréder 1967) is given by where   is the Thompson cross-section and  is given by The optical depth for pair creation of a -ray photon emitted by a source at a redshift  via interacting with the EBL field (Dwek and Krennrich 2013) is given by where l is the proper distance such that l (5) assuming a flat (Ω  = 0) and matter-dominated is the specific comoving number density of the EBL.(1 + ) 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 .Calculating the pair opacity of -rays requires the knowledge of the EBL comoving specific photon number density n  (, ) as a function of redshift.The conversion of specific comoving intensity to comoving specific photon number density is given by where  is the energy of EBL photons and the coefficient in the second term is calculated for  in eV, n  in cm −3 eV −1 , and   in nW m −2 sr −1 .
The pair opacity leads to the steepening of the observed VHE spectrum from the intrinsic VHE -ray emitted from the source due to EBL absorption, and the observed spectrum is related to intrinsic one as: 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 (  ; Equation ( 2)) can change the observed VHE spectrum from a simple power-law to more complex with steepening followed by a hardening depending on   , and thus, a correct understanding of the EBL is essential to explore the VHE and TeV emission mechanisms (e.g., (Mazin et al. 2013)).

OBSERVATIONAL CONSTRAINTS
VERITAS detected an almost persistent, marginally variable VHE -ray emission from OJ 287 over a relatively extended duration from December 2016-March 2017 at ∼10 above the background in the energy range 100 GeV-560 GeV (O'Brien 2017).The total time-averaged VHE spectrum is consistent with a powerlaw photon spectral index (Γ    ) of 3.49 ± 0.28.The total time-averaged flux above 150 GeV for the entire dataset is (4.61 ± 0.62) × 10 −12 cm −2 s −1 , which corresponds to ∼(1.3 ± 0.2)% of the flux of the Crab Nebula above the same threshold (Hillas et al. 1998).
The reported VHE activity followed in the declining phase, after the brightest-ever-recorded X-ray flare of the source (Kushwaha et al. 2018;Komossa et al. 2020).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) (Kushwaha et al. 2018;Kushwaha 2022) 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 Γ 2   and the photon energy increases by a factor of Γ   , where Γ   is the bulk Lorentz factor.Hence, the Thomson limit is    < 1.9eV × (   /Γ   ).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:

GeV
(10) and for IR radiation, the observed threshold is 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 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 (Γ    ) cannot be harder than the corresponding synchrotron counterpart, i.e., the -ray spectral index here (Γ −  ).A fit of a power-law 2-left) results in a satisfactory fit ( 2 /.. = 3.15/5) with flux normalization  0 = (3.32 ± 1.02) × 10 −6 MeV −1 cm −2 s −1 and spectral index Γ −  = 2.66 ± 0.05 at normalization energy  0 = 0.2 TeV, providing the maximum hardness the intrinsic VHE spectrum can attain, 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 (Kushwaha et al. 2018(Kushwaha et al. , 2021;;Kushwaha 2022), 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 (Singh et al. 2022).The NuS-TAR observation during an intermediate spectral state also supports a simple power-law spectrum up to 30 keV.Beyond this, the spectrum is background-dominated (Singh et al. 2022).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 For OJ 287, the observed VHE range of 100 GeV to 560 GeV corresponds to an EBL wavelength range of 0.8-2 μm.

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 (Foreman-Mackey et al. 2013), a Python implementation of an affine invariant MCMC ensemble sampler (Goodman and Weare 2010), 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, ) for the data (hereafter, 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: where (    )  is the model flux,  is the number of data points, and (    ) , and  , are the observed flux and the variance, respectively.Following Baye's theorem, the posterior probability is given by: Due to the simplicity of the models, we used uninformative (flat) priors (()) 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.Since the EBL range implied by the OJ 287 VHE spectrum is very small, we adopted the power-law (PL) shape (de Jager et al. 1994).However, the first EBL peak (Figure 1) resides within the constraint range, and thus, complex models parabola (PB, with axis  0 = 1.4 μ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 () 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 ( 2 ) method is sufficient for comparing models, and the model with the least  2 value implies that it has a more acceptable fit with the data than the others.However, a reduced chi-squared ( 2   ) is used extensively in testing the goodness-of-fit, so we have considered  2   for a better interpretation of the result, which is defined as follows: where  is the number of data points,  is the number of free parameters for the model, and for Gaussian errors,  2  = −2((L (|D))  ).When the variance of the measurement error is  2   ≈ 1, this implies that the measurement error is in accord with the error variance.Having the value of  2   ≫ 1 indicates a poor model fit.A  2   ≪ 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,  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 (1973) and the Bayesian Information Criterion (BIC) by Schwarz (1978), defined as: Now, after considering a reference model (denoted by ), for any model , the difference ΔAIC= |AIC  − AIC  | will get us the following conclusions.As explained in Burnham and Anderson (2004), (i) if ΔAIC≤ 2, then there is significant support for the   ℎ model, and   ℎ is probably the proper description of the data, (ii) 4 ≤ ΔAIC ≤ 7 implies less support for the   ℎ model over the reference model, and, (iii) ΔAIC≥ 10 means that the   ℎ model essentially has no support, in principle, this model is of no use.For BIC model selection criteria, following the guidelines of Raftery (1995), if the difference in the BICs is 0 < ΔBIC ≤ 2, this implies weak evidence, for 2 ≤ ΔBIC ≤ 6, this implies positive evidence, and 6 ≤ ΔBIC ≤ 10 is strong evidence, whereas more than 10 is very strong evidence against the model with a higher BIC.

RESULTS
Figure 3 shows the best-fit EBL curves for different assumed phenomenological models listed in Table 1 along with the 1 credible range, extracted from the posterior.The corresponding best-fit parameters are given in Table 2.
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 -ray absorption at     ∼ 100 GeV is negligible, hence fixing the normalization of the intrinsic VHE spectrum to the first VHE observed data at     = 0.1187 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 range is the median, the 68% credible intervals (analogous 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 ∼(0.8-2)μm, which is in agreement with the concave-like SED of the EBL at around ∼1.4 μm (see Figure 3d).The 1 contour lies in between the constraints derived from the galaxy counts and the direct measurements with a peak amplitude at 1.4 μm of   = 17.700 Wm −2 sr −1 .
The 1 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).

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., (Kushwaha et al. 2013;Arsioli and Chang 2018;Roychowdhury et al. 2022)) 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 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 (   ), 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 -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 (Mirzoyan 2016;Manganaro et al. 2021) and AP Librae (Hofmann 2010).However, the OT 081 broadband SED associated with the VHE activity is unavailable to date (Manganaro et al. 2021), while the AP Librae VHE SED is argued to be due to IC scattering of CMB radiation from the extended jet (Sanchez et al. 2015), while MeV-GeV (Fermi-LAT) from the IC of the IR field has been recently observed (Roychowdhury et al. 2022).The VHE emission seems to be the continuity of the gamma-ray spectrum from the Fermi-LAT (Roychowdhury et al. 2022).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.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 -rays will continue to improve as more sources over a greater range of redshifts are observed.Next-generation TeV -ray telescopes, such as the Cherenkov Telescope Array (CTA) (Acharya et al. 2013;Cherenkov Telescope Array Consortium et al. 2019), 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., (Mazin et al. 2013)) 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 (Genaro et al. 2024), 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 (Fermi- LAT Collaboration et al. 2018)).

Figure 1 .
Figure 1.EBL SED from direct and indirect measurement.Data points are from direct measurements (Mazin and Raue 2007), and the curves are the theoretical models 2 predicting limits on the EBL.
than the traditional spectral class of the sourceAbdo et al. (2010);Kushwaha et al. (2013).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 250 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(Kushwaha et al. 2013) in the LBL spectral state of the source while the 2015 gamma-ray spectral change was explained by the IC-BLR(Kushwaha et al. 2017;Kushwaha 2020).The photons emitted from the broad-line region (BLR) and the infrared (IR) emission from the dusty torus have typical photon energies   ≃ 10 eV and    ≃ 0.3 eV (∼ 1200 K), respectively.For the IC scattering of soft photons     , the Lorentz factor of ultra-relativistic electrons that can produce the VHE emission is  , ≃ [(1 + )   , /(       )] 1/2 , where    is the Doppler factor and    , is the typical VHE photon energy taken to be 100 GeV for our purposes.The scattering occurs in the Thomson regime if  ,     /(   2 ) < 1, which modifies to

For
our study, we have taken bulk Lorentz factor Γ   = 16.5±4.0and Doppler factor    = 18.9±6.4computed by Jorstad et al. (2005) 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.
with a new radio feature at parsec scales (∼10 pc,Lico et al. (2022))-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.,Kushwaha (2022)) show a hardening of the optical-UV spectrum and a much softer X-ray, implying a peak at UV energies (e.g.,Kushwaha et al. (2018);Singh et al. (2022);Kushwaha (2020Kushwaha ( , 2022)); O'Brien (2017)).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 sequenceAbdo et al. (2010) establishes that such spectral features are characteristics of HBL/HSP emission as stated above and argued in previous works as well (e.g.,Kushwaha et al. (2018);Singh et al. (2022)).

Figure 2 .
Figure 2. Left: Swift-XRT measurements of OJ 287.The line is the best fit of a power-law (Γ = 2.66 ± 0.05).Right: Photon number density () within the BLR corresponding to different optical depths () for 100 GeV photons.

Figure 4 .
Figure 4.The constructed EBL-corrected VHE spectrum (blue dashed, PL (ii)) along with the 1 confidence band (striped red region).The X-ray spectrum is plotted for reference (normalization scaled to match the VHE data)

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 (O'Brien 2017; Mukherjee and VERITAS Collaboration 2017).The activity was coincident with the observation of a new jet feature in the parsec-scale jet (Lico et al. 2022) at ∼10 pc, indicating the new HBL-like emission component being associated with it and supporting the IC-IR origin of gamma-rays.
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  > 0.01.Being a fully Bayesian framework enabled improvement in the typical resolution of gamma-ray measurements of the EBL to around 15%(Gréaux and Biteau 2023).Another recent work by Genaro et al. (2024) 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.

Table 1 .
Functions describing the different EBL forms.