The HAWC real-time flare monitor for rapid detection of transient events

TL;DR

This paper introduces a real-time flare monitor for the HAWC observatory, capable of detecting very high energy transient events from blazars on short timescales, enabling rapid multiwavelength and multimessenger follow-up.

Contribution

It presents the development and operational performance of a novel real-time VHE flare detection system for HAWC, expanding survey capabilities for blazar transient events.

Findings

Detection of three high-confidence VHE events from blazars

Demonstrated sensitivity of the flare monitor

Operational since January 2017 with real-time capabilities

Abstract

We present the development of a real-time flare monitor for the High Altitude Water Cherenkov (HAWC) observatory. The flare monitor has been fully operational since 2017 January and is designed to detect very high energy (VHE; $E\gtrsim100$ GeV) transient events from blazars on time scales lasting from 2 minutes to 10 hours in order to facilitate multiwavelength and multimessenger studies. These flares provide information for investigations into the mechanisms that power the blazars' relativistic jets and accelerate particles within them, and they may also serve as probes of the populations of particles and fields in intergalactic space. To date, the detection of blazar flares in the VHE range has relied primarily on pointed observations by imaging atmospheric Cherenkov telescopes. The recently completed HAWC observatory offers the opportunity to study VHE flares in survey mode, scanning 2/3 of the entire sky every day with a field of view of $\sim$1.8 steradians. In this work, we report on the sensitivity of the HAWC real-time flare monitor and demonstrate its capabilities via the detection of three high-confidence VHE events in the blazars Markarian 421 and Markarian 501.

Tables

Table 1. Table 1: Definition of Analysis Bins

	Low	High	Smoothing
Bin	Fraction	Fraction	Radius [∘]
0	0.044	0.067	1.20
1	0.067	0.105	1.20
2	0.105	0.162	0.75
3	0.162	0.247	0.60
4	0.247	0.356	0.40
5	0.356	0.485	0.40
6	0.485	0.618	0.30
7	0.618	0.740	0.30
8	0.740	0.840	0.30
9	0.840	1	0.25

Table 2. Table 2: Definition of Target Categories

Category	Targets	Threshold Parameter $\gamma$
Mrk 421 & Mrk 501	2	$1.6\times {10}^{-7}$
Other TeVCat Extragalactic	44	$4.2\times {10}^{-9}$
Nearby 2FHL Extragalactic ()	22	$9.3\times {10}^{-9}$
Other 2FHL Extragalactic & Unknown	119	$1.3\times {10}^{-9}$
Total	187

Table 3. Table 3: High-Confidence Markarian Triggers Identified by the Archival Search

	False Alarm		Trigger	Time to	Search	Bayesian Block	Approximate
Event	Rate [${\mathrm{century}}^{-1}$]	Target	MJD	Detection [min.]	Threshold	Duration [min.]	Flux [Crab units]
1	0.68	Mrk 501	57252.068	34	1%	74	3.8 $\pm$ 1.1
2	2.1	Mrk 421	57020.357	112	5%	158	3.9 $\pm$ 0.8
3	2.7	Mrk 501	57032.614	80	5%	148	2.1 $\pm$ 0.6

Equations

$$P(\mathbf{n},\mathbf{m}|\mathbf{\nu },\mathbf{\mu })=\prod _{i=j}^k\frac{e^{-{\nu }_i}{\nu }_i^{n_i}}{n_i!}\frac{e^{-{\mu }_i}{\mu }_i^{m_i}}{m_i!}=L_{b;j,k},$$

$$L_b^{\mathrm{tot}}=\prod _{i=1}^B{L}_{b;l_i,u_i}.$$

$$\begin{array}{rl}\ln L_{b;j,k}=& \sum _{i=j}^k[-\alpha {\mu }_i{\rho }_i+n_i\ln (\alpha {\mu }_i{\rho }_i)-{\mu }_i\\ & +m_i\ln {\mu }_i-S_i],\end{array}$$

$${\rho }_{j,k}=\frac{N_{j,k}}{\alpha M_{j,k}};$$

$$N_{j,k}\equiv \sum _{i=j}^k{n}_i;$$

$$\begin{array}{rl}\ln L_{b;j,k}^{\max }=& -N_{j,k}-N_{j,k}\ln \left(1+\frac{M_{j,k}}{N_{j,k}}\right)\\ & -M_{j,k}-M_{j,k}\ln \left(1+\frac{N_{j,k}}{M_{j,k}}\right)\\ & +\sum _{i=j}^k\left[(n_i+m_i)\ln (n_i+m_i)-S_i\right].\end{array}$$

$$D_b=\Delta \ln L_b^{\max }=\sum _{i=1}^B\ln L_{b;l_i,u_i}^{\max }-\ln L_{b;1,T}^{\max }.$$

$$\begin{array}{rl}{D}_b=& N_{1,T}\ln \left(1+\frac{M_{1,T}}{N_{1,T}}\right)+M_{1,T}\ln \left(1+\frac{N_{1,T}}{M_{1,T}}\right)\\ & -\sum _{i=1}^B{N}_{l_i,u_i}\ln \left(1+\frac{M_{l_i,u_i}}{N_{l_i,u_i}}\right)\\ & -\sum _{i=1}^B{M}_{l_i,u_i}\ln \left(1+\frac{N_{l_i,u_i}}{M_{l_i,u_i}}\right).\end{array}$$

$$D=\sum _b{D}_b.$$

$$\begin{array}{rl}{D}_b=& N\ln \left(1+\frac{M}{N}\right)+M\ln \left(1+\frac{N}{M}\right)\\ & -N_1\ln \left(1+\frac{M_1}{N_1}\right)-M_1\ln \left(1+\frac{N_1}{M_1}\right)\\ & -N_2\ln \left(1+\frac{M_2}{N_2}\right)-M_2\ln \left(1+\frac{N_2}{M_2}\right).\end{array}$$

$$D=\sum _b\max \left[0,\mathrm{sign}({\rho }_2-{\rho }_1)D_b\right],$$

$$D^{\max }>-\ln \gamma ,$$

$$r_{\mathrm{target}}=\frac{n_{\mathrm{target}}-\alpha m_{\mathrm{target}}}{\alpha m_{\mathrm{target}}}.$$

$$D^{\max }>-\ln \gamma +k_{\mathrm{target}}.$$

Full text

astroparticle physics — gamma rays: galaxies — galaxies: BL Lacertae objects: individual (Mrk 421, Mrk 501) — methods: data analysis

\AuthorCollaborationLimit

=300

1 Blazar Transients at Very High Energies

Active galactic nuclei (AGNs) comprise a class of extremely luminous galactic cores that are widely understood to be powered by accretion onto a supermassive black hole (Begelman et al., 1984). Many AGNs exhibit highly collimated jets, first observed by Curtis (1918) in the object M87, that efficiently transport plasma via bulk relativistic outflows (Boettcher et al., 2012). In a small number of AGNs, the jet is serendipitously aligned along the line of sight to Earth. Under the standard AGN unification scheme (Urry & Padovani, 1995), such objects are referred to as blazars, and their enhanced emission due to Lorentz boosting serves as a probe into the production and dynamics of relativistic jets.

Perhaps the most striking feature of blazars is their extreme variability. Long-term observations in the VHE band reveal that the flux may reach levels several times higher than the quiescent emission and persist in such states for years (Tluczykont et al., 2007; Acciari et al., 2014). Extreme VHE flares surpassing the quiescent emission by factors of 100 and with observed variability time scales as short as minutes have also been observed in some sources (e.g. Aharonian et al., 2007). These flares are of particular interest because they facilitate studies of the mechanisms powering the blazar central engine (Katarzyński et al., 2005; Madejski & Sikora, 2016) and they place stringent constraints on the bulk jet Lorentz factors (Spada et al., 2001; Krawczynski et al., 2001; Begelman et al., 2008). They also constrain the location of the region where electromagnetic energy in the initially Poynting-dominated jet flux (Blandford & Znajek, 1977; Blandford & Payne, 1982) is dissipated into accelerated particles, producing flares (Nalewajko et al., 2014).

In general, multiwavelength observations providing spectral and, where possible, morphological or polarization information during and after a flare have yielded the most complete pictures of the dissipation process (e.g. Marscher et al., 2008; Acciari et al., 2009; Abdo et al., 2010). Developing an understanding of dissipation in jets will in turn help explain jet composition and structure, accretion, jet formation, and AGN feedback (Bykov et al., 2012). Furthermore, the increased fluxes and short variability time scales defining extreme VHE flares render them useful for studying Lorentz invariance violation (Biller et al., 1999; Albert et al., 2008; Abramowski et al., 2011), the extragalactic background light (EBL; e.g. Mazin & Raue, 2007; Biteau & Williams, 2015), or the intergalactic magnetic field (Aharonian et al., 1994; Plaga, 1995; Neronov & Semikoz, 2009).

Although many informative studies on extreme VHE flares have been performed, until the present time no dedicated VHE survey instrument has been available for monitoring large numbers of blazars simultaneously. Imaging atmospheric Cherenkov telescopes (IACTs) are the most sensitive VHE instruments, but their narrow fields of view prohibit unbiased and continuous monitoring of more than a few sources. The Large Area Telescope (LAT) on board the Fermi Gamma-Ray Space Telescope (hereafter Fermi-LAT) surveys the sky in the energy range from 100 MeV to 300 GeV (Atwood et al., 2009), but flares detected by the Fermi-LAT are not known a priori also to produce VHE emission. Consequently, the advent of the High Altitude Water Cherenkov (HAWC) observatory, the most sensitive VHE survey instrument yet constructed, opens a new window of opportunity for identifying extreme VHE flares. The rapid identification of these events with HAWC permits follow-up observations with the more sensitive IACTs—VERITAS, HESS, MAGIC, and FACT—as well as with instruments sensitive to other energies or other particle types, such as Swift, Fermi, NuSTAR, Chandra, and IceCube.

Many VHE flares exhibit soft spectra with photon indices $\Gamma$ in the range from 3–4 (e.g. Arlen et al., 2013; Abeysekara et al., 2017a) for power-law spectra of the form $dN/dE\propto E^{-\Gamma}$. However, some flares have been observed with harder indices, in the range from 2–3, either over the entire VHE range (Archambault et al., 2015; Buson et al., 2016) or as a hard low-energy component up to several hundred GeV (Aharonian et al., 2007). This makes them good candidates for detection by a dedicated survey instrument such as HAWC. Since the most luminous flares (e.g. those that exceed the flux from the Crab Nebula) tend to have short durations, obtaining a multiwavelength picture of the brightest transients requires a rapid alert system capable of issuing alerts before the transient activity has finished.

In this work, we present the development and deployment of a real-time flare monitor designed to identify extreme VHE flares with the HAWC observatory and issue alerts rapidly to other instruments. Following a brief summary of HAWC operations in Section 2, we present the mechanics of the flare monitor in Section 3, describe the tuning of the false alarm rate in Section 4, and characterize the flare monitor sensitivity in Section 5. Section 6 describes the selection of flare monitor targets. In Section 7, we apply the flare monitor to archival data on the two nearest blazars, Markarian 421 and Markarian 501 (hereafter Mrk 421 and Mrk 501), located respectively at redshifts of $z=0.031$ and $z=0.033$ (Mao, 2011). The flare monitor identifies three high-confidence events in these data. We discuss our future plans for the flare monitor in Section 8.

2 The HAWC Observatory

The HAWC observatory consists of a close-packed array of 300 optically isolated water Cherenkov detectors, situated at 4100 m above sea level in the state of Puebla, Mexico. Located at $19^{\circ}$ north latitude, HAWC surveys 2/3 of the entire sky with an uptime in excess of $90\%$, making it an ideal instrument for searching for VHE transient phenomena. Each water Cherenkov detector contains four upward-facing photomultiplier tubes (PMTs) that sample the energy deposition at ground level from air showers initiated by the interactions of primary particles with the atmosphere. HAWC detects events in the energy range from $\sim$100 GeV to $>$100 TeV. Further details concerning the performance of HAWC may be found in Abeysekara et al. (2017b).

Over $99.9\%$ of the events that HAWC detects are cosmic rays that must be separated efficiently from any gamma-ray signal present in the data. Due to the large transverse momentum carried by the products of hadronic interactions in the air showers, the distribution of energy is much less smooth for showers initiated by cosmic rays than it is for gamma-ray showers. To remove a large fraction of the cosmic rays, we apply gamma-hadron separation cuts in two parameters designed to characterize the smoothness of the energy distribution. As described in Abeysekara et al. (2017b), one parameter tests for a single localized excess far from the shower core, while another tests for rotational symmetry about the shower axis.

We apply the gamma-hadron separation cuts separately to analysis bins defined by the fractional number of operational PMTs that participate in the event. The present work includes the nine analysis bins used in Abeysekara et al. (2017b) but also includes events that trigger between $4.4\%$ and $6.7\%$ of the detector in an additional bin. The bins are numbered from bin 0, in which the smallest fraction of PMTs participate, to bin 9, in which nearly every PMT in the detector participates. Table 1 shows the high and low fractions of operational PMTs that define each of the ten bins.

Each event that triggers at least $4.4\%$ of the operational PMTs falls into exactly one analysis bin. When determining the events that are associated with a given source, we smooth the events with a top-hat function with a radius that depends on the analysis bin. The radius that maximizes the significance of the Crab Nebula, the strongest steady source of VHE gamma rays, appears in the fourth column of Table 1.

While other recent HAWC publications (e.g. Abeysekara et al., 2017b, c) report results in terms of source fluxes, the real-time flare monitor uses only the observed event counts, not fluxes, to decide whether to issue an alert. It is therefore less affected by an existing discrepancy between observed event counts and Monte Carlo simulation in bin 0 that prevents the inclusion of bin 0 events in the flux estimation. Efforts are currently ongoing to resolve this discrepancy. Bin 0 generally contains the lowest energy gamma rays that are of considerable interest to the study of blazars. We have tested the effect of including bin 0 by performing simulations as described in Section 5 with and without the bin 0 data included. We find that the inclusion or exclusion of bin 0 does not affect the sensitivity of the real-time flare monitor, and therefore we include it in the present analysis.

In general, events with larger energies and smaller zenith angles tend to fall into the higher analysis bins. However, there is substantial overlap in the distribution of energies within the bins, so the analysis bins cannot be interpreted as strict energy bins. Instead, the analysis bins correspond more directly to the amount of information available about a given event, with the higher analysis bins providing more information and better angular resolution. This renders the analysis bins suitable for the likelihood method described in the following section. As was done for the smoothing radius, we tune the cuts in each analysis bin to maximize the sensitivity to the Crab Nebula.

While many AGNs exhibit spectra harder than that of the Crab at lower energies, these are softened at higher energies due to EBL absorption. Typical EBL models (e.g. Franceschini et al., 2008; Domínguez et al., 2010; Gilmore et al., 2012) indicate that the attenuation of gamma rays via pair production becomes of order unity around a redshift of $z\approx 1$ for an observed energy of 100 GeV, or a redshift of $z\approx 0.2$ for gamma rays with energies of $\sim$400 GeV. Furthermore, in contrast to that of the IACTs, HAWC’s sensitivity varies slowly with energy, reaching its maximum only in the multi-TeV regime (see Figure 14 of Abeysekara et al., 2017b). As a result of these effects, HAWC is most sensitive to hard flares from nearby blazars. Since several known VHE blazars lie at redshifts $z\lesssim 0.3$, and intrinsically hard ($\Gamma\lesssim 2$) VHE spectra are likely to be rendered softer than the Crab, we expect that the sensitivity reported here is somewhat optimistic for the detection of realistic blazar flares. At present, we reserve an in-depth investigation into the precise sensitivity to different flare properties for future work.

3 Flare Monitor Search Method

The HAWC real-time flare monitor is a fully automated system capable of rapid transient detection on minute to hour time scales. It complements other ongoing HAWC search efforts that are tuned to longer time scales or focused on gamma-ray bursts. As described in Section 6, we choose a number of targets to monitor based on their high probability to produce VHE gamma rays. We collect the data from each monitored target into a buffer containing 300 observations. Each observation has a duration of 2 minutes, which is consequently the minimum time resolution for the method. The buffer depth of 300 observations is chosen as a compromise between computer memory utilization and providing a sufficient number of observations to detect flares lasting up to a few hours. We record data in the buffer exclusively at times for which the target is located at zenith angles less than $45^{\circ}$. As a result, the 10 hours of buffer time correspond to a much larger amount of real time that includes at least 4 hours of time from the previous day’s observations. For targets with declinations far from the latitude of HAWC, the buffer contains several days of observations.

Following the language of Li & Ma (1983), we use an off-source region to estimate the background contribution to the event count in the on-source region. The off-source count is derived from the direct integration technique (Atkins et al., 2003) with an integration time of 2 hours. Because the flare monitor must operate in real time, the result of the direct integration from the previous 2 hours is used when deriving the background estimation for the current observations. Most of the time, the response of the detector to air showers is sufficiently stable that the distribution of reconstructed events in horizontal (Earth-fixed) coordinates does not change significantly over the course of 4 hours. This permits us to use the previous direct integration result. We address occasional changes in the detector response in Section 4.3.

3.1 Likelihood Ratio Test for Flare Detection

Each time the buffer is updated with a new observation, we perform a likelihood ratio test to search for an increase in the observed flux. In this test, the null hypothesis is that the data can be represented by a constant flux for the entire duration of the buffer. For the alternative hypothesis, we assume that the flux is different at different points in the buffer. In what follows, we first derive the likelihood ratio for the case when the buffer is divided into any number of distinct flux states. We then specialize to the case for which there are only two states: an initial low state followed by a high state. This is the case of interest for the detection of a flare.

We begin by considering a single analysis bin $b$. For time bin $i$, we represent the observed number of on-source events in the analysis bin as $n_{i}$ and the observed number of off-source events as $m_{i}$. We assume that the on-source count $n_{i}$ is derived from a Poisson distribution with mean $\nu_{i}$, and that the off-source count $m_{i}$ is derived from a Poisson distribution with mean $\mu_{i}$. Following Li & Ma (1983), we designate the ratio of on-source to off-source exposure as $\alpha$. The value of $\alpha$ depends on the analysis bin, ranging from $0.082$ for bin 0 down to $0.017$ for bin 9. The probability to observe the measured counts given the model parameters for time bins running from $j$ to $k$ inclusive is then

[TABLE]

for analysis bin $b$. Here, $L_{b;j,k}$ is the likelihood of the model specified by $\bm{\nu}$ and $\bm{\mu}$ giving rise to the observations in the specified range of time bins. For clarity, we omit the subscript $b$ on $n_{i}$, $m_{i}$, $\nu_{i}$, $\mu_{i}$, and $\alpha$.

Let us now assume that we have $T$ time bins in the buffer, and that the observations are partitioned into $B$ distinct flux states, with $1\leq B\leq T$. The null hypothesis is thus formulated as $B=1$ and the alternative hypothesis as $B>1$. Let the position of the lower edge of flux state $j$ in the buffer be denoted by $l_{j}$ and the upper edge by $u_{j}$, such that the flux state runs from time bins $l_{j}$ to $u_{j}$, inclusive. Clearly, $l_{1}=1$ and $u_{B}=T$, and since the flux states are not overlapping, $l_{j+1}=1+u_{j}$. The total likelihood of the model for bin $b$ is then the product of the likelihoods for each flux state:

[TABLE]

Although we have assumed a constant flux in each state, Equation 2 can accommodate any flare shape with a time resolution equal to the time bin size of 2 minutes.

It is convenient to re-cast Equation 1 in terms of the ratio $\rho_{i}\equiv\nu_{i}/(\alpha\mu_{i})$. In this case, the log-likelihood takes the form

[TABLE]

where $S_{i}=\ln(n_{i}!)+\ln(m_{i}!)$.

Adopting a frequentist approach, we maximize the likelihood with respect to the parameters $\mu_{i}$ and $\rho_{i}$. If we assume that the flux is constant from time bin $j$ through time bin $k$, then $\rho_{i}$ takes on a single, constant value $\rho_{j,k}$. Since $\rho_{i}$ is constructed as a ratio of on-source to off-source counts, the dependence of the detector response on the zenith angle as the target moves across the sky is canceled. We allow the values of $\mu_{i}$ to vary from time bin to time bin because we expect these to be affected by the detector response, which depends on the instantaneous zenith angle of the target. The solutions that maximize the likelihood are

[TABLE]

where

[TABLE]

With these solutions, the expression for the maximum log-likelihood becomes

[TABLE]

The log of the ratio of the alternative hypothesis likelihood to the null hypothesis likelihood is given by

[TABLE]

Here, $D_{b}$ is a test statistic that is equal to half the usual test statistic for the likelihood ratio test, $2\Delta\ln L_{b}$. By inspection of Equation 6, we see that the sums drop out of Equation 7, leaving us with

[TABLE]

The total test statistic is simply the sum of the test statistics for the individual analysis bins:

[TABLE]

Thus, $D$ is the log of the full likelihood ratio for any alternative hypothesis that divides the buffer into states of constant flux, compared to the null hypothesis of a single uniform flux state.

3.2 Flare Monitor Trigger Condition

The flare monitor searches for a single increased flux state in the buffer, which can be represented as an alternative hypothesis with $B=2$. We allow the increased flux state to begin at any point $C$ in the buffer. For convenience let us define $N\equiv N_{1,T}$, $M\equiv M_{1,T}$, $N_{1}\equiv N_{1,C-1}$, $M_{1}\equiv M_{1,C-1}$, $N_{2}\equiv N_{C,T}$, and $M_{2}\equiv M_{C,T}$. That is, $N_{1}$ and $M_{1}$ are the sums of the on-source and off-source counts prior to the proposed increased flux state, and $N_{2}$ and $M_{2}$ are the respective sums during the proposed increased flux state. Under these assumptions, Equation 8 takes the form

[TABLE]

Since we are interested in issuing alerts only for flux increases, we require that our alternative hypothesis show an increase in the ratio $\rho$. That is, if we define $\rho_{1}\equiv N_{1}/(\alpha M_{1})$ and $\rho_{2}\equiv N_{2}/(\alpha M_{2})$, then we replace Equation 9 with

[TABLE]

where as before we have omitted the subscript $b$ on $\rho_{1}$ and $\rho_{2}$ for clarity.

The value of $D$ given by Equation 11 depends on $C$, the starting point of the hypothesized flare in the buffer. The flare monitor determines the value of $D$ for each possible value of $C$, selecting the largest value, $D^{\mathrm{max}}$. Due to the correlations involved in this search, we cannot easily derive the analytic distribution of $D^{\mathrm{max}}$ in the absence of any flares. We elect instead to set a trigger threshold condition on the value of $D^{\mathrm{max}}$ and characterize the false alarm rate of the method as a function of this threshold. Our trigger condition takes the form

[TABLE]

where $\gamma$ is a threshold parameter that fulfills $0<\gamma<1$.

4 Tuning the False Alarm Rate

The threshold parameter $\gamma$ allows us to control the false alarm rate of the flare detection method. Under ideal circumstances, we would be able to calculate the false alarm rate analytically. However, because the search is performed on a sliding time window, the large degree of correlation between searches prohibits such calculation. We therefore estimate the false alarm rate by applying the search algorithm to simulated data. To do this, we first determine the false alarm rate for sources located at the same declination as the Crab Nebula, which culminates at a zenith angle of $3^{\circ}$. We then extend the simulation to other declinations in order to derive a declination-dependent correction factor for $\gamma$. This correction brings the relationship between false alarm rate and threshold parameter into agreement at all declinations. Following this correction, we identify data quality cuts necessary to resolve discrepancies between the simulated false alarm rate and the false alarm rate derived from data collected from regions of the sky that are at least $2^{\circ}$ away from any known VHE source or monitored target. Finally, we discover that an additional correction to the threshold parameter is required for targets with strong steady VHE emission, on account of the small number of events in several of the analysis bins.

4.1 Simulating the False Alarm Rate

We derive our simulated observations based on the observed off-source counts in the data. The simulation thus mimics the rise and fall of the rates as the target rises and sets. The only parameter that we require from the data is the observed off-source count $m_{i}$ in time bin $i$, for each analysis bin $b$. After obtaining $m_{i}$ from the data, we sample a simulated background count $m_{i}^{\prime}$ from a Poisson distribution with mean $m_{i}$ and a simulated source count $n_{i}^{\prime}$ from a Poisson distribution with mean $\alpha m_{i}$. This procedure enables us to simulate many thousands of years of data in order to understand precisely the relationship between the parameter $\gamma$ and the false alarm rate, while maintaining a realistic distribution of off-source counts as the target changes zenith angles while crossing the sky.

Since we sample $m_{i}^{\prime}$ and $n_{i}^{\prime}$ from the off-source distribution, the observed on-source distribution does not affect the simulation, so we can use data collected from regions with strong VHE sources as inputs to the simulation. Figure 1 shows the relationship between $\gamma$ and the false alarm rate derived from this method, using simulation sampled from Crab Nebula data. We fit the points in Figure 1 with a third-order spline in order to convert false alarm rate to threshold parameter. The threshold parameters for false alarm rates of 1 event per target per year and 1 event per target per century are $\gamma=1.2\times 10^{-7}$ and $\gamma=5.5\times 10^{-10}$, respectively.

4.2 False Alarm Rate Declination Dependence

When a target is located at a large zenith angle, the on-source and off-source event counts are reduced compared to when the target is at a smaller zenith angle. This is a result of the extensive air showers traversing a greater amount of atmosphere to reach the detector. These smaller event counts make the data more likely to be fit well by a constant value, shifting the distribution of $D^{\mathrm{max}}$ to smaller values. Since targets with declinations far from $+19^{\circ}$ spend more time at large zenith angles, we expect that the value of $\gamma$ required to meet the trigger condition of Equation 12 for a fixed false alarm rate will increase for these targets. In this section, we derive a correction to the threshold parameter $\gamma$ that accounts for this effect.

Following the simulation procedure from Section 4.1, we simulate results from nine different declinations, ranging from $-21^{\circ}$ to $+59^{\circ}$, in order to derive the declination dependence of the relationship between $\gamma$ and the false alarm rate. The upper panel of Figure 2 shows the value of $\gamma$ required to obtain a fixed false alarm of 1 event per year and 1 event per century, as a function of declination. The increase in $\gamma$ for declinations below $+4^{\circ}$ and above $+34^{\circ}$ demonstrates that our expectations are borne out by the simulation.

The lower panel of Figure 2 reveals that the ratio between the threshold parameters required to obtain these two false alarm rates does not change as a function of declination. A fit of this ratio to a constant yields a value of 214.2, which describes the data well ($\chi^{2}$ of 7.2 for 8 degrees of freedom). This independence of the ratio on declination allows us to apply a single declination-dependent correction to the threshold parameter that is constant for all false alarm rates. We derive this correction from the results for a false alarm rate of 1 event per year. After this correction is applied, the relationship between $\gamma$ and the false alarm rate follows the curve depicted in Figure 1 for all targets, independent of the target declination.

4.3 Data Quality Cuts

Our implementation of the flare monitor algorithm is sensitive to changes that affect the detector response on time scales shorter than 4 hours, due to our use of the previous direct integration period to provide the off-source counts. Although most detector maintenance leaves the response unchanged, occasionally the response shifts enough to affect the algorithm. PMT rate spikes associated with lightning and similar random events can also affect the detector response. Removing these unstable periods from the data processing is necessary for achieving the same false alarm rate in data that we obtain in simulation. Since the flare monitor operates in real time, this requires a set of data quality cuts that can be applied automatically.

In order to assess the data quality, every 10 seconds we determine the overall event rate, the distribution of event zenith angles, and the two-dimensional angular distribution of events on the sky. We compare the current distributions to the previous distributions, triggering a 2-hour pause in the accumulation of observations into the buffer if the rate changes by $\geq 5$%, or if the two-sample Kolmogorov-Smirnov probability for the angular distributions to be drawn from the same underlying distribution falls below $10^{-5}$. The time interval of 2 hours is chosen in order to correspond to the direct integration time duration. We find that applying these cuts retains 96% of the data, showing that the detector is stable most of the time.

To verify the false alarm rates after applying the data quality cuts, we select a set of 150 regions at the declination of the Crab Nebula. Each region is located at least $2^{\circ}$ away from any known VHE source or monitored target, in order to avoid possible contamination from actual flares. We run the flare monitor algorithm over the data from these regions and compare the resulting false alarm rate to our simulation.

Figure 3 shows the comparison between the false alarm rates in data and simulation. Applying the data quality cuts removes events that would trigger with small threshold parameters. These are events that appear to be highly significant, but actually arise due to detector instabilities. As shown in Figure 3, all events that would trigger with a threshold parameter less than $10^{-10}$ are removed, and the observed false alarm rate matches the rate expected from simulation. In contrast, in the absence of any data quality cuts, the false alarm rate derived from these 150 regions fails to match the simulated false alarm rate, becoming several orders of magnitude too large for small values of $\gamma$. This mismatch reveals that our use of the previous direct integration period does indeed affect the data, necessitating our use of the data quality cuts. Caution should be used in interpreting the error bars in Figure 3 because the results presented are cumulative and therefore the points are correlated.

4.4 Threshold Correction for Strong VHE Targets

The higher analysis bins suffer from low count values. Simulations reveal that this causes the distribution of the test statistic constructed from Equations 10 and 11 to have a weak dependence on the presence of steady VHE emission. Targets with stronger steady emission tend to produce higher values of $D$ than targets with weaker or no steady emission in gamma rays, even in the absence of any flares. Accounting for this effect is necessary to cause the false alarm rates for strong steady VHE emitters to match our expectations from Figure 1.

In order to quantify this effect, we perform simulations following the procedure outlined in Section 4.1. We account for the presence of steady emission by including an additional term in the mean of the Poisson distribution from which the on-source counts are drawn. Instead of using a mean of $\alpha m_{i}$ for this distribution, which models a source region that has only background counts, we use a mean of $\alpha m_{i}+r_{\textnormal{target}}\alpha m_{i}$, where $r_{\mathrm{target}}$ is the long-term relative excess observed at the target location:

[TABLE]

Here, $n_{\mathrm{target}}$ and $m_{\mathrm{target}}$ are the total on-source and off-source counts collected from the vicinity of the strong target over the entire data set, yielding a highly accurate estimate of $r_{\mathrm{target}}$. Thus, the additional term of the relative excess $r_{\mathrm{target}}$ multiplied by $\alpha m_{i}$ represents the excess expected from the source. The values of $n_{\mathrm{target}}$, $m_{\mathrm{target}}$, and $r_{\mathrm{target}}$ depend on the analysis bin.

Figure 4 depicts the cumulative ratio $R_{\mathrm{Crab,sim}}/R_{\mathrm{0,sim}}$ as a function of $\gamma$, where $R_{\mathrm{Crab,sim}}$ is the simulated false alarm rate for a source with steady VHE emission equal to that observed from the Crab Nebula, and $R_{\mathrm{0,sim}}$ is the simulated false alarm rate in the absence of any VHE source as shown in Figure 1. Thus, $R_{\mathrm{Crab,sim}}(\gamma^{\prime})$ and $R_{\mathrm{0,sim}}(\gamma^{\prime})$ are the rates at which false alarms occur when the threshold parameter is set to $\gamma=\gamma^{\prime}$. The open black circles in Figure 4 show that the false alarm rate for a Crab-like source is $\sim$2.5 times higher than what we would expect in the absence of any source. We combat this effect by introducing a target-dependent correction factor $k_{\mathrm{target}}$ to the trigger condition expressed by Equation 12, which now takes the form

[TABLE]

For $k_{\mathrm{target}}>0$, this requires a larger value of $D^{\mathrm{max}}$ in order for the algorithm to trigger than would have been required in the absence of any steady VHE emission. Using a correction factor of $k_{\mathrm{Crab}}=1.2$ for a source with Crab-like VHE emission results in the filled red squares shown in Figure 4. These points demonstrate that a single correction factor suffices to bring the false alarm rate into agreement with simulation for all values of the threshold parameter $\gamma$.

The value of $k_{\mathrm{target}}$ is different for every VHE source. As outlined in Section 6, the only strong sources in our current target list are the blazars Mrk 421 and Mrk 501. These sources are generally weaker than the Crab, and the above simulation procedure produces a value of $k_{\mathrm{Mrk}}=0.2$ for both Mrk 421 and Mrk 501. For all other targets, we use a value of $k_{\mathrm{target}}=0$.

5 Sensitivity of the Flare Monitor

We next investigate the sensitivity of the flare monitor to flares with various properties. The flux, duration, structure, and source declination of the flare can all affect both the probability of detection and the speed with which the flare is detected.

5.1 Flare Simulation Procedure

We characterize the sensitivity of the flare monitor by injecting fiducial flares represented by a constant increase in flux for a duration of 1 hour with the same spectral shape as that of the Crab Nebula as measured by HAWC. We adapt the procedure outlined in Section 4.4, sampling the simulated off-source count $m_{i}^{\prime}$ from a Poisson distribution with mean $m_{i}$ and the simulated on-source count $n_{i}^{\prime}$ before and after the flare from a Poisson distribution with mean $\alpha m_{i}$. During the flare, we adjust the mean of the on-source count distribution to $(1+fr_{\mathrm{Crab}})\alpha m_{i}$, where $r_{\mathrm{Crab}}$ is the observed relative excess of events collected from the vicinity of the Crab Nebula as defined by Equation 13 and $f$ is the flare flux in Crab units. Although the relative excess $r_{\mathrm{Crab}}$ depends on the analysis bin, in our fiducial case the flare flux $f$ does not. The spectrum of the simulated flares therefore differs from that of the Crab Nebula only by the overall constant $f$. We assume no underlying steady state emission for the source, so that the analysis presented here is applicable to targets with little or no VHE emission. The results are equally applicable to strong VHE emitters, provided that correction factors are derived according to the procedures outlined in Sections 4.2 and 4.4.

Similarly to the tuning of the analysis bins referred to in Section 2, we use the measured spectrum of the Crab Nebula as the input flare spectrum. Blazar flare spectra are likely to be harder than the Crab, but they are also affected by EBL absorption. A thorough assessment of the sensitivity of the flare monitor to a range of different intrinsic spectra and source redshifts is beyond the scope of this work. Further studies that address the effects of more realistic flare spectra are currently in progress.

We inject the simulated flares uniformly randomly in time, independent of whether the target spends any time within the field of view throughout the duration of the flare. We consider a flare to be detectable in principle if any part of the flare occurs while the source is at zenith angles smaller than the $45^{\circ}$ zenith angle cut on the events. As a result, we do not expect the probability of detection to converge to 1, since some flares will pass this cut only for a small fraction of their duration.

The flares are separated by an amount of time sufficient to fill the entire 10-hour buffer in order to guarantee that only one flare exists in the buffer at any given time. This procedure gives us the sensitivity to isolated flares. We run the algorithm detailed in Section 3 over the simulated data set. Each time the buffer advances forward in time, we evaluate Equations 10 and 11 and apply the trigger condition expressed by Equation 12 to determine whether the flare is detected. We generate a sufficient number of simulated data, 300 observations of 2 minutes’ duration each, to fill the buffer before the flare starts, and we continue to generate data until the flare has completely exited the buffer. This procedure insures that the probability to detect the flare is unaffected by the presence of other flares.

5.2 Sensitivity to Flare Strength

Figure 5 presents the sensitivity of the flare monitor to fiducial flares of varying fluxes. The upper panel of the figure shows the detection probability as a function of flare flux. The data points are derived from observations at the declination of the Crab Nebula, and the false alarm rate is set to 10 events per year from the entire target collection defined in Section 6. Since the flares are injected at uniformly random times, the detection probability for strong flares converges to a value somewhat less than 1 because flares that occur at the edges of the transits are counted as detectable, even though their lower event counts make their detection less likely. The probability of detecting a given flare reaches $50\%$ at approximately 5.5 Crab units.

The lower panel of Figure 5 shows the time to detection—the amount of time between the start of the flare and the time at which the trigger is issued—as a function of flare flux. The points in the lower panel of the figure show the median time to detection, while the error bars represent the central $68\%$ of the distribution. This distribution is wide partly because flares occurring near target culmination are detected more quickly than flares occurring when the target is at larger zenith angles. Of the flares that are detected, the majority are identified in less than their duration of 1 hour, demonstrating that the method can trigger sufficiently rapidly to enable follow-up observations.

5.3 Sensitivity to Flare Duration

Figure 6 expands the results of Figure 5 to flares of durations other than 1 hour, for flares with fluxes equal to 3, 5, 7, and 10 times the Crab Nebula flux. As before, the upper panel shows the detection probability as a function of flare duration, while the lower panel shows the time to detection. An increase in flare duration improves the detection probability for a fixed flare flux. For instance, whereas Figure 5 shows that the probability of detecting a 3-Crab flare that lasts for 1 hour is $\sim$10%, Figure 6 reveals that this probability rises to $\sim$40% if the duration increases to 5 hours. Most importantly, the lower panel of Figure 6 shows that the time to detection rises more slowly than the flare duration. Thus, long-duration flares provide a larger window of time for other instruments to conduct follow-up observations before the flare terminates.

5.4 Sensitivity to Flare Structure

While our fiducial flare comprises no internal structure aside from a constant elevated flux, many detected VHE flares exhibit a large degree of structure in their light curves. We have tested three different light curve shapes, shown in Figure 7, to ascertain how the flare monitor sensitivity may change in the presence of such structure. In addition to the fiducial square wave flare, we have also tested a flare that rises linearly for half of the flare duration and falls linearly for the other half (hereafter referred to as the linear case), as well as a flare that rises exponentially over the first fifth of the flare duration and then falls exponentially over the remainder (hereafter referred to as the exponential case).

Figure 8 shows the detection probability for the fiducial flare as a function of both flare flux and flare duration. In the figure, the flare flux is the average flux over the entire flare duration, which in the fiducial case is equivalent to the constant flux of the flare. As expected, increasing either the flux or the duration improves the detection probability.

The detection probability for the linear flare structure appears in Figure 9. As in the fiducial case, the flare flux is the average over the duration of the flare, which in the linear case is half of the peak flux. Points in Figure 9 therefore represent points with the same fluence—total flux integrated over the duration of the flare—as the corresponding points in Figure 8.

In the linear case, the probability tends to peak for flares of around 1.5 hours’ duration and gradually decreases for flares that last longer than this. This trend is most apparent for the strongest flares, and can be understood as an effect of the source transit. Since we inject flares uniformly randomly in time, a fraction of flares will occur near the edges of the transit, as the source enters or leaves the field of view. We therefore sample only the beginning or end of these truncated flares. For such flares, a smaller fraction of the total fluence from linear flares appears in the field of view than in the fiducial case. When the flare duration is small compared to the transit duration of $\sim$6 hours, the fraction of truncated flares is negligible. However, as the flare duration increases, a larger fraction of flares are subjected to this effect, which manifests itself as a reduction in the probability of detection for a fixed total fluence. The white contour in Figure 9 shows the point at which the probability of flare detection in the linear case drops below that of the fiducial case as a result of this effect. Overall, the linear structure does not reduce the detection probability by more than $10\%$, and for weak or short flares the structure increases the detection probability by as much as a factor of 2. This suggests that light-curve structure, which is expected to be present in most flares, will enhance the detection probability in many cases.

We have also tested the sensitivity of our method to flares with the exponential profile from Figure 7. The results are very similar to those for the linear case and we consequently do not discuss them here.

5.5 Sensitivity to Flare Source Declination

The sensitivity depends strongly on the declination of the target, primarily due to the reduced number of events at large zenith angles, which renders discrimination of elevated flux states from the background rate increasingly difficult. Figure 10 depicts the sensitivity to flares from targets at a range of declinations. The figure is constructed for the same false alarm rate that was used to construct Figure 5. This rate accounts for the declination-dependent threshold parameter correction shown in Figure 2. The asymmetry apparent in Figure 10 arises due to HAWC being located in the northern hemisphere. At declinations of $-16^{\circ}$ and $+54^{\circ}$, located $35^{\circ}$ away from the latitude of the HAWC observatory, the probability of detecting a flare of 10 Crab units is small but non-negligible. We monitor targets at a maximum of $40^{\circ}$ away from the HAWC latitude in order to provide sensitivity to these extreme events. Since we cut events with zenith angles above $45^{\circ}$, these sources are monitored for much less time per day than sources located nearer to the HAWC latitude.

6 Target Selection

In the ideal case, we would be able to tile the sky with targets for the search algorithm and perform all-sky monitoring for VHE transients. However, achieving an acceptably low false alarm rate under this approach requires us to lower the value of the threshold parameter $\gamma$, in turn lowering the sensitivity of the search. For this reason, we concentrate on identifying a specific list of the most promising candidates, reserving an all-sky search for future work.

We select targets from two source catalogs: the online catalog of known VHE sources, TeVCat111See http://tevcat.uchicago.edu. (Wakely & Horan, 2008); and the second catalog of hard-spectrum Fermi-LAT sources, or 2FHL (Ackermann et al., 2016). The former catalog comprises sources that have been observed to produce VHE emission, while the latter represents sources strongly detected by the Fermi-LAT above 50 GeV that consequently have a large probability to produce VHE-band emission. Although some Galactic sources, such as binaries, are known to flare in the VHE band (Aharonian et al., 2005; Albert et al., 2006), for this work we restrict our attention to extragalactic objects because we regard these as being far more likely to produce strong VHE transients.

Both catalogs also contain several unidentified sources with no known or obvious counterpart at other wavelengths. Since the TeVCat unidentified sources cluster along the Galactic Plane, we exclude them from our target list as being likely Galactic in origin. In contrast, the contribution of extragalactic sources to the 2FHL unidentified source population is likely $\sim$80% (Ackermann et al., 2016), so we include these objects in our analysis. We also exclude targets culminating at a zenith angle of greater than $40^{\circ}$ in order to limit events to the region where HAWC is sensitive and the gamma-hadron response is stable. For the latitude of HAWC, we therefore accept target locations between $-21^{\circ}$ and $+59^{\circ}$ in declination.

After these cuts, our final selection comprises a total of 187 targets: 46 TeVCat sources and 141 additional 2FHL sources that are not in TeVCat. To enhance the algorithm’s sensitivity to targets that are more likely to produce transients, we divide the targets into four categories, defined in Table 2. We split the TeVCat sources into two categories: one containing the very strong VHE emitters Mrk 421 and Mrk 501, and another containing everything else. We also split the 2FHL sources into two categories: one comprising sources with measured redshifts less than 0.3, because these sources are less likely to be affected by gamma-ray attenuation by the EBL; and a separate category covering more distant sources, sources without a measured redshift, and the unidentified sources. We then choose the threshold parameter $\gamma$ independently for each category, such that each category contributes 1/4 of the total expected false alarm rate. Thus, we choose a less restrictive value of $\gamma$, for instance, for the nearby 2FHL targets, reflecting our opinion that these targets are more likely to produce flares than the rest of the 2FHL sources, and their smaller number allows us to increase our sensitivity to them. The specific value of the threshold parameter chosen for each category appears in the third column of Table 2. Of course, we could have chosen to divide the false alarm rate contributions differently among the target categories, but since this is degenerate with the definitions of the target categories themselves, our arbitrary choice of an even division is unimportant.

7 Verification of the Method with Data

As a verification of the search method, we process archival HAWC data collected between 2014 November 26 and 2016 June 1 for Mrk 421 and Mrk 501, which are known to produce frequent strong flares. In this data set, we identify six events surpassing the threshold that would have been set to achieve a false alarm rate of 10 events per year for the entire source selection described in Section 6. This section first takes a detailed look at the flare monitor’s view of the strongest of these events, which is associated with Mrk 501. We then summarize the rest of the triggers, reserving a thorough investigation of these events and the remainder of the targets for future work.

7.1 Markarian 501 Flare on 2015 August 18

The strongest of the six events identified by the flare monitor in the Mrk 421 & Mrk 501 target class is associated with Mrk 501. It occurred on 2015 August 18, at modified Julian day (MJD) 57252.062. This event crossed the trigger threshold for the Markarian target class after 34 minutes. Figure 11 shows the equivalent false alarm rate of this event as a function of the age of the flare, which is the amount of time after the trigger time MJD 57252.062. The minimum equivalent false alarm rate occurs 54 minutes after the identified start time of the event, achieving a value of $6.8\times 10^{-3}$ events per year. In other words, in the absence of any flares, we would expect to see an event at least this extreme only once every $\sim$150 years from a single source, or once every $\sim$75 years from the Mrk 421 & Mrk 501 target class.

Figure 12 shows the light curve for the flare in terms of excess counts above background for each of the ten analysis bins. The start of the flare as determined by the flare monitor search algorithm appears as a dashed red vertical line in the figure. Analysis bins 3, 4, 5, 7, and 9 dominate the contributions to the test statistic $D$ from Equations 10 and 11, together accounting for 96% of the total. This assessment of $D$ uses all observations in the buffer, including data from the previous transit of the source. The inclusion of these data produces a more precise estimate of the background event rate than if we used only data collected over a shorter period. This procedure improves the algorithm’s determination of the importance of the higher analysis bins where the events are sparse. Our confidence in this event being a real flare increases because the contributions to $D$ come from several analysis bins and are not dominated by a single bin.

As a consistency check, we also run the Bayesian block algorithm of Scargle et al. (2013), which among other applications is commonly used to identify significantly distinct states in astrophysical light curves (e.g. Barrière et al., 2014; Aartsen et al., 2015; Ackermann et al., 2016). The advantages of this algorithm include insensitivity to gaps in the data taking and, for a given set of $N$ observations, the capacity to identify the optimum partition of the data in $N^{2}$ steps. Because of this scaling, it is possible to run the Bayesian block algorithm over the few thousand 2-minute points that fit into several days, but it is computationally prohibitive to run the algorithm over the complete data set of $\sim$1.5 years with 2-minute time resolution.

Following Section 3 of Scargle et al. (2013), we use the full likelihood ratios constructed by Equations 8 and 9 to identify regions of the light curve where the flux from the source is well represented by a constant value. Our threshold parameter $\gamma$ from Equation 12 naturally becomes the prior parameter required by the Bayesian block algorithm. For a given set of observations, $\gamma$ controls the probability that the Bayesian block algorithm will produce a single false trigger.

We select four separate time windows centered on the identified start point of the flare, with durations of 2, 4, 6, and 8 days. For each separate time window, we run the Bayesian block algorithm with a value of $\gamma$ such that the probability for a false trigger is 5%. We choose a relatively high probability because our goal here is to identify the duration of a flare that has already been detected by an independent method, not to quantify the significance of our observation.

All of the windows produce a strong detection of a flare beginning 88 $\pm$ 2 minutes after the start of MJD 57252 and ending 74 $\pm$ 2 minutes later. We also run the Bayesian block algorithm with a false trigger rate of 1%. In this more restrictive search, only the 6-day and 8-day time windows detect the flare. It is likely that the shorter time windows do not contain sufficient information in the region before and after the flare to overcome the stricter trigger condition.

The consistency between the various windows gives us confidence that the trigger is not due to a peculiarity of the data on MJD 57252, but is a real flare that does not depend on where the edges of the data selection occur. The gray band in Figure 12 shows the block that is selected by the Bayesian block algorithm for the 8-day time window. The start of this region lines up well with the starting point identified by the faster but more limited algorithm used by the flare monitor.

Figure 13 shows the relative excess as a function of analysis bin for the data collected from Mrk 501 during the flare identified by the Bayesian block algorithm. The long-term relative excess from the Crab Nebula taken over the entire data set is shown for comparison, and the long-term relative excess from Mrk 501 over the entire data set also appears in the figure. These long-term relative excesses have been computed from a weighted average that matches the zenith angle distribution to the observed distribution from Mrk 501 during the flare. Furthermore, the background used to compute the relative excess during the flare has been derived from zenith-matched data collected 2.5 days before and after the flare. This procedure was necessary to insure that a sufficient number of background events were collected to produce an accurate determination of the relative excess in the higher analysis bins. In addition to deriving the background counts from 2.5 days before and after the flare, we have also checked the results using 10 days before and after the flare and the full data set. We have also derived the background from the time of the flare alone. We find good agreement between the relative excesses for all analysis bins in each case.

As discussed in Section 2, the analysis bins do not correspond directly to energy, but rather to the amount of information available about the events. Deriving a spectrum requires the detailed analysis techniques presented in Abeysekara et al. (2017b) that include the HAWC detector response and account for the energy overlap and zenith angle response of the analysis bins. Consequently, the temptation to interpret the data in Figure 13 in terms of energy spectra should be avoided, as such interpretations will not produce reliable results. Rather, we include Figure 13 to demonstrate that the detected flare is substantially stronger than both the average Mrk 501 flux and the Crab Nebula flux.

Since the flare monitor is sensitive to multi-Crab level flares lasting on the scale of hours, we expect the detected flare to surpass both the flux from the Crab and the long-term average flux from Mrk 501. We obtain a rough estimate of the flux of the flare in Crab units by fitting a constant value to the ratio of the relative excess from Mrk 501 during the flare to that from the Crab Nebula over the full data set. This estimate is only approximate because the analysis bins overlap substantially in energy and the spectrum of the flare may differ from the Crab spectrum. Nevertheless, the errors on the relative excess ratio are large enough to admit a good fit with a constant value, as shown in Figure 14. The flare flux according to the fit is 3.8 $\pm$ 1.1 Crab units, in accord with expectations based on the sensitivity of the method presented in Figure 5.

7.2 Summary of Mrk 421 and Mrk 501 Triggers

In the data set considered in this work, we expect 3.8 events surpassing our trigger threshold in the Mrk 421 & Mrk 501 target class. We observe a total of six events: four from Mrk 421 and two from Mrk 501. Three of these events have large equivalent false alarm rates and appear consistent with background fluctuations. Here, we consider the three other events, which all have false alarm rates much lower than we would expect from background alone. We designate these as high-confidence triggers from the Markarian source class.

For each of the three high-confidence triggers, we run the Bayesian block algorithm with an 8-day window in order to estimate the flux following the methods of Section 7.1. If the flare is not identified in the Bayesian block algorithm with a 1% false trigger rate, we run the algorithm with a more permissive false trigger rate of 5%. Since our goal is to identify the properties of events that have already been detected by an independent algorithm and not to quantify the significance of our detections, we are not concerned with statistical trials in this case.

Table 3 summarizes the properties of the three high-confidence events. The columns in Table 3 show: (1) a convenient reference number for the event, (2) the equivalent false alarm rate for the event at its most significant detection, (3) whether the event is associated with Mrk 421 or Mrk 501, (4) the MJD of the trigger time, (5) the time to detection for the event, (6) whether the event is identified in the Bayesian block search run with an 8-day window at the 1% false alarm or 5% false alarm level, (7) the duration of the flare as identified by the Bayesian block search, and (8) the approximate flux of the flare as determined by the methods of Section 7.1.

Event 1 in Table 3 is the Mrk 501 flare from Section 7.1 and has the lowest equivalent false alarm rate. Events 2 and 3 are associated with Mrk 421 and Mrk 501 respectively. These events have false alarm rates in the range of 2 to 3 events per century, and their estimated durations and fluxes in Crab units are in agreement with the expectations from Section 5. Furthermore, their times to detection are smaller than their durations. Thus all three of these events would have been detected by the real-time flare monitor sufficiently rapidly to issue alerts for follow-up observations, had the flare monitor been operational at the time.

The HAWC collaboration recently reported on a separate search of Mrk 421 and Mrk 501 for flares using the standard HAWC analysis with a base time scale of 1 transit (Albert et al., 2017). This search is sensitive to weaker flares but lacks the rapid response and flexible time scale search of the real-time flare monitor. In Albert et al. (2017), several periods of significant enhanced activity appear in both sources. The three high-confidence events identified by the present work all correspond to regions of elevated flux according to the standard HAWC analysis, providing further support for the detections presented in this work. It is plausible that the real-time flare monitor, which is dedicated to finding short, strong flares, fails to detect the other flares identified by the single-transit light curves because they are weaker in flux. Flares that rise slowly over the course of several days would also elude the present analysis due to its focus on short time scales. Further assessment of the consistency of the two methods is ongoing.

Event 1 in Table 3 occurred during quasi-simultaneous observations from the FACT experiment (Anderhub et al., 2013). The publicly available light curves from FACT222See the quick-look analysis at http://fact-project.org/monitoring. show a clear increase in the excess rate from Mrk 501 on the same night as event 1. This excess is significantly larger than the long-term FACT excess measured from Mrk 501, as well as the FACT excess measured on the surrounding nights. The FACT observations took place between MJD 57251.87 and MJD 57251.96, approximately 2.6 hours before the flare detected by HAWC at MJD 57252.068. Owing to the short delay between the flare detections by the two instruments, it is likely that the increased activity evident in the FACT observations is related to the subsequent outburst detected by HAWC a few hours later.

Near the time of event 2, which is associated with Mrk 421, FACT observations are available on MJD 57015 and 57021, with the latter spanning the range from MJD 57021.14 to 57021.29. These observations occurred approximately 17 hours after the end of event 2 as identified by HAWC. Nevertheless, the available FACT monitoring data show that Mrk 421 was in a relatively high state compared to its long-term average during these dates, in agreement with the daily HAWC observations reported by Albert et al. (2017). The HAWC event detected on MJD 57020 may be related to this high state.

Event 3 occurred when the source was not visible to FACT during nighttime observations. As a result, FACT monitoring data are unavailable for event 3.

8 Flare Monitor Outlook

The HAWC real-time flare monitor has been fully operational since 2017 January, when it began sending alerts to the VERITAS experiment. Alerts to the Fermi-LAT began in 2017 March, and additional instruments are currently being added. Ultimately, we aim to provide public alerts via the Gamma-ray Coordinates Network (Barthelmy et al., 1998).

This work has explored the capability of the HAWC real-time flare monitor to detect VHE flares at the level of several times the flux from the Crab Nebula in under 1 hour. We have demonstrated this capability via the detection in archival HAWC data of three high-confidence flares from Mrk 421 and Mrk 501. These observations are consistent with fluxes several times that from the Crab Nebula and with durations ranging from one to a few hours. We have also shown that the false alarm rate of the method matches the expectation from simulation, allowing us confidently to assign an equivalent false alarm rate for each detection.

The HAWC real-time flare monitor will provide a unique opportunity to study an unbiased collection of extreme VHE flares in order to unravel further the mysteries surrounding blazars, their VHE emission, and their jets. Of particular interest are the nature and location of the dissipation region, where the flares are produced. Understanding the particle populations and acceleration mechanisms that produce the flares, investigating the possible correlation of X-ray and gamma-ray fluxes during the flares, and determining the causes of apparent orphan VHE flares all proceed from the detection of these events. By helping to identify the most extreme flares and the shortest time scales on which they occur, the HAWC real-time flare monitor will facilitate coordinated observations seeking new clues for addressing these and other important topics of interest.

We acknowledge the support from: the US National Science Foundation (NSF); the US Department of Energy Office of High-Energy Physics; the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory; Consejo Nacional de Ciencia y Tecnología (CONACyT), México (grants 271051, 232656, 260378, 179588, 239762, 254964, 271737, 258865, 243290, 132197), Laboratorio Nacional HAWC de rayos gamma; L’OREAL Fellowship for Women in Science 2014; Red HAWC, México; DGAPA-UNAM (grants IG100317, IN111315, IN111716-3, IA102715, 109916, IA102917); VIEP-BUAP; PIFI 2012, 2013, PROFOCIE 2014, 2015; the University of Wisconsin Alumni Research Foundation; the Institute of Geophysics, Planetary Physics, and Signatures at Los Alamos National Laboratory; Polish Science Centre grant DEC-2014/13/B/ST9/945; Coordinación de la Investigación Científica de la Universidad Michoacana. Thanks to Luciano Díaz and Eduardo Murrieta for technical support. HAWC