Probing the extragalactic mid-infrared background with HAWC

The extragalactic background light (EBL) contains all the radiation emitted by nuclear and accretion processes in stars and compact objects since the epoch of recombination. Measuring the EBL density directly is challenging, especially in the near- to far-infrared waveband, mainly due to the zodiacal light foreground. Instead, gamma-ray astronomy offers the possibility to indirectly set limits on the EBL by studying the effects of gamma-ray absorption in the very high energy (VHE:$>$100 GeV) spectra of distant blazars. The High Altitude Water Cherenkov gamma ray observatory (HAWC) is one of the few instruments sensitive to gamma rays with energies above 10 TeV. This offers the opportunity to probe the EBL in the near/mid IR region: $\lambda$ = 1 $\mu$m - 100 $\mu$m. In this study, we fit physically motivated emission models to Fermi Large Area Telescope (LAT) GeV data to extrapolate the intrinsic TeV spectra of blazars. We then simulate a large number of absorbed spectra for different randomly generated EBL model shapes and calculate Bayesian credible bands in the EBL intensity space by comparing and testing the agreement between the absorbed spectra and HAWC extragalactic observations of two blazars. The resulting bands are in agreement with current EBL lower and upper limits, showing a downward trend towards higher wavelength values $\lambda>10\, \mu$m also observed in previous measurements.


INTRODUCTION
The extragalactic background light (EBL) comprises all radiation released by nuclear and accretion processes since the epoch of recombination. It consists of all emitted radiation from stars and compact object surroundings, including that absorbed/re-emitted by dust and accumulated over all redshifts. Measuring and constraining this background radiation is crucial to understand star formation processes and galaxy evolution models. Our current knowledge of the EBL is limited, its direct measurements are challenging due to foreground contamination coming, mainly, from the zodiacal light. However, upper and lower limits have been established using various methods, e.g. integrated galaxy counts from optical observations with the Hubble Space Telescope (Gardner et al. 2000;Madau & Pozzetti 2000) and infrared (IR) observations using the Spitzer Space Telescope (Fazio et al. 2004;Papovich et al. 2004) and the Infrared Space Observatory (Elbaz et al. 2002). An extensive discussion on EBL related matters can be found in Elbaz et al. (2002), Orr et al. (2011) and Cooray (2016).
Over the past two decades, VHE gamma-ray observations have been used to help constrain the spectral properties of the EBL, particularly with observations from blazars, a sub-type of Active Galactic Nuclei (AGN) with ultra-relativistic jets oriented close to the observer's line of sight. VHE gamma rays coming from blazars interact via pair production with the EBL photons (Gerasimova et al. 1962;Gould & Schréder 1967) producing imprints in the observed energy spectra of distant sources. These imprints, along with intrinsic spectral assumptions, can be used to derive limits on the EBL near/mid IR range using VHE observations of blazars (e.g. Aharonian et al. 2007;Mazin & Raue 2007;Orr et al. 2011;Biteau & Williams 2015;Abdalla et al. 2017; Abeysekara & Archer 2019; Acciari & Ansoldi 2019). The gamma-ray horizon establishes the energy at which the intensity of radiation is diminished by a factor 1/e with respect to the emitted intensity. This energy value depends on the distance of the source and it needs to be taken into account when selecting suitable candidate sources for this type of study. If the source is too close, the absorption effect will only be measurable at higher energies, where, depending on the detector and the energy of the gamma ray, the sensitivity might be too low. On the other hand, distant sources will be dimmer, precisely due to the EBL absorption, so these are not ideal candidates to work with either (Franceschini et al. 2019). Equation 1 approximately relates the energy of a gamma ray (E γ ) with the wavelength range of the EBL radiation (λ EBL ) involved in the pair production interaction: where z is the redshift of the source emitting the gamma ray. This equation is useful to estimate the EBL probing power when considering a specific source observed with a specific instrument. The High Altitude Water Cherenkov gamma ray observatory (HAWC) is a water Cherenkov detector that has been operational since 2015 (further technical details can be found in section 2), and that has detected extragalactic sources significantly up to 10 TeV (Albert et al. 2021a). This energy is close to the one established by the gamma-ray horizon for sources like the two blazars Markarian 501 and Markarian 421 (from now on referred to as Markarians), putting HAWC in an advantageous position to potentially probe the mid-IR region of the EBL using observations from these sources.
In this study, physically motivated emission models and data from the Fermi Large Area Telescope (Fermi-LAT) are used to construct an intrinsic spectrum for each of the Markarians (see section 3.1). A large number of randomly generated EBL model shapes are used to apply the EBL attenuation effect to these intrinsic spectra to compare with HAWC data (see section 3.2). The comparison is performed using threeML 1 , a software package for likelihood fitting, in a way that each EBL model can be assigned a likelihood value that expresses the agreement between that particular spectral realization and HAWC data (see section 3.3). This method has the advantage of being independent of any particular EBL shape and of assuming physically motivated intrinsic spectral properties for the sources. Finally, weights are applied to each model in accordance to their calculated likelihood value, and credible intervals in the EBL spectral energy distribution (SED) space are derived (section 4).

