Measurement of the Crab Nebula Spectrum Past 100 TeV with HAWC

The GP algorithm is based primarily on the charge density at a fixed optimal distance from the shower axis. As discussed in the introduction, this is a robust estimator of the energy reaching the ground.

The radius at which the uncertainty in the shower energy density is minimized (hereafter known as the "optimal radius") must be far from the air shower axis owing to the presence of large shower-to-shower fluctuations that make energy estimation difficult, but it also must be close enough to the shower axis that the measured PMT signal is large enough that threshold effects in the electronics are not a concern.

To determine this optimal radius, the LDF is fit to a modified version of the Nishimura-Kamata-Greisen (NKG) function. The canonical NKG function measures particle density. The signals in water Cherenkov detectors scale with energy deposited in the water. Therefore, a factor of 1/r is introduced to measure energy density rather than particle density. The lateral distribution of the energy is steeper by a factor of 1/r because the highest-energy particles are less likely to be scattered away from the air shower axis (Kamata & Nishimura 1958). Note that this technique is similar to the method used by the Tibet air shower array (Kawata et al. 2017). The main difference is that the Tibet implementation measures particle density, while this method measures energy density.

The LDF fit function, which gives the PMT signal (charge, hereafter called sig_r) as a function of several variables, is

$$\begin{eqnarray}\begin{array}{rcl}{\mathrm{log}}_{10}\left({\mathrm{sig}}_{r}\right) & = & A+s\left[{\mathrm{log}}_{10}\displaystyle \frac{r}{{r}_{m}}+{\mathrm{log}}_{10}\left(1+\displaystyle \frac{r}{{r}_{m}}\right)\right]\\ & & -3{\mathrm{log}}_{10}\displaystyle \frac{r}{{r}_{m}}-4.5{\mathrm{log}}_{10}\left(1+\displaystyle \frac{r}{{r}_{m}}\right).\end{array}\end{eqnarray} \tag{ 1 }$$

Here A is the logarithm of the amplitude of the fit, s is related to the shower age, and r_m is the Moliére radius, which is ∼124 m at the HAWC site. Since this is the NKG/r function instead of the NKG function, the amplitude of the fit absorbs some of the dependence on the shower age, and a conversion is needed to calculate the true shower age. Parameter r is the distance from the PMT to the air shower axis. A and s are the only two free parameters in the fit. See Figure 2 for a depiction of an LDF fit to this NKG/r function. When doing this fit, PMTs with zero signal are ignored.

Zoom In Zoom Out Reset image size

After obtaining the best fit, s is varied by ±10%, and additional fits are performed, leaving the normalization free. This creates a band of fits (see Figure 2). The location where the width of the band is the smallest is the point where the uncertainty in the signal is minimized. This distance is known as the optimal radius. For HAWC, this value is found to be ∼40 m from the shower axis irrespective of zenith angle or primary-particle energy. This mirrors the findings of Newton et al. (2007), who note that the optimal radius is almost solely a function of array geometry.

Equation (1) is evaluated at r = 40 m. This value is known as log₁₀sig₄₀; it is then translated to energy. For a fixed primary energy, the signal measured on the ground varies strongly with zenith angle owing to the differing amount of atmosphere that air showers entering at different angles must travel through. The formula must be parameterized by zenith angle: ${\mathrm{log}}_{10}\hat{E}=f({{sig}}_{40},\theta )$. The exact functional form of the fit is chosen empirically to provide a good match to simulated events:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\hat{E}=m\left(\theta \right){\mathrm{log}}_{10}{{sig}}_{40}+c\left(\theta \right).\end{eqnarray} \tag{ 2 }$$

Here m(θ) is chosen to be a piecewise linear function and c(θ) is chosen to be a piecewise quadratic function.

It is important to note that 40 m is only the mean optimal radius and that for a given event the true optimal radius may be higher or lower. To quantify the effect of using one optimal radius for all events, the procedure above was repeated with the optimal radii set to both 30 and 50 m. No systematic shifts in the assigned energy were observed.

