Long-term Spectra of the Blazars Mrk 421 and Mrk 501 at TeV Energies Seen by HAWC

Blazars are a particular class of radio-loud active galactic nuclei (AGN), characterized by ultra-relativistic jets escaping from a supermassive black hole that are oriented very close to the observer's line of sight (Urry & Padovani 1995). The spectral energy distributions (SEDs) of blazars are characterized by two peaks: the first at low to medium energies (radio to X-ray), and the second at high energies (gamma-ray; Fossati et al. 1998). The first peak is produced by synchrotron emission from ultra-relativistic charged particles embedded in a magnetic field within the plasma jet. The high-energy peak is thought to be produced by inverse Compton (IC) scattering of low-energy photons. These low-energy photons can be assumed to come from the same charged particle population that generates the synchrotron emission (the synchrotron-self Compton model, SSC; Jones et al. 1974) and/or from an external region, like the AGN accretion disk or the torus of dust that surrounds it (the external Compton model; Sikora et al. 1994). Also, there are multiple models that use different assumptions about the nature of the charged particles that generate the electromagnetic radiation that we see from blazars. On the one hand, there are the lepton models that assume that the accelerated particles are mainly electrons; and, on the other hand, those that assume a population mainly of hadronic particles. There are also lepto-hadronic models that consider multiple physical processes that could take place in the blazar emission zone. The relativistic jets are key to understanding the nature of black holes and their environment; however, little is known about their composition, production, and how gravitational energy is transported to the dissipation zone, where radiation is generated. Therefore, to constrain the physical parameters of the jet, such as the size of the emission region and the magnetic field, emission models have to be applied to gamma-ray observations. The very high-energy (VHE; ≳ 0.1 TeV) observations of blazars are often motivated by flaring activity, and are mainly performed by Imaging Atmospheric Cherenkov Telescopes (IACTs). However, the average activity levels of blazars have been poorly characterized. Although IACTs have performed observation campaigns to characterize the properties of the average emission of blazars, the best instruments for getting averaged spectra over long periods of time are satellites and extensive air shower (EAS) experiments, since their observation times are not biased by flaring activities due to their large duty cycle.

The blazars Mrk 421 (z = 0.031) and Mrk 501 (z = 0.034) were first detected in the VHE range by the Whipple Observatory: Mrk 421 in 1992, above 0.5 TeV (Punch et al. 1992), and Mrk 501 in 1996, above 0.3 TeV (Quinn et al. 1996). Their gamma-ray flux has been measured extensively with IACTs during different activity periods in the 0.1–10 TeV energy range. Also, both sources are continuously detected by the Large Area Telescope on board NASA's Fermi Gamma-ray Space Telescope (Fermi-LAT) in the 50 MeV to 1 TeV energy range (Abdollahi et al. 2020). In the literature, the quasi-contemporaneous multiwavelength spectra have been measured mainly through observation campaigns.

Mrk 421 has been extensively studied in the VHE range. Between 2004 and 2005, MAGIC performed a total of 25.6 hr of observations during a period where Mrk 421 had low-energy flux, identifying an energy cutoff and setting upper limits for energies at and above 4 TeV (Albert et al. 2007a). With the same telescope, in 2006 the average spectra of Mrk 421 was observed during a high-energy flux state for eight nights (12.7 hr in total), where the highest energy detection was reported at 3.3 TeV (Aleksić et al. 2010). Between 2006 and 2008, observations of the source were carried out with the VERITAS telescope, where different periods of activity were identified, having a total of 35.19 hr of observations, of which 9% correspond to observations made when the source had a high-activity state (Acciari et al. 2011a). The ARGO-YBJ experiment used 676 days of observation to report four different averaged spectra of Mrk 421 between 2007 and 2010, where simultaneous X-ray data were used to find correlations between the different states of activity of the source; the highest energy signal was reported at 10 TeV (Bartoli et al. 2011). In 2009, the MAGIC telescope observed Mrk 421 for 27.7 hr as part of the multiwavelength campaign organized by the Fermi collaboration, detecting the source at 4 TeV; emission models were fit to the SED and the results differed from those obtained in previous works, since they were based on observations when the source was in a state of high activity (Abdo et al. 2011a). In Bartoli et al. (2016), using 4.5 yr of data, the ARGO-YBJ experiment reported the averaged spectrum of Mrk 421 for observations performed between 2008 and 2013; the highest energy signal detected was reported at 4.5 TeV. Recently, FACT, an IACT whose core program focuses on long-term observations of Mrk 421 and Mrk 501 (Dorner et al. 2019), reported observations of Mrk 421 for 279 hr between several hundreds of GeV and 10 TeV (MAGIC Collaboration et al. 2021).

