Observation of the Crab Nebula with the HAWC Gamma-Ray Observatory

The HAWC instrument is modeled using a combination of community-standard simulation packages and custom software. The CORSIKA package (v7.4000) is used for simulation of air showers, propagating the primary particles through the atmosphere to the ground (Heck et al. 1998). At ground level, a GEANT 4 (Agostinelli et al. 2003) simulation (v4.10.00) of the shower particles is used to propagate the ground-level particles through the HAWC tanks and to track the Cherenkov photons to the faces of the PMTs.

The response of the PMTs and the calibration are approximated with a custom simulation that assumes that recorded light is faithfully detected with some efficiency and an uncertainty in the logarithm of the total charge recorded. Decorrelated single PE noise is added. The absolute PMT efficiency for detecting Cherenkov photons is established by scaling the simulated PMT response to vertical muons to match data. Most muons passing through HAWC are minimum ionizing with nearly constant energy loss. Vertical muons, therefore, are a nearly constant light source and are convenient for establishing to total PMT efficiency. Simulated events are subsequently reconstructed by the same procedure as experimental data to study the performance of the algorithms.

Table 1 summarizes the steps in reconstruction of HAWC events. The first step in the reconstruction process is calibration, the processes by which true time and light level in each PMT are estimated from the TDC-measured threshold-crossing times of each PMT (Lauer 2013; Solares et al. 2016).

Table 1.  Steps in the HAWC Event Reconstruction

Step	Description	Hit Selection
1	Calibration	⋯
2	Hit selection	⋯
3	Center-of-mass core reconstruction	Selected hits
4	SFCF core—first pass	Selected hits
5	Direction—first pass	Selected hits
6	SFCF core—second pass	Selected hits within 50 ns of first-pass plane
7	Direction—second pass	Selected hits within 50 ns of first-pass plane
8	Compactness, ${ \mathcal C }$	Selected hits within 20 ns of second-pass plane
9	PINCness, ${ \mathcal P }$	Selected hits within 20 ns of second-pass plane

Note. After calibration (Section 2.2) and the selection of hits (Section 2.3), the reconstruction is applied. The best shower core (Section 2.4) and direction (Section 2.5) reconstructions are applied with gradual narrowing of the hits used (Section 2.7). After the final core and direction determination, the photon/hadron discrimination algorithms, compactness, and PINCness are applied (Section 2.6).

Download table as:  ASCIITypeset image

The calibration associates the measured ToT in each PMT with the true number of PEs. To give a sense of scale, the ToT for a single PE crossing the low-threshold discriminator (about 1/4 PE) is ∼100 ns. Above a few PEs, the higher-threshold (about 4 PEs) ToT is used for charge assignment, and a high-threshold ToT of 400 ns roughly corresponds to a charge of 10⁴ PEs. The timescale for these ToTs is determined by the shaping of the front-end electronics and is chosen to be longer than the characteristic arrival time distribution of PEs during an air shower so as to integrate the whole air shower arrival into one PMT hit.

In addition to the PE measurement, the calibration system accounts for electronic slewing of the PMT waveform: lower-PE waveforms cross the threshold later than contemporaneous high-PE waveforms.

Subsequent reconstruction algorithms treat all PMTs, the 8 inch and 10 inch, as identical, despite the larger size and greater efficiency of the 10-inch PMTs. To accommodate this, an effective charge Q_eff is defined. For Q_eff, PE values from the central 10-inch PMTs are scaled by a factor of 0.46 to place them on par with the 8-inch PMTs.

Finally, each PMT has a single calibrated timing offset that accounts for the different cable lengths and any other timing delays that may differ from PMT to PMT. These delays are measured to within a few nanoseconds by the calibration system and are refined to subnanosecond precision by iterated fits to hadronic air showers. Since the hadronic background is isotropic, the point of maximum cosmic-ray density is overhead and the PMT timing pedestals are chosen to ensure that this is true. A final small (∼02) rotation of all events is performed to ensure that the Crab Nebula appears in its known location.

Throughout the analysis, the Crab is assumed to be at a location of 8363 right ascension and 2201 declination, in the J2000.0 epoch, taken from Aharonian et al. (2004). While the pulsar position is known more precisely (e.g., Comella et al. 1969), this precision is sufficient for use in HAWC.

As described in Section 1, the HAWC Data Acquisition System (DAQ) records 1.5 μs of data from all PMTs that have a hit during an air shower event. A subset of these hits are selected for the air shower fit. To be used for the air shower fit, hits must be found between −150 and +400 ns around the trigger time. Hits are removed if they occur shortly after a high-charge hit under the assumption that these hits are likely contaminated with afterpulses. Additionally, hits are removed if they have a pattern of TDC crossings that is not characteristic of real light; they cannot be calibrated accurately. Finally, each channel has an individual maximum calibrated charge, typically a few thousand PEs, but no more than 10⁴ PEs, above which the PMTs are not used. Above ∼10⁴ PEs, corresponding to a ToT of ∼400 ns, prompt afterpulsing in the PMTs can artificially lengthen the ToT measurement, giving a false measurement. Channels are considered available for reconstruction if they have a live PMT taking data that have not been removed by one of these cuts.