DATA: HAWC & FERMI-LAT
HAWC is an array of 300 water Cherenkov detector tanks located in Sierra Negra, Mexico, at an altitude of 4100 m above the sea level and covering a total area of 22000 m 2 . Each tank has four photo-multiplier tubes (PMTs) facing upwards that can detect the Cherenkov light in the water from the transit of secondary particles, which are produced by gamma rays and cosmic rays interacting with the atmosphere. HAWC triggers with a rate of 25 kHz and has a duty cycle of >95%. The observatory continuously monitors 2/3 of the sky, detecting gamma rays in the energy range between 100 GeV and several hundred TeV. More information on the HAWC observatory operation, performance and the way air shower event data are reconstructed can be found in Abeysekara & Albert (2017). For this analysis, specific Fermi-LAT and HAWC data from the blazars Markarian 421 (z = 0.031) and Markarian 501 (z = 0.034) were selected (redshift sourced from NED 2 ). These are two extensively studied extragalactic sources and the brightest in the TeV band. Both sources have been significantly detected by HAWC above 300 GeV up to 10 TeV (Albert et al. 2021a) and by Fermi-LAT between 100 MeV and 1 TeV . For the HAWC dataset, 1343 days of data were used, taken between June 2015 and June 2019 . The analysis was performed using maps created with a special algorithm for reconstructing and determining the energy of the primary gamma rays. This algorithm, denominated ground parameter algorithm, is based primarily on the charge density at a fixed optimal distance from the shower axis and it divides the energy range into quarter-decade bins in log10(E), beginning at log(E/TeV)=-0.5 (0.316 TeV) and ending at log(E/TeV)=2.5 (316 TeV). The ground parameter energy estimator was chosen because it is optimal for higher energies, between 10 TeV and 316 TeV , where the instrument could potentially probe the EBL mid-IR region. For the Fermi-LAT dataset, a time period corresponding to the HAWC dataset was selected between 57180 -58640 MJD. This dataset is within the energy range of 100 MeV -316 GeV, where the absorption is, at most, around 5% for the Markarians redshift 3 . The Fermi-LAT analysis is described in more detail in section 3.1.

EBL ANALYSIS
The approach adopted in this study consists of calculating EBL intensity limits by comparing the expected absorbed flux with actual HAWC observations. The method is similar to previous EBL studies: Mazin & Raue (2007), Orr et al. (2011), Biteau & Williams (2015) and, in particular, to Abeysekara & Archer (2019), presented by the VERITAS collaboration. In the latter analysis, a large number of random generated EBL models is used to calculate the corresponding de-absorbed spectra of selected blazars, and then each model is weighted using criteria based on intrinsic assumptions for these sources to finally derive limits to the EBL intensity. In the present analysis, the EBL models are used to simulate the absorption effect which is then applied over the assumed intrinsic spectra. The resulting absorbed spectra are then compared to the real HAWC data by calculating a likelihood value. Finally, each EBL model is weighted according to their corresponding likelihood value.