Mrk 501 is also a very well-studied VHE source. Shortly after its detection at VHE, during a period of intense high activity in 1997, the HEGRA collaboration reported the averaged spectra from this source for 140 and 85 hr observations performed with the Cherenkov telescopes CT1 and CT2, respectively, with the highest energy detection above 10 TeV (Aharonian et al. 1999). In 2005, the MAGIC telescope made observations of Mrk 501, where it reported variability in its energy flux, identifying three states of activity; the average spectrum during the period of low activity was obtained with 17.2 hr of observations, the intermediate with 11 hr, and the high with 1.52 hr, with the highest energy detection at approximately 4.5 TeV (Albert et al. 2007b). Between 2005 and 2006, the TACTIC telescope observed Mrk 501 for a total of 112.5 hr, reporting an average spectrum in the energy range of 400 GeV to 6 TeV (Godambe et al. 2008). A multiwavelength campaign to study the low-activity state of Mrk 501 was carried out in 2005, where the MAGIC telescope observed the source for 9.1 hr, detecting it at up to 2 TeV; the obtained SED, using Fermi and Suzaku data in X-rays, was fitted with an SSC model that suggested that the high-activity states could be due to variations in the electron population (Anderhub et al. 2009). A few years later, in 2009, the Fermi collaboration organized a multiwavelength campaign to study the Mrk 501 SED, which was built with data from radio to gamma-rays; the MAGIC telescope participated with 16.2 hr of observation and the VERITAS array telescopes with 9.7 hr, of which 2.4 hr were during a high-activity state of the source; an SSC model was fit to the data, showing that relativistic proton shock waves are related to the bulk of the energy dissipation within the source emission zone (Abdo et al. 2011b). Also, between 2012 December and 2018 April, FACT monitored Mrk 501 at TeV energies for a total of 1344 hr, distributed over 633 nights (Arbet-Engels et al. 2021).

The average spectrum of Mrk 501 was also reported by the ARGO-YBJ experiment using 1179.6 days of data taken between 2008 and 2012. This was done for different periods of activity; it was deduced that there was a correlation between the X-ray and gamma-ray fluxes, as well as a hardening of the spectra when the flux increased, thus favoring the SSC model in explaining the physical processes of this source (Bartoli et al. 2012).

We remark that most of the highest gamma-ray photons from Mrk 421 and Mrk 501 were detected in periods of high activity. Also, from the works cited above, we can summarize a range of physical parameters that describe the jet of each source, such as the emission zone, which is assumed to be spherical with a radius R that moves at a relativistic speed v = β c throughout the jet with a bulk Lorentz factor Γ, so that the observed radiation is amplified by a Doppler factor $\delta ={({\rm{\Gamma }}(1-\beta \cos \theta ))}^{-1}$, where θ is the jet inclination angle. For Mrk 421, the Doppler factor value, which accounts for the relativistic effects and depends on the speed of the emission zone and the pitch angle of the jet, ranges from δ = 15–50, the magnetic field varies between B = 39–200 mG, and the size of the emission zone ranges from R = (0.25–5.2) × 10¹⁶ cm. For Mrk 501, the value for the Doppler factor varies between δ = 12–25, the range of the magnetic field goes from B = 15–310 mG, and the size of the emission zone varies from R = (0.1–13) × 10¹⁶ cm. It should be mentioned that the values for the highest δ and B and the lowest R values, for both sources, correspond to models fitted to observations carried out over short periods of time, mainly by IACTs.