For the performance of the GP on simulation, see Section 3.3.

The NN energy reconstruction algorithm employs an artificial NN to estimate primary energies of photon events based on several quantities that are computed as part of HAWC's event reconstruction. The Toolkit for Multivariate Analysis NN implementation, described in Hoecker et al. (2007), is used.

The NN energy estimator uses a multilayer-perceptron architecture with two hidden layers and a logistic activation function. The first and second hidden layers have 15 and 14 nodes, respectively.

The values of the 479 NN weights are chosen to minimize the error function

$$\begin{eqnarray}&&D\left({\boldsymbol{w}}\right)\equiv \displaystyle \frac{1}{2}\sum _{i=1}^{n}{u}_{i}{\left[{\mathrm{log}}_{10}\hat{E}\left({{\boldsymbol{x}}}_{i};{\boldsymbol{w}}\right)-{\mathrm{log}}_{10}{E}_{i}\right]}^{2}.\end{eqnarray} \tag{ 3 }$$

This is evaluated using Monte Carlo events, where ${\boldsymbol{w}}$ is the vector of NN weights, n is the number of events, u_i is the relative importance of the ith event, ${{\boldsymbol{x}}}_{i}$ is the vector of input variables for the ith event, $\hat{E}$ is the function returning an energy estimate for a given vector of inputs and vector of weights, and E_i is the simulated energy of the ith event. The values of u_i are chosen to resemble an E⁻² power law, which was selected because an NN trained on such a spectrum was found to produce a relatively constant rms error between 1 and 100 TeV, as shown in Figure 7. The minimization of the error function is performed via the Broyden–Fletcher–Goldfarb–Shanno algorithm, described in Hoecker et al. (2007).

For the performance of the NN on simulation, see Section 3.3.

3.2.1. Input Variables

The mixing matrices, which compare the energy estimate to the simulated energy, can be seen in Figure 3. Each plot is normalized as a joint distribution in true and reconstructed energy, so that its 2D integral is 1. This figure assumes an isotropic E⁻² spectrum of gamma-rays; this assures that there are sufficient events at high energy to evaluate the performance. Several data quality cuts have been applied here: only simulated gamma-ray events whose shower core is successfully reconstructed on the HAWC array, that have PMT signal in more than 6.7% of the active detectors in the array (corresponding to ${ \mathcal B }$ bins 1 and above), and that have a zenith angle of <45° are used. Additionally, the events are selected to have a reconstructed zenith angle less than 075 from the true Monte Carlo value. To more accurately show what this figure would look like for data, gamma/hadron separation cuts have been applied to the simulated gamma-ray events.

Zoom In Zoom Out Reset image size

Figure 4 shows an event-by-event comparison of the two estimators. All of the quality cuts described in the preceding paragraph are also used here. The optimal gamma/hadron separation cuts are different for each estimator. Only events passing both sets of gamma/hadron cuts are shown.

Zoom In Zoom Out Reset image size

Figure 5 shows the difference between the two energy estimators as a function of E. A systematic difference can be seen at low energies, with the NN returning, on average, a lower estimate than the GP. At high energies, there is almost no systematic difference.

Zoom In Zoom Out Reset image size

Two quantities are used to evaluate the energy-dependent performance of the estimators. The first is the resolution: the standard deviation of the energy estimate in log-energy space. The second is the bias, defined as the average difference between the reconstructed and true energies in log space:

$$\begin{eqnarray}&&b\equiv \left\langle {\mathrm{log}}_{10}\hat{E}-{\mathrm{log}}_{10}E\right\rangle .\end{eqnarray} \tag{ 4 }$$

The bias and resolution for both estimators can be seen in Figure 6. Both the NN and GP have a large bias below 1 TeV. This is due to the event selection requirement that a minimum of 6.7% of the array be hit, which removes the vast majority of events below this energy. The only events left are from air showers with upward fluctuations in the number of PMTs hit. Due to the substantial bias and poor resolution below 1 TeV, events with reconstructed energies below this threshold are excluded from the spectral fit. HAWC is not as sensitive to GeV energies as it is to TeV energies, so this choice does not affect the fit.

