Daily Monitoring of TeV Gamma-Ray Emission from Mrk 421, Mrk 501, and the Crab Nebula with HAWC

HAWC can record extensive air showers from all directions visible above the horizon. Due to the increasing absorption of secondary particles in the atmosphere, the actual effective area for gamma-rays is a function of the zenith angle of the primary particle and the contribution of events from outside a cone with an opening angle of ∼45° around zenith is usually small. The field of view thus spans a solid angle of ∼2 steradians (sr) and HAWC is most sensitive to sources between declinations −26° and +64°. With the rotation of the Earth, any location in this declination range passes over HAWC once every sidereal day. In the following, a transit is defined by visibility over HAWC at zenith angles $\theta \lt 45$° and lasts approximately 6 hr for the sources discussed in this paper. The detection efficiency is not uniform during the transit and Figure 1 shows the expected fraction of signal as a function of time relative to culmination. Approximately 90% of the signal events arrive within the central ∼4 hours of a transit for a source modeled on the Crab Nebula (photon index ${\rm{\Gamma }}=2.63$ and declination 22°). While the shape of this event distribution as a function of transit time can, in principle, change for different spectra and declinations, it is not significantly altered for the sources discussed in this paper, which culminate within $\leqslant 20$° of zenith. For the flux measurement over a full transit, the zenith dependence of HAWC's sensitivity simplifies to a dependence on the source's declination that determines the expected excess and energy distribution.

Zoom In Zoom Out Reset image size

In order to process the data in units that do not contain more than one full transit for any source, all reconstructed events are sorted into sidereal days, starting at midnight local sidereal time at the HAWC site.³³ For each sidereal day and each of the nine analysis bins, a sky map of event counts is produced by populating pixels on a HEALPix grid (Gorski et al. 2005) with an average spacing of $\sim 0\buildrel{\circ}\over{.} \,06$. These maps are still dominated by hadronic background events and we use direct integration (Atkins et al. 2003) to obtain a background estimate. In this procedure, a local efficiency map is created by averaging counts in a strip of pixels over two hours in right ascension around any location. We smooth this efficiency map via a spline fit to compensate for the limited statistics in higher analysis bins. The pixels near the strongest known sources and the galactic plane are excluded during the averaging in order not to bias the result by counting gamma-ray events as background. The estimated background counts in each pixel are stored in a second map with the same grid structure.

A quality selection is applied before including data in the maps. First, monitoring of the stability of the angular distributions of reconstructed background events is used to exclude data taken during unstable conditions, for example, related to maintenance. In order to control rate stability during a sidereal day, we then fit the detector rate with a function that follows tidal effects of the atmosphere and reject short periods of data that significantly deviate from this fit. To ensure a uniform detector response, sidereal days with partial coverage are not included for a given source if the lost signal fraction is expected to exceed 50%. This expected coverage fraction is calculated by integrating the signal distribution from Figure 1 only over those sections of the transit that are included in the data, assuming a uniform flux during 6 hr. The different right ascensions of the three sources lead to slightly different exposures, which are reported in Table 1. The total observation time is calculated based on 6 hr of effective HAWC observations for an uninterrupted transit and is corrected for gaps in the case of partial coverage. For the period of 513 sidereal days included in this analysis, on average, 92% of transits or 22% of actual time per source are covered.

Table 1.  Observation Time Per Source after Quality Cuts

Source	Included Transits	Time At Zenith Angles $\lt 45$°
		(hr)
Crab	472	2700
Mrk 421	471	2665
Mrk 501	479	2750

Download table as:  ASCIITypeset image

For this light-curve analysis, the standard HAWC maximum-likelihood method (Younk et al. 2016) is applied to the sidereal day maps in order to fit photon fluxes for each transit of selected source locations. The two extragalactic sources discussed here are modeled as gamma-ray point sources with differential flux energy spectra described by a power law with normalization F at 1 TeV, photon index Γ and an optional exponential cut-off E₀:

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{\mathrm{ph}}}{{dE}}=F{\left(\displaystyle \frac{E}{1\mathrm{TeV}}\right)}^{-{\rm{\Gamma }}}\exp \left(-\displaystyle \frac{E}{{E}_{0}}\right).\end{eqnarray} \tag{ 1 }$$

In the HAWC likelihood analysis framework, this input flux is convolved with a detector response function that includes the point-spread function and efficiency of triggers and cuts, depending on primary energy and incident angle. For one source transit over HAWC, the signal hypothesis contributions as a function of zenith angle are summed and yield the expected number of events ${S}_{b,p}$ per analysis bin b (ranging from 1 to 9) and pixel p (for all pixels within a radius of 3° around the source).

In cases where the coverage of a source transit is interrupted, for example, due to detector down time, the lost signal fraction compared to a full transit is calculated by excluding the gap period from the integration over zenith angles (see Figure 1) and the expected event count is reduced accordingly. For the source hypothesis defined by $\{F,{\rm{\Gamma }},{E}_{0}\}$ and the observation ${\boldsymbol{N}}$ of numbers of events in all bins and pixels, we express the likelihood as

$$\begin{eqnarray}&&{{ \mathcal L }}_{S}({\boldsymbol{N}},\{F,{\rm{\Gamma }},{E}_{0}\})=\displaystyle \prod _{b}\displaystyle \prod _{p}P({N}_{b,p},{\lambda }_{b,p}),\end{eqnarray} \tag{ 2 }$$