Both Mrk 421 and Mrk 501 have been detected by the High Altitude Water Cherenkov (HAWC) Gamma-Ray Observatory above 300 GeV. Each source is observed over 6.2 hr per day, on average, within the field of view of the observatory. The HAWC collaboration first published the detection of these sources in Abeysekara et al. (2017a), which reported on the monitoring carried out over a period of 513 days of observation. The spectral analysis did not consider the extragalactic background light (EBL) attenuation, and was performed by taking into account only the size of the air shower, which is weakly related to the energy of the primary gamma-rays. In the second HAWC catalog of VHE sources (Abeysekara et al. 2017b), which comprised 507 days of observations, the detection of both sources was also reported. The energy spectra of Mrk 421 and Mrk 501 were recently reported using a preliminary implementation of the energy estimators developed by the HAWC collaboration, according to which the spectrum of each source was divided into energy bins for the first time (Coutiño de Leon et al. 2019). That analysis comprised a period of 837 days, and was carried out by using a different framework than the one used in this work (see Section 2.1).

In this work, we report the observations of Mrk 421 and Mrk 501 with HAWC during ∼1038 days of exposure, with the spectral analysis above 0.5 TeV, and we include the systematic uncertainties calculation. We built a contemporaneous SED with Fermi-LAT data with the aim of modeling it with an SSC model in order to improve the constraints on the physical parameters of the blazar jets. These results can be used to study the characteristics of secondary gamma-rays that are produced by IC scattering between the cosmic microwave background and the electron–positron pairs (caused by the interaction between primary gamma-rays and the extragalactic background light), and thus to restrict/constrain the properties of the intergalactic magnetic field (Arlen et al. 2014).

The paper is organized as follows. In Section 2, the HAWC observatory is described, along with the method of measuring the energy flux. In Section 3, the energy spectra of each source are presented. In Section 4, a comparison with previous results is made. In Section 5, the results are used to build a multiwavelength SED to test a one-zone leptonic blazar emission model; and finally a summary and future work are presented in Section 6.

HAWC is located at latitude + 19°N and at an altitude of 4100 m in Puebla, Mexico. It consists of 300 water Cherenkov detectors (WCDs), of 7.3 m diameter and 4.5 m height, spread over an area larger than 22,000 m². Each WCD is filled with 180,000 liters of water and instrumented with four photomultiplier tubes (PMTs) that measure the arrival times and directions of cosmic and gamma-ray primaries, mostly above 300 GeV, within its ∼2 sr field of view.

The data are divided into nine analysis bins (f_hit), according to the fraction of PMTs that are triggered in each shower event (Abeysekara et al. 2017c). Recently developed energy estimators are used to divide the these f_hit bins into 12 quarter-decade energy bins, covering the 0.316–316 TeV range. In this paper, we use the ground parameter method presented in Abeysekara et al. (2019), which uses the measured charge 40 m from the air shower axis, along with the zenith angle of the air shower, to estimate the primary gamma-ray energy. The binning scheme is identical to that in Abeysekara et al. (2019). The data used for this analysis are from 2015 June to 2018 July.

2.1. Fitting Technique

A forward-folding method is performed to fit the spectral shapes of the sources using a maximum-likelihood technique, maximizing the test statistics (TS) so the input parameters have the highest likelihood of providing a good description of the observed data. Assuming a point-source model, the TS is defined as follows:

$$\begin{eqnarray}&&{TS}\equiv 2\mathrm{ln}\displaystyle \frac{{ \mathcal L }({H}_{1})}{{ \mathcal L }({H}_{0})},\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal L }$ is the likelihood function, H₀ is the background hypothesis, and H₁ is the signal plus background hypothesis, which depends on the spectral parameters assumed to describe the sources. The significance maps (Figure 1) are obtained, as explained in Abeysekara et al. (2017b), at the positions of the sources, maximizing the TS value in each pixel of an N_size = 1024 HEALPix grid (Górski et al. 2005). The likelihood calculation is performed as in Younk et al. (2016), using the multi-mission maximum likelihood framework (Vianello et al. 2015) along with the HAWC accelerated likelihood (HAL) ²⁹ plugin.

Zoom In Zoom Out Reset image size