The angular error and the ability to distinguish photon events from hadron events are strongly dependent on the energy and size of events on the ground. We adopt analysis cuts and an angular resolution description that depends on this measured size. The data are divided into nine size bins, ${ \mathcal B }$, as outlined in Table 2. The size of the event is defined as the ratio of the number of PMT hits used by the event reconstruction to the total number of PMTs available for reconstruction, f_hit. This definition allows for relative stability of the binning when PMTs are occasionally taken out of service.

Table 2.  Cuts Used for the Analysis

${ \mathcal B }$	fhit	${\psi }_{68}(\deg )$	${ \mathcal P }$ Maximum	${ \mathcal C }$ Minimum	Crab Excess per Transit
1	6.7%–10.5%	1.03	<2.2	>7.0	68.4 ± 5.0
2	10.5%–16.2%	0.69	3.0	9.0	51.7 ± 1.9
3	16.2%–24.7%	0.50	2.3	11.0	27.9 ± 0.8
4	24.7%–35.6%	0.39	1.9	15.0	10.58 ± 0.26
5	35.6%–48.5%	0.30	1.9	18.0	4.62 ± 0.13
6	48.5%–61.8%	0.28	1.7	17.0	1.783 ± 0.072
7	61.8%–74.0%	0.22	1.8	15.0	1.024 ± 0.053
8	74.0%–84.0%	0.20	1.8	15.0	0.433 ± 0.033
9	84.0%–100.0%	0.17	1.6	3.0	0.407 ± 0.032

Note. The definition of the size bin ${ \mathcal B }$ is given by the fraction of available PMTs, f_hit, that record light during the event. Larger events are reconstructed better, and ψ₆₈, the angular bin that contains 68% of the events, reduces dramatically for larger events. The parameters ${ \mathcal P }$ and ${ \mathcal C }$ (Section 2.6) characterize the charge topology and are used to remove hadronic air shower events. Events with a ${ \mathcal P }$ less than indicated and a ${ \mathcal C }$ greater than indicated are considered photon candidates. The cuts are established by optimizing the statistical significance of the crab and trend toward harder cuts at larger size events. The number of excess events from the crab in each ${ \mathcal B }$ bin per transit is shown as well.

Download table as:  ASCIITypeset image

For this analysis, events are only used if they have more than 6.7% of the available PMTs' seeing light. Since typically 1000 PMTs are available, typically a minimum of 70 PMTs is needed for an event. This is substantially higher than the trigger threshold. The data between the trigger threshold and the threshold for ${ \mathcal B }=1$ in this analysis consist of real air showers, and techniques to recover these events and lower the energy threshold, beyond what is presented here, are under study.

Figure 2 shows the distribution of true energies as a function of the ${ \mathcal B }$ of the events. The distribution of energies naturally depends heavily on the source itself, both its spectrum and the angle at which it culminates during its transit. A pure power-law spectrum with a shape of ${E}^{-2.63}$ and a declination of 20° was assumed for this figure. As ${ \mathcal B }$ is a simple variable—containing no correction for zenith angle, impact position, or light level in the event—the energy distribution of ${ \mathcal B }$ bins is wide. Section 5.3 discusses planned improvements to this event parameter that will measure the energy of astrophysical gamma rays better.

Zoom In Zoom Out Reset image size

Bin ${ \mathcal B }=9$ bears particular attention. It is an "overflow" bin containing events that have between 84% and 100% of the PMTs in the detector seeing light. Typically, a 10 TeV photon will hit nearly every sensor, and the ${ \mathcal B }$ variable has no dynamic range above this energy. This limit is not intrinsic to HAWC, and variables that utilize the light level seen in PMTs on the ground, similar to what was used in the original sensitivity study (Abeysekara et al. 2013), have a dynamic range above 100 TeV. These variables, not used in this analysis, will improve the identification of high-energy events. This is discussed further in Section 5.3.

In an air shower, the concentration of secondary particles is highest along the trajectory of the original primary particle, termed the air shower core. Determining the position of the core on the ground is key to reconstructing the direction of the primary particle. In the sample event, Figure 3, the air shower core is evident in Figure 3(a). The image is an overhead view of the HAWC detector, with circles indicating the WCD location and the PMTs within the WCDs. The colors indicate the amount of light (measured in units of PEs) seen in each PMT. The air shower core is evident as the point of maximum PE density.

Zoom In Zoom Out Reset image size

The PE distribution on the ground is fit with a function that decreases monotonically with the distance from the shower core. The signal in the ith PMT, S_i, is presumed to be