Intrinsic Spectra: Fermi-LAT & Naima
A possible and reasonable assumption is to consider the intrinsic TeV spectrum of a given source as an extension of a physical emission model that is in agreement with GeV observations. This assumption relies on the fact that gamma rays in the high energy regime (HE: 100 MeV -100 GeV) and relatively low redshifts (z 0.1) are not significantly affected by EBL absorption, and therefore the observed spectrum can be safely considered to be practically the same as the intrinsic one (e.g. Abdo et al. 2010;Orr et al. 2011;Furniss et al. 2013). In the present study, a Fermi-LAT standard fitting analysis was performed to obtain spectral flux points corresponding to Markarian 421 and Markarian 501, using observations from a similar time period as the one observed by HAWC (for details see section 2). The Fermi-LAT analysis was carried out using the instrument response function P8R2 SOURCE V6, and the spectral parameters are estimated by the binned maximum likelihood method using the Fermipy v.0.17.3 package. The analysis was also performed with a zenith cut of 90 o and, as mentioned in section 2, within the energy range of 100 MeV -316 GeV to minimize possible absorption effects while including as much data as possible. The events were extracted within a 10 degree region of interest (RoI) centered on each source position. The background model includes sources from the Preliminary Fermi-LAT 8-year point sources catalog 4 , Galactic diffuse emission glliemv06.fits and the isotropic diffuse emission isoP8R2SOURCEV6v06.txt models.
To avoid the overestimation of the intrinsic flux, each source lowest energy fluxpoint observed by HAWC (Albert et al. 2021a) is used as a guide; observed Mrk 421 and Mrk 501 fluxpoints at an energy of 830 GeV and 750 GeV, respectively, are de-asbsorbed according to Franceschini et al. (2008) model. These de-absorbed fluxpoints are included in the set of flux points obtained from the Fermi-LAT analysis described above, and altogether are used to fit a physically motivated emission model with Naima, a Python software software that calculates the non-thermal emission from a leptonic or hadronic population of particles (Zabalza 2015).
In this case, the dataset is fit with a Synchrotron-Self-Compton (SSC) model, a scenario in which synchrotron radiation is produced by electrons moving at relativistic velocities in randomly oriented magnetic fields and upscattered by inverse Compton into higher energies by the same electron population. The leptonic distribution is chosen to follow an exponentially cut-off powerlaw (ECPL) distribution: where A is a normalization factor, E 0 is a reference energy set at 1 TeV, α is the power law index and E c the cut-off energy. The decision for this particular leptonic 4 https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/ distribution is motivated by other studies performed on the Markarians by HAWC (Albert et al. 2021a) and others (Zhu et al. 2016;. Based on the results reported in Albert et al. (2021a), the magnetic field and emission radius values are set at 0.03 G and 4 ×10 16 cm for Mrk 421 and 0.02 G and 1 ×10 17 cm for Mrk 501, respectively. Table 1 shows the best fit parameters for the electron distributions of each source. The best fit models are then extrapolated into TeV energies (see figure 2) and used as templates for the intrinsic spectra of the sources (see section 3.3).

EBL models
To generate different EBL model shapes, a grid of 12 λ-values across the wavelength of the EBL photons between 1 µm and 100 µm is set. For each λ-value,an intensity (nW m −2 sr −1 ) value is assigned by randomly generating a number between 1 nW m −2 sr −1 and 50 nW m −2 sr −1 , resulting in a flat initial distribution in EBL intensity; the range is chosen to contain the EBL upper and lower limits discussed in section 1. A particular EBL shape is defined by interpolating each intensity point corresponding to each wavelength, using second order splines. To avoid sharp model shapes, intensity at consecutive grid points is not allowed to change by more than a factor of 2. For this reason and, since the interpolation is performed in increasing order in λ, intensities for higher λ values are more conditionally sampled than lower λ values. This bias, is then corrected by weighting the intensities in each grid point to recover the initial flat intensity distribution. Avoiding sharp models that are nonphysical has the caveat of neglecting possible features coming from polycyclic aromatic hydrocarbons, known to be present in the mid-IR range (Lagache et al. 2005). However, for relatively small redshifts, like those considered in this analysis, these sharp features would be smoothed out to a bump (Domínguez et al. 2011b). Figure 1 shows examples of EBL shapes generated in this way. A total of 30000 models are generated at redshift z = 0. The EBL intensity can be then used to compute the optical depth τ (E, z), which quantifies the absorption effect of gamma rays with energy (E) and traveling a given distance associated with a redshift value z. The theoretical framework behind these calculations can be found in e.g. Dwek & Krennrich (2013); Biteau & Williams (2015); Abdalla et al. (2017). In this analysis, the τ values are computed using the ebltable Python package 5 to read in and interpolate tables for EBL density and to calculate the resulting opac-  The calculated opacity for each redshift-energy is then used to simulate the absorption process by applying the attenuation factor to the assumed intrinsic differential flux in the following way: where dN dE obs is the observed differential flux and dN dE int the assumed intrinsic differential flux. Figure 2 shows Mrk 421 Fermi-LAT flux points with its corresponding Naima-SSC fit extrapolated into TeV energies. An example of the absorbed spectrum according to an EBL model along with HAWC data is also shown.

threeML framework, likelihood and limit extraction
The comparison between the expected absorbed model and HAWC data, is performed adopting a Bayesian approach and using the Multi-Mission Maximum Likelihood framework . This analysis pipeline is capable of handling data from a wide variety of astrophysical detectors. In this particular study the HAWC plugin (HAL 6 ) is used. The HAL fitting technique is based on a forward-folding method that assumes a spectral model shape for the source. In this case, a SSC source emission model is built using the astromodels package, a useful framework to define models for likelihood or Bayesian analysis of astrophysical data 7 . This emission model is customized for each source by plugging in the fit parameters obtained from the Naima fit described in section 3.1, to create an intrinsic spectrum that serves as an input for the threeML fitting pipeline. The EBL absorption is applied to the emission model while their spectral parameters are kept fixed. Finally, the forward-folding method, including detector response effects, is performed to fit the absorbed spectrum of the source using a maximum likelihood technique  to calculate a likelihood value for the fit (L). The process is repeated for each of the 30000 EBL models generated as explained in section 3.2. In this way, instead of optimizing the parameters of the source to find the ones that give the maximum likelihood, these are kept fixed, and only the EBL models are evaluated by their ability to reproduce the data. After the test, each EBL model is assigned a likelihood value that quantifies the agreement between that particular emission + absorption model and the actual data. The maximum likelihood value corresponds to the EBL model whose simulated absorption best reproduces HAWC observations. Starting off from a prior flat intensity distribution for each wavelength, resulting from the uniform sampling described in section 3.2, each model intensity is then assigned a weight as follows: where L i is the likelihood value corresponding to the "ith" EBL model. The intensities corresponding to the maximum-likelihood model are then assigned a weight of 1 and the rest of the model intensities are weighted by a factor of L i /L max , disfavoring EBL shapes that differ from the EBL model that best agrees with data. This approach is based on the idea of relative likelihood, which is an estimate of the probability of a model to reproduce data, normalized by the maximum likelihood value (Kalbfleisch 1985). When combining results from both sources, each EBL model is weighted by multiplying the individual source weights previously calculated: W mrk421 ×W mrk501 . For each λ-value, the maximumlikelihood EBL models from each source are used to calculate an average intensity value for the combined result. A projection of the EBL intensity probability for each value of λ can be obtained from the histogram of weighted intensities, which represents the posterior distribution of intensities at that λ given the observed 7 https://astromodels.readthedocs.io/en/latest/ Figure 3. Weighted intensity distributions for λ = 3.5µm and λ = 12.3µm. The red hatching corresponds to the 68% containment region. Black dashed line represents the intensity corresponding to the maximum likelihood model at that particular value of λ. data. Figure 3 shows an example of these distributions for two values of λ. The distributions are initially (pre-weighting) flat on average, but fluctuate in accordance with the random sampling described in section 3.2. This explains why the intensity corresponding to the maximum-likelihood model (black-dashed line in figure 3), in spite of having the greatest weight, does not necessarily correspond with the maximum of the distribution.
The credible intervals (analogous to confidence intervals in frequentist statistics) are obtained by integrating the 0.68 and 0.95 containment regions from the λ value corresponding to the maximum-likelihood model outwards. In this way, the 68% and 95% containment regions are defined around the model which best reproduces the data. For some wavelength values, the corresponding distributions are such that the limits fall outside considered range of intensity. For these cases, the containment region is reported without a lower/upper boundary. To test the sensitivity of the method to the prior assumptions, the procedure is repeated varying the smoothing condition between consecutive knots and alternatively sampling a flat distribution in logarithm of the intensities instead of the original linear sampling (see section 3.2). No significant changes are observed throughout the process with respect to the original prior conditions.

RESULTS AND DISCUSSION
The containment regions are calculated assuming SSC emission models from ECPL-leptonic distributions for both Mrk 421 and Mrk 501 intrinsic spectral models, as explained in section 3.1. Figure 4 shows the resulting 68% and 95% containment regions for each of the considered sources as a function of wavelength. The in-tensities corresponding to the model with the highest likelihood are shown in black circles, along with lower limits from galaxy counts and upper limits from direct measurements shown as upward and downward triangles, respectively. Figure 5 shows the resulting 68% and 95% containment regions when combining the results from both sources. The combined intensities from each source's highest likelihood model are shown in black circles along with results from other experiments (Abeysekara & Archer 2019;Acciari & Ansoldi 2019; Abdalla et al. 2017). Everything else in the plot is identical to what is shown in figure 4. From figure 5 it is clear that for some λ-values, namely 30µm < λ for Mrk 421 and λ < 4µm ∪ 10µm < λ for Mrk 501, boundaries cannot be established (see section 3.3), limits in these cases are not reported (these are shown as a dash in table 2).
The containment bands are in good agreement with the region delimited between the lower and upper limits from galaxy counts and direct measurements respectively. However, for wavelengths lower than λ=∼5µm, the containment region is broader than the limits, so results are non-constraining in this range. Combined results show a general downward tendency in the EBL measurement for λ values >5µm, following the lower limits from galaxy counts. This downward trend has been observed also by VERITAS (Abeysekara & Archer 2019) and MAGIC (Acciari & Ansoldi 2019). In the case of HAWC, this could be explained by the fact that the highest observed data points (> 7 TeV), especially for Mrk 501, show an upward trend (see Albert et al. (2021a)), leading to lower EBL density values for this wavelength range. HAWC observed spectrum of Mrk 421 shows a cutoff energy around 5 TeV, this could, in principle, be related to the bump around 20µm seen in Mrk 421 results (see figure 4).

Systematic uncertainties
There are many sources of systematic uncertainties that affect the flux estimation from HAWC observations, mostly coming from discrepancies between data and simulations, arising from mis-modeling of the detector. A complete treatment of HAWC possible systematics is presented in , in which they found that the dominant systematic uncertainties for the spectral flux come from mis-modeling the late light in the air shower and the uncertainty in the modeling of the PMT efficiencies and the charge measured by the PMTs. To estimate the potential effect of these systematics in the overall results, the analysis is repeated simulating extreme detector responses considering these dominant systematics. Results in each case are com-pared with nominal results. Figure 6 shows the relative difference for each wavelength between each considered systematic and the nominal results quantified in the following manner: where I nom ν and I sys ν are the EBL density values, at each wavelength, resulting from the nominal model and the considered extreme systematic model respectively; σ 2 nom and σ 2 sys are the corresponding errors, at each wavelength, for the nominal and systematic results respectively.
No significant (∆ > 3) difference is observed between the results from the nominal model and those coming when considering extreme systematics. It is important to note that the fact systematics are not significant is, in this case, mostly due to the relatively large errors associated with the method, rather than the instrumental uncertainty being small.

Intrinsic model: potential bias
The selection of a particular emission model, magnetic field and emission radius (see section 3.1) introduces a bias that will eventually impact the calculated EBL limits. As mentioned in section 3.1, the choice on the parameter values and model is based on previous studies performed on the markarians (Albert et al. 2021a;Zhu et al. 2016;. To estimate the potential bias coming from the selection of the emission model, the analysis was repeated using an electron distribution following a broken powerlaw (BPL) function instead of the exponentially cut-off powerlaw as the lepton distribution for the emission model. In addition, two models following the contours of the ±1σ confidence band (see figure 2) given by the fit described in section 3.1, are also considered. These latter models are obtained by sampling the fit parameters of the nominal model within their corresponding ±1σ errors and retrieving the minimum and maximum flux values for different values of energy. These flux points are then refitted using Naima to get the corresponding parameters of the ±1σ models to be used in the analysis in an analogous way as performed with the nominal model.
The analysis described in section 3 is repeated for these ±1σ models and the corresponding results are used to quantify the overall effect of the intrinsic uncertainties in the final containment bands. Figure 7 shows the relative difference in the derived limits when considering these models with respect to the nominal model and computed as in section 4.1 by using equation 5. No significant change is observed between the results of the  nominal model and the ones resulting from assuming the ±1σ models and assuming the BPL lepton distribution. This suggests that potential biases introduced by the assumptions in the intrinsic emission model don't have a significant impact in the final results within this wavelength range.

CONCLUSIONS
After years of operation, the HAWC observatory is able to significantly detect extragalactic sources, like the Markarians, up to energies of around 10 TeV. This motivates an analysis that could, in principle, probe the mid-IR region of the EBL, a region that is inaccessible to other gamma-ray instruments. In this study, physically motivated emission models and Fermi-LAT data are used to construct an intrinsic spectrum for each of the Markarians, using the Naima python package. A large number of EBL model shapes are randomly generated and used to obtain "observed" spectra to compare with HAWC data of the Markarians. The com-  parison is performed using threeML, in a way that each EBL model can be assigned a likelihood value that expresses the agreement between that particular spectral realization and HAWC data. The present method has the advantage of being independent of any particular EBL shape and assuming physically motivated emission models as intrinsic spectra. The EBL intensity measurements and containment bands are calculated from 1µm to 100µm, probing higher wavelength values than previous measurements performed with VERITAS by using a similar method (Abeysekara & Archer 2019). Results for both sources are in agreement with current upper limits from direct IR observations and lower limits from galaxy counts. Limits are, in general, less constraining for wavelengths below λ=∼5µm and there is a downward trend when moving to higher λ values, roughly following the lower limits. This trend is clearer in the case of Mrk 501. A bump around 20µm is observed for Mrk 421, possibly due to an observed cutoff drop present in the source observed spectrum around ∼5 TeV Albert et al. (2021a). Results are also in agreement with other instruments EBL intensity estimations reported by HESS (Abdalla et al. 2017) and MAGIC (Acciari & Ansoldi 2019). These results would surely improve their constraining power if more adequate sources were included. HAWC keeps collecting data from extragalactic sources that could soon be used to derive EBL limits by applying the present method. The radio galaxy M87 was pointed out by Franceschini et al. (2019) as a good candidate to perform EBL related studies. At the time of writing, this source was detected at too low significance to be reliable for the analysis, therefore it is not included. However, M87 data are being accumulated (Albert et al. 2021b) and will soon provide a significant detection that will make the source useful for this type of analysis.