Since the sources are of extragalactic origin, the attenuation due to the EBL is taken into account. The input spectral model is assumed to be the intrinsic one, and it is then attenuated using an EBL model. The resulting spectrum is then convolved with the detector response to be compared with the observed counts. This way, the output parameters correspond to the intrinsic ones. The EBL model used to perform the fits in this work is from Gilmore et al. (2012).

The spectral shapes tested were a simple power law (PL; Equation (2)) and a PL with an exponential energy cutoff (PL+CO; Equation (3)):

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}={N}_{0}{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-\alpha }\times \exp \left[-\tau (E,z)\right],\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}={N}_{0}{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-\alpha }\times \exp \left(\displaystyle \frac{-E}{{E}_{c}}\right)\times \exp \left[-\tau (E,z)\right],\end{eqnarray} \tag{ 3 }$$

where N₀ is the normalization flux [TeV⁻¹cm⁻² s⁻¹], E₀ is the pivot energy fixed at 1 TeV, α is the spectral index, E_c is the energy cutoff [TeV], and τ is the opacity value given by the EBL model, which is an increasing function of E and the source redshift, z.

Depending on the TS values in the global fit using all the available energy bins, a preferred spectral shape is chosen. The flux points are estimated as in Abeysekara et al. (2019), by fitting N₀ in each energy bin, with α and E_c fixed using the resulting values from the global fit. If a fit from an individual energy bin has a TS < 4, an upper limit at a 95% confidence interval is set.

2.2. Energy Range

2.3. Systematic Uncertainties

3.1. Mrk 421

All the best-fit parameters are quoted with their respective statistical and systematic errors. Above 0.5 TeV, the intrinsic spectrum of Mrk 421 is better described by a PL+ CO with ${N}_{0}=\left[4.0\pm {\left(0.3\right)}_{\mathrm{stat}}{\left({}_{-0.2}^{+0.9}\right)}_{\mathrm{sys}}\right]\times {10}^{-11}\,{{\rm{TeV}}}^{-1}{{\rm{cm}}}^{-2}\,{{\rm{s}}}^{-1},$ $\alpha \,=2.26\,\pm \,{\left(0.12\right)}_{\mathrm{stat}}{\left({}_{-0.39}^{+0.17}\right)}_{\mathrm{sys}}$, and ${E}_{c}=5.10\pm \,{\left(1.60\right)}_{\mathrm{stat}}{\left({}_{-2.5}^{+1.4}\right)}_{\mathrm{sys}}$TeV(Table 1). The intrinsic and observed differential energy spectra are shown in Figure 2. After integrating Equation (3) above 0.5 TeV, the observed integrated photon flux is ${N}_{\mathrm{obs}}\left(\gt 0.5\,{\rm{TeV}}\right)$ $=\left[4.3\pm {\left(0.5\right)}_{\mathrm{stat}}{\left({}_{-0.6}^{+0.1}\right)}_{\mathrm{sys}}\right]\times {10}^{-11}\,{\text{ph}\,\text{cm}}^{-2}\,{{\rm{s}}}^{-1}$, and the intrinsic energy flux is f_E ( > 0.5 TeV) $=(25.7\pm {\left(4.2\right)}_{\mathrm{stat}}{\left({}_{-2.6}^{+3.5}\right)}_{\mathrm{sys}})\times {10}^{-12}\,{\text{erg cm}}^{-2}\,{{\rm{s}}}^{-1}$. The maximum energy at which the source is detected is 9 TeV at a 2σ level.

Zoom In Zoom Out Reset image size

Table 1. Best-fit Spectral Parameters for Mrk 421 and Mrk 501, Following the Method Described in Section 2.1

	Mrk 421	Mrk 501