$$\begin{eqnarray}{S}_{i} & = & S(A,{\boldsymbol{x}},{{\boldsymbol{x}}}_{i})\\ & = & A\left(\displaystyle \frac{1}{2\pi {\sigma }^{2}}{e}^{-| {{\boldsymbol{x}}}_{i}-{\boldsymbol{x}}{| }^{2}/2{\sigma }^{2}}+\displaystyle \frac{N}{{(0.5+| {{\boldsymbol{x}}}_{i}-{\boldsymbol{x}}| /{R}_{m})}^{3}}\right),\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{x}}$ is the core location, ${{\boldsymbol{x}}}_{i}$ is the location of the measurement, R_m is the Molière radius of the atmosphere, approximately 120 m at HAWC altitude, σ is the width of the Gaussian, and N is the normalization of the tail. Fixed values of σ = 10 m and N = 5 × 10⁻⁵ are used. This leaves three free parameters, the core location and overall amplitude, A.

The functional form used in this algorithm, termed the Super Fast Core Fit (SFCF), is a simplification of a modified Nishimura–Kamata–Greisen (NKG) function (Greisen 1960) and is chosen for rapid fitting of air shower cores. The NKG function has an additional free parameter, the shower age, and involves computationally intensive power-law and gamma function evaluation. The SFCF hypothesis in Equation (1) is similar, but numerical minimization can converge faster because the function is simpler, the derivatives are computed analytically, and there is a lack of a pole at the core location. This fit is seeded with a simple center-of-mass algorithm that computes the center of mass from the recorded charge in each PMT.

Figure 3(b) shows the recorded charge in each PMT as a function of the PMT's distance along the ground to the reconstructed shower core. While the full NKG function would describe the lateral distribution better, the SFCF form allows rapid identification of the center of showers, and this is sufficient for the present analysis. Cores can be localized to a median error of ∼2 m for large events (${ \mathcal B }=8$) and ∼4 m for small events (${ \mathcal B }=3$) for gamma-ray events with a core that lands on the HAWC detector. The error in reconstructing the shower core increases as the core moves farther from the array. For example, a shower with a core that is 50 m from the edge of the array will have an error in the location of the core of ∼35 m.

To first order, the air shower particles arrive on a plane defined by the speed of light and direction of the primary particle. In fact, the shower front is curved, since the shower particles originate from a common interaction point. Particles with paths along the shower axis have shorter travel distances compared with those that stray from it owing simply to their shorter path. In practice, this means that particles detected near the shower core are observed to arrive earlier than particles far from the core.

Multiple scattering of shower particles also leads to variation of arrival times at the observation level. Generally, the shower front is several nanoseconds thick near the core and wider far from the core. Since the HAWC electronics only measures the arrival time of the earliest arriving particle, there is an unavoidable pulse-amplitude-dependent timing correction. If only a single particle is detected, its most probable time is at the center of the shower front, but when many particles are detected, the measured time is the time of the first particle, which is likely near the leading edge of the shower front. Since particle density is higher near the core than in the shower periphery, there is a resulting position-dependent shower timing offset. We refer to this phenomenon as sampling. Generally, it is hard to disentangle curvature and sampling, so in the reconstruction they are regarded as a single effect.

The curved shape of the air shower front can be easily seen in data with the sample event in Figure 3. Figure 3(c) shows, for each PMT in the sample event, the calibrated time each PMT saw light. The color trend is due to the inclined direction of the air shower. Taking out this inclination, we can see the curved shower front. The event in Figure 3 was chosen for display because—taken from a high signal-to-background sample of data—it is very likely a photon from the Crab. If we assert that the origin of this particle is the Crab Nebula, we know the air shower plane precisely and can make Figure 3(d): we adjust the times of each PMT hit assuming a pure planar air shower from the direction of the Crab Nebula. Figure 3(d) shows the plane-corrected time as a function of the PMT distance from the core of the air shower along the ground. A pure planar shower would correspond to a horizontal line on this figure. The delay of particles far from the core is evident. Note that the deviation as a function of distance is about 5–10 ns/100 m, more than a degree of deviation from a plane. Understanding and correcting for shower sampling and curvature is the main challenge in achieving the best possible angular resolution of 0.15 deg.

A combined curvature/sampling correction—a function of the distance of hits from the shower core and the total charge recorded in the PMT—is utilized to account for the deviation of the front shape compared to a plane. The curvature/sampling correction is based on a combination of simulation and Crab observations. The functional form of the correction was derived from studies of simulated showers, but the parameters of the function were optimized using a subset of our Crab data set and then tested using a separate subset. We use data instead of simulations to optimize the shower front shape, since the subtitles of the timing calibration methodology and PMT and electronics response are not well reproduced by our simulation.

After correcting for the sampling and curvature, the angular fit is a simple χ² planar fit and has been described before (Atkins et al. 2003).

Hadronic cosmic rays are the most abundant particles producing air showers in HAWC and constitute the chief background to high-energy photon observation. The air showers produced by high-energy cosmic rays and gamma rays differ: gamma-ray showers are pure electromagnetic showers with few muons or pions. Conversely, hadronic cosmic rays produce hadronic showers rich with pions, muons, other hadronic secondaries, and structure due to the showering of daughter particles created with high transverse momentum. In HAWC, these two types of showers appear quite different, particularly for showers above several TeV.