where $P({N}_{b,p},{\lambda }_{b,p})$ is the Poisson distribution for a mean expectation ${\lambda }_{b,p}={S}_{b,p}+{B}_{b,p}$, the sum of the expected signal (S) and the number of background (B) events estimated from data for analysis bin b and pixel p. In the likelihood ratio test, the result of Equation (2) is compared to the likelihood value ${{ \mathcal L }}_{B}$ for a background-only assumption (${S}_{b,p}=0$). We express this ratio as the difference of the logarithms of the two likelihood values and define the standard test statistic as

$$\begin{eqnarray}&&\mathrm{TS}=2{\rm{\Delta }}\mathrm{ln}{ \mathcal L }=2(\mathrm{ln}({{ \mathcal L }}_{S})-\mathrm{ln}({{ \mathcal L }}_{B})).\end{eqnarray} \tag{ 3 }$$

TS is then numerically maximized by iteratively changing the input parameters, yielding those values that have the highest likelihood of describing the observed data for the point-source model assumption.

For the analysis in this paper, the normalization F, the spectral index Γ, and the cut-off value E₀ in Equation (1) were allowed to vary when fitting the spectral shape with the time-integrated data of the whole period. For the light-curve measurements, the spectral parameters Γ and E₀ were kept constant and only the normalization F was left free to vary in the likelihood maximization, since the counts during a single transit are often not sufficient for a multi-parameter fit to converge.

In the light curves shown in the results section, we include all flux measurements and their uncertainties (1 standard deviation), even if they do not constitute a significant detection by themselves. The likelihood-maximization procedure can produce negative flux normalizations. These are obviously non-physical as gamma-ray flux measurements but occur when low statistics lead to an underfluctuation of the event count compared to the background estimate in a sufficient number of analysis bins.

3.3.1. Likelihood Variability Test

The maximum-likelihood approach is also used to test if the daily flux measurements in a light curve are consistent with a source flux that is constant in time over the whole period under consideration. We consider the likelihood ${{ \mathcal L }}_{i}(M)$ for the observation in time interval i under two different assumptions for the hypothesis M.

• 1.  
  ${{ \mathcal L }}_{i}({F}_{i})$, where F_i is the best-fit flux value for time interval i, as obtained from a likelihood maximization with only this flux as a free parameter. The light curves show these flux values F_i.
• 2.  
  ${{ \mathcal L }}_{i}({F}_{\mathrm{const}})$, where ${F}_{\mathrm{const}}$ is the best-fit flux value for the time-integrated data set, as obtained from a likelihood maximization with only this flux as a free parameter.

These definitions allow us to compare the likelihood of individual flux measurements with that of a constant flux. Similar to Section 3.6 of Nolan et al. (2012), we define a test statistic as twice the differences between the logarithms of these likelihood values, summed over all intervals:

$$\begin{eqnarray}&&{\mathrm{TS}}_{\mathrm{var}}=2\displaystyle \sum _{i}(\mathrm{ln}{{ \mathcal L }}_{i}({F}_{i})-\mathrm{ln}{{ \mathcal L }}_{i}({F}_{\mathrm{const}})).\end{eqnarray} \tag{ 4 }$$

If the null hypothesis of a constant flux is true, then the distribution of TS_var can be approximated as ${\chi }^{2}({n}_{\mathrm{dof}}-1)$, according to Wilks' theorem (Wilks 1938). By applying this variability test to light curves of empty sky locations, we found that ${\chi }^{2}(n-1)$ indeed matches the distribution of TS_var for random fluctuations around zero if we use an effective $n=1.06{n}_{\mathrm{dof}}$, where ${n}_{\mathrm{dof}}$ is the number of degrees of freedom in the light curve. We calculate the probability for a given source to be consistent with the constant flux hypothesis by integrating the ${\chi }^{2}(n-1)$ distribution above the TS_var value obtained for the light curve of that source.

3.3.2. Bayesian Blocks

A detailed comparison of simultaneous multiwavelength data for features observed in the HAWC light curves is beyond the scope of this paper. Instead, we present a first look at multi-instrument comparisons of unbiased, long-term monitoring that, like the HAWC data, provide daily binning and are not affected by seasonal visibility or weather-related gaps. Public data with comparable sampling and duty cycle for observations of Mrk 421 and Mrk 501, in particular, no gaps larger than a few days, are currently only available from very few other monitoring instruments. We checked the lower energy gamma-ray light curves with daily binning from the Fermi Large Area Telescope (LAT) Monitored Source List,³⁴ with an energy coverage from 100 MeV to 300 GeV. For neither Mrk 421 nor Mrk 501 strong flares were detected on a 1 day timescale and none of the daily averaged integral fluxes exceeded a typical Fermi-LAT alert threshold of 10⁻⁶ ph cm⁻² s⁻¹. These results were generated by an automated analysis pipeline and are not suitable for detailed comparisons of absolute fluxes. A dedicated analysis and the study of the correlation between the high energy emission detected by Fermi-LAT and the HAWC TeV data will be presented in a forthcoming publication.

In the X-ray band, we can compare our data to the daily light curves provided by the Swift/Burst Alert Telescope (BAT; Krimm et al. 2013). This instrument covers energies between 15 and 50 keV and, for catalog sources like those discussed here, has a median exposure of 1.7 hr per day that can vary throughout the year but stays $\lt 5.4$ hr for 95% of the days. The Swift-BAT light curves are sampled with one data point per day, based on Modified Julian Dates (MJD), and are thus not perfectly aligned with the binning in local sidereal days that was chosen as a natural frequency for the HAWC data. The observations are also not necessarily exactly simultaneous on the scale of hours.

A detailed analysis of systematic uncertainties of gamma-ray fluxes measured with HAWC in Abeysekara et al. (2017b) concludes with estimating a ±50% uncertainty in the flux normalization. This uncertainty affects all per-transit flux measurements only as a common change in absolute scaling and thus does not impact the relative magnitude of daily flux measurements. The results of variability studies and change point identification can only be affected by systematic uncertainties that change between individual sidereal days or time periods. The calibration is monitored and is very stable. Updates of calibration parameters were only performed to accommodate hardware changes, for example, additions of PMTs. The remaining hardware-related potential source of variability is removal or replacement of individual PMTs due to maintenance and repairs. This has been found to affect flux measurements by less than ±5%. A higher level systematic uncertainty that could, in principle, affect individual fluxes is due to the possibility that the blazar spectra vary with time or flux state (see, e.g., Krennrich et al. 2002). We simulated gamma-ray fluxes with different spectral parameters that are allowed within the uncertainties of our data and analyzed them with the fixed parameters used in the light-curve analysis. For each source, we determined an optimal threshold above which we perform the analytical flux integration by requiring that the difference between the photon fluxes for different spectral hypotheses is minimal. These threshold values are used when quoting the photon fluxes in the results sections: 1 TeV for the Crab Nebula, 2 TeV for Mrk 421, and 3 TeV for Mrk 501. The resulting uncertainty on individual flux values under spectral hardening or softening is ±5%. The combination of these two potentially time-dependent systematic uncertainties is significantly smaller than the statistical uncertainty of the per-transit flux values and thus marginal with respect to the analysis of variability features.

We have performed further tests of the robustness of the likelihood variability estimation. Using the Crab Nebula as a reference, we found that changes in the analysis procedure with respect to background estimation (different smoothing procedures within the direct integration), data selection (excluding the three highest bins with an average $\leqslant 1$ photon per day for our sources), and spectral model (power law and log parabola) only changed the resulting TS_var value of the variability by $\leqslant 0.2$ standard deviations of the null hypothesis. We therefore find no indication that the map-making and likelihood analysis can introduce significant variability features. We conclude that our analysis of flux variations and identification of flaring states in this paper is not limited by these systematic uncertainties.

For the interpretation of absolute flux values, it is helpful to compare our measurements to a gamma-ray reference flux. We therefore convert fluxes to multiples of the HAWC-measured Crab Nebula flux ($1.89\,\times \,{10}^{-11}$ ph cm⁻² s⁻¹, see the detailed analysis in Abeysekara et al. 2017b) as Crab Units (CU) with a common threshold of 1 TeV. This threshold was chosen to provide easier comparisons between the different sources and the literature. We still have to consider the flux uncertainty introduced by the choice of a fixed spectral assumption for which the analytical integration above 1 TeV is performed. We used the time-integrated HAWC data to fit the flux normalization of each of the three sources discussed here with a number of different power-law indices and cut-off values that cover the individual statistical and systematic uncertainty range. We find that the maximum uncertainty of the photon fluxes in CU due to the spectral assumption is ±25%.

The Crab Nebula is the brightest galactic TeV point source. A detailed analysis of time-integrated HAWC data for this source is presented in Abeysekara et al. (2017b). In Figure 2, we show the results of applying the likelihood analysis to the sidereal day maps at the location of the Crab Nebula. We use a fixed spectrum with index ${\rm{\Gamma }}=2.63$ in Equation (1) and no exponential cut-off, ${E}_{0}\to \infty$, based on the best-fit value obtained in the HAWC catalog (Abeysekara et al. 2017a).

Zoom In Zoom Out Reset image size

The left-hand y-axis in Figure 2 indicates the photon flux (ph cm⁻² s⁻¹) after analytically integrating the spectrum above 1 TeV for the best-fit normalization. The right-hand axis shows the CU defined by dividing the flux by the time-averaged HAWC measurement of the Crab flux, also indicated as a dashed line in the figure. We use MJD for labeling the time axes and highlight the duration of HAWC measurements (6 sidereal hours) through horizontal bars.

We applied the variability test outlined in Section 3.3 to the light curve with 1-transit intervals and found a ${\mathrm{TS}}_{\mathrm{var}}=517.9$, with a probability of 0.292 (1.1 standard deviations) of measuring the same or a larger TS value for a constant flux hypothesis. An analysis of the light curve with the Bayesian blocks algorithm with a false positive probability of 5% reveals no change points. HAWC daily flux measurements thus show no indication of variability in data from the Crab Nebula.

In Figure 3, we show a histogram of $({F}_{i}-\bar{F})/{\sigma }_{i}$, where F_i and ${\sigma }_{i}$ are the fluxes and uncertainties from Figure 2 and $\bar{F}$ is the best-fit value for a constant flux. A fit to a Gaussian function yields a center at 0.035 ± 0.050 and a width of 1.033 ± 0.036, confirming that the observed flux distribution is consistent with arising from a constant source flux.

Zoom In Zoom Out Reset image size

Based on measurements by other instruments, the Crab is generally believed to be a steady source³⁵ at TeV energies (Abramowski et al. 2014; Aliu et al. 2014; Bartoli et al. 2015). The non-detection of variability in TeV emission from the Crab Nebula with HAWC is in agreement with these results. We conclude that HAWC daily light-curve measurements and the likelihood variability check are a robust test of the steady gamma-ray source hypothesis and that any systematic uncertainties in the HAWC data are very unlikely to mimic significant variability in this analysis.

Given that the Crab Nebula is known to flare in lower energy bands, we can use the unique daily TeV light-curve data to constrain any TeV flux enhancement during such episodes. During the 17 months included here, the Fermi-LAT collaboration reported an increased gamma-ray flux for energies $\gt 100\,\mathrm{MeV}$ between 2015 December 28 and 2016 January 9, reaching up to ∼1.7 times the average flux (Buehler & Ciprini 2016). The maximum HAWC 1-transit flux during this period was $(2.14\pm 0.54)\,\times \,{10}^{-11}$ ph cm⁻² s⁻¹ above 1 TeV on 2016 January 7, only 0.46 standard deviations above the average flux. This corresponds to an upper limit at the 95% confidence level of $3.04\,\times \,{10}^{-11}$ ph cm⁻² s⁻¹ above 1 TeV, 1.6 times the average flux. When we combine the HAWC measurements over the 12 transits³⁶ included in this period, we obtain a flux measurement of $(1.42\pm 0.15)\,\times \,{10}^{-11}$ ph cm⁻² s⁻¹ above 1 TeV, 0.75 times the average flux and consistent with a random fluctuation. We conclude that we observe no significant change in the TeV flux during this MeV flare period.

Mrk 421 is a BL Lacertae type blazar with a redshift of z = 0.031 (Mao 2011). It was the first extragalactic object discovered in the TeV band (Punch et al. 1992) and has been extensively studied by many TeV gamma-ray observatories. Mrk 421 is known to exhibit a high degree of variability in its emission and yearly average fluxes are known to vary between a few tenths and ∼1.9 times the flux of the Crab Nebula (Acciari et al. 2014). Variability has been observed down to timescales of hours or less and its spectral shape is known to vary with its brightness (Krennrich et al. 2002).

By using the time-integrated HAWC data for Markarian 421, we fit the spectral shape with the likelihood methods discussed in Section 3.2, leaving the normalization F, the photon index Γ, and the exponential cut-off E₀ free. The resulting best-fit values are $F=(2.82\pm {0.19}_{\mathrm{stat}}\pm {1.41}_{\mathrm{sys}})\,\times \,{10}^{-11}$ TeV⁻¹ cm⁻² s⁻¹ for the time-averaged normalization at 1 TeV, a photon index ${\rm{\Gamma }}\,=2.21\pm {0.14}_{\mathrm{stat}}\pm {0.20}_{\mathrm{sys}}$ and an exponential cut-off at ${E}_{0}\,=5.4\pm {1.1}_{\mathrm{stat}}\pm {1.0}_{\mathrm{sys}}$ TeV. The significance of this description over the background-only hypothesis is $\mathrm{TS}=1232.47$ or 35.1 standard deviations. When compared to a pure power-law hypothesis (${E}_{0}\to \infty$), the fit with a cut-off is clearly the better description, preferred at ${\rm{\Delta }}\mathrm{TS}=64.8$ or 8 standard deviations. We use these values for index and cut-off as a fixed set of parameters when fitting the flux normalization for each sidereal day, reported in the following section.

The flux light curve for Mrk 421 with 1-transit intervals is shown in Figure 4. Photon flux units (left y-axis) are based on analytical integration of the fixed spectral shape above a threshold of 2 TeV that minimizes flux uncertainties due to spectral variations. The conversion to CU (right axis) is based on average HAWC Crab Nebula measurements, see Section 4, for a common threshold of 1 TeV in order to allow comparisons between the different sources. The average flux for the 17 month period is determined via a fit of the combined data under a constant flux assumption and yields $(4.53\pm 0.14)\,\times \,{10}^{-12}$ ph cm⁻² s⁻¹ above 2 TeV.

Zoom In Zoom Out Reset image size

Applying the likelihood variability test to this light curve yields TS${}_{\mathrm{var}}=1154.9$, which corresponds to a p-value $4.40\,\times \,{10}^{-54}$ based on the expected ${\chi }^{2}$ distribution for constant flux models and clearly shows the variable nature of the TeV emission from Mrk 421. The highest per-transit flux value, $(2.94\pm 0.37)\,\times {10}^{-11}$ ph cm⁻² s⁻¹ above 2 TeV was measured for MJD 57238.74–57238.99 (2015 August 04 UTC 17:40–23:40), with a pre-trial significance of 9.3 standard deviations compared to the null hypothesis. A flux that is only slightly lower, $(2.91\,\pm 0.38)\,\times \,{10}^{-11}$ ph cm⁻² s⁻¹, was observed during MJD 57020.33–57020.58 (2014 December 29 UTC 8:00–14:00) and highlights that the variability occurs on timescales of less than one day, since the flux value for the day before and after this maximum are a factor of ∼3 and ∼4 lower, respectively.

The Bayesian blocks algorithm with a prior corresponding to a false positive probability of 5% identifies 18 change points in the light curve shown in Figure 4. The flux amplitudes for the periods between two change points that are consistent with a constant flux are included as blue lines with a shaded region for the statistical uncertainty of one standard deviation. These block positions and amplitudes are listed in Table 2.

Table 2.  HAWC Bayesian Blocks for Mrk 421

MJD Start	MJD Stop	Duration	Flux $\gt \,2$ TeV
		(days)	(ph cm−2 s−1)
56988.38	56999.64	12.01	$(1.1\pm 0.1)\,\times \,{10}^{-11}$
57000.39	57005.63	5.98	$(3.9\pm 1.3)\,\times \,{10}^{-12}$
57006.37	57009.61	3.99	$(1.6\pm 0.2)\,\times \,{10}^{-11}$
57010.36	57019.59	9.97	$(6.4\pm 1.0)\,\times \,{10}^{-12}$
57020.33	57020.58	1.00	$(2.9\pm 0.4)\,\times \,{10}^{-11}$
57021.33	57045.52	24.93	$(7.2\pm 0.7)\,\times \,{10}^{-12}$
57046.26	57086.40	40.89	$(3.8\pm 0.5)\,\times \,{10}^{-12}$
57087.15	57090.39	3.99	$(1.4\pm 0.2)\,\times \,{10}^{-11}$
57091.14	57143.25	52.85	$(6.9\pm 0.5)\,\times \,{10}^{-12}$
57144.00	57236.99	93.74	$(4.0\pm 0.4)\,\times \,{10}^{-12}$
57237.74	57239.98	2.99	$(2.2\pm 0.2)\,\times \,{10}^{-11}$
57240.73	57254.95	14.96	$(8.1\pm 0.9)\,\times \,{10}^{-12}$
57255.69	57275.89	20.94	$(1.3\pm 0.7)\,\times \,{10}^{-12}$
57276.63	57319.76	43.88	$(4.7\pm 0.5)\,\times \,{10}^{-12}$
57320.51	57368.63	48.87	$(4.8\pm 4.2)\,\times \,{10}^{-13}$
57369.38	57382.59	13.96	$(5.0\pm 0.8)\,\times \,{10}^{-12}$
57383.34	57387.58	5.98	$(1.5\pm 0.2)\,\times \,{10}^{-11}$
57389.32	57411.51	22.94	$(5.3\pm 0.6)\,\times \,{10}^{-12}$
57412.26	57496.28	83.77	$(1.5\pm 0.3)\,\times \,{10}^{-12}$

Download table as:  ASCIITypeset image

The spectral fit results, ${\rm{\Gamma }}=2.21\pm {0.14}_{\mathrm{stat}}\pm {0.20}_{\mathrm{sys}}$ and ${E}_{0}=5.4\pm {1.1}_{\mathrm{stat}}\pm {1.0}_{\mathrm{sys}}$ TeV, are consistent with spectral shapes previously observed (see, e.g., Albert et al. 2007b). If we compare to the range of VERITAS spectral fits as a function of flux state in Acciari et al. (2011a), we find that the average HAWC spectrum is closest to the parameters for the Mid-state level, ${{\rm{\Gamma }}}^{{\rm{M}}}=2.278\pm 0.037$ and ${E}_{0}^{{\rm{M}}}=4.36\pm 0.58$ TeV. A more detailed analysis of the HAWC spectral fits and a discussion of the absorption features of the EBL is beyond the scope of this paper. We will revisit this in a separate paper and take advantage of better energy estimation techniques for HAWC data that are currently under development to enhance the sensitivity to the curvature at the highest energies.

HAWC data can confirm and track the variability of Mrk 421 via daily flux measurements. By applying the Bayesian blocks algorithm, we identified 19 distinct flux states. The apparent substructure within some blocks in Figure 4 is likely to be due to flux variations on shorter timescales, which cannot be resolved by the present analysis, given the predetermined 5% false positive probability and the uncertainties of our measurements. As a stability check, we lowered the false positive condition from 5% to 10%, which led to the identification of only one additional block around MJD 57065.

In Figure 5, a histogram of all flux measurements from Mrk 421 highlights the spread of the observed flux states. We can compare this distribution to a function from Tluczykont et al. (2010) that was derived as a good fit to archival data, composed of the sum of a normal distribution (${f}_{\mathrm{Gauss}}$) around a low flux peak and a log-normal part (${f}_{\mathrm{LogN}}$) describing a tail to higher fluxes. The number of observations as a function of the flux x is given by

$$\begin{eqnarray}{f}_{T}(x) & = & {f}_{\mathrm{Gauss}}+{f}_{\mathrm{LogN}}\\ & = & \displaystyle \frac{{n}_{\mathrm{Gauss}}}{{\sigma }_{\mathrm{Gauss}}\sqrt{2\pi }}\exp \left(-\displaystyle \frac{{(x-{\mu }_{\mathrm{Gauss}})}^{2}}{2{\sigma }_{\ \mathrm{Gauss}}^{2}}\right)\\ & & +\displaystyle \frac{{n}_{\mathrm{LogN}}}{x{\sigma }_{\mathrm{LogN}}\sqrt{2\pi }}\exp \left(-\displaystyle \frac{{(\mathrm{ln}(x)-{\mu }_{\mathrm{LogN}})}^{2}}{2{\sigma }_{\mathrm{LogN}}^{2}}\right),\end{eqnarray} \tag{ 5 }$$

with the best-fit parameters ${n}_{\mathrm{Gauss}}=48.08$, ${\mu }_{\mathrm{Gauss}}=0.3285$, ${\sigma }_{\mathrm{Gauss}}=0.1137$, ${n}_{\mathrm{LogN}}=45.55$, ${\mu }_{\mathrm{LogN}}=0.1025$, and ${\sigma }_{\mathrm{LogN}}\,=1.022$. Here we follow the convention from the reference of measuring x in CU above 1 TeV but setting its unit to 1 in the formula. In order to account for the HAWC measurement uncertainties, we use a two-step process to define samples that each have 471 flux values, matching the size of the data set. First, we draw 471 random values according to the distribution in Equation (5) and use these value as centers and the standard deviation values from data as widths to define 471 normal distributions. Then, we draw one random value from each of these normal distributions to obtain a set of fluxes that reflects the uncertainties of the HAWC data. We average over 10,000 such samples to obtain the prediction for an HAWC measurement of this flux distribution. It is included in Figure 5 (blue, dashed line). We compare this expectation with the histogram of HAWC data via a Kolmogorov–Smirnov (KS) test and find a probability of 0.0008 that they arise from the same distribution. This value is stable under changes of the histogram binning and rescaling the HAWC fluxes within the systematic uncertainty of ±25% leads to a maximum KS probability value of 0.0046. Since we account for HAWC flux uncertainties in our sampling procedure, we also tested reducing ${\sigma }_{\mathrm{Gauss}}$ from Tluczykont et al. (2010) under the assumption that it is mostly reflecting measurement uncertainties in the fitted data, but found only smaller KS probabilities. Since Equation (5) is based on the fit to a compilation of measurements from many different instruments, it is hard to assess the systematic uncertainties of this parametrization. We have to consider that the large gaps in time coverage and a likely bias due to observations triggered by multiwavelength alerts for the public data in Tluczykont et al. (2010) can lead to a fit that does not well represent the average daily flux distribution for Mrk 421. We compare this function here for the first time with data from an unbiased, regular monitoring with a single, stable detector and conclude that these 17 months of HAWC observations cannot be well described by Equation (5). The current level of statistical uncertainties of the HAWC data prevents us from obtaining a stable fit of the six parameters from Equation (5) or testing if the tail of higher fluxes is indeed best described with a log-normal distribution, which could indicate an origin of variability from multiplicative processes. Increased statistics and the coverage of more high flux states with new HAWC data will provide the basis for obtaining a better functional description of Mrk 421 flux states.

Zoom In Zoom Out Reset image size

Since 2016 February, the daily flux measurements for both Mrk 421 and Mrk 501 have been automatically performed at the HAWC site immediately after the end of each transit. The preliminary analysis is based on the so-called online reconstruction, performed with only a few seconds' time lag on all recorded events and a preliminary calibration and data quality selection. Our threshold for issuing alerts about flaring states for both sources is a flux value equivalent to 3 CU in a single transit, which corresponds to a detection at ∼5 standard deviations for Mrk 421. In the 17 months of data included here, Mrk 421 surpassed this threshold during 11 transits, 2.3% of the time.

In Figure 6, we compare the HAWC TeV measurements to light curves from the Swift-BAT (15–50 keV)³⁷ monitoring instrument that provides very similar sampling and instrument duty cycle. The Swift-BAT data allow us to apply the same Bayesian blocks algorithm as used for HAWC data. The ratio of average error to mean flux value is larger for Swift-BAT light curves than for the HAWC data, but simulations show that the same Bayesian prior value (${\mathrm{ncp}}_{\mathrm{prior}}=6$) will guarantee the same false positive probability of 5%. This analysis identifies eight change points, i.e., nine blocks, in the X-ray data. The only Swift-BAT flux state that matches one of the HAWC flux states with less than 10 days' difference in start and end times is the lowest one.³⁸ None of the highest HAWC-measured flaring states are mirrored in the blocks for the X-ray light curve. We cannot exclude that the size of the statistical uncertainties hides correlated features at the day scale, considering that at least the Swift-BAT energy band seems to cover mostly a steeply falling part of the spectral energy distribution observed in the past (Abdo et al. 2011b). It is possible that this is due to insufficient overlap between the instruments' exposures during one day, since we established significant TeV variability within less than one sidereal day in Section 5.2. Similar cases of missing correlations for bright TeV flares have been observed before (see, e.g., Blazejowski et al. 2005; Acciari et al. 2011a). On the other hand, when we compare all daily averaged fluxes by calculating the Spearman rank correlation coefficient for the HAWC and the Swift-BAT data, we find a positive correlation of 0.341 ± 0.030. The probability for this result to occur for uncorrelated data sets of the same size is $3\,\times \,{10}^{-13}$ and we can thus qualitatively confirm previous observations of TeV-to-keV correlations based on (partially biased) IACT data (Blazejowski et al. 2005; Albert et al. 2007b; Fossati et al. 2008; Horan et al. 2009; Tluczykont et al. 2010).

Zoom In Zoom Out Reset image size

The only public notification about flaring activity for Mrk 421 that was sent during the period under investigation is an Astronomer's Telegram (Kapanadze 2015) about an increased X-ray flux observed with Swift-XRT (0.3–10 keV) between 2015 June 8 and June 16. Our Bayesian blocks analysis does not identify any significant flux state changes within 1.2 months around these dates.

Mrk 501 is a BL Lacertae type blazar that is similar to Mrk 421, given its distance of z = 0.033 (Mao 2011) and its classification as a high-peaked BL Lac object. It is the second extragalactic object that was discovered at TeV energies (Quinn et al. 1996). Various studies at TeV energies have shown different features of low flux state emission and extreme outbursts, for example, in Acciari et al. (2011b).

Our initial fit of the spectral shape uses the integrated 17 months of HAWC data. When we do not allow curvature in the spectral model, ${E}_{0}\to \infty$ in Equation (1), we obtain the best-fit values $F=(4.50\pm {0.28}_{\mathrm{stat}}\pm {2.25}_{\mathrm{sys}})$ × 10⁻¹² TeV⁻¹ cm⁻² s⁻¹ for normalization at 1 TeV and a photon index ${\rm{\Gamma }}=2.84\,\pm {0.04}_{\mathrm{stat}}\pm {0.20}_{\mathrm{sys}}$. This is consistent with results reported in Acciari et al. (2011b) and Abdo et al. (2011a). Leaving also the exponential cut-off free yields a normalization $F=(4.40\,\pm {0.60}_{\mathrm{stat}}\pm {2.20}_{\mathrm{sys}})\,\times \,{10}^{-12}$ TeV⁻¹ cm⁻² s⁻¹, a photon index ${\rm{\Gamma }}=1.60\pm {0.30}_{\mathrm{stat}}\pm {0.20}_{\mathrm{sys}}$, and an exponential cut-off value of ${E}_{0}=5.7\pm {1.6}_{\mathrm{stat}}\pm {1.0}_{\mathrm{sys}}$ TeV. The latter result is clearly preferred by ${\rm{\Delta }}\mathrm{TS}=48.64$ or 7.0 standard deviations over the power-law fit without a cut-off. Its significance compared to the background-only hypothesis is $\mathrm{TS}=610.49$ or 24.7 standard deviations. For the flux normalization fits performed to construct the flux light curve, we use the description with the cut-off and keep the index and cut-off parameters fixed at the HAWC-measured values.

The Mrk 501 flux light curve for all 1-transit intervals that have $\gt 50$% coverage with HAWC is shown in Figure 7. The photon flux is calculated as the analytical integration above 3 TeV, the optimal threshold value for Mrk 501 in order to minimize the systematic uncertainties of the flux measurement due to the fixed spectral assumption. For the 17 month period included here, we find an average flux of $(1.74\pm 0.08)\,\times {10}^{-12}$ ph cm⁻² s⁻¹ above 3 TeV.

Zoom In Zoom Out Reset image size

The result of the likelihood variability calculation for Mrk 501 is TS${}_{\mathrm{var}}=1115.4$, corresponding to a p-value $9.18\,\times {10}^{-48}$, and thus clearly establishes variability of the TeV emission measured here. The highest daily flux, $(1.67\pm 0.23)\,\times {10}^{-11}$ ph cm⁻² s⁻¹ above 3 TeV, was observed during the transit MJD 57251.94–57252.19 (2015 August 17 UTC 22:40 to 2015 August 18 UTC 4:40) with a pre-trial significance of 9.5 standard deviations compared to the null hypothesis. This is approximately a factor of 10 higher than the constant flux fit average and shows a variability timescale of less than one day, since the flux is higher by a factor of ∼4 compared to the previous transit and by a factor of ∼8 compared to the next transit.

In order to find significant flux state changes in this light curve, we applied the Bayesian blocks algorithm. Given the prior for 5% false positive probability, the algorithm identifies 13 change points. The amplitudes of periods between these change points are consistent with a constant flux and are shown as blue lines in Figure 7, with shaded bands indicating one standard deviation around the mean amplitude. These block positions and amplitudes are listed in Table 3.

Table 3.  HAWC Bayesian Blocks for Mrk 501

MJD Start	MJD Stop	Duration	Flux $\gt 3$ TeV
		(days)	(ph cm−2 s−1)
56989.66	57024.82	35.90	$(1.1\pm 0.3)\,\times \,{10}^{-12}$
57025.56	57033.79	8.98	$(7.0\pm 0.7)\,\times \,{10}^{-12}$
57034.54	57062.67	28.92	$(1.7\pm 0.3)\,\times \,{10}^{-12}$
57063.46	57076.67	13.96	$(5.7\pm 0.5)\,\times \,{10}^{-12}$
57077.42	57127.53	50.86	$(1.5\pm 0.3)\,\times \,{10}^{-12}$
57128.28	57129.53	1.99	$(1.4\pm 0.2)\,\times \,{10}^{-11}$
57130.28	57133.52	3.99	$(5.4\pm 1.0)\,\times \,{10}^{-12}$
57134.27	57251.20	117.68	$(7.6\pm 1.5)\,\times \,{10}^{-13}$
57251.94	57252.19	1.00	$(1.7\pm 0.2)\,\times \,{10}^{-11}$
57252.94	57394.80	142.61	$(4.5\pm 1.4)\,\times \,{10}^{-13}$
57395.55	57407.77	12.96	$(3.2\pm 0.5)\,\times \,{10}^{-12}$
57408.51	57483.56	75.79	$(1.1\pm 0.2)\,\times \,{10}^{-12}$
57484.31	57485.55	1.99	$(9.8\pm 1.5)\,\times \,{10}^{-12}$
57486.30	57497.52	10.97	$(2.1\pm 0.5)\,\times \,{10}^{-12}$

Download table as:  ASCIITypeset image

The integrated HAWC data for Mrk 501 is best described via a curved spectrum that we model with a photon index ${\rm{\Gamma }}=1.6$ and an exponential cut-off at E₀ = 5.7 TeV. In Figure 8, we compare this result to spectra measured by MAGIC (Ahnen et al. 2016), VERITAS and Whipple (Aliu et al. 2016), as well as HEGRA (Aharonian et al. 2001). We include only the values designated as low-state by these observatories since the HAWC light curve shows long periods of low activity that dominate our averaged measurement. The analyses of spectra during flaring states in the same publications show generally harder photon indices. The HAWC spectrum averaged over 17 months is consistent with these measurements within the statistical and systematic uncertainties as shown in Figure 8. The spectral curvature in our measurement manifests itself primarily outside the energy range covered by the IACT measurements. The HAWC energy range was determined from simulation as the central interval containing 90% of expected signal events for the best-fit spectrum. While the strong curvature in the spectrum of Mrk 501 prevents us from constraining the spectral shape above an energy of ∼15 TeV with the current analysis, HAWC is generally sensitive to gamma-ray energies up to ∼100 TeV, as discussed in Abeysekara et al. (2017a). The HEGRA analysis also obtains a better fit with an exponential cut-off than a pure power law. The HEGRA cut-off value, $5.1{({}_{-2.3}^{+7.8})}_{\mathrm{stat}}$ TeV, is consistent with the HAWC value. A cut-off in the energy spectrum can arise from gamma-ray absorption through the EBL (see, e.g., Domínguez et al. 2011) or originate in processes intrinsic to the source, for example, a limit to the energies of injected particles, changes in the Klein–Nishina scattering cross-section (Hillas 1999), or absorption through photon fields in the lower jet (Dermer & Schlickeiser 1994). The best-fit photon index ${\rm{\Gamma }}=1.6\pm {0.30}_{\mathrm{stat}}\pm {0.20}_{\mathrm{stat}}$ that we measure is hard compared to, for example, Mrk 421 but still greater than the lower limit of ∼1.5 for Fermi-acceleration in shocks (Malkov & Drury 2001). At lower energies, up to $\sim 300\,\mathrm{GeV}$, spectral hardening with photon indices down to ∼1.0 has been observed for Mrk 501 by Fermi-LAT (see, e.g., Shukla et al. 2016). We will provide a more detailed study of the spectral energy distribution and EBL absorption for Mrk 501 with HAWC data in a separate publication.

Zoom In Zoom Out Reset image size

The TeV light curve for Mrk 501 shows various flaring periods and we find 14 Bayesian blocks defining distinct flux states under a false detection prior of 5%. The three highest per-transit fluxes exceed the level of three times the Crab Nebula flux, corresponding to 0.6% of all observations included here. The last of these flares (2015 August 17) was also captured in tests of the HAWC real-time fast transient monitor (Weisgarber 2017) that resolves a sub-transit light curve of event rates. A histogram of all flux measurements is shown in Figure 9. The shape is generally similar to Figure 5, though with the peak and maximum fluxes shifted to lower CU values. As in the case of Mrk 421, the limited statistics currently prevent us from distinguishing between different functional descriptions of this distribution.

Zoom In Zoom Out Reset image size

The automated daily light-curve monitor that performs the per-transit light-curve analysis at the HAWC site, which has been operational since 2016 February and identified the increased flux state of Mrk 501 on 2016 April 6, is visible on the right side of Figure 7 at a level of ∼2.4 CU. We reported this flare in Sandoval et al. (2016), including the fact that the following day still showed a higher than average flux before returning to a lower state. The gamma-ray excess rates measured by FACT³⁹ above 750 GeV also show a rising trend before the HAWC alert.

The Swift-BAT (15–50 keV) data⁴⁰ for Mrk 501 have a very high ratio of average error to weighted mean flux (2.5), but we can apply the same Bayesian blocks analysis with prior value ${\mathrm{ncp}}_{\mathrm{prior}}=6$ as for HAWC data, keeping a 5% false positive probability. Only two change points, i.e., three distinct flux states are found. The resulting comparison between the three light curves is shown in Figure 10. The X-ray data reveal no day-scale light-curve features that are correlated with the TeV flares observed with HAWC. With the Bayesian blocks analysis, we find no short flaring periods in the Swift-BAT light curves that mirror the activity observed at TeV energies, but the large relative uncertainties (average error corresponds to 2.5 times the mean rate value) precludes us from obtaining a quantitative limit for this absence of correlation. We can calculate the Spearman rank correlation coefficient for all daily averaged fluxes and obtain 0.164 ± 0.032, with a probability of 10⁻³ to occur for an uncorrelated system. This positive correlation is qualitatively similar to results obtained previously (e.g., Krawczynski et al. 2000; Tluczykont et al. 2010) but is less significant than that observed for Mrk 421 in Section 5.4.

Zoom In Zoom Out Reset image size