$\sqrt{{TS}}$	48	12
N0 [TeV−1cm−2 s−1]	$\left[4.0\pm {\left(0.3\right)}_{\mathrm{stat}}{\left({}_{-0.2}^{+0.9}\right)}_{\mathrm{sys}}\right]\times {10}^{-11}$	$\left[6.6\pm {(0.9)}_{\mathrm{stat}}{\left({}_{-0.6}^{+0.9}\right)}_{\mathrm{sys}}\right]\times {10}^{-12}$
α	$2.26\pm {\left(0.12\right)}_{\mathrm{stat}}{\left({}_{-0.39}^{+0.17}\right)}_{\mathrm{sys}}$	$2.61\pm {\left(0.11\right)}_{\mathrm{stat}}{\left({}_{-0.07}^{+0.01}\right)}_{\mathrm{sys}}$
Ec [TeV]	$5.1\pm {\left(1.6\right)}_{\mathrm{stat}}{\left({}_{-2.5}^{+1.4}\right)}_{\mathrm{sys}}$	∞
${E}_{{\rm{\max }}}$ [TeV]	9	12

Note. The photon ${E}_{{\rm{\max }}}$ value corresponds to the maximum energy at which the signal is detected at a 2σ level.

Download table as:  ASCIITypeset image

3.2. Mrk 501

For Mrk 501, the best spectral model that describes the data above 0.5 TeV is a single PL with ${N}_{0}=\left[6.6\pm {\left(0.9\right)}_{\mathrm{stat}}{\left({}_{-0.6}^{+0.9}\right)}_{\mathrm{sys}}\right]$ ×10⁻¹² TeV⁻¹cm⁻² s⁻¹ and $\alpha =2.61\pm {\left(0.11\right)}_{\mathrm{stat}}{\left({}_{-0.07}^{+0.01}\right)}_{\mathrm{sys}}$. The maximum energy at which Mrk 501 is detected is 12 TeV at a 2σ level (Table 1). Integrating Equation (2) above 0.5 TeV, the observed integrated photon flux is N_obs( > 0.5 TeV) $=\left[9.1\pm {\left(1.2\right)}_{\mathrm{stat}}{\left({}_{-0.8}^{+1.3}\right)}_{\mathrm{sys}}\right]\times {10}^{-12}\,{\text{ph cm}}^{-2}\,{{\rm{s}}}^{-1}$, and the intrinsic energy flux is ${f}_{E}(\gt 0.5\,{\rm{TeV}})=\left[25.7\pm {\left(4.2\right)}_{\mathrm{stat}}{\left({}_{-2.6}^{+3.5}\right)}_{\mathrm{sys}}\right]$ ×10⁻¹² ergcm⁻² s⁻¹. The intrinsic and observed spectra are shown in Figure 3.

Zoom In Zoom Out Reset image size

4.1. Mrk 421