Figure 4 shows the lateral distributions for two showers, an obvious cosmic ray (left) and a strong photon candidate (right) from the Crab Nebula. The effective light level Q_eff falls off for hits farther from the shower core in both showers, but in the hadronic shower there are sporadic high-charge hits far from the air shower's center. This clumpiness is characteristic of hadronic showers and arises from a combination of penetrating particles (primarily muons) and hadronic subshowers that are largely absent in photon-induced showers.

Zoom In Zoom Out Reset image size

Two parameters are used to identify cosmic-ray events. The first parameter, compactness, was used in the sensitivity study (Abeysekara et al. 2013). The variable CxPE₄₀ is the effective charge measured in the PMT with the largest effective charge outside a radius of 40 m from the shower core. We then define the compactness, ${ \mathcal C }$, as

$$\begin{eqnarray}&&{ \mathcal C }=\displaystyle \frac{{N}_{\mathrm{hit}}}{{\mathrm{CxPE}}_{40}},\end{eqnarray} \tag{ 2 }$$

where N_hit is the number of PMT hits during the air shower. CxPE₄₀ is typically large for a hadronic event, so ${ \mathcal C }$ is small.

In addition to the largest hit outside the core, the "clumpiness" of the air shower is quantified with a parameter ${ \mathcal P }$, termed the PINCness of an event (short for Parameter for Identifying Nuclear Cosmic rays). ${ \mathcal P }$ is defined using the lateral distribution function of the air shower, seen in Figure 4. Each of the PMT hits, i, has a measured effective charge ${Q}_{\mathrm{eff},i}$. ${ \mathcal P }$ is computed using the logarithm of this charge ${\zeta }_{i}={\mathrm{log}}_{10}({Q}_{\mathrm{eff},i})$. For each hit, an expectation is assigned $\langle {\zeta }_{i}\rangle$ by averaging the ζ_i in all PMTs contained in an annulus containing the hit, with a width of 5 m, centered at the core of the air shower.

${ \mathcal P }$ is then calculated using the χ² formula:

$$\begin{eqnarray}&&{ \mathcal P }=\displaystyle \frac{1}{N}\displaystyle \sum _{i=0}^{N}\,\displaystyle \frac{{({\zeta }_{i}-\langle {\zeta }_{i}\rangle )}^{2}}{{{\sigma }_{{\zeta }_{i}}}^{2}}.\end{eqnarray} \tag{ 3 }$$

The errors ${\sigma }_{{\zeta }_{i}}$ are assigned from a study of a sample of strong gamma-ray candidates in the vicinity of the Crab.

The ${ \mathcal P }$ variable essentially requires axial smoothness. Figure 4 shows the moving average $\langle {\zeta }_{i}\rangle$ for two sample events. The hadronic event in Figure 4 is "clumpy" and has several hits that differ sharply from the moving average, yielding a large ${ \mathcal P }.$

The ${ \mathcal C }$ and ${ \mathcal P }$ variables are well modeled in simulation. Figure 5 shows the measured distribution of ${ \mathcal C }$ in the vicinity of the Crab Nebula (the Crab region) and in an annular reference region around the Crab (the background region). The background region is scaled to have the same solid angle as the Crab region. The distributions in the vicinity of the Crab are made of a combination of hadronic cosmic rays and true photons from the Crab. Figures 5(c) and (d) show these distributions with the background distribution subtracted. The subtraction yields the data-measured distribution of ${ \mathcal C }$ for gamma rays from the Crab. Figure 6 is a comparable figure for ${ \mathcal P }$. Figures 5 and 6 are compared to a simulation prediction from the final fitted flux from Section 4; the simulation agreement is evident.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

HAWC's outer 8-inch PMTs individually trigger at some 20–30 kHz, and the 10-inch central PMTs trigger at 40–50 kHz. Of this random noise, roughly half (with large uncertainties) is believed to be due to real shower fragments and roughly half due to nonshower sources like radioactive contaminants in the PMT glass or PMT afterpulsing. Approximately 1 kHz of the noise is "high-PE" noise from individual accidental air shower muons (10–200 PEs) and can be correlated between the PMTs within a WCD that the muon hits.

This noise can bias the air shower reconstruction, and the muon noise has the potential, if not removed, to confuse the identification of gamma-ray showers because a single muon is enough to indicate that a shower is of hadronic origin.

In order to achieve the best direction reconstruction and to avoid falsely rejecting true photons, the SFCF core reconstruction and the plane fit are each performed twice, as outlined in Table 1. During the first pass, all selected hits are used to locate the core and initial direction. After this first "rough" fit, hits that are more than ±50 ns from the curvature/sampling-corrected air shower plane are removed, and the shower is fit a second time.

The computation of photon/hadron separation variables is done with hits that are within ±20 ns of the curvature/sampling-corrected air shower plane. With these cuts, only ∼4% of gamma-ray events will have an accidental muon contributing to the photon/hadron variable computation and risk being falsely rejected.

We consider here data taken by HAWC between 2014 November 26 and 2016 June 2, a total elapsed time of 553 days. The detector was not taking data for a cumulative time of 40 days for various operational reasons. Additionally, a further 7 days of data were rejected owing to trigger rate instability. This yields a total livetime of 506.7 days, an average fractional livetime of 92%. Figure 7 shows the fractional livetime achieved in blocks of 10 days. With the exception of one period of extended downtime (due to a failure of the power transformer at the site), during no single 10-day block was the detector live less than 75% of the time.

Zoom In Zoom Out Reset image size

The occasional downtime does not heavily bias the exposure. Figure 7 shows the relative exposure of the instrument, measured as the fractional deviation of the number of shower events observed as a function of reconstructed right ascension. The anisotropy in cosmic-ray arrival direction is subdominant at ∼10⁻³ (Abdo et al. 2009), and we can conclude that the exposure is flat to within ±2%. During the 507-day livetime, the Crab is visible above a zenith angle of 20° for ∼1400 hr and above a zenith angle of 45° for ∼3200 hr.

HAWC is operated 90% of the time. The sensitivity depends very weakly on the season or on local weather. Our exposure in R.A. is uniform to the few-percent level, and the background rate for the analysis bins, a measure of the sensitivity, has a day-to-day variability of only 2%, which is much lower than other systematic errors discussed later. For this reason, we treat the sky exposure for this data set as uniform in R.A. and do not attempt to correct for weather or seasonal variability. With these assumptions, the exposure is completely described by a source's declination and the total number of operational days in the data set.

Even with strict photon/hadron discrimination cuts, the data are still dominated by hadronic cosmic-ray events. Fortunately, the directions of hadronic cosmic rays are randomized by magnetic deflection in transit from their sources, and the population of cosmic rays is isotropic to a few parts in 10³ (Abdo et al. 2009). Gamma-ray sources, by comparison, appear as localized "bumps" on this smooth background. In order to identify gamma-ray sources, this background contamination must be estimated.

Figure 8 shows the background computed in the declination of the Crab across the entire sky in the example bins ${ \mathcal B }=3$ and ${ \mathcal B }=8$. The background is estimated using an algorithm developed for analysis of Milagro data (Abdo et al. 2012) called direct integration. The algorithm assumes that, for 2 hr blocks of time, the all-sky rate is independent of the spatial arrival distribution of events in the detector. This assumption is accurate to a few parts in 10³, the level of the cosmic-ray anisotropy. The background at any R.A., decl. point in the sky, during the 2 hr interval, is determined by convolving the all-sky rate with the spatial distribution of events in detector coordinates. The convolution integral output is averaged over 05 patches of sky to account for sparse events in higher-${ \mathcal B }$ bins. Events from known strong sources, the Galactic plane, Mrk 421 and Mrk 501, and the Crab Nebula, are excluded from the background computation.

Zoom In Zoom Out Reset image size

This convolution naturally handles sporadic downtime in the detector, natural rate drifts as the atmosphere changes, and has been shown to accurately determine the background to a few parts in 10³. Figure 8(a) shows the slight underexposure of events from R.A. ≈ 100. The background estimate accurately captures this feature.

Finally, being completely determined by the data themselves, the background estimate is wholly independent of any simulation of the background. Any remaining systematic error is subdominant and is neglected in the eventual systematic error tabulation.

The point-spread function (PSF), ρ(ψ), describes how accurately the directions of gamma-ray events are reconstructed. Here ψ is the space-angle difference between the true photon arrival direction and the reconstructed direction. To a good approximation, the PSF of HAWC is the sum of two 2D Gaussians with different widths

$$\begin{eqnarray}&&\rho (\psi )=\alpha {G}_{1}(\psi )+(1-\alpha ){G}_{2}(\psi ),\end{eqnarray} \tag{ 4 }$$

where G_i is a Gaussian distribution with width σ_i,

$$\begin{eqnarray}&&{G}_{i}(\psi )=\displaystyle \frac{1}{2\pi {\sigma }_{i}^{2}}{e}^{-\displaystyle \frac{{\psi }^{2}}{{2{\sigma }_{i}}^{2}}}\end{eqnarray} \tag{ 5 }$$

which is normalized to unity across the unit sphere.

Figure 9 exhibits the measured angular resolution in HAWC data in two size bins ${ \mathcal B }=3$ and ${ \mathcal B }=8$ obtained by assuming that the Crab Nebula is a point source and that all the angular spread observed in HAWC is due to instrumental precision. The solid-angle density of recorded events dN/dΩ in the vicinity of the Crab is shown as a function of ψ². Bins of ψ² have constant solid angle (in the small-angle approximation), so any remaining cosmic-ray background shows up as a flat component and the gamma rays are evident as a peak near ψ² = 0. The improvement in angular resolution for larger events is clear.

Zoom In Zoom Out Reset image size

Fits to this functional form of Equation (4) can have highly coupled parameters. It is more useful and traditional to quantify the resulting fits with the 68% containment radius, ψ₆₈, the angular radius around the true photon direction in which 68% of events are reconstructed. Figure 10 shows ψ₆₈, for each ${ \mathcal B }$ of the analysis, measured on the Crab and predicted from simulation. At best, events are localized to within 017, the best angular resolution achieved for a wide-field ground array detector directly sampling the particle cascade on the ground.

Zoom In Zoom Out Reset image size

Knowing the angular resolution is critical to subsequent steps of the analysis. Figure 9 indicates that the simulated angular resolution is in good agreement with measurements of the Crab Nebula. This is important because the angular resolution of HAWC for objects at declinations above and below the Crab will differ. While the measured PSF at the position of the Crab cannot be easily extrapolated to other declinations, the simulation can be used to predict the shape of the PSF at any declination. Therefore, the data–simulation agreement shown in Figure 9 is an important verification step. The angular resolution is weakly dependent on the assumed spectral index, and the impact on the analysis is described later in Section 4.3.

The two parameters described in Section 2.6, the compactness ${ \mathcal C }$ and PINCness ${ \mathcal P }$, are used to remove hadrons and keep gamma rays. Events are removed using simple cuts on these variables, and the cuts depend on the size bin, ${ \mathcal B }$, of the event. The cuts are chosen to maximize the statistical significance with which the Crab is detected in the first 337 days of the 507-day data set. Concerns of using the data themselves for optimizing the cuts are minimal with a source as significant as the Crab.

Table 2 shows the cuts chosen for each ${ \mathcal B }$ bin. The rates of events across the entire sky going into the 9 bins, after hadron rejection cuts, vary dramatically, from ∼500 Hz for ${ \mathcal B }=1$ to ∼0.05 Hz for ${ \mathcal B }=9.$ Figure 11 shows the predicted efficiency for gamma rays (from simulation) along with the measured efficiency for hadronic background under these cuts. The efficiency of photons is universally greater than 30% while keeping, at best, only 2 in 10³ hadrons. The efficacy of the cuts is a strong function of the event size, primarily because larger cosmic-ray events produce many more muons than gamma-ray events of a similar size.

Zoom In Zoom Out Reset image size

The limiting rejection at high energies is better than predicted in the sensitivity design study (Abeysekara et al. 2013). The original study was conservative in estimating the rejection power that HAWC would ultimately achieve. With more than a year of data, we now know the hadron rejection of the cuts and can accurately compute the background efficiency.

The HAWC data are fit using the maximum likelihood approach to find the physical flux of photons from the Crab (Wilks 1938; Younk et al. 2016). In this approach, the likelihood of observations is found under two "nested" hypotheses where some number of free parameters are fixed in one model. We form the likelihood of our observations under the null hypothesis, ${{ \mathcal L }}_{\mathrm{null}}$, and an alternative hypothesis, ${{ \mathcal L }}_{\mathrm{alt}}$, with additional free parameters that are not in the null-hypothesis model. If the null model is true, the Test Statistic, $\mathrm{TS}=-2\cdot \mathrm{ln}({L}_{\mathrm{null}}/{L}_{\mathrm{alt}})$, will be distributed as a χ² distribution with a number of degrees of freedom equal to the number of additional free parameters in the alternative model, allowing a quantitative description of how much improvement the additional parameters provide.

The likelihood function is formed over the small (on the scale of the angular resolution) spatial pixels within 2° of the Crab. Each pixel p has an expected number of background events of B_p and, for a specific flux model, an expected number of true photons ${S}_{p}({\boldsymbol{a}})$, where ${\boldsymbol{a}}$ denotes the parameters of our spectral model of the Crab. The predicted photon counts fall off from the source according to the assumed PSF. The likelihood ${ \mathcal L }({\boldsymbol{a}})$ is then the simple Poisson probability of obtaining the measured events in each pixel, M_p, under the assumption of the flux given by ${\boldsymbol{a}}$. The ${ \mathcal B }$ dependence of each term in Equation (6) is suppressed,

$$\begin{eqnarray}&&\mathrm{ln}({ \mathcal L }({\boldsymbol{a}}))=\displaystyle \sum _{{ \mathcal B }=1}^{9}\,\displaystyle \sum _{p=1}^{N}\mathrm{ln}\left(\displaystyle \frac{{({B}_{p}+{S}_{p}({\boldsymbol{a}}))}^{{M}_{p}}{e}^{-{B}_{p}-{S}_{p}({\boldsymbol{a}})}}{{M}_{p}!}\right).\end{eqnarray} \tag{ 6 }$$

Specifically, we fit a differential photon flux ϕ(E) of the log parabola (LP) form:

$$\begin{eqnarray}&&\phi (E)={\phi }_{0}{(E/{E}_{0})}^{-\alpha -\beta \cdot \mathrm{ln}(E/{E}_{0})}.\end{eqnarray} \tag{ 7 }$$

Here ϕ₀ is the flux at E₀, α is the primary spectral index, and β is a second spectral index that governs the changing spectral power change of spectral shape across the energy range of the fit. In this formulation, E₀ is not fitted but is chosen to minimize correlations between the free parameters in the fit. When fit to an LP function, E₀ = 7 TeV produces good results.

We find the parameters for ${\boldsymbol{a}}$ that maximize the likelihood function under signal and background hypotheses and quantify the error region of ${\boldsymbol{a}}$ using Wilks's theorem (Wilks 1938). Figure 12 shows the corresponding regions of α, β, and ϕ₀ for the LP fit that are consistent with HAWC data at 1σ and 2σ. The maximum likelihood occurs at α = 2.63 ± 0.03, β = 0.15 ± 0.03, and ${\mathrm{log}}_{10}({\phi }_{0}\,{\mathrm{cm}}^{2}\,{\rm{s}}\,\mathrm{TeV})=-12.60\pm 0.02$. At this best flux the TS, compared to the background-only hypothesis, is 11,225, a more than 100σ detection.

Zoom In Zoom Out Reset image size

Figure 13 confirms that the fit is a faithful representation of the HAWC data. We show the excess number of events above background in each ${ \mathcal B }$ bin. These excesses come from a fit to the total number of events in the Crab in each ${ \mathcal B }$ bin. The measured number of events in each bin agrees well with the prediction from simulation at the fitted spectrum for the Crab.

Zoom In Zoom Out Reset image size

The TS between an unbroken power-law hypothesis (with β = 0) and the full LP fit is 142, so the spectrum is inconsistent with an unbroken power law at 12σ. The best-fit hypothesis with an assumption that β = 0 is α = 2.57 ± 0.02 and ${\mathrm{log}}_{10}({\phi }_{0}\,{\mathrm{cm}}^{2}\,{\rm{s}}\,\mathrm{TeV})=-12.72\pm 0.01.$

We quantify the energy range of this fit in two ways. First, we take the spectral fit solution and compute the lower 10% quantile of true energy for ${ \mathcal B }=1$ and the upper 90% quantile of true energy for ${ \mathcal B }=9.$ These are 375 GeV and 85 TeV, respectively. This is the energy range over which, under the fitted hypothesis, most of the measured photons are expected to lie.

This energy range is representative of the true energies over which we believe that HAWC can detect the Crab. However, the energy analysis is rudimentary. Though we certainly detect the Crab Nebula for a ${ \mathcal B }$ of 1 through 9, the energies these bins correspond to can be uncertain by nearly an order of magnitude, depending on the spectral shape of the source. Hence, while 375 GeV–85 TeV represents the energy range over which the Crab is detected under the assumed functional form, we cannot definitively say that we detect either 375 GeV or 85 TeV photons with the current analysis.

A more conservative approach can be made focusing on the lowest and highest energies where HAWC data can definitively claim detection of photons. To do this, we separately fit functions of the form