Zoom In Zoom Out Reset image size

Note that both estimators have very good resolution (less than the bin width of log₁₀(E/TeV) = 0.25) and almost no bias in the high-energy regime (between 10 and 316 TeV). The NN has a more favorable bias below ∼32 TeV, while the GP performs better above this energy.

The log rms error is defined as

$$\begin{eqnarray}&&\rho \equiv \sqrt{\left\langle {\left({\mathrm{log}}_{10}\hat{E}-{\mathrm{log}}_{10}E\right)}^{2}\right\rangle }.\end{eqnarray} \tag{ 5 }$$

This is the bias and resolution added in quadrature. Figure 7 shows the log rms error for both energy estimators. The NN performs better using this metric.

Zoom In Zoom Out Reset image size

The data used in this analysis were collected between 2015 June and 2017 December. The total live time is ∼837 days. The detector had >90% uptime during this period. The loss of live time comes from days where the detector was off for maintenance or as a result of operational difficulties. Additionally, a small amount of data were removed owing to large variances in the zenith angle distribution, which is an indication that the detector was unstable during that period. These instabilities make the background estimation method unusable, so such data are removed. This background estimation technique is described in Section 4.3.

The spectral fit is performed using a binned likelihood technique. This forward-folding method accounts for bias and resolution in the energy estimate. We use a 2D binning scheme based on ${ \mathcal B }$ (described in Section 2) and the estimated energy. This 2D binning scheme was chosen instead of binning solely in energy because the gamma/hadron separation parameters, as well as the angular resolution, depend on both the energy and size of the event. We use nine ${ \mathcal B }$ bins each subdivided into 12 energy bins, for a total of 108 bins. The energy bins are quarter-decade bins in ${\mathrm{log}}_{10}(\hat{E})$, beginning at ${\mathrm{log}}_{10}(\hat{E}/\mathrm{TeV})=-0.5$ (0.316 TeV) and ending at ${\mathrm{log}}_{10}(\hat{E})=2.5$ (316 TeV). See Tables 1 and 2 for the bin definitions. For example, bin 9k would be all of the events with >84% of the array hit and energies between 100 TeV and 177 TeV.

In practice, not all 108 bins are used. Some of these bins have little to no probability that events will populate them; for example, there are no low-energy events where the entire array is hit. Additionally, some of these bins contain so few events that they are not modeled well in the Monte Carlo simulation and are also excluded from the fit. The bins used in the fit were chosen a priori by looking at the distribution of estimated energies across each simulated ${ \mathcal B }$ bin and keeping the central 99% of the events. This removes empty bins, as well as the tails of the distribution, where statistics are low and there are more likely to be mismodeled events and a data/Monte Carlo simulation discrepancy. The effect of this exclusion is discussed in the systematic uncertainty section (Section 4.5). In this analysis, 40 bins are used in the GP fit and 37 bins in the NN fit.

Several improvements have been made from Abeysekara et al. (2017a) to strengthen the analysis. A requirement that the shower core be reconstructed on the HAWC array has been added. This improves the angular and energy resolutions, as events with cores on the array are typically better reconstructed.

Gamma/hadron separation has also been improved (through refinements to the simulation that have improved the data/Monte Carlo simulation agreement), although the gamma/hadron separation variables used in this analysis are unchanged from Abeysekara et al. (2017a). Compactness, first described in Abeysekara et al. (2013), is effective at identifying air showers containing muons. Muons, dominantly present in hadronic (background) showers, appear as localized charge depositions far from the shower core. The second parameter is known as PINCness (Abeysekara et al. 2017a) and measures the smoothness of the LDF. Gamma-ray showers have smoother profiles than hadronic air showers.

Although the gamma/hadron separation variables are unchanged, the actual cut values have been reoptimized in each 2D ${ \mathcal B }$/energy bin. This allows for better identification of the highest-energy events as compared to Abeysekara et al. (2017a), where nearly everything above 30 TeV was included in one analysis bin and had the same gamma/hadron cuts. The cut values are determined a priori using simulated Crab signals and background data, with a requirement that each bin has at least 50% gamma-ray efficiency. The efficiency to gamma-rays in a given bin ranges from 50% to nearly 100%. The gamma/hadron cuts are optimized separately for each estimator. Additionally, the data quality cuts described in Section 3.3 have also been applied to the data.

Another improvement is the point-spread function (PSF). As before, this is modeled as a linear combination of two 2D Gaussians, determined from simulated events. Better modeling of this PSF is one of the significant changes from the previous Crab analysis. The radius required to contain 68% of the photons has a strong dependence on both the event size and energy, so the 2D binning scheme used here allows for a more precise determination of the PSF (see Figure 8). For example, all events from ${ \mathcal B }$ bin 1 in Abeysekara et al. (2017a) had a 68% containment of ∼1°. Here these events have a 68% containment between ∼027 and ∼075, depending on the energy of the shower.

Zoom In Zoom Out Reset image size

Lastly, note that the definition of ${ \mathcal B }$ has changed slightly from Abeysekara et al. (2017a): there ${ \mathcal B }$ was defined as the number of PMTs detecting light divided by the total number of PMTs that were operational at the time. Here the numerator is changed to the fraction of PMTs detecting light within 20 ns of the shower front. This change reduces the number of noise hits contributing to the size of the event.

For small showers, hadronic cosmic rays dominate over gamma-rays even after gamma/hadron separation cuts have been applied. An estimate of this cosmic-ray background is performed individually in each analysis bin. For the lower-energy bins, where there are many events, the standard HAWC background estimation technique is applied. This is known as "direct integration." This algorithm was originally developed by the Milagro Collaboration (Atkins et al. 2003) and has become the standard HAWC background estimation algorithm. As described in Abeysekara et al. (2017a), the background estimate from direct integration is smoothed by an additional 05 to compensate for the sparseness of the background.

In the highest-energy bins, the statistics are too low to give a spatially smooth background estimate using the 2 hr chunks of data that are the backbone of direct integration. A different algorithm known as "background randomization," similar to the one in Alexandreas et al. (1991), is used to average over the entire data set and give a spatially smooth background estimate for these low-background bins. For each bin where the all-sky rate is less than 500 events per day, a 2D distribution of the local coordinates (zenith and azimuth) is constructed. A random (zenith, azimuth) pair is drawn from this distribution for each event and used with the time of the event to calculate an R.A. and decl., which is added to the background map. This process is repeated 10,000 times for each event; the background map is then normalized to the number of events in the map. This produces a background estimate much smoother than given by direct integration. Direct integration is still used for higher-statistics bins, as it is less computationally intensive and is needed to correctly incorporate the cosmic-ray anisotropy into the background estimate.

The background estimation technique described above has the potential to be systematically biased if the local coordinate distributions are not stable in time. The zenith and azimuthal angle distributions have been checked and found to have the required stability.

The functional form assumed for the forward-folded fit is a log parabola:

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}={\phi }_{0}{\left(E/{E}_{0}\right)}^{-\alpha -\beta \mathrm{ln}\left(E/{E}_{0}\right)}.\end{eqnarray} \tag{ 6 }$$

Previous measurements indicate that a log parabola is likely to be a good fit to the Crab Nebula spectrum. The pivot energy, E₀, was chosen to be 7 TeV to minimize correlations with the other parameters. The other parameters are free in the fit, which is performed using the HAWC plug-in to the Multi-Mission Maximum Likelihood framework (Vianello et al. 2015; Younk et al. 2015), an analysis pipeline that is capable of handling data from a wide variety of astrophysical detectors. The spectral parameters ϕ₀, α, and β are chosen to maximize the test statistic

$$\begin{eqnarray}&&{TS}\equiv 2\mathrm{ln}\displaystyle \frac{{L}_{{\rm{S}}+{\rm{B}}}\left({\phi }_{0},\alpha ,\beta \right)}{{L}_{{\rm{B}}}},\end{eqnarray} \tag{ 7 }$$

where L_S+B is the likelihood for the signal-plus-background hypothesis and L_B is the likelihood for the background-only hypothesis.

Although the Crab Nebula is slightly extended at TeV energies (Holler et al. 2017), it is modeled as a point source here. HAWC lacks the angular resolution to measure the extent.

The spectra of the Crab Nebula obtained using the two energy estimators can be seen in Figure 9, and the global best-fit parameters over the HAWC energy range can be seen in Table 3. Uncertainties quoted in the table are statistical only. Systematic uncertainties are discussed in Section 4.5. The test statistic is 17,995 for the GP and 19,402 for the NN. The chi square per degree of freedom (χ²/NDF) is approximately 1.7 for both the GP and NN log-parabola fits, dominated by the low-energy (high-statistics) bins. Adding a systematic uncertainty of 1%–2% in quadrature with statistical uncertainties reduces the χ²/NDF to 1. This value was computed using 2 times the optimal top-hat radius in each 2D bin.

Zoom In Zoom Out Reset image size

Alternative spectral models were also considered. For a power law, the test statistic is 17,865 (19,347) for the GP (NN). For a power law with an exponential cutoff, the test statistic is 17,979 (19,395) for the GP (NN). We report the log parabola as the nominal spectrum because it offers the most improvement over a power law for both energy estimation methods over the HAWC energy range.

Flux points are calculated by holding α and β constant from the global fit and fitting the normalization (ϕ₀) individually for each group of ${ \mathcal B }$ bins that contribute to a given reconstructed energy bin, similar to Yuan et al. (2013). While this is not a full unfolding prescription, it allows one to see whether any energy bins are inconsistent with the fitted log-parabola spectrum. Figure 9 shows the Crab Nebula spectrum computed in this manner. The error bars are statistical only. Systematic uncertainties (discussed in Section 4.5) are shown as a band over the global forward-folded fit. Points are shown for each reconstructed energy bin where the statistical significance in the fit is above 2σ. The flux points are plotted at the median simulated energy in each bin, as determined from the Monte Carlo simulation.

The test statistics for each flux point in each of the two spectra (Figure 9) are listed in Table 4. Since the test statistic in the last bin (>177 TeV) is only 0.33 for the GP and 0.14 for the NN, upper limits are set in this bin using a 95% upper confidence interval following Feldman & Cousins (1998).

Table 4.  Test Statistic as a Function of Energy and Flux Points

Bin	$\hat{E}$ Energy Range	GP TS	GP Median	GP Flux	NN TS	NN Median	NN Flux
	(TeV)		Energy (TeV)	(TeV cm−2 s−1)		Energy (TeV)	(TeV cm−2 s−1)
c	1–1.78	3896	0.932	(3.73 ± 0.07) ×  10−11	2734	1.04	(3.63 ± 0.08) ×  10−11
d	1.78–3.16	3754	1.46	(3.11 ± 0.07) ×  10−11	4112	1.83	(2.67 ± 0.05) ×  10−11
e	3.16–5.62	3543	2.68	(2.37 ± 0.06) ×  10−11	4678	3.24	(1.92 ± 0.04) ×  10−11
f	5.62–10.0	3481	5.41	(1.37 ± 0.04) ×  10−11	3683	5.84	(1.24 ± 0.03) ×  10−11
g	10.0–17.8	1864	9.82	(8.26 ± 0.33) ×  10−12	2259	10.66	(8.15 ± 0.31) ×  10−12
h	17.8–31.6	975	18.4	(5.04 ± 0.31) ×  10−12	1237	19.6	(5.23 ± 0.29) ×  10−12
i	31.6–56.2	365	33.9	(2.47 ± 0.27) ×  10−12	572	36.1	(3.26 ± 0.28) ×  10−12
j	56.2–100	107	59.3	(1.26 ± 0.25) ×  10−12	105	66.8	(1.23 ± 0.24) ×  10−12
k	100–177	19.9	102	(6.79 ± 2.70) ×  10−13	28.8	118	(8.37 ± 2.91) ×  10−13
l	177–316	0.33	174	<5.92 × 10−13	0.14	204	<8.14 × 10−13

Note. The test statistic for each energy bin, corresponding to the flux points in Figure 9. The "$\hat{E}$ Energy Range" column gives the range in reconstructed energy for each bin, while the columns labeled "GP Median Energy" and "NN Median Energy" give the median energy from simulation for the GP and NN, respectively, assuming that the fitted log-parabola spectra are the true spectra. Some median energies fall outside the reconstructed energy range because the Crab Nebula spectrum is steep, so that there are more photons with lower energy than higher that are reconstructed at a given $\hat{E}$. The flux gives statistical uncertainties only and is reported at the median energy in each bin. The last bin is a 95% upper limit following Feldman & Cousins (1998).

Download table as:  ASCIITypeset image

The measured excess per transit, along with the expected value from simulation, can be seen in Figure 10. Assuming that the true spectrum is the measured HAWC spectrum, the simulation predicts 57.87 gamma-rays with a reconstructed energy above 1 TeV per day from the Crab Nebula using the GP analysis chain and 48.36 using the NN analysis chain. The values are different because the gamma/hadron cuts were optimized separately for the two techniques and therefore have different efficiencies to gamma-rays. In the data, we observe an excess of 60.85 ± 2.10 gamma-rays with the GP and 47.72 ± 1.28 with the NN, consistent with expectations. All values are computed using a 2° radius centered on the Crab Nebula location.

Zoom In Zoom Out Reset image size

4.4.1. Bin Contamination in Spectral Fits

Bin purity measures the contamination of a reconstructed energy bin by mis-reconstructed events. It is defined here as the fraction of events in a quarter-decade reconstructed energy bin whose simulated energy is also within that bin:

$$\begin{eqnarray}&&p\equiv P\left(E\in B| \hat{E}\in B\right),\end{eqnarray} \tag{ 8 }$$

where B is a quarter-decade energy bin. Both bias and energy resolution can affect the bin purity.

Because astrophysical sources emit following roughly power-law spectra with negative spectral indices, there are many more lower-energy gamma-rays than higher-energy ones. If even a small percentage of these low-energy gamma-rays are mis-reconstructed with a higher energy, the bin purity can be adversely affected. Thus, this parameter is a function of spectral assumption. A softer spectrum will have worse bin purity. Figure 11 shows the bin purity for both estimators. For a power-law spectrum with index between 2 and 3, bin purity is ≳50% above 100 TeV. Bin purity can worsen if the observed gamma-ray spectrum has a cutoff or curvature.

Zoom In Zoom Out Reset image size

Note that bin contamination is not a concern in the likelihood fit described above; since the fit is forward-folded, biases and energy resolution in the energy assignments are taken in account. However, bin purity is a concern if one wants to make a claim about the emission at the highest energies.

4.4.2. Significance of the Highest-energy Detection

Detection of the Crab Nebula above ∼75 TeV would be the highest-energy detection of any astrophysical source to date. Figure 12 gives a significance map of the region for $\hat{E}\gt 56\,\mathrm{TeV}$ for both estimators. Figure 9 provides flux estimates above 100 TeV. While at these energies the cosmic-ray background is significantly suppressed, the possibility of events that might have their energy overestimated and are statistical upward fluctuations from lower-energy bins becomes a concern. We investigate this possibility by fitting the Crab Nebula to the product of a log parabola and a step function, which effectively introduces the null hypothesis of no events above a certain hard-energy cutoff value. All of the parameters of the log parabola are left free.

Zoom In Zoom Out Reset image size

We find that the conventional log-parabola fit is significantly preferred over the log parabola convolved with a hard cutoff at 56 TeV for both estimators (5.12σ for the GP and 6.99σ for the NN, respectively). Moving the hard cutoff to 100 TeV, the conventional log-parabola fit is preferred over the cutoff by 0.2σ for the GP and 2.4σ for the NN. The differences between the two methods can be explained by a combination of differences in the gamma/hadron cuts (which causes differences in gamma-ray efficiency) and statistical fluctuations. The NN has a higher efficiency to gamma-rays above 100 TeV. We interpret this as evidence for emission up to at least 100 TeV from the Crab Nebula. This forward-folding procedure accounts for the energy resolution and bias but ignores systematic uncertainties on the energy scale. This should be taken as a conservative approach to the maximum energy that emission from the Crab Nebula is detected at.

Assuming that the measured energy spectrum of the Crab Nebula extends significantly past 100 TeV, we can use the procedure outlined in the following paragraphs to estimate the highest energy of the photons actually detected by HAWC.

Table 4 gives the median energy from simulation for a source transiting at the Crab decl. and with the best-fit spectra for each energy estimator. This number takes into account events that may have their energies overestimated and are upward fluctuations from a lower-energy bin. For both estimators, the last bin with a significant detection has a median energy above 100 TeV. The median energy is 102 TeV for the GP and 118 TeV for the NN. The somewhat large difference in median energies between the estimators can be explained by differing bin purities that stem from differences in energy resolution (see Figures 6 and 11). This calculation assumes that the true spectrum of the Crab Nebula is the fitted log parabola.

We can expect roughly half of the ∼11 events in the 100–177 TeV bin to be above the median energy. From the binomial distribution, the probability of seeing zero events above the median is simply (0.5)¹¹, or 0.000488. This corresponds to a 3.3σ detection of gamma-rays above the median energy (102 TeV for the GP and 118 TeV for the NN).

If instead 2σ is used as the threshold in the binomial calculation (which is the same threshold chosen for plotting flux points vs. setting an upper limit), the spectrum is detected to 121 TeV for the GP and 137 TeV for the NN. The spectra shown in Figure 9 are plotted up to the 2σ numbers from this binomial calculation.

The main sources of systematic uncertainties with HAWC come from discrepancies between the data and the simulated Monte Carlo events that stem from uncertainties in the modeling of the detector. The systematic uncertainties described in Section 4.3 of HAWC's previous Crab analysis (Abeysekara et al. 2017a) are present here, although improved detector modeling and constraints on the simulation parameters based on low-level data distributions have decreased the size of these uncertainties. Rather than quoting one number for the systematic uncertainty on the flux, all of the uncertainties are treated in an energy-dependent manner for the first time. This is an improvement over Abeysekara et al. (2017a), where the systematic uncertainty was quoted at ±50% across HAWC's entire energy range.

We have looked for correlations between the sources of systematic uncertainty and have not found any. Therefore, the effect of each source of systematic uncertainty can be added in quadrature to the others. The systematic uncertainties on each of the fit parameters in the log-parabola likelihood fit can be seen in Table 5.

The major sources of systematic uncertainty are described below. Figure 13 shows the shift due to systematics in E²dN/dE as a function of energy for each estimator.

Zoom In Zoom Out Reset image size

4.5.1. Angular Resolution Discrepancy

4.5.2. Late Light Simulation

4.5.3. Charge Uncertainty

4.5.4. Absolute PMT Efficiency/Time Dependence

4.5.5. PMT Threshold

4.5.6. Bin Selection

4.5.7. Additional Sources of Systematic Uncertainty

The two energy estimators take very different approaches in deriving an estimate of the gamma-ray energy and give compatible results. There is a small discrepancy above ∼90 TeV, where the GP and NN forward-folded global best fits do not agree within systematic uncertainties. However, the flux points, which show by how much the data agree with the forward-folded fit at a given energy, agree within statistical uncertainties.

The disagreement in the forward-folded fit is likely a combination of two effects. First, potentially mismodeled parameters in the Monte Carlo simulation may affect the two energy estimators differently. Second, due to the gamma/hadron cuts chosen, the two estimators have different efficiencies to gamma-rays above 100 TeV, with the NN's being slightly higher. Given the small number of events above this energy, the inclusion or exclusion of a single event can easily account for the difference in significance. This disagreement does not affect the significance of the observed high-energy emission from the Crab Nebula. Regardless of which analysis technique is used, the Crab Nebula is seen past 100 TeV.

The energy resolution is lognormal with a linear equivalent of 40% (NN) to 55% (GP) at 1 TeV and 23% (NN) to 30% (GP) at 50 TeV. This is a significant improvement over the previously published HAWC analysis (see Figure 2 of Abeysekara et al. 2017a) For comparison, IACTs typically have a resolution of ∼8–15% at 1 TeV and ∼15–35% at 50 TeV (Aleksić et al. 2012; Parsons & Hinton 2014; Park 2015). Starting around 50 TeV, the techniques presented here give comparable energy resolution to what IACTs achieve.

There is good agreement between the spectra presented here and results from other experiments, as can be seen in Figure 9. This is true regardless of which energy estimator is chosen. In particular, improved detector modeling has eliminated the tension at the low-energy end (∼1 TeV) between the original HAWC Crab fit presented in Abeysekara et al. (2017a) and measurements from IACTs.

Compared to the inverse Compton model in Meyer et al. (2010), which has been used as a reference spectrum to compare the energy scale of IACTs, both methods presented here have a 20% higher flux at 7 TeV. When applying a scaling of 0.94 on the energy scale, a deviation less than 10% from the IC model is achieved below 20 and 100 TeV for the GP and NN methods, respectively. The more curved spectrum of the GP method tends more toward the recent publication by MAGIC (Aleksić et al. 2015).

This detection is the highest-energy observation of the Crab Nebula to date. Additionally, the development of two methods to identify gamma-rays above 100 TeV lays the foundation for future high-energy analyses across the entire HAWC field of view. Extending the measured energy range of previously discovered sources up to 100 TeV or higher may allow us to distinguish between leptonic and hadronic gamma-ray emission mechanisms, as they have different signatures. This, in turn, may help us determine whether any Galactic gamma-ray sources are good candidates to be the source of the astrophysical neutrinos discovered by IceCube (Aartsen et al. 2013). Due to gamma-ray attenuation, it is expected that >50 TeV gamma-rays will only arrive at Earth from nearby sources (<100 Mpc), excluding nearly all active galactic nuclei and cosmological sources (Hoffman 2009). Extending the spectra to high energies may also identify PeVatron candidates and give insight into the origins of cosmic rays (Gabici & Aharonian 2007).

Additionally, high-energy observations also naturally lead to studies of Lorentz invariance violation. Particle physics models that add Lorentz-invariance-violating terms to the electromagnetic part of the Standard Model Lagrangian allow photon decay to electron/positron pairs above some energy. Since the decay probability is very nearly 1 for photons propagating across astrophysical distance scales, observations of high-energy photons constrain the energy at which such decay becomes allowed (Martínez-Huerta 2017). The measurements presented in this paper do not by themselves imply a limit on this energy scale; rather, it must be shown that there is a statistically significant excess of events above some reconstructed energy compared to the event rate expected owing to hadronic events and lower-energy photons whose energies are overestimated. Both the uncertainty on the true spectrum of the source and the systematic uncertainties of the HAWC detector must be considered. Such an analysis will be carried out in a future paper.

HAWC recently obtained a boost in high-energy sensitivity with the completion of an upgrade. This sparsely populated "outrigger" array allows for better reconstruction of the largest, most energetic events (Joshi & Jardin-Blicq 2017). Data from the outrigger array are not used here but will be used in future analyses.