The value of the highest energy flux point in the HAWC spectrum for Mrk 421 is at 8.8 TeV. This value is one of the highest energy detections for a long-term time-averaged spectrum reported to date. The intrinsic spectrum of Mrk 421 is in good agreement, within the statistical and systematic errors, with the spectra previously reported by HAWC (Coutiño de Leon et al. 2019), and with the averaged spectra reported in Albert et al. (2007a) and Bartoli et al. (2011). As mentioned before, most of the IACT observations are biased to the high-activity state of the source, showing not only a higher flux, but also a harder spectrum, following a "brighter–harder" relation. This can be seen in Figure 4, where the spectral index (left panel) and the energy cutoff (right panel) are plotted against the normalization flux for different intrinsic values reported in the literature. The fitted values in this work for Mrk 421 lie between the values from short- and long-term observations. The observed spectra of Mrk 421 reported in the literature that best coincide with our flux points are those reported by VERITAS for a very low activity state (Acciari et al. 2011a); the spectrum measured by MAGIC during the Fermi multiwavelength campaign, which was fitted to a single log-parabola (Abdo et al. 2011a); and the average spectrum reported by the ARGO-YBJ experiment (Bartoli et al. 2016), fitted with a single PL (see Figure 5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2. Mrk 501

For Mrk 501, the intrinsic spectrum is in agreement with the previous results obtained by HAWC (Coutiño de Leon et al. 2019); it is also in agreement with the intrinsic spectrum of the ARGO-YBJ experiment for a 1179.6 day observation period, where α = 2.59 ± 0.27 (Bartoli et al. 2012). The trend of having a harder spectrum when the source is in a high-activity state is not as noticeable as with Mrk 421, as shown in Figure 6, where the spectral index is plotted versus the normalization flux, from the intrinsic spectra reported in the literature.

Zoom In Zoom Out Reset image size

The energy of the last flux point bin is at 10.90 TeV, which is also one of the highest energy detections for time-averaged spectra to date. The flux points obtained with the HAWC data, compared to previous observations made with IACTs and ARGO-YBJ, are shown in Figure 7, where the HAWC flux points at higher energies are below previous observations by a factor of ∼6–7. This difference can be explained in terms of the Mrk 501 activity state during those observations, such as those reported in Abdo et al. (2011b), where a high-energy state was detected.

Zoom In Zoom Out Reset image size

We note that the intrinsic spectral parameters depend on the choice of EBL model. For the redshifts of Mrk 421 and Mrk 501 (z = 0.031 and z = 0.034, respectively), and in the detection energy range (0.5 < E < 10 TeV), the opacity values of the different EBL models that were used in the works cited in Sections 4.1 and 4.2 are not significantly different.

Mrk 421 and Mrk 501 are blazars classified as high-synchrotron peaked BL Lacs (HBLs; Padovani & Giommi 1995; Fossati et al. 1998). These types of objects are characterized as emitting most of their power in the UV and X-ray range. HBL blazars are also characterized as having a low luminosity, so it is thought that the only seed photons that are scattered at very high energies are synchrotron photons; that is, there is no contribution from external photons from the broad-line region, torus, or accretion disk (Tavecchio & Ghisellini 2008; Madejski & Sikora 2016). For this reason, we assume the SSC scenario for this work. In addition, SSC models have been widely used to describe the long-term averaged spectra (Abdo et al. 2011a, 2011b; Bartoli et al. 2011, 2012, 2016; Fraija et al. 2017), the short-term averaged spectra (Albert et al. 2007a, 2007b; Anderhub et al. 2009; Aleksić et al. 2010; Acciari et al. 2011a), and SEDs and flares (Acciari et al. 2011b; Bartoli et al. 2012, 2016; Fraija et al. 2017); and, in general, all describe the observational data well. However, it should be noted that this family of models produces parameter degeneration (Ahnen et al. 2017), so it is important to be able to restrict some of the models, as far as possible.

The electron population within the emission zone has a total energy given by

$$\begin{eqnarray}&&{W}_{e}={\int }_{{E}_{\min }}^{{E}_{\max }}{E}_{e}\displaystyle \frac{{{dN}}_{e}}{{{dE}}_{e}}{{dE}}_{e},\end{eqnarray} \tag{ 4 }$$

where E_min and E_max are the minimum and maximum electron energies, respectively, and ${{dN}}_{e}/{{dE}}_{e}$ is the particle distribution embedded in a magnetic field B. The electron population is accelerated to relativistic velocities by the magnetic field, producing synchrotron radiation, which is then used as a seed photon field for the IC scattering.

The total radiative energy output of the jet L_jet = L_e + L_p + L_B is the sum of the radiative output carried by the electrons, protons (under the assumption of one cold proton per emitting relativistic electron) and magnetic field, which are defined as

$$\begin{eqnarray}&&{L}_{i}\simeq \pi {R}^{2}c\,{{\rm{\Gamma }}}^{2}{U}_{i},\quad i=e,p,B,\end{eqnarray} \tag{ 5 }$$

with U_e = 3W_e /4π R³, U_p = m_p N_p , and U_B = B²/8π being the energy densities, and m_p being the proton mass (Celotti & Ghisellini 2008).

5.1. Multifrequency Data

5.2. SED Modeling

5.3. SED Modeling Results

5.3.1. Mrk 421

For Mrk 421, the best fit corresponds to a Doppler factor of δ = 25 ± 1 for an electron energy distribution that follows a broken PL with energy break E_break = 112 ± 41 GeV and spectral indexes before and after the break of α₁ = 2.21 ± 0.01 and α₂ = 5.4 ± 0.1, respectively. The minimum and maximum electron energies are best fits to ${E}_{\min }=425\pm 8\,\mathrm{MeV}$ and ${E}_{\max }=51\pm 13$ TeV. The magnetic field results in a value of B = 24 ± 6 mG. Figure 8 shows the multifrequency SED of Mrk 421 and the best SSC model (red curve). Table 2 shows a summary of the SED fitting results, and Table 3 shows the comparison of the results from Table 2 with previous analyses.

Zoom In Zoom Out Reset image size

Table 2. Fitted Parameters in the SSC Leptonic Model for Mrk 421

Parameter	Symbol	Mrk 421
Doppler factor	δ	25 ± 1
Magnetic field	B [mG]	24 ± 6
Spectral index before break	α1	2.21 ± 0.01
Spectral index after break	α2	5.4 ± 0.1
Energy break	Ebreak [GeV]	112 ± 41
Minimum electron energy	Emin [MeV]	425 ± 8
Maximum electron energy	Emax [TeV]	51 ± 13
Jet power in electrons	Le × 1044 [erg s−1]	2.3
Jet power in protons	Lp × 1044 [erg s−1]	4.0
Jet power in magnetic field	LB × 1042 [erg s−1]	6.8

Download table as:  ASCIITypeset image

Table 3. Comparison between the Previous SSC Parameters, δ, B, and R, and the Results in This Work for Mrk 421

δ	B	R	Flux State	Reference
	[mG]	× 1016 [cm]
15	200	1.1	Low	Albert et al. (2007a)
15	150	5	High	Bartoli et al. (2011)
16	80	5	Low
40	200	0.25	Low–Mid–High	Acciari et al. (2011a)
21	38	5.2	Long-term averaged	Abdo et al. (2011a)
${38}_{-4}^{+6}$	48 ± 0.012	1	Long-term averaged	Bartoli et al. (2016)
25 ± 1	24 ± 6	5	Long-term averaged	This work

Download table as:  ASCIITypeset image

5.3.2. Mrk 501

For Mrk 501, the SED is better described with a Doppler factor value of δ = 13 ± 0.7. The electron energy distribution that results in a better fit is a broken PL with energy break E_break = 470 ± 20 GeV and spectral indexes before and after the break of α₁ = 2.1 ± 0.01 and α₂ = 4.6 ± 0.1, respectively. The best-fit values of the minimum and maximum electron energy are ${E}_{\min }=166\pm 8\,\mathrm{MeV}$ and ${E}_{\max }=19\pm 0.5$ TeV. According to the IC scattering process, the energy of the VHE photons must not exceed that of the electrons, so accounting for the Doppler boosting, the E_max value agrees with our observations. The magnetic field is fitted to a value of B = 20 ± 5 mG. This model is shown in Figure 9 (red curve), and a summary of these results is given in Table 4.

Zoom In Zoom Out Reset image size

Table 4. Fitted Parameters in the SSC Leptonic Model for Mrk 501

Parameter	Symbol	Mrk 501
Doppler factor	δ	13 ± 0.7
Magnetic field	B [mG]	20 ± 5
Spectral index before break	α1	2.1 ± 0.01
Spectral after before break	α2	4.6 ± 0.1
Energy break	Ebreak [GeV]	470 ± 20
Minimum electron energy	Emin [MeV]	166 ± 8
Maximum electron energy	Emax [TeV]	19 ± 0.5
Jet power in electrons	Le × 1044 [erg s−1]	2.7
Jet power in protons	Lp × 1044 [erg s−1]	4.1
Jet power in magnetic field	LB × 1042 [erg s−1]	4.1

Download table as:  ASCIITypeset image

In Tables 3 and 5, we provide a comparison between our results and previous analyses using an SSC model for Mrk 421 and Mrk 501, respectively. As can be seen, the most notable difference lies in the value of the magnetic field, which is up to an order of magnitude larger than our results for analyses that were carried out using VHE data averaged over short periods or flares. This can also be noted for the size of the emission zone of Mrk 501, whose value is up to two orders of magnitude larger than for the long-term averaged spectra. The small differences between our results and those in the literature that include VHE data averaged over long periods of time show a good agreement between the outcomes of the newly analyzed HAWC data presented here and those previously reported from other experiments.

Table 5. Comparison between Previous SSC Parameters, δ, B, and R, and the Results in This Work for Mrk 501

δ	B	R	Flux State	Reference
	[mG]	× 1016 [cm]
25	310	0.1	Low	Albert et al. (2007b)
20	313	0.103	Low	Anderhub et al. (2009)
12	15	13	Long-term averaged	Abdo et al. (2011b)
12	70	3	Long-term averaged	Bartoli et al. (2012)
10	100	3	High
13 ± 0.7	20 ± 5	10	Long-term averaged	This work

Download table as:  ASCIITypeset image

The jet power in the electrons for both sources is comparable to that of the protons, L_e ∼ L_p , and both are larger than the jet power carried by the magnetic field, L_e > L_B , thus the Poynting flux does not contribute significantly to the total radiation of the jet. The total radiative energy output of the Mrk 421 jet is then L_jet−421 = 6.5 × 10⁴⁴ erg s⁻¹, which corresponds to ∼4% of the Eddington luminosity, and for Mrk 501 it is L_jet = 6.1 × 10⁴⁴ erg s⁻¹, which represents ∼0.3% of the Eddington luminosity.

We report the detection of Mrk 421 and Mrk 501, above 0.5 TeV, with the High Altitude Water Cherenkov (HAWC) Gamma-Ray Observatory, using 1038 days of exposure comprising the period between 2015 June and 2018 July.

• 1.  
  For Mrk 421, the VHE intrinsic spectrum is well described by a PL+CO. For Mrk 501, the intrinsic VHE spectrum is described by a single PL.
• 2.  
  These results are in good agreement with those previously obtained with HAWC for both sources, once the EBL attenuation is taken into account. Additionally, the reported values for the intrinsic spectra in this work are compatible with those in the previous averaged spectra reported by IACTs and EAS experiments, setting a baseline energy spectrum for each source. It is also important to mention that the obtained flux points in this work are in good agreement with the observed spectra reported in the literature for Mrk 421; however, for Mrk 501, the HAWC flux points lie below the observed spectra reported in the literature, which could be related to the activity state of the source when it was observed in the past.
• 3.  
  Compared to the previously published results using HAWC data, this is the first time that we have estimated the highest energies of the detected signals, with 9 TeV for Mrk 421 and 12 TeV for Mrk 501 at 2σ confidence levels; for each source, these individual values are some of the highest energy detections reported to date, for spectra averaged over long periods of time. This contributes to the restriction of the energy detection limits for both sources.
• 4.  
  The SEDs built using contemporaneous data from HAWC, Fermi-LAT, Swift-XRT, and Swift-UVOT, together with the previously published data in the radio to optical energy range, were modeled using one-zone SSC scenarios.
• 5.  
  The estimated physical parameters from the jets are in general agreement with the values found in the literature for long-term observations, confirming that both sources are intrinsically different, assuming that the same physical processes take place.

To characterize the spectra at VHE of Mrk 421 and Mrk 501 in greater detail, it is important to identify the periods of variability of both sources and thus carry out the spectral analyses for each of them; this way, the physical processes that give rise to these energy flux variations can be constrained. To achieve this, a time-resolved analysis of the data set used in this work is necessary, and will be addressed in future publications.

We acknowledge support from: the US National Science Foundation (NSF); the US Department of Energy, Office of High-Energy Physics; the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory; Consejo Nacional de Ciencia y Tecnología (CONACyT), México, grants 271051, 232656, 260378, 179588, 254964, 258865, 243290, 132197, A1-S-46288, A1-S-22784, cátedras 873, 1563, 341, and 323; Red HAWC, México; DGAPA-UNAM grants IG101320, IN111315, IN111716-3, IN111419, IA102019, and IN112218; VIEP-BUAP; PIFI 2012, 2013; PROFOCIE 2014, 2015; the University of Wisconsin Alumni Research Foundation; the Institute of Geophysics, Planetary Physics, and Signatures at Los Alamos National Laboratory; Polish Science Centre grant DEC-2017/27/B/ST9/02272; Coordinación de la Investigación Científica de la Universidad Michoacana; Royal Society—Newton Advanced Fellowship 180385; Generalitat Valenciana, grant CIDEGENT/2018/034; and Chulalongkorn University's CUniverse (CUAASC) grant. Thanks to Scott Delay, Luciano Díaz, and Eduardo Murrieta for technical support. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester.