A Survey of Active Galaxies at TeV Photon Energies with the HAWC Gamma-Ray Observatory

HAWC is an extensive air shower (EAS) array, sampling in detail secondary particles produced by primary cosmic rays in the upper atmosphere. HAWC data analysis can distinguish between hadronic and γ-ray-induced cascades through their different charge distributions at the ground. The observatory consists of a dense array of 300 large water Cherenkov detectors (WCDs) covering collectively a physical area larger than 22,000 m². Each WCD is a cylindrical tank of 7.3 m diameter and 5 m height, filled with 180 m³ of water and instrumented with four upward-facing photomultiplier tubes (PMTs) at its base. The signals from the 1200 PMT channels are brought to the data acquisition system, located near the center of the array, to be processed in real time. HAWC has been in full operation since its inauguration on 2015 March 20, after two years of gathering data with a partial array configuration (Abeysekara et al. 2016). Further details about the observatory can be found in Abeysekara et al. (2017b).

The HAWC array records about 25,000 events per second, the vast majority caused by hadronic cosmic rays. Each data record contains the particle arrival timing and deposited charge on each of the 1200 PMTs, which are used to locate the event in the sky and to perform photon/hadron discrimination. The analysis presented here follows the validation observation of the Crab Nebula, both in the γ-hadron cuts used and in partitioning the data in nine bins according to the fraction of channels hit, as indicated in Table 2 of Abeysekara et al. (2017b). The bin number provides a coarse measure of the primary energy, with an important overlap in the energy distributions of different bins due to the fluctuations inherent in the development of particle cascades. Lower bins relate to lower energies, and the spatial resolution improves with increasing bin number, with the detailed detector response depending on the spectrum of the source and its declination. The lowest bin used in this analysis, ${ \mathcal B }=1$, has peak sensitivity around 0.5 TeV for a source with a power-law spectrum of index 2.63 culminating at the zenith, as the Crab Nebula.

The HAWC sky surveys have a mean photon energy of about 7 TeV and one-year sensitivity between 50 and 100 mCrab (Abeysekara et al. 2017c). HAWC data analysis is based on computing the likelihood ratio of a source+background to a background-only model, given by the test statistic,

$$\begin{eqnarray}&&\mathrm{TS}=2\mathrm{ln}\left\{\displaystyle \frac{{ \mathcal L }(S+B)}{{ \mathcal L }(B)}\right\},\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal L }(A)$ is the likelihood of model A, given by the product of the probability density function computed at each point of the region of interest. The test statistics (Equation (1)) refers to the comparison of a background-only model (B) and a background+source model (S + B). Given a TS value, its statistical significance can be approximated by $s=\pm \sqrt{\mathrm{TS}}$, with the sign indicating an excess or deficit relative to the background. For all-sky surveys, like 2HWC, the analysis is performed by optimizing TS on every pixel of a N_side = 1024 HEALPix grid model of the observable sky (Górski et al. 2005). The source model generally consists of either a point-source or an extended-source hypothesis following a simple power-law spectrum of fixed spectral index, with free normalization. Joint normalization and spectral index optimizations are then performed to further characterize detected sources. The analysis presented here is performed similarly to that of the 2HWC, although on predefined sky locations and including the attenuation of TeV photons caused by their interaction with extragalactic background light. We use 1523 days of live data acquired between 2014 November 26 and 2019 June 3. The live data comprises 92.3% of the total time span. The data deficit is due mostly to quality cuts and run losses during bad weather.

The comparison of HAWC and Fermi-LAT data requires the consideration of the respective systematic uncertainties of each experiment. HAWC fluxes presented here have an estimated 15% systematic uncertainty (Section 2.2); those in the 3FHL catalog are quoted to have uncertainties of 9% in the 150–500 GeV band, and 15% in 0.5–2.0 TeV (Ajello et al. 2017).

The 3FHL catalog contains 1556 objects detected at photon energies between 10 GeV and 2 TeV in the first seven years of Fermi operations, from 2008 August 4 to 2015 August 2 (Ajello et al. 2017). Of the 3FHL entries, 79% are identified or associated with extragalactic objects, mostly BL Lac objects and flat-spectrum radio quasars (48% and 11% of the 3FHL, respectively). As defined in the different LAT catalogs, an association refers to the positional coincidence of the HE γ-ray source with an object having suitable properties, while an identification requires measuring correlated variability between the γ-ray source and its associated counterpart. These criteria result in 9% of the sources in 3FHL listed as identified and 78% as associated, while the remaining 13% are unassociated or unclassified. It is customary in Fermi catalogs to distinguish between identifications using uppercase letters, such as RDG for identified radio galaxies, and associations using lowercase letters, such as rdg for an association with a radio galaxy.

The 3FHL catalog partitions its nominal wide 10 GeV–2 TeV spectral interval into five bands. The majority of AGNs are detected at TS > 10 in the two lower-energy bands, 10–20 and 20–50 GeV. Spectral cutoffs at energies ≲50 GeV are common, as can be noticed in our sample: of the 138 AGNs, 50 (14) are detected above 5σ in the LAT 50–150 GeV (150–500 GeV) intermediate band(s). Furthermore, Mkn 421 and Mkn 501 are the only two AGNs studied here with TS > 25 in the 0.5–2.0 TeV LAT band, the spectral intersection with HAWC.

The 3FHL catalog assigns a flag to sources as follows: TeV = P, when reported at VHE energies; TeV = C for candidates for TeV detection; and TeV = N for nonreported and not candidates. TeV candidates are by definition sources undetected with VHE ground-based instruments whose LAT data satisfy three conditions: (i) significance s > 3 above 50 GeV; (ii) spectral index <3; (iii) integrated photon flux N(>50 GeV) > 10⁻¹¹ cm⁻² s⁻¹. In addition, the 3FHL catalog assesses source variability through the V_bayes parameter, the number of Bayes blocks needed to model the light curve. A source with V_bayes = 1 is consistent with a constant flux (Ajello et al. 2017).

We select 3FHL catalog sources identified or associated with AGNs with redshifts z ≤ 0.3 that culminate within 40° of the zenith as viewed from the HAWC site. For the selection of the follow-up sample we used the current version of the 3FHL catalog available at the Fermi Science Support Center.³³ The sample of 138 objects is summarized in Table 1. It contains 32 objects flagged TeV = P (positive VHE detections), and 32 TeV = C (candidates). The sample is grouped in five source classes, defined in the 3FHL, of distinct properties (Table 1):

• 1.  
  Starburst galaxies are the nearest and apparently less luminous AGNs in Fermi-LAT. While prone to host active nuclei, the prevailing γ-ray emission is dominated by cosmic rays produced in star formation processes. NGC 1068 is the only 3FHL starburst inside the decl. range of our selection, and the lowest redshift AGN in our sample (Figure 1). Intriguingly, it has recently been associated with one of the hot spots in the neutrino sky, as observed by the IceCube observatory (Aartsen et al. 2020).
• 2.  
  Radio galaxies (RDG) are the nearest extragalactic GeV sources dominated by an active nucleus. They appear up to three orders of magnitude more luminous than starbursts (Figure 1). Radio galaxies are attractive targets for HAWC because the relatively close distance translates in a reduced photon–photon attenuation, potentially allowing sampling of the far-infrared portion of the EBL through observations above ∼10–30 TeV. Six 3FHL catalog radio galaxies transit through the field of view of HAWC, with redshifts between z = 0.0042, for M87, and z = 0.029, for NGC 1218. Four of them have been claimed as VHE sources by IACT collaborations³⁴  (Rieger & Levinson 2018).
• 3.  
  BL Lacertae objects (BLL) constitute the majority of known GeV and VHE γ-ray emitters. They dominate the 3FHL catalog, in particular for redshifts z ≲ 0.7. As expected, BL Lac objects completely dominate our sample with 6 identifications and 111 associations. They span most of the z ≤ 0.3 range, led by Mkn 421 at z = 0.031, the nearest and brightest 3FHL AGN.
• 4.  
  Flat-spectrum radio quasars (FSRQ) are the most distant and seemingly luminous blazars. Our sample includes six such sources, with redshifts ranging from z = 0.158 to z = 0.222. Of these, only PKS 0736+017 has been reported as a TeV source, the nearest FSRQ claimed in the VHE range so far (H.E.S.S. Collaboration et al. 2020).
• 5.  
  Blazars candidates of uncertain type (bcu) are AGNs poorly characterized across the electromagnetic spectrum. The 3FHL catalog contains 290 bcus, a good fraction of them in the Southern sky, all with radio-loud associations, and 90% of them lacking redshift measurements. Only eight such objects satisfy our selection criteria. They are relatively near objects, between z = 0.036 and z = 0.128 (Figure 1). None is a VHE source, but three of them are flagged as TeV candidates.

Figure 1 shows the luminosity per solid angle unit, d_L(z) f_e, with f_e being the (10 GeV–1 TeV) energy flux from the 3FHL catalog, and d_L(z) the luminosity distance. The actual luminosity of each source depends on the unknown solid angle of emission.

Zoom In Zoom Out Reset image size

Prior to the analysis, we identified five objects in our sample close in the sky to bright 2HWC sources that could affect their analysis. We set a conservative distance threshold of 5^◦, equivalent to five times the 68% containment radius of ${ \mathcal B }=1$, the bin where contamination by a bright source nearby is more likely. As our sample excludes by construction low Galactic latitudes, with all selected sources located at $| b| \gt 5^\circ$, the only concerns were for:

• 1.  
  3FHL J0521.7+2112, located at 3.07 degrees from the Crab Nebula;
• 2.  
  3FHL J1041.7+3900, 3FHL J1100.3+4020, and 3FHL J1105.8+3944, located respectively at 4.51, 2.28, and 1.57 degrees from Markarian 421;
• 3.  
  3FHL J1652.7+4024, located at 0.68 degrees from Markarian 501.

Only 3FHL J0521.7+2112 is associated with a well-known VHE source, the BL Lac object VER J0521+211. The other four sources are flagged as TeV = N. These five objects were analyzed and tested later for contamination from the bright neighbor source, as detailed for VER J0521+211 in Section 5.3.2. From the respective tests we decided to exclude 3FHL J1105.8+3944 and 3FHL J1652.7+4024 from our sample.

We also revised the redshift measurements listed in the 3FHL catalog. Our sample is dominated by BL Lac objects, often with questionable redshifts. Redshifts from the 3FHL catalog were systematically collated with the SIMBAD and NED databases, with the Sloan Digital Sky Survey (SDSS) database consulted in some particular cases. We used as references the redshift surveys of Shaw et al. (2013) and Healey et al. (2008), and in particular the dedicated survey of TeV sources performed with the 10.4 m Gran Telescopio Canarias (GTC) by Paiano et al. (2017). While we decided to systematically use the redshifts listed in the 3FHL catalog, we list in Table 2 sources with disputed values, for future reference.

AGNs have been extensively studied with IACTs, with evidence of VHE γ-ray emission in almost a hundred of them (Wakely & Horan 2008; Madejski & Sikora 2016). These observations have shown that TeV flaring is an intrinsic characteristic of blazars and radio galaxies. IACT deep observations, able to reach down a few percent of the flux of the Crab Nebula in a single run, have identified high, medium, and low states of activity in several sources. However, continuous long-term coverage of the AGNs population cannot be performed with IACTs. The definition of the base level of AGN VHE emission is a pending task. Despite their lower instantaneous sensitivity, EAS arrays with efficient γ/hadron discrimination can quantify better time-averaged fluxes, integrated over long periods of time, while monitoring for flaring activity, as they continuously drift over large portions of the sky.

The first significant AGN detections at TeV energies were those of Mkn 421 and Mkn 501 (Punch et al. 1992; Quinn et al. 1996). These two Markarian galaxies have been extensively studied in the TeV range for more than two decades. They are highly variable but remain bright enough over long periods of time to have been detected by EAS arrays like the Tibet Air Shower Array and MILAGRO, both providing first unbiased views of their TeV emission on timescales of years (Bartoli et al. 2011; Abdo et al. 2014). In the last five years, HAWC has been performing an increasingly deeper monitoring of these two BL Lac objects (Abeysekara et al. 2017a). Measurements of their long-term averaged emission are presented in this paper.

IACTs have been able to go deeper and beyond these two well-known blazars by implementing sophisticated technologies to achieve lower energy thresholds, around or below 100 GeV, with improved sensitivities. This has permitted us to peer through the EBL horizon up to redshifts z ≲ 0.9 (Abeysekara et al. 2015; Ahnen et al. 2015), while sampling different types of extragalactic sources, such as:

• 1.  
  The starbursts galaxies M82 and NGC 253, detected by VERITAS and H.E.S.S. respectively, with fluxes below 1% of the Crab (Ohm 2016).
• 2.  
  Four Faranoff-Riley type I radio galaxies: Centaurus A (too Southern for HAWC); NGC 1275, the most prominent galaxy of the massive Perseus cluster; 3C 264, newly reported as a VHE source; and M87, the massive central elliptical galaxy in the nearby Virgo cluster. We note the ambiguous classification of IC 310 and PKS 0625–35, referred as RDG in the 3FHL, and as unknown type of AGN in TeVCat (Wakely & Horan 2008) and by Rieger & Levinson (2018).
• 3.  
  At least 60 BL Lacertae objects as VHE γ-ray sources (Wakely & Horan 2008; Ajello et al. 2017), mostly high-frequency peaked BL Lac objects (HBL; 50 sources), plus a few intermediate-frequency peaked BLL (IBL; 8 sources, including BL Lac, W Comae, and VER J0521+211, in our sample); and only two low-frequency peak BL Lacs (LBL). While most of the IACT observations have been reported with sub-TeV thresholds, we note the observations of H 1426+428 by VERITAS and HEGRA at energies between 0.25 and 2.5 TeV, reporting fluxes from 3% to 10% of the Crab (Aharonian et al. 2002; Petry et al. 2002).
• 4.  
  Seven FSRQs, of which PKS 0736+017 is the only FSRQ detected at VHE with z < 0.3 and entering the field of view of HAWC. This object was found in a flaring state at a flux level of 100 mCrab between 100 and 300 GeV, and showing a steep spectrum (H.E.S.S. Collaboration et al. 2020).

We calculated the expected integral photon flux above 0.5 TeV by extrapolating the spectral models and parameters listed in the 3FHL. Although in HAWC a distinction is made between the intrinsic spectrum of a source (as emitted) and the related observable parameters (attenuated by the EBL), Fermi-LAT spectral models do not need to make this distinction. The difference is minor in most of the LAT energy regime, with only a moderate increase in spectral indices for sources with z ≳ 1 (Ajello et al. 2017). Therefore, we added the attenuation effect of the EBL to the LAT spectral models. The extrapolated fluxes as a function of redshift are shown in Figure 2. Note the uncertainties in the flux extrapolations due to the uncertainties in the fit parameters, specially as they are propagated more than one order of magnitude above the 3FHL pivot energies. The horizontal red line corresponds to a photon flux of 3% of the Crab Nebula, indicative of potential HAWC detectability within the current data.  Mkn 421 and Mkn 501 stand clearly at about an order of magnitude above the red line, followed by 1ES 2344+514 and I Zw 187 (1ES1727+502), two of the nearest BL Lac objects known. The radio galaxies M87 and IC 301 are also above the flux level of interest, although with fairly large uncertainties in the extrapolation. Of the sources potentially detectable, only TXS 0210+515 is tagged as undetected in the VHE range in the 3FHL. However, we note its recent detection with MAGIC (Acciari et al. 2020).

Zoom In Zoom Out Reset image size

We also computed photon flux extrapolations of VHE sources within TeVCat, using simple power laws and approximating the EBL attenuation by a direct exponential in redshift (Section 3). The predictions from this second set are also rather uncertain, because IACTs observations are by nature short and sparse, reflecting different activity states, and spectra are not always fitted above 1 TeV. Five candidates stood out with extrapolated fluxes N(>0.5 TeV) ≥ 20 mCrab, four of them in common with the 3FHL extrapolation: M87, Mkn 421, Mkn 501, H 1426+428, and 1ES 2344+514. Results for the sources named in this section are shown and discussed in Section 5.

The default spectral model for the AGN search was an index α = 2.5, as used in the 3HWC (Albert et al. 2020). However, note that here the index refers to the intrinsic spectrum and that the observed spectra will be softer. While we know from previous analyses that this index can represent fairly well the data from both Mkn 421 and Mkn 501 (Coutiño de Leon et al. 2019), photon indices in extragalactic 3FHL sources peak in the interval between 2.0 and 2.5 (Ajello et al. 2017).

5.1.1. Overall Results

The results for the complete AGN sample are shown in Table 3. The left-hand side of Figure 3 shows the histogram of significances for the sample, also displayed as the first entry of Table 4. Mkn 421 is detected with a significance $\sqrt{\mathrm{TS}}=+64.6$, and Mkn 501 with $\sqrt{\mathrm{TS}}=+16.6$. These two high-significance detections drive up the statistics of the complete sample to an overall 9σ deviation from the null hypothesis. When removing Mkn 421 and Mkn 501, the joint significance drops to a p-value of 2.2% (Figure 3, right-hand side).

Zoom In Zoom Out Reset image size

Table 3.  HAWC Power-law Fits and Significances for the Revised Sample of 136 AGNs from the 3FHL Catalog

3FHL Source	Counterpart	Class	Redshift	TS	$\pm \sqrt{\mathrm{TS}}$	K ± ΔK	K2σ
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
J0007.9+4711	RX J0007.9+4711	bll	0.2800	−4.49	−2.12	−31.9 ± 15.1	9.16
J0013.8–1855	RBS 0030	bll	0.0949	−0.10	−0.32	−2.20 ± 6.81	11.5
J0018.6+2946	RBS 0042	bll	0.1000	+0.01	+0.11	+0.10 ± 0.93	1.97
J0037.8+1239	NVSS J003750+123818	bll	0.0890	+0.47	+0.69	+0.52 ± 0.76	2.03
J0047.9+3947	B3 0045+395	bll	0.2520	−0.18	−0.42	−2.68 ± 6.38	10.3
J0056.3–0936	TXS 0053–098	bll	0.1031	−0.21	−0.46	−1.33 ± 2.88	4.53
J0059.3–0152	1RXS J005916.3–015030	bll	0.1440	+0.95	+0.98	+2.61 ± 2.67	7.93
J0112.1+2245	S2 0109+22	BLL	0.2650	+0.33	+0.57	+1.99 ± 3.47	8.88
J0123.0+3422	1ES 0120+340	bll	0.2720	+6.53	+2.56	+13.1 ± 5.13	23.4
J0131.1+5546	TXS 0128+554	bcu	0.0365	−0.85	−0.92	−1.61 ± 1.75	2.10
J0152.6+0147	PMN J0152+0146	bll	0.0800	−0.89	−0.94	−0.89 ± 0.94	1.10
J0159.5+1047	RX J0159.5+1047	bll	0.1950	+1.09	+1.04	+2.52 ± 2.41	7.31
J0211.2+1051	MG1 J021114+1051	BLL	0.2000	−0.58	−0.76	−1.92 ± 2.51	3.26
J0214.5+5145	TXS 0210+515	bll	0.0490	+2.63	+1.62	+2.43 ± 1.51	5.43
J0216.4+2315	RBS 0298	bll	0.2880	−0.08	−0.29	−1.14 ± 3.97	6.74
J0217.1+0836	ZS 0214+083	bll	0.0850	+4.58	+2.14	+1.70 ± 0.80	3.31
J0219.1–1723	1RXS J021905.8–172503	bll	0.1287	+0.01	+0.10	+1.0 ± 9.9	20.9
J0232.8+2017	1ES 0229+200	bll	0.1400	−0.28	−0.53	−0.70 ± 1.32	1.98
J0242.7–0002	NGC 1068	sbg	0.0038	+1.45	+1.21	+0.24 ± 0.20	0.65
J0308.4+0408	NGC 1218	rdg	0.0288	+0.24	+0.49	+0.17 ± 0.35	0.87
J0312.8+3614	V Zw 326	bll	0.0710	+0.48	+0.69	+0.53 ± 0.77	2.07
J0316.6+4120	IC 310	RDG	0.0189	+0.74	+0.86	+0.31 ± 0.36	1.01
J0319.8+1845	RBS 0413	bll	0.1900	+0.05	+0.22	+0.45 ± 2.08	4.64
J0319.8+4130	NGC 1275	RDG	0.0176	+0.36	+0.60	+0.21 ± 0.35	0.91
J0325.6–1646	RBS 0421	bll	0.2910	+0.07	+0.27	+12.6 ± 46.4	107
J0326.3+0226	1H 0323+022	bll	0.1470	−0.22	−0.47	−1.02 ± 2.16	3.35
J0334.3+3920	4C+39.12	rdg	0.0203	−0.91	−0.96	−0.33 ± 0.34	0.40
J0336.4–0348	1RXS J033623.3–034727	bll	0.1618	+0.19	+0.44	+1.64 ± 3.75	9.10
J0339.2–1736	PKS 0336–177	bcu	0.0656	+0.27	+0.52	+1.83 ± 3.55	8.97
J0349.3–1159	1ES 0347–121	bll	0.1850	+0.74	+0.86	+9.00 ± 10.5	30.0
J0416.8+0105	1ES 0414+009	bll	0.2870	+1.71	+1.31	+9.16 ± 7.00	22.9
J0424.7+0036	PKS 0422+00	bll	0.2680	+1.60	+1.27	+8.17 ± 6.46	21.1
J0521.7+2112	TXS 0518+211	bll	0.1080	+9.49	+3.08	+2.85 ± 0.93	4.72
J0602.0+5316	GB6 J0601+5315	bcu	0.0520	−0.17	−0.42	−0.79 ± 1.88	3.01
J0617.6–1715	TXS 0615–172	bll	0.0980	+1.51	+1.23	+7.46 ± 6.07	19.7
J0648.7+1517	RX J0648.7+1516	bll	0.1790	+0.26	+0.51	+1.00 ± 1.97	4.92
J0650.7+2503	1ES 0647+250	bll	0.2030	+0.01	+0.10	+0.24 ± 2.43	5.12
J0656.2+4235	4C +42.22	bll	0.0590	−0.65	−0.81	−0.71 ± 0.88	1.13
J0725.8–0056	PKS 0723–008	bcu	0.1270	−2.76	−1.66	−3.48 ± 2.09	1.58
J0730.4+3307	1RXS J073026.0+330727	bll	0.1120	+0.96	+0.98	+1.21 ± 1.23	3.66
J0739.3+0137	PKS 0736+01	fsrq	0.1910	−0.96	−0.98	−3.40 ± 3.47	3.93
J0753.1+5354	4C +54.15	bll	0.2000	+0.78	+0.88	+15.2 ± 17.3	49.6
J0757.1+0957	PKS 0754+100	bll	0.2660	−0.12	−0.34	−1.39 ± 4.04	6.80
J0809.7+3457	B2 0806+35	bll	0.0830	+0.09	+0.30	+0.26 ± 0.87	2.03
J0809.8+5218	1ES 0806+524	bll	0.1371	+0.27	+0.52	+3.59 ± 6.9	17.3
J0816.4–1311	PMN J0816–1311	bll	0.0460	−4.67	−2.16	−3.23 ± 1.49	0.89
J0816.4+5739	SBS 0812+578	bll	0.2940	−0.55	−0.74	−2.32 ± 3.14	4.16
J0816.9+2050	SDSS J081649.78+205106.4	bll	0.0583	−2.72	−1.65	−0.77 ± 0.47	0.36
J0828.3+4153	GB6 B0824+4203	bll	0.2262	+0.84	+0.92	+5.80 ± 6.33	18.4
J0831.8+0429	PKS 0829+046	bll	0.1738	+2.00	+1.41	+3.62 ± 2.56	8.72
J0847.2+1134	RX J0847.1+1133	bll	0.1982	+0.09	+0.29	+0.72 ± 2.45	5.65
J0850.6+3454	RX J0850.5+3455	bll	0.1450	+1.43	+1.20	+2.33 ± 1.95	6.20
J0908.9+2311	RX J0908.9+2311	bll	0.2230	−2.15	−1.47	−3.98 ± 2.71	2.32
J0912.4+1555	SDSS J091230.61+155528.0	bll	0.2120	−0.47	−0.69	−1.73 ± 2.52	3.54
J0930.4+4952	1ES 0927+500	bll	0.1867	−4.58	−2.14	−20.4 ± 9.53	5.60
J1015.0+4926	1H 1013+498	bll	0.2120	+0.03	+0.19	+2.14 ± 11.5	25.1
J1027.0–1749	1RXS J102658.5–174905	bll	0.2670	+0.10	+0.32	+14.6 ± 45.5	105
J1041.7+3900	B3 1038+392	bll	0.2084	−0.51	−0.71	−3.15 ± 4.43	6.01
J1053.6+4930	GB6 J1053+4930	bll	0.1404	+0.31	+0.56	+2.97 ± 5.33	13.8
J1058.6+5628	TXS 1055+567	BLL	0.1433	−1.62	−1.27	−15.5 ± 12.2	11.7
J1100.3+4020	RX J1100.3+4019	bll	0.2250	+2.67	+1.63	+9.08 ± 5.56	20.3
J1104.4+3812	Mkn 421	BLL	0.0310	+4166.97	+64.55	+29.5 ± 0.5	30.6
J1117.0+2014	RBS 0958	bll	0.1380	+5.28	+2.30	+3.00 ± 1.31	5.63
J1120.8+4212	RBS 0970	bll	0.1240	+1.42	+1.19	+2.73 ± 2.29	7.29
J1125.9–0743	1RXS J112551.6–074219	bll	0.2790	−2.74	−1.65	−23.2 ± 14.0	10.8
J1136.8+2549	RX J1136.8+2551	bll	0.1560	−1.13	−1.06	−1.74 ± 1.64	1.75
J1140.5+1528	NVSS J114023+152808	bll	0.2443	−0.59	−0.77	−2.41 ± 3.13	4.00
J1142.0+1546	MG1 J114208+1547	bll	0.2990	−0.93	−0.96	−4.08 ± 4.24	4.94
J1145.0+1935	3C 264	rdg	0.0216	+3.61	+1.90	+0.44 ± 0.24	0.92
J1150.3+2418	OM 280	bll	0.2000	+7.15	+2.67	+6.24 ± 2.33	10.9
J1154.1–0010	1RXS J115404.9–001008	bll	0.2535	−0.28	−0.53	−3.27 ± 6.19	9.46
J1204.2–0709	1RXS J120417.0–070959	bll	0.1850	−2.30	−1.52	−9.58 ± 6.31	5.09
J1217.9+3006	1ES 1215+303	bll	0.1300	+11.36	+3.37	+4.64 ± 1.38	7.39
J1219.7–0312	1RXS J121946.0–031419	bll	0.2988	−0.94	−0.97	−10.0 ± 10.3	12.4
J1221.3+3010	PG 1218+304	bll	0.1837	+5.02	+2.24	+5.23 ± 2.34	9.92
J1221.5+2813	W Comae	bll	0.1029	+6.03	+2.45	+2.33 ± 0.95	4.24
J1224.4+2436	MS 1221.8+2452	bll	0.2187	+0.55	+0.74	+1.99 ± 2.68	7.32
J1229.2+0201	3C 273	FSRQ	0.1583	−1.97	−1.40	−3.48 ± 2.48	2.15
J1230.2+2517	ON 246	bll	0.1350	+0.43	+0.66	+0.86 ± 1.31	3.46
J1230.8+1223	M87	rdg	0.0042	+12.93	+3.60	+0.56 ± 0.16	0.88
J1231.4+1422	GB6 J1231+1421	bll	0.2559	−0.09	−0.30	−1.01 ± 3.40	5.82
J1231.7+2847	B2 1229+29	bll	0.2360	+0.00	+0.05	+0.18 ± 3.30	6.75
J1253.7+0328	MG1 J125348+0326	bll	0.0657	−3.80	−1.95	−1.36 ± 0.70	0.46
J1256.2–1146	PMN J1256–1146	bcu	0.0579	+0.26	+0.51	+0.84 ± 1.63	4.11
J1310.3–1158	TXS 1307–117	bll	0.1400	−3.50	−1.87	−11.4 ± 6.10	4.11
J1341.2+3959	RBS 1302	bll	0.1715	−0.97	−0.99	−3.35 ± 3.40	3.92
J1402.6+1559	MC 1400+162	bll	0.2440	+1.78	+1.33	+4.14 ± 3.10	10.3
J1411.8+5249	SBS 1410+530	bcu	0.0765	+0.24	+0.49	+1.40 ± 2.84	7.12
J1418.0+2543	1E 1415.6+2557	bll	0.2363	+0.27	+0.52	+1.60 ± 3.07	7.78
J1419.4+0444	SDSS J141927.49+044513.7	bll	0.1430	+0.16	+0.40	+0.74 ± 1.85	4.50
J1419.7+5423	OQ 530	bll	0.1525	+0.30	+0.55	+5.94 ± 10.8	27.6
J1428.5+4240	H 1426+428	bll	0.1292	+1.58	+1.26	+3.18 ± 2.54	8.18
J1436.9+5639	RBS 1409	bll	0.1500	+2.38	+1.54	+21.0 ± 13.6	48.1
J1442.8+1200	1ES 1440+122	bll	0.1631	+3.53	+1.88	+3.37 ± 1.80	6.95
J1449.5+2745	B2.2 1447+27	bll	0.2272	+0.24	+0.49	+1.49 ± 3.03	7.46
J1500.9+2238	MS 1458.8+2249	bll	0.2350	+3.06	+1.75	+5.10 ± 2.91	11.0
J1508.7+2708	RBS 1467	bll	0.2700	+0.15	+0.38	+1.48 ± 3.89	9.29
J1512.2+0203	PKS 1509+022	fsrq	0.2195	−1.84	−1.36	−5.73 ± 4.23	3.81
J1518.5+4044	GB6 J1518+4045	bll	0.0652	+0.91	+0.96	+0.83 ± 0.87	2.57
J1531.9+3016	RX J1531.9+3016	bll	0.0653	+0.03	+0.18	+0.10 ± 0.57	1.25
J1543.6+0452	CGCG 050–083	bcu	0.0400	−5.82	−2.41	−1.04 ± 0.43	0.24
J1554.2+2010	1ES 1552+203	bll	0.2223	+1.76	+1.33	+3.50 ± 2.64	8.81
J1603.8+1103	MG1 J160340+1106	bll	0.1430	−0.17	−0.41	−0.61 ± 1.50	2.41
J1615.4+4711	TXS 1614+473	fsrq	0.1987	−1.08	−1.04	−8.29 ± 7.99	8.82
J1643.5–0646	NVSS J164328–064619	bcu	0.0820	−0.00	−0.01	−0.02 ± 1.61	3.23
J1647.6+4950	SBS 1646+499	bll	0.0475	+0.10	+0.32	+0.38 ± 1.21	2.82
J1653.8+3945	Mkn 501	BLL	0.0330	+276.97	+16.64	+7.74 ± 0.49	8.72
J1719.2+1745	PKS 1717+17	bll	0.1370	−0.32	−0.56	−0.72 ± 1.28	1.90
J1725.4+5851	7C 1724+5854	bll	0.2970	−1.94	−1.39	−99.0 ± 71.0	63.5
J1728.3+5013	I Zw 187	bll	0.0550	−0.05	−0.22	−0.32 ± 1.44	2.57
J1730.8+3715	GB6 J1730+3714	bll	0.2040	−0.21	−0.45	−1.71 ± 3.78	5.86
J1744.0+1935	S3 1741+19	bll	0.0830	+1.62	+1.27	+0.83 ± 0.66	2.17
J1745.6+3950	B2 1743+39C	bll	0.2670	+0.04	+0.21	+1.47 ± 7.11	15.5
J1813.5+3144	B2 1811+31	bll	0.1170	−1.68	−1.30	−1.59 ± 1.23	1.14
J1917.7–1921	1H 1914–194	bll	0.1370	+0.01	+0.12	+1.64 ± 13.6	28.9
J2000.4–1327	NVSS J200042–132532	fsrq	0.2220	−0.76	−0.87	−15.2 ± 17.4	21.1
J2014.4–0047	PMN J2014–0047	bll	0.2310	+1.23	+1.11	+6.12 ± 5.51	17.1
J2039.4+5219	1ES 2037+521	bll	0.0540	−0.94	−0.97	−1.70 ± 1.75	2.04
J2042.0+2428	MG2 J204208+2426	bll	0.1040	+0.58	+0.76	+0.67 ± 0.89	2.44
J2055.0+0014	RGB J2054+002	bll	0.1508	−0.26	−0.51	−1.27 ± 2.51	3.78
J2108.8–0251	TXS 2106–030	bll	0.1490	−0.55	−0.74	−2.23 ± 3.01	4.03
J2143.5+1742	OX 169	fsrq	0.2110	+0.18	+0.42	+1.03 ± 2.43	5.82
J2145.8+0718	MS 2143.4+0704	bll	0.2350	−1.02	−1.01	−3.65 ± 3.61	4.14
J2150.2–1412	TXS 2147–144	bll	0.2290	+0.01	+0.10	+2.02 ± 20.2	42.7
J2202.7+4216	BL Lacertae	BLL	0.0690	−0.25	−0.50	−0.50 ± 1.0	1.53
J2250.0+3825	B3 2247+381	bll	0.1190	−0.27	−0.52	−0.86 ± 1.67	2.54
J2252.0+4031	MITG J2252+4030	bll	0.2290	−0.78	−0.88	−5.03 ± 5.71	6.92
J2314.0+1445	RGB J2313+147	bll	0.1625	+1.37	+1.17	+1.97 ± 1.68	5.29
J2322.6+3436	TXS 2320+343	bll	0.0980	+0.01	+0.09	+0.09 ± 1.06	2.23
J2323.8+4210	1ES 2321+419	bll	0.0590	−3.79	−1.95	−1.64 ± 0.84	0.57
J2329.2+3755	NVSS J232914+375414	bll	0.2640	+0.04	+0.19	+1.17 ± 6.03	13.3
J2338.9+2123	RX J2338.8+2124	bll	0.2910	+2.91	+1.71	+6.76 ± 3.97	14.7
J2346.6+0705	TXS 2344+068	bll	0.1720	+2.24	+1.50	+3.33 ± 2.23	7.82
J2347.0+5142	1ES 2344+514	bll	0.0440	+2.09	+1.45	+1.93 ± 1.34	4.62
J2356.2+4035	NVSS J235612+403648	bll	0.1310	+0.01	+0.10	+0.22 ± 2.22	4.64
J2359.3–2049	TXS 2356–210	bll	0.0960	−0.06	−0.23	−1.97 ± 8.42	15.0

Note. Column (3) uses the identification (uppercase) or association (lowercase) class as given in the 3FHL (described in Section 4.1). Column (4) shows the redshift entry in the 3FHL, used in our analysis. Columns (5) and (6) show the test statistics (TS) and significances ($\pm \sqrt{\mathrm{TS}}$) estimated, allowing fluxes to be positive or negative. Column (7) shows that K is the power-law normalization at 1 TeV in units of 10⁻¹² TeV⁻¹ cm⁻² s⁻¹, and column (8), K_2σ, its corresponding 2σ upper limit in the same units, computed following Feldman & Cousins (1998).

Download table as:  ASCIITypeset images: 1 2 3

Table 4.  Average Significances for Different Groups of Sources

Source	Number	Significances		p-value
Group	of Objects	Mean	Std. dev
	N	μ(s)	σs	$P\left[\mu (s)\gt x\right]$
(1)	(2)	(3)	(4)	(5)
All	136	+0.769	5.762	5.58 × 10−19
All–Mkn	134	+0.175	1.209	0.022
Starburst	1	+1.21	⋯	0.113
Radio galaxies	6	+1.082	1.403	4.03 × 10−3
BL Lacs–Mkn	113	+0.220	1.179	9.75 × 10−3
FSRQ + fsrq	6	−0.872	0.609	0.984
bcu	8	−0.487	1.028	0.916
TeV = P (clean)	30	+0.768	1.228	1.29 × 10−5
TeV = C (clean)	32	+0.124	1.152	0.241
TeV = N (clean)	72	−0.050	1.142	0.665

Download table as:  ASCIITypeset image

Three more objects showed marginally significant test statistic TS > 9: M87 (TS = 12.9); 1ES 1215+303 (TS =11.4); and VER J0521+211 (TS = 9.5). We note that the random probability of having three TS > 9 values in 134 trials, once Mkn 421 and Mkn 501 are excluded, is 8.5 × 10⁻⁴. These three tentative sources have favorable declinations for HAWC, $\left|\delta -19^\circ \right|\lesssim 10^\circ$, allowing for transits of more than 6 hr. Spectral fits for the five TS > 9 sources are presented in Section 5.3.

5.1.2. Upper Limits and Sensitivity

The K_2σ upper limits span a large range of values: from 2.4 × 10⁻¹³ TeV⁻¹ cm⁻² s⁻¹ (3FHL J1543.6+0452 at z =0.040) to 10⁻¹⁰ TeV⁻¹ cm⁻² s⁻¹ (3FHL J0325.6–1646 at z =0.291). Upper limits are sensitive to both the declination and redshift of the source, as shown in Figure 4. We performed a fit to the set of K_2σ values that assumes a Gaussian dependence in decl. and exponential in redshift. Three parameters were computed: the normalization limit at the reference point δ = 19° and z = 0, the angular width in decl., and the redshift exponential scale. The fit obtained is

$$\begin{eqnarray}\,\begin{array}{rcl}{\mathrm{log}}_{10}{K}_{2\sigma } & = & -12.33+\displaystyle \frac{1}{2\mathrm{ln}10}{\left(\displaystyle \frac{\delta -19^\circ }{18\buildrel{\circ}\over{.} 93}\right)}^{2}\\ & & +\ \displaystyle \frac{1}{\mathrm{ln}10}\left(\displaystyle \frac{z}{0.089}\right).\end{array}\,\end{eqnarray} \tag{ 6 }$$

This fit excludes the five TS > 9 sources. Its correlation coefficient is r = +0.870 and the dispersion with the data is 0.241 dex. The value for the fit of K_2σ at z = 0, δ = 19° corresponds to 14 mCrab, a factor of two lower than the predefined depth of the survey. This would apply to near and optimally located sources only. The sensitivity of the survey degrades with a Gaussian angle of 18.9 degrees, leading to a 11% response at 40° from zenith. The dependence of the upper-limit normalizations with redshift gives z_h = 0.089, matching that of an observed power-law spectrum integrated from E₀ = 0.63 TeV, close to the initial assumption E₀ = 0.5 TeV (Section 3).

Zoom In Zoom Out Reset image size

5.1.3. Comparison with LAT Extrapolations

In Figure 5 we compare photon fluxes, N(>0.5 TeV), computed from the power-law fit (Equation (5)) with the extrapolations of LAT spectra. HAWC photon fluxes relate to the normalizations K through the integration of the differential spectra (Equations (4) and (5)),

$$\begin{eqnarray}&&{N}_{\mathrm{obs}}(\gt 0.5\,\mathrm{TeV})=\displaystyle \frac{4\sqrt{2}}{3}\,K\cdot 1\,\mathrm{TeV}\cdot \,{e}^{-z/{z}_{h}}.\end{eqnarray} \tag{ 7 }$$

The green (TeV = P) dots shown in Figure 5 are for the objects with TS > 9, while the rest are shown with the respective HAWC upper limit, using K_2σ in Equation (7).  Thirteen objects have HAWC limits below the LAT extrapolations, although compatible within the uncertainties. These include IC 310, I Zw 187 (1ES 1727+502), B3 2247+381, and 1ES 2344+514, identified in Figure 2.

Zoom In Zoom Out Reset image size

5.1.4. Source Classes

Table 4 presents average significances for distinct groups of sources. While GeV emission is known to occur in these sources, 80% of the objects in our sample have not been detected in the VHE range. Furthermore, HAWC is testing long-term average emission, in contrast to the relatively brief IACTs observations. The group statistics, reinforced by the HAWC measurements as a function of redshift for different source classes (Figures 6 and 7), are briefly summarized as follows:

• 1.  
  Except for the one starburst, radio galaxies constitute the nearest class of sources in our sample (Figure 6). All six are nearer than z = 0.03, and happen to be at favorable declinations, from δ = +4° (NGC 1278) to δ = +41° (IC 310 and NGC 1275). Four of them are known VHE sources, including 3C 264, reported after the publication of the 3FHL, which shows a +1.9σ excess in the HAWC data. As shown in Figure 6, their bounds are at the N(>0.5 TeV) ∼ 10⁻¹² cm⁻² s⁻¹ level. The mean significance for this group has a p-value of 0.4%.
• 2.  
  The blazars candidates of uncertain type studied here are at redshifts intermediate between those of radio galaxies and FSRQs. None of the bcus studied here has been reported in the VHE regime. The HAWC upper limits shown in Figure 6 are not particularly constraining, owing to the unfavorable declinations of these sources, seven out of eight outside +0° ≤ δ ≤ +50°.
• 3.  
  Flat-spectrum radio quasars constitute the most distant class of source in our sample. Five out of six have not been detected as VHE sources, and do not comply with the LAT requirement for the TeV = C flag. They appear mostly as underfluctuations in the HAWC data. Figure 6 shows the upper limits on photon fluxes, most of them below the 30 mCrab level.
• 4.  
  BL Lacertae objects constitute the clear majority of our sample, and of VHE sources. Still, 89 of the 117 BL Lac objects studied here have not been reported in the VHE regime. Figure 7 shows the HAWC fluxes for Mkn 421, Mkn 501, VER J0521+211, and 1ES 1215+303, at about (2–3) × 10⁻¹² cm⁻² s⁻¹ for the two weaker cases.

The TeV = P flagged sources appear with a p-value at the 10⁻⁵ level, excluding Mkn 421 and Mkn 501, as expected for a sub-threshold persistent TeV emission in these known sub-TeV emitters. The TeV = C candidates and TeV = N groups do not provide any collective hint of emission. The p-values quoted do not account for the number of trials used in testing different groups; they are quoted as indicative of the potential presence of persistent TeV emission at levels ≲10⁻¹² cm⁻² s⁻¹.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We computed optimized spectra for the TS > 9 sources, fitting together normalizations and spectral indices. The fits are summarized in Figure 8, together with the systematic uncertainties. These have been quantified as 15% in K, the 1 TeV normalization, and 5% in α, the spectral index. We discuss here the two high-significance sources, Mkn 421 and Mkn 501.

Zoom In Zoom Out Reset image size

5.2.1. Markarian 421

Markarian 421 is the brightest persistent extragalactic object in the TeV sky. For the default α = 2.5 search, we computed TS = 4167 for K = (29.5 ± 0.5_stat ± 4.4_syst) ×10⁻¹² TeV⁻¹ cm⁻² s⁻¹. The optimized power-law fit for the intrinsic spectrum resulted in

$$\begin{eqnarray}&&\,\begin{array}{l}\displaystyle \frac{{dN}}{{dE}}=(33.0\pm {0.6}_{\mathrm{stat}}\pm {4.9}_{\mathrm{syst}})\times {10}^{-12}\\ \ \ \times \ {\left(\displaystyle \frac{E}{1\mathrm{TeV}}\right)}^{-2.63\pm {0.02}_{\mathrm{stat}}\pm {0.13}_{\mathrm{syst}}}\,{\mathrm{TeV}}^{-1}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1},\end{array}\,\end{eqnarray} \tag{ 8 }$$

with TS = 4193, i.e., an increase in the test statistic ΔTS = 26 relative to α fixed at 2.5. The integration to the observed photon flux, using z_h = 0.116 for α = 2.63 in Equation (4), gives N_obs = (48 ± 8) × 10⁻¹² cm⁻² s⁻¹, somewhat larger than the α = 2.5 estimate shown in Figure 5. The intrinsic energy flux, f_E = (132 ± 22) × 10⁻¹² erg cm⁻² s⁻¹, translates into a luminosity per solid angle of f_Ed_L(z)² = L(>0.5 TeV)/ΔΩ = (2.1 ± 0.4) × 10⁴³ erg s⁻¹ sr⁻¹, about one-third (32%) of the (10 GeV–1 TeV) luminosity per steradian inferred from the Fermi-LAT measurements (Figure 1). The HAWC spectrum is consistent with IACT observations made by MAGIC between 2007 and 2009 (Ahnen et al. 2016). Using an energy threshold of 400 GeV, these authors found photon fluxes varying from N_min (>0.5 TeV) = 9.33 × 10⁻¹² cm⁻² s⁻¹ to N_max(>0.5 TeV) = 2.22 × 10⁻¹⁰ cm⁻² s⁻¹, weakly dependent on their assumed differential index of 2.5.

The 3FHL catalog and HAWC spectra for Mkn 421 are shown together in Figure 9. There is a fairly good match between the two fits, with the curves intersecting at about 0.66 TeV. The local LAT spectral index at 1 TeV, 2.5 ± 0.2, is consistent within uncertainties with the HAWC spectral index. IACT observations between 100 GeV and 5 TeV, performed around 2005 with MAGIC, resulted in a spectral index of (2.20 ± 0.08), with indications of a cutoff in the intrinsic spectrum (Albert et al. 2007a). A detailed spectral analysis of Mkn 421 and Mkn 501 with HAWC data will be presented in a separate publication. Preliminary results can be found in Coutiño de Leon et al. (2019).

Zoom In Zoom Out Reset image size

5.2.2. Markarian 501

The TeV emission from Markarian 501 observed by HAWC is not as steady as that of Mkn 421 (Abeysekara et al. 2017a). In fact, the statistical significance of the time-averaged TeV emission of Mkn 501 has decreased with increased HAWC exposure. The α = 2.5 search resulted in a test statistic TS = 276.97 for K = (7.74 ± 0.49_stat ± 1.16_syst) × 10⁻¹² TeV⁻¹ cm⁻² s⁻¹. The optimized power-law fit is

$$\begin{eqnarray}&&\,\begin{array}{l}\displaystyle \frac{{dN}}{{dE}}=(6.21\pm {0.69}_{\mathrm{stat}}\pm {0.93}_{\mathrm{syst}})\times {10}^{-12}\\ \ \ \times \ {\left(\displaystyle \frac{E}{1\mathrm{TeV}}\right)}^{-2.31\pm {0.08}_{\mathrm{stat}}\pm {0.12}_{\mathrm{syst}}}\,{\mathrm{TeV}}^{-1}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1},\end{array}\,\end{eqnarray} \tag{ 9 }$$

with TS = 280.28, representing a moderate increase in the test statistic ΔTS = 3.33. The integrated observed photon flux is N_obs = (8.5 ± 2.6) × 10⁻¹² cm⁻² s⁻¹, accounting for EBL attenuation. The intrinsic energy flux, f_E = (40 ± 16) × 10⁻¹² erg cm⁻² s⁻¹, translates into a luminosity per solid angle of d_L(z)²f_E = L(>0.5 TeV)/ΔΩ = (7.3 ± 2.8) × 10⁴² erg s⁻¹ sr⁻¹, which is about 25% of the (10 GeV–1 TeV) luminosity per steradian measured by Fermi-LAT (Figure 1). This fraction is similar to that observed for Mkn 421.

The 3FHL catalog and HAWC spectra for Mkn 501 are shown together in Figure 10. The agreement is not as good as for Mkn 421: the Mkn 501 spectrum is harder and lies below the measurement by the LAT. The local LAT spectral index at 1 TeV is 2.58 ± 0.35, just consistent with HAWC when accounting for the propagation of the uncertainty in the curvature parameter β of the 3FHL fit. We note that Mkn 501 has a variability index V_bayes = 4 in the 3FHL catalog.

Zoom In Zoom Out Reset image size

The literature on the VHE characteristics of Mkn 501 is extensive. Its high TeV variability was noted shortly after its 1996 discovery, when the HEGRA group reported flux variations of an order of magnitude observed in mid-1997 (Quinn et al. 1996; Aharonian et al. 1997). Contemporaneous 10 m Whipple data confirmed the high state, adding indications of a curved spectrum favored over a simple power law (Samuelson et al. 1998). Further monitoring showed strong variations in flux, although with stable spectra best-fitted by a power law plus an intrinsic cutoff at around 6 TeV (Aharonian et al. 1999). Observations in 1998–1999 showed lower activity, with evidence for spectral curvature and steepening (Aharonian et al. 2001). Mkn 501 observations renewed in 2005, when the MAGIC telescope measured strong and very fast variability, with spectral indices ranging from ∼2.0 to 2.7 (Albert et al. 2007b). Data taken the following year (2006) by MAGIC permitted the characterization of a low-activity state, with fluxes similar to those reported by VERITAS and MAGIC from data taken in the 2009 joint observations with the Fermi-LAT (Anderhub et al. 2009; Acciari et al. 2011). The contemporaneous LAT data showed spectral variability also present in the GeV range, with index variations Δα ∼ 1 during the first 480 days of Fermi observations (Abdo et al. 2011).

In the last decade, EAS arrays have been able to reach the sensitivities needed to perform long-term monitoring of Mkn 501: the ARGO collaboration reported variations of a factor of six in flux observed between 2011 October and 2012 April (Bartoli et al. 2012). The HAWC collaboration presented light curves for the Crab, Mkn 421, and Mkn 501 from its first 17 months of observations (Abeysekara et al. 2017a). The HAWC light curve of Mkn 501 showed a low flux baseline with a handful of very short strong flares. The simple power-law fit (of index 2.84) was disfavored against a power law with exponential cutoff.

The variability of Mkn 501 has prevented a baseline characterization of this object. In Coutiño de Leon et al. (2019) we presented the HAWC spectrum of Mkn 501 using data acquired between 2015 June and 2017 December. Even though spectra with exponential cuts were tested, the Mkn 501 data proved consistent with a pure power law of index 2.40 ± 0.06. While there might be a slight decrease in the spectral index (2.31 ± 0.08 ± 0.20), the flux reported here is about half that reported in Coutiño de Leon et al. (2019). The four year average TeV flux observed by HAWC is a factor of two above the lowest activity observed so far (Anderhub et al. 2009; Acciari et al. 2011).

5.3.1. M87

M87 is the central galaxy of the Virgo cluster, a giant elliptical at a distance of just (16.4 ± 0.5) Mpc, as measured independently of its redshift, z = 0.0042 (Bird et al. 2010). With an optical magnitude V = 8.6 and an angular diameter of about 8', M87 has been imaged in detail for more than a century, showing a distinct bright active nucleus and a single optical jet (Curtis 1918; Tsvetanov et al. 1998). The mass of the central black hole has been measured to be (6.5 ± 0.7) × 10⁹ M_⊙, through its imaging with the Event Horizon Telescope (Event Horizon Telescope Collaboration et al. 2019).

Also known as Virgo A, it is a bright object all throughout the electromagnetic spectrum. As 3FHL J1230.8+1223, it is a 12.1σ detection above 10 GeV, as reported in the 3FHL catalog, with 4.8σ in 150–500 GeV, and a 2σ upper limit in the top (0.5–2.0) TeV band. M87 has been observed frequently in the VHE regime since its 2003 discovery by the HEGRA collaboration (Aharonian et al. 2003). The temporal behavior of M87 in the high-energy and VHE bands has been reviewed by Ait Benkhali et al. (2019). Fluxes can vary by a factor of 1o between low, mid, and high states, with indications of the spectral index varying from 2.6 in low states to 2.2 in high states. In addition to the low and high states, time variability on single-day timescales during high states has been reported by the H.E.S.S. and MAGIC collaborations (Beilicke et al. 2007; Albert et al. 2008).

The default HAWC search at the M87 location gave a test statistic of TS = 12.93 for a normalization K = (0.56 ± 0.16_stat ± 0.08_syst) × 10⁻¹² TeV⁻¹ cm⁻² s⁻¹. The optimized power-law fit for the intrinsic spectrum is

$$\begin{eqnarray}&&\,\begin{array}{l}\displaystyle \frac{{dN}}{{dE}}=(0.69\pm {0.22}_{\mathrm{stat}}\pm {0.10}_{\mathrm{syst}})\times {10}^{-12}\\ \ \ \times \ {\left(\displaystyle \frac{E}{1\mathrm{TeV}}\right)}^{-2.63\pm {0.18}_{\mathrm{stat}}\pm {0.13}_{\mathrm{syst}}}\,{\mathrm{TeV}}^{-1}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1},\end{array}\,\end{eqnarray} \tag{ 10 }$$

with TS = 13.19, i.e., a nonsignificant increase in the test statistic, ΔTS = 0.26, showing α = 2.5 to be an acceptable solution within the statistics of the optimized fit. The integrated photon flux is N_obs(>0.5 TeV) = (1.3 ± 0.9) × 10⁻¹² cm⁻² s⁻¹, with less than 4% of attenuation by the EBL. The energy flux, f_E = (2.7 ± 2.4) × 10⁻¹² erg cm⁻² s⁻¹, translates into a luminosity per solid angle of L(>0.5 TeV)/ΔΩ = (6.9 ± 6.3) × 10³⁹ erg s⁻¹ sr⁻¹, which is about 25% of that in the 10 GeV–1 TeV range, from the respective 3FHL parameters. The lower apparent power with respect to the two nearest BL Lacs is attributed to the off-axis viewing of the jet. Still, the relative power when compared with the LAT regime appears to be similar.

From Ait Benkhali et al. (2019), we computed photon fluxes from M87 observations by IACTs, in different activity states:

• 1.  
  High state: N(>0.5 TeV) = (5.74 ^+ 1.14_−1.47) × 10⁻¹² cm⁻² s⁻¹,
• 2.  
  Mid state: N(>0.5 TeV) = (1.39 ^+ 0.47_−0.43) × 10⁻¹² cm⁻² s⁻¹,
• 3.  
  Low state: N(>0.5 TeV) = (2.85 ^+ 2.07_−1.53) × 10⁻¹³ cm⁻² s⁻¹.

The 4.5 yr averaged emission, as indicated by the HAWC data, matches the mid state, with a relatively steep spectral index. The comparison with the 3FHL data, shown in Figure 11, points to a steepening in the spectrum. While the LAT and HAWC data are not contemporaneous, with a variability index V_bayes = 1, M87 does not stand as a variable in the 3FHL catalog. Still, an analysis of joint contemporaneous LAT-HAWC data is desirable.

Zoom In Zoom Out Reset image size

5.3.2. VER J0521+211

VER J0521+211 was discovered as a TeV γ-ray source by the VERITAS collaboration, during observations following up its detection above 30 GeV by LAT (Ong 2009). 3FHL J0521.7+2112 itself is associated with the radio-loud BL Lac object TXS 0518+211, and the X-ray ROSAT source RX J0521.7+2112. The redshift survey of Shaw et al. (2013) assigned z = 0.108 to the optical counterpart. This is the value listed in the 3FHL catalog and assumed for this analysis. Archambault et al. (2013) did not confirm this redshift, and neither did Paiano et al. (2017), who only set a lower limit z > 0.18 using deep spectroscopy with the GTC 10.4 m. The VHE emission from this object has been measured at least up to 1 TeV, with no apparent decline in the spectrum (Prokoph et al. 2015). This makes the high-energy characterization of VER J0521+211 of particular interest.

The α = 2.5 search resulted in TS = 9.49 for K =(2.85 ± 0.93_stat ± 0.43_syst) × 10⁻¹² TeV⁻¹ cm⁻² s⁻¹. The best power-law fit for the intrinsic spectrum is

$$\begin{eqnarray}&&\,\begin{array}{l}\displaystyle \frac{{dN}}{{dE}}=(2.39\pm {0.89}_{\mathrm{stat}}\pm {0.36}_{\mathrm{syst}})\times {10}^{-12}\\ \ \ \times \ {\left(\displaystyle \frac{E}{1\mathrm{TeV}}\right)}^{-2.01\pm {0.38}_{\mathrm{stat}}\pm {0.10}_{\mathrm{syst}}}\,{\mathrm{TeV}}^{-1}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1},\end{array}\,\end{eqnarray} \tag{ 11 }$$

with TS = 10.34, representing a modest increase in the test statistic ΔTS = 0.85 with respect to the default search. This is the hardest AGN spectrum for the TS > 9 HAWC sample. The integrated observed photon flux is N_obs(>0.5 TeV) =(1.5 ± 1.1) × 10⁻¹² cm⁻² s⁻¹, attenuated by a factor of 2/3 due to the EBL. The energy flux and luminosity per solid angle cannot be accurately determined with the spectral index so close to 2.0. If we take E²dN/dE at 1 TeV as indicative, we get an estimated f_E ∼ (3.8 ± 1.4) × 10⁻¹² erg cm⁻² s⁻¹ and d_L²f_E ∼ (8.4 ± 3.1) × 10⁴² erg s⁻¹ sr⁻¹. Notwithstanding the uncertain distance, this source appears particularly luminous in the GeV regime, with L(10 GeV − 1 TeV)/ΔΩ ∼ 1.4 ×10⁴⁴ erg s⁻¹ sr⁻¹, standing above the two Markarians in Figure 1.

The 3FHL data, shown together with the HAWC spectral fit in Figure 12, have strong detections up to the 150–500 GeV band, with a highest-energy photon of 370 GeV. The 3FHL fit is a log-parabola, transiting from a hard spectrum at about 30 GeV to a very steep local spectral index of 3.7 ± 0.7 at 1 TeV (4.3 when attenuating the 3FHL spectrum). The attenuated HAWC spectrum corresponds to an observed spectral index of ∼2.7, mostly inconsistent with the LAT fit. The data may be reconciled through an intrinsic hardening at about 200 or 300 GeV. The HAWC quasi-differential bound in the 0.5–2.0 TeV band is a factor of 4.6 lower than the 3FHL limit in the same band. We note that the quasi-differential analysis gave TS = 9.1 for the 2.0–8.0 TeV band, optimal in terms of the HAWC instrumental response. With the assumed redshift, τ is due to range from 1.6 to 3.2 in that energy interval, pointing to an intrinsically hard spectrum. The analysis gave TS = 0 for 8.0–32.0 TeV, expected to be heavily attenuated for the assumed redshift. Additional HAWC data analysis should allow the further constraint of the shape of the TeV spectrum VER J0521+211, testing the redshift assumption. We also note that VER J0521+211 has a variability index V_bayes = 4 in the 3FHL catalog, indicating that a joint analysis of contemporaneous LAT and HAWC data would be relevant.

Zoom In Zoom Out Reset image size

As mentioned in Section 4.1, VER J0521+211 is located 307 from the Crab Nebula, the brightest source in the 2HWC and 3HWC catalogs. This angular distance corresponds to three times the 68% containment angle ψ₆₈ for ${ \mathcal B }=1$, and >6ψ₆₈ for ${ \mathcal B }\gt 2$. We tested for potential contamination, repeating the maximum-likelihood test with α = 2.5 at five locations equidistant from the Crab Nebula, forming together with VER J0521+211 a hexagon around the Crab. These provided test statistics TS between −4.46 and +1.54, in contrast with TS = +9.49 at the location of VER J0521+211.

5.3.3. 1ES 1215+303

1ES 1215+303 is one of six γ-ray emitting BL Lac objects located in the Northern part of the Coma Berenices constellation, five of them known to be VHE sources. The redshift of this HBL is now confirmed to be z = 0.130, discarding the early measurement z = 0.237 by Lanzetta et al. (1993). The lower value was confirmed through optical spectroscopy at the GTC and, more recently, through the direct identification of a Ly α emission line (Paiano et al. 2017; Furniss et al. 2019). Cataloged as 3FHL J1217.9+3006, this object is well detected up to 500 GeV, and modeled with a power law of index α = 2.3 ± 0.1. It has a variability index V_bayes = 2.

The default HAWC search gave TS = 11.36 for K =(4.64 ± 1.38_stat ± 0.70_syst) × 10⁻¹² TeV⁻¹ cm⁻² s⁻¹ as the intrinsic normalization. The optimized power-law fit for the intrinsic spectrum is rather softer,

$$\begin{eqnarray}&&\,\begin{array}{l}\displaystyle \frac{{dN}}{{dE}}=(3.78\pm {1.36}_{\mathrm{stat}}\pm {0.57}_{\mathrm{syst}})\times {10}^{-12}\\ \ \ \times \ {\left(\displaystyle \frac{E}{1\mathrm{TeV}}\right)}^{-3.07\pm {0.37}_{\mathrm{stat}}\pm {0.15}_{\mathrm{syst}}}\,{\mathrm{TeV}}^{-1}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1},\end{array}\,\end{eqnarray} \tag{ 12 }$$

with TS = 12.80, representing a test-statistic increase ΔTS =1.44. The attenuated photon flux is N_obs(>0.5 TeV) = (2.53 ± 1.35) × 10⁻¹² cm⁻² s⁻¹. The integrated energy flux results in a luminosity per solid angle of L(>0.5 TeV)/ΔΩ =(3.9 ± 2.7) × 10⁴³ erg s⁻¹ sr⁻¹, given the luminosity distance of 585 Mpc. This is about 30% of the corresponding value between 10 GeV and 1 TeV. The 3FHL and HAWC spectra, compatible within the uncertainties, are shown in Figure 13.

Zoom In Zoom Out Reset image size

1ES 1215+303 was first detected as a VHE source by MAGIC in 2011 (Lombardi et al. 2011; Aleksić et al. 2012). It is usually observed together with PG 1218+304, located just 0.88 degrees away in the sky, for which we obtained a test statistic TS = +2.24. Given that both sources are known VHE emitters, the distance is not large enough to be certain that there is no overlap between both sources. With a redshift z = 0.184, PG 1218+304 is more distant and prone to be heavily attenuated above 1 TeV. A dedicated study with improved HAWC analysis tools is now pending.

Long-term monitoring of this source by VERITAS was presented by Aliu et al. (2013), prior to the report of a single short and intense flare seen in 2014 with VERITAS and Fermi-LAT (Abeysekara et al. 2017d). During this episode the VHE flux of this HBL reached 2.4 times that of the Crab Nebula, with a variability timescale ≲3.6 hr. The MAGIC spectrum of Aleksić et al. (2012), ranging from 70 GeV to 1.8 TeV, had an intrinsic spectral index of 2.96, deattenuated with the model of Domínguez et al. (2011), used here. The MAGIC and HAWC spectral indices are in good agreement. The long-term joint monitoring of this source by LAT and VERITAS, spanning more than 10 years, confirms the spectral index and points to four strong and short flares that occurred during full HAWC operations by Valverde et al. (2020). We checked the dates of these four episodes and did not find evidence for them in the HAWC data. An optimized time-dependent analysis of this source with HAWC, beyond the scope of this paper, will be the subject of future work.

The extrapolation of 3FHL catalog spectra (Figure 2) and of IACT spectral fits allowed us to identify AGN targets for potential HAWC detection (Section 4.3). We presented evidence of persistent TeV emission for three of those targets—Mkn 421, Mkn 501, and at a weaker level, M87—and, also marginally, for two sources below the 30 mCrab reference limit (marked as a red line in Figure 2), VER J0521+211 and 1ES 1215+303. On the other hand, sources like IC 310, 1ES 2344+514, TXS 0210+515, 1ES 1727+502, B3 2247+381, and H 1426+428 were not detected. As shown in Figure 4, the HAWC sensitivity is dependent on source decl. We note that while the sources with $\sqrt{\mathrm{TS}}\gt 9$ are all in the range +12° < δ < +40°, five of the undetected candidates are either North of δ = +50°, or farther than z = 0.1. Decl. is of particular relevance here, as the response of EAS arrays to low-energy events is compromised at large zenith angles. Still, as shown in Table 5 and in Section 5.1.3, four of the corresponding upper limits are below the extrapolation of the corresponding 3FHL spectrum.

Table 5.  Upper Limits on 3FHL AGN Candidate Targets

3FHL Entry	Source	z	$\pm \sqrt{\mathrm{TS}}$	K2σ	${N}_{0.5}^{x}$	${N}_{0.5}^{\mathrm{UL}}$
3FHL J0214.5+5145	TXS 0210+515	0.049	+1.62	5.43	2.40	6.50
3FHL J0316.6+4120	IC 310	0.019	+0.86	1.01	4.37	1sd.60
3FHL J1428.5+4240	H 1426+428	0.129	+1.26	8.18	1.21	4.67
3FHL J1728.3+5013	I Zw 187	0.055	−0.22	2.57	3.25	2.91
3FHL J2250.0+3825	B3 2247+381	0.119	−0.27	2.54	2.37	1.59
3FHL J2347.0+5142	1ES 2344+514	0.044	+1.45	4.62	7.12	5.80

Note. K_2σ, the spectrum normalization at 1 TeV, is in units of 10⁻¹² TeV⁻¹ cm⁻² s⁻¹. N_0.5 are extrapolated photon fluxes above 0.5 TeV and the corresponding upper limit, in units of 10⁻¹² cm⁻² s⁻¹.

Download table as:  ASCIITypeset image

IC 310 is a target that is both nearby and culminates at an adequate zenith angle, but remained undetected in this HAWC analysis. The upper limit set in the photon flux is close to a third of the rather uncertain LAT extrapolation, and the HAWC upper limit in the 0.5–2.0 TeV band is 2.5 lower than the one in the 3FHL catalog. We note that while the 3FHL catalog reports a relatively low variability index, V_bayes = 2, this source is known to display extreme variability in the VHE range, on timescales as low as five minutes, challenging models and severely constraining the emission region to scales smaller than its event horizon (Aleksić et al. 2014).

In addition to the preselected targets, we point here to two additional sources of intrinsic interest: 3C 264 and NGC 1068. 3C 264 is a radio galaxy hosted by the elliptical galaxy NGC 3862, at a distance of about 90 Mpc. A compact radio source powers a relativistic jet imaged in radio and in the optical (Crane et al. 1993; Lara et al. 2004). 3C 264 shows in the LAT data up to a highest-energy photon of 97 GeV. It was listed as a TeV = N source in the 3FHL catalog, and considered as such in this analysis. However, it was later detected in the VHE range by the VERITAS Collaboration (Mukherjee 2018). The photon flux measured by VERITAS above 300 GeV indicates that this object should be too faint for HAWC. The HAWC data show a +1.9σ excess, statistically consistent with such a low flux and providing an upper limit $N(\gt 0.5\,\mathrm{TeV})\lt 1.42\times {10}^{-12}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$.

NGC 1068 has become a source of interest due to its coincidence with a hotspot in the IceCube all-sky map (Aartsen et al. 2020). While two other starburst galaxies, NGC 253 and M82, have been detected in the VHE range as faint sources, with fluxes 1% of the Crab Nebula (Ohm 2016), NGC 1068 remains undetected and an unlikely candidate for HAWC detection. It is a weak 5.3σ detection in the 3FHL catalog, with practically no signal above 20 GeV and a rather steep spectral index α = 3.8 ± 1.0. The HAWC data have a +1.2σ excess at the location of NGC 1068, for an upper-limit normalization K_2σ = 6.46 × 10⁻¹³ TeV⁻¹ cm⁻² s⁻¹ at 1 TeV. The quasi-differential HAWC limit in the common energy band, N(0.5–2.0 TeV) ≤ 2.32 × 10⁻¹² cm⁻² s⁻¹, is a factor of eight lower than the respective LAT limit.