$$\begin{eqnarray*}\phi (x)=\left\{\begin{array}{ll}0 & \mathrm{if}\,E\geqslant {E}_{\mathrm{high}}\\ {\phi }_{0}{E}^{-\alpha }, & \mathrm{otherwise}\end{array}\right.\end{eqnarray*}$$

and

$$\begin{eqnarray*}\phi (x)=\left\{\begin{array}{ll}0 & \mathrm{if}\,E\leqslant {E}_{\mathrm{low}}\\ {\phi }_{0}{E}^{-\alpha }, & \mathrm{otherwise}\end{array}\right.\end{eqnarray*}$$

to find the highest E_low and lowest E_high that are, at 1σ, inconsistent with the HAWC observation. This is a contrived, simplified test with the sole objective of making conservative estimates of the maximum and minimum energy of photons contributing to the Crab result. With this approach, we believe that we have positive detection of photons from the Crab between 1 and 37 TeV. That is, the ${ \mathcal B }=1$ detection means photons of at least 1 TeV, and the ${ \mathcal B }=9$ detection means photons of no lower than 37 TeV. This is not to say that higher- or lower-energy photons cannot be a part of the HAWC observation, but using the event size ${ \mathcal B }$ to measure the energy of photons limits the dynamic range of the observation. Other sources at other declinations may yield different high and low energies.

Table 3 summarizes the major systematic errors contributing to the measurement of the Crab spectrum with HAWC. These systematic errors have been investigated by computing the spectrum from the Crab under varying assumptions to study the stability of the results under perturbation of the assumptions.

For spectral measurements, a systematic error in three quantities is shown: the overall flux, the spectral index measured, and the energy scale. The errors are summed in quadrature to arrive at a total systematic error. In addition to these systematic errors, a systematic error in the absolute pointing of the instrument has been studied.

4.3.1. Charge Resolution and Relative Quantum Efficiency

4.3.2. PMT Absolute Quantum Efficiency

4.3.3. Time Dependence, PMT Layout, and Crab Optimization

4.3.4. Angular Resolution

4.3.5. Late Light Simulation

4.3.6. Absolute Pointing

Figure 14 shows the Crab spectrum measured with HAWC between 1 and 37 TeV compared to the spectrum reported by other experiments. It is consistent with prior measurements within the systematic errors of the HAWC measurement.

Zoom In Zoom Out Reset image size

A number of improvements to the HAWC measurement—most notably the inclusion of a proper energy reconstruction—will reduce the systematic errors and increase the dynamic range of our measurements. The observation of the Crab validates this analysis for subsequent application to the sources across the rest of the sky.

Figure 15 shows the effective area of HAWC, in this analysis, to photons arriving within 13° from overhead. The effective area is defined as the geometrical area over which events are detected, convolved with the efficiency for detecting events, and is determined from simulations. The exact conditions for a photon to be considered detected are complicated for this analysis because, since we are performing a likelihood fit of all pixels in the vicinity of the Crab, photons that are poorly reconstructed play some role in the analysis. In order to have a well-formed effective area, we consider only photons defined within the 68% containment radius, from Table 1. For comparison, Figure 15 includes the progression of cuts, from the effective area without any photon/hadron discrimination without a strong angular accuracy cut to the full analysis cuts. The effective area can exceed the geometrical area of the instrument (about 2 × 10⁴ m²) because events with a core location off the detector occasionally pass the imposed cuts.

Zoom In Zoom Out Reset image size

The effective area computed in this work is smaller than in Abeysekara et al. (2012). For this analysis, developed for steady, multi-TeV sources, we employ tight angular cuts and have a higher energy threshold than used in the initial design study. Furthermore, many small events are excluded by demanding that they are larger than ${ \mathcal B }=1$. Excluding these events limits the effective area below 1 TeV.

Figure 16 shows the computed differential sensitivity of HAWC to sources at the declination of the Crab utilizing the procedure of Abeysekara et al. (2013) with the analysis presented here. A point source of differential photon spectrum ${E}^{-2.63}$ is simulated and fitted using the full likelihood fit. The flux required to be detected at 5σ 50% of the time is shown for each bin, ${ \mathcal B }$. The lines for each ${ \mathcal B }$ are shown with a width corresponding to the width required to contain 68% of the events under the ${E}^{-2.63}$ hypothesis. A correction is made to adjust the ${ \mathcal B }$ separation to a quarter decade in true energy, and the result is fitted. The sensitivity prediction from Abeysekara et al. (2013) suffered from uncertain background at the highest energies. Now, with more than a year of data, we know the background precisely and can set the cuts appropriately. This has resulted in a more accurate (and more sensitive) analysis above 10 TeV. Below about 1 TeV, for a number of reasons, the sensitivity is somewhat worse than predicted in the original study. The background is larger than the original simulation-only prediction. Furthermore, in the current analysis we employ a relatively high cut (defined by ${ \mathcal B }=1$) so that improperly modeled noise can be neglected.

Zoom In Zoom Out Reset image size

The main limitation of this analysis is the reliance on the number of PMTs, used for the definition of ${ \mathcal B }$, to simultaneously constrain the energy of photons, angular resolution, and photon/hadron efficiency. Figure 2 shows that this is a poor energy estimation, with each bin ${ \mathcal B }$ spanning roughly an order of magnitude of energy. More critically, an overhead ∼10 TeV photon can trigger nearly every PMT in HAWC if the core lands near the center of the detector. Consequently, ${ \mathcal B }=9$ is an overflow bin of everything above 10 TeV. These limitations can be removed with an event parameter that accounts for the light level in the event and the specific geometry and inclination angle of events. Approaches like this are under development. The planned deployment of a sparse "outrigger" array should further increase the sensitivity to photons above 10 TeV (Sandoval 2016).

Additionally, the principal systematic error (the modeling of late light) is conservatively estimated here and is being studied using the calibration system. It is likely that the effects of late light will be better modeled in the future.

Finally, the threshold for this analysis is established by including only events where more than 6.7% of the PMTs detect light. The typical number of live, calibrated PMTs is ∼1000, corresponding to a threshold of ∼70 PMTs. Events with 20–30 PMTs could be reconstructed if the noise could be confidently identified. A relatively high event size threshold is used in this analysis to reduce its dependence on the modeling of noise hits. Planned improvements in the modeling should lower the energy threshold of the spectrum analysis in future studies.

The HAWC instrument is performing well, with survey sensitivity exceeding current-generation instruments above 10 TeV, sensitivity that HAWC maintains across much of its field of view. The all-sky survey Abeysekara et al. 2017 conducted by HAWC probes unique flux space and reveals the highest-energy photon sources in the northern sky. Understanding the Crab gives confidence in the survey results.