PeV $\gamma$-rays from the Crab

The Crab Nebula is a remnant of the explosion of a massive star, formed on Aug 24, 1054 A.D., as recorded in Chinese chronicles as a ‘guest star’  (?). It is the brightest pulsar wind nebula, an extended nonthermal structure powered by the ultrarelativistic electron-positron wind from the central neutron star (the Crab Pulsar). This makes the Crab Nebula one of the brightest $\gamma$-ray sources in the sky, which has been observed for decades at TeV energies (?). At the transition from GeV to TeV energies, the spectral energy distribution (SED; $E^{2}{\rm d}N/{\rm d}E$, where ${\rm d}N/{\rm d}E$ is the differential flux of radiation at the photon energy $E$) achieves a maximum around 100 GeV  (?, ?). Previous observations reported a detection of $\approx 80$ TeV (?), and measured the spectrum up to 300 TeV  (?, ?, ?). The angular size of the $\gamma$-ray source at TeV energies has been reported $\approx 50$ arcsec (?). The Crab Nebula was among the first sources detected at energies well beyond 100 TeV using Large High Altitude Air Shower Observatory (LHAASO)  (?).

LHAASO is a dual-purpose complex of particle detectors designed for the study of cosmic rays (CRs) and $\gamma$-rays in the sub-TeV to 1000 PeV energy range (?, ?). When a very high energy extraterrestrial particle enters Earth’s atmosphere, it initiates a cascade consisting of secondary hadrons, leptons, and photons, known as an air shower. The LHAASO detectors record different components of air showers which are used to reconstruct the type, energy and arrival direction of the primary particles.

LHAASO consists of three arrays  (?). The largest is the square kilometer array (KM2A) composing surface counters and subsurface muon detectors. KM2A is designed to detect cosmic rays (CR) from 10 TeV to 100 PeV and identify $\gamma$-ray photons among the CRs using the muon detector array. The muon content in an air shower event measured by later can be used to effectively reject hadronic showers initiated by CRs. The surface array of KM2A is used to measure the energy and arrival direction of primary particles. Thus, KM2A serves as an Ultra High Energy (UHE, $E_{\gamma}>0.1$ PeV) $\gamma$-ray telescope (?) with energy resolution of $\leq 20\%$ angular resolution and 0.25^∘, respectively. For $\gamma$-rays above 50 TeV, the sensitivity of KM2A reaches the flux level of $10^{-14}\rm erg~{}cm^{-2}s^{-1}$ for a single source like the Crab Nebula per year (?).

The Water Cherenkov Detector Array (WCDA) consists of three ponds with total area of $0.08\ \rm km^{2}$, sensitive to $\gamma$-rays down to 0.1 TeV with an angular resolution $\delta\psi\approx 0.2^{\circ}$ (?). WCDA is designed to perform deep surveys of very high energy $\gamma$-ray sources, including the Crab for both pulsed and unpulsed signals. The sensitivity reaches the flux level of $10^{-12}\rm erg~{}cm^{-2}s^{-1}$ for a single source like the Crab Nebula per year at energies above 2 TeV (?) . At energies above 0.1 PeV, WCDA also serves as a large muon detector to enhance the capability of separation between electromagnetic and hadronic showers detected by KM2A.

KM2A and WCDA are complemented by the Wide-Field-of-view Cherenkov Telescope Array (WFCTA) (?) consisting of 18 telescopes designed for detection of the Cherenkov radiation emitted by air showers induced by CRs with energy ranging from 0.1 to 1000 PeV. Cherenkov light in showers initiated by $\gamma$-rays at energies above 0.1 PeV is recorded by the currently operating fourteen WFCTA telescopes which cover a $32^{\circ}\times 112^{\circ}$ patch of the northern sky through which the Crab passed every day.

On 2020 January 11 at 17:59:18 Coordinated Universal Time (UTC), a giant air shower was recorded by all three LHAASO detectors. The shower arrived from the direction of Crab and landed in the western part of KM2A (see Fig.1). WCDA was triggered, despite being 150 m away from the shower core. The event occurred after local midnight when eight WFCTA telescopes were operational. The Cherenkov radiation of the air shower appeared in the Field of View (FoV) of Telescope No.10 and triggered it (also see Fig.1).

We identify the event as a $\gamma$-ray induced shower based on 4996 particles (electrons, photons, muons, hadrons) recorded by 395 surface detectors and 15 muons recorded by 11 under-surface detectors of KM2A. The chance probability of this event to be misidentified is estimated as 0.1%. Two independent estimates of the shower energy were derived from the KM2A and WFCTA data are $0.88\pm 0.11$ PeV and $0.92^{+0.28}_{-0.20}$ PeV, respectively. This value has been reported as the maximum energy of $\gamma$-rays detected from the Crab direction by LHAASO elsewhere (?).

Approximately one year later, 2021 January 4 at 16:45:06, another shower at even higher energy, i.e. 1.12$\pm 0.09$ PeV, was registered by KM2A at zenith angle $12.9^{\circ}$ which is more vertically arrived and better measured than the previous one by KM2A. Unfortunately, the primary photon arrived one hour before the Crab entered the FoV of the WFCTA telescopes. While the number of detected secondary particles (5094) exceeded those in the previous event, the number of muons (14) was fewer, we also identified this event as a $\gamma$-ray induced shower with a 0.03% chance probability of mis-identification.

KM2A operated for 314 days at half of its design capacity and a further 87 days at three quarters of its capacity. In that time, a total of 89 UHE $\gamma$-rays with energy exceeding 0.1 PeV were detected from the Crab. Fig.2 shows the integrated $\gamma$-ray detection rate, re-normalized for the nominal $1\ \rm km^{2}$ array and for one hour exposure time, assuming the Crab is within the FoV of KM2A (approximately 7.4 hours per day with zenith angle less than 50^∘). Above 0.1 PeV, we find about 0.05 events per hour, equivalent to 135 events per year. Fig.2 also shows the detection rate of the CR induced showers within a cone of $1^{\circ}$ centred on the Crab, both before and after applying the CR shower rejection based on the muon content in the showers, namely the ‘muon cut’ filter requiring that the number of muons detected by KM2A in the shower must be less than 1/230 of the number of particles detected by the KM2A surface counters. This cuts the cosmic ray background by the factors of 1,000 and 500,000 at 50 TeV and 1 PeV, respectively. At energies above 0.1 PeV, the detection rate of $\gamma$-rays from the Crab exceeds the CR induced background by an order of magnitude. Because the Point Spread Function (PSF) of KM2A is $\delta\psi\approx 0.25^{\circ}$, the CR background could be lower by an additional factor of $(1^{\circ}/\delta\psi)^{2}\approx 16$.

The $\gamma$-ray fluxes in the energy range from 0.5 TeV to 13  TeV are measured using the firstly built pond of WCDA. From 2019 September to 2020 October, the total exposure was 343.5 transits of the Crab. The measurement using KM2A in the observation period reported above covers the higher energy range from 10 TeV to 1.6 PeV. The SED of the Crab are shown in Fig.3. The two independent measurements of $\gamma$-ray flux using two detector components of LHAASO fit a simple SED functional form ${\rm d}N/{\rm d}E=(8.2\pm 0.2)\times 10^{-14}(E/10\ \rm TeV)^{-\Gamma}$ cm^-2s^-1TeV^-1, where $N$ is the number of $\gamma$-rays, $E$ is the $\gamma$-ray energy and $\Gamma$ is the energy dependent spectral index. The two measurements are smoothly connected in the small overlapping region around 12.5 TeV. In this energy bin, the discrepancy between flux measured by KM2A and WCDA is 1.3$\sigma$. The overall $\chi^{2}/dof$ is 9.3/14, where $dof$ refers to degrees of freedom. No systematic deviation between the two segments of the SED is observed. Function $\Gamma=(2.90\pm 0.01)+(0.19\pm 0.02)\log_{10}(E/10\ \rm TeV)$ implies a gradual steepening of the spectrum characterised by the local index $\Gamma$, from $\approx 2.5$ at 1 TeV to 3.7 at 1 PeV.

Fig.3 also shows previous measurements using atmospheric Cherenkov telescopes  (?, ?, ?, ?) and the air shower arrays  (?, ?, ?). These are consistent with the WCDA and KM2A data from sub-TeV to multi-TeV energies, including the flux at the highest energy $E_{\gamma}\approx 300$ TeV previously reported (?).

Photons from the Crab Nebula have been detected over 22 decades of energy, from MHz radio to UHE $\gamma$-rays, and consists of several pulsed and unpulsed components. $\gamma$-rays can be produced in three physically distinct sites - in the pulsar’s magnetosphere, ultrarelativistic electron-positron (hereafter, simply electron) wind and nebula. UHE $\gamma$-rays are absorbed in the strong magnetic field of the pulsar, thus a pulsed component has been expected only from the magnetosphere in the form of MeV-to-GeV $\gamma$-rays. However, the surprise detection of pulsed TeV $\gamma$-rays from the Crab  (?) initiated a new concept allowing the location of the pulsed $\gamma$-ray production in the wind  (?, ?, ?). Although the extension of this component to the UHE domain seems rather problematic, we plan to search for a pulsed component in the LHAASO data after adequate photon statistics is accumulated. Below we assume that the entire flux detected by LHAASO consists of an unpulsed component and is produced in the nebula.

The broad-band nonthermal emission of the Crab Nebula is dominated by two mechanisms - synchrotron radiation and inverse Compton (IC) scattering of relativistic electrons interacting with the ambient magnetic and radiation fields  (?, ?). In the standard paradigm, the acceleration of electrons is initiated by the termination of the wind by a standing reverse shock at a distance $R\approx 0.1$ parsec (pc) from the pulsar  (?, ?). Although the details of the acceleration remain unknown, the detection of PeV photons allows us to estimate the accelerator size $l$, the magnetic field strength $B$, and the minimum acceleration rate $\eta$.

In the Crab Nebula, several radiation fields supply target photons for the IC scattering of electrons. However, at energies above 100 TeV, the 2.7 K Cosmic Microwave Background radiation (CMBR) dominates the $\gamma$-ray production (?, ?). Because the CMBR is well quantified, the $\gamma$-ray data provides direct information about the parent electrons. $\gamma$-ray production above 100 TeV proceeds in the Klein-Nishina regime, where the energies of the upscattered photon $E_{\gamma}$ and the parent electron $E_{\rm e}$ are linked through the simple relation $E_{\gamma}=0.37(E_{\rm e}/1\,\rm PeV)^{1.3}\,PeV$, which over two energy decades, from 30 TeV to 3 PeV, provides accuracy better than 10% (see Fig.S7), or equivalently

Thus, for the 1.1 PeV photon, the energy of the parent electron is 2.3 PeV. Correspondingly, the mean energies of the synchrotron ($\varepsilon_{\rm syn}$) and IC ($E_{\gamma}$) photons produced by the same electron in the ambient magnetic field $B$, are related by

Simultaneous modelling  (?) of the synchrotron and IC components constrains the magnetic field strength within a narrow interval, $B\simeq 112^{+15}_{-13}\mu$G (see Fig.4 and discussion below). The upper is set by the requirement for the synchrotron radiation of the same electrons responsible for the production of 1 PeV photons to not overshoot the measured MeV flux  (?). The lower limit is set by requiring the electron energy $E_{\rm e}$ and the accelerator’s linear size, $l$, to meet the condition that the electron gyroradius $R_{\rm g}=E_{\rm e}/eB$ ($e$ is the charge of electron) cannot exceed $l$. Using Eq.(1), we find

Magnetohydrodynamic (MHD) models of the Crab Nebula (?, ?) postulate that electrons are accelerated at the termination of a electron wind, then advected into the nebula through the MHD outflow. X-ray imaging of the inner parts of the nebula  (?) provides support for the location of the acceleration site being close to the termination shock, at the distance $R\approx 0.1-0.14$ pc from the pulsar  (?). The linear size of the accelerator should exceed the gyroradius, $l\geq R$; this imposes a lower limit on the magnetic field from Eq.(3), $B\geq 20\,\mu$G. On the other hand, the standard one-zone model, which assumes that both the synchrotron and IC components of radiation are produced in the same region by the same electron population, gives $B\simeq 112\,\mu$G (see Fig.4). Then, from Eq.(3) we find that $l\geq 0.02$ pc. These constraints are inconsistent with estimates of the characteristic size and magnetic field in the region(s) where the flares of the MeV/GeV $\gamma$-ray emission (“Crab flares”)  (?) originate. The variation of $\gamma$-ray fluxes on timescales of days are interpreted by fast acceleration and synchrotron cooling of PeV electrons in compact ($R\leq 0.01$ pc) highly magnetized ($B\geq 1$ mG) regions (see  (?) for a review). In the presence of such a large magnetic field, the IC $\gamma$-ray component is suppressed; however this does not exclude an indirect link between the PeV electrons responsible for the UHE $\gamma$-ray emission and the synchrotron MeV/GeV flares.

The detection of $\sim 1$ PeV photons implies an acceleration rate which overcomes the synchrotron losses of the parent electrons up to PeV energies. The acceleration rate of electrons is $\dot{E}=e\mathcal{E}c=\eta eBc$, where $\eta$ is the ratio of the projection of the electric field $\mathcal{E}$, averaged over the particle trajectory, to the magnetic field, $\eta=\mathcal{E}/B$. This parameter characterizing the acceleration efficiency is always smaller than 1; in objects where $\eta\to 1$, the accelerator proceeds at the maximum rate allowed by the classical electrodynamics and ideal MHD  (?). The maximum energy of electrons then is determined by the balance between the acceleration and the energy loss rates: $E_{\rm e,max}\approx 5.8\eta^{1/2}(B/100\,\mu\rm G)^{-1/2}\ \rm PeV$. Using the relation between $E_{\gamma}$ and $E_{\rm e}$ given by Eq.(1), we find

Thus, for the detected $E_{\gamma}=1.1\ \rm PeV$ and the magnetic field $B\simeq 112\,\mu$G derived from the one-zone model (Fig.4), we find that the acceleration proceeds at the rate $\eta\approx 0.16$. For comparison, at the diffusive shock acceleration in young supernova remnants (?), $\eta$ is smaller by at least 3 orders of magnitude.

For the distance to the Crab, $d=2$ kpc  (?), the luminosity in PeV $\gamma$-rays is estimated $L_{\rm\gamma,PeV}=4\pi d^{2}F_{\gamma}\approx 5\times 10^{31}\ \rm erg\ s^{-1}$, where $F_{\gamma}\approx 10^{-13}\ \rm erg\ s^{-1}$ is the energy flux of 0.5 to 1.1 PeV $\gamma$-rays (Fig.4). In the Klein-Nishina limit, the IC cooling time of electrons in 2.7 K CMBR is $t_{\rm IC}\simeq 5\times 10^{11}(E_{\rm e}/1\ \rm PeV)^{0.7}\,$s  (?) or, for the given $E_{\gamma}$, using Eq. 1 we have $t_{\rm IC}\simeq 8\times 10^{11}(E_{\gamma}/1\,{\rm PeV})^{0.54}\,$s. Thus, the total energy held by the $\geq 1$ PeV electrons responsible for production of $\geq 0.5$ PeV $\gamma$-rays is estimated as $W_{\rm e,PeV}=L_{\gamma,\rm PeV}t_{\rm IC}\approx 4\times 10^{43}\ {\rm erg}$. Because the overall cooling of PeV electrons is dominated by synchrotron losses ($t_{\rm syn}\ll t_{\rm IC}$), the injection rate of PeV electrons $\dot{W}_{\rm e,PeV}=W_{\rm e,PeV}/t_{\rm syn}\approx 2\times 10^{36}(B/100\mu\rm G)^{2}\ \rm erg\ s^{-1}$. Thus, within the framework of the standard one-zone model with $B=112\mu$G, the acceleration power of PeV electrons constitutes $\approx 0.5\,\%$ fraction of the pulsar’s spin-down luminosity, $L_{\rm SD}\approx 5\times 10^{38}\ \rm erg\ s^{-1}$ (?).

We consider whether the detection of PeV photons agrees with predictions for the standard MHD paradigm of the Crab Nebula which assumes that nonthermal emission from X-rays to multi-TeV $\gamma$-rays  (?, ?, ?, ?, ?) is produced by electrons accelerated at the termination of the pulsar wind. We modelled the Crab’s multi-wavelength radiation within the idealized Synchrotron-IC one-zone model, assuming a homogeneous spatial distribution of the magnetic field and electrons (Fig. 4). For $E_{\gamma}\geq 100$ TeV $\gamma$-rays the dominant target for IC scattering is the 2.7 K CMBR, with properties that are known more precisely than the targets at lower energies. For a steady-state electron energy distribution, above 1 TeV, we assumed a power-law function terminated by a super-exponential cutoff at the high-energy end. Fig. 4 demonstrates that using three free parameters - the power-law slope $\alpha=3.42$, cutoff energy $E_{0}=2.15\,$PeV and magnetic field $B=112\,\mu$G, we reproduce the SED fit of reasonable accuracy from the X-rays to multi-MeV $\gamma-$rays with synchrotron radiation, and the TeV to PeV $\gamma$-rays with IC radiation. Below 1 TeV, the electron spectrum must undergo a break to avoid a conflict with the synchrotron radiation at optical to radio frequencies, and to provide a smooth transition of the IC radiation from TeV to GeV energies  (?).

The one-zone model in Fig.4 is tightly constrained; the 3-sigma variation of the index $\alpha$ (3.37 to 3.47) is determined by the uncertainty of the X-ray data. Because $\geq 10$ TeV electrons are cooled on short timescales, the index of the initial (acceleration) spectrum must be close to 2.4. The magnetic field is determined by the flux ratio of the synchrotron and IC components, with 3 sigma range from 100 $\mu$G to 130 $\mu$G. For the given $\alpha$ and $B$, the cutoff energy $E_{0}$ is determined by the synchrotron and IC spectra above 10 MeV and 100 TeV, respectively. The derived ranges of magnetic field and the power-law index of electrons are consistent with the previous studies based on the Synchrotron-IC one-zone model  (?, ?) as well as the MHD flow  (?, ?, ?) models in which the magnetic field’s radial distribution is determined by the wind-magnetization parameter $\sigma_{B}$. The magnetic field $B\simeq 112$ G derived for the production region of multi-TeV to PeV $\gamma$-rays, is a factor of 2-3 smaller than the average nebular magnetic field, consistent with the MHD flow model (?). The latter predicts reduced magnetic field at the termination shock for a broad range of the magnetization parameter $\sigma_{B}$ between 0.001 and 0.01  (?).

Within the one-zone model, the IC $\gamma$-ray spectrum is precisely predicted. While the KM2A spectral points from 10 TeV to 1 PeV do agree, within the statistical uncertainties, with the one-zone fit, there are possible deviations from its predictions. Between 60 TeV and 500 TeV, two differ with a significance of $4\sigma$ indicating a steeper spectrum than the one-zone model predictions. The possible excess around one PeV indicates an opposite tendency - a hardening of the spectrum. A hardening of the electron spectrum is difficult to accommodate with plausible assumptions. The problem of suppression of the one-zone spectrum at 1 PeV can be circumvented by introducing a second population of PeV electrons. This could also explain the inconsistency of the synchrotron part of the SED with the one-zone model (Fig.S4) by decoupling the highest energy Synchrotron and IC components, assuming that the MeV synchrotron radiation is predominantly produced in compact highly magnetized regions, while the PeV IC photons originate from regions with $B\leq 100\mu\rm G$. A second electron component could extend the SED to a few PeV but not much further. From Eq.(4) follows that, even for $\eta=1$ and minimum allowed magnetic field, $B\geq 20\mu$G, the the maximum energy of photons cannot exceed 4 PeV. Any detection of $\gamma$-rays well beyond 1 PeV would require non-leptonic origin of the extra component of radiation, i.e. the presence of multi-PeV protons and nuclei in the nebula.

Because of the limited energy budget available for acceleration of protons, hadronic interactions cannot be responsible for the overall broad-band $\gamma$-ray luminosity. Yet, the contribution of hadronic interactions at PeV energies could be noticeable. The $\gamma$-ray production efficiency of hadronic interactions, i.e. the fraction of the proton kinetic energy converted to $\gamma$-rays, is determined by the ratio $\kappa=t_{\rm esc}/t_{\rm pp}$, where $t_{\rm esc}$ is the confinement time of protons inside the nebula, and $t_{\rm pp}\approx 1.5\times 10^{8}n^{-1}\ \rm yr$ is the cooling time of protons through the production and decay of the secondary $\pi^{0}$-mesons. For the average density of the nebular gas $n\approx 10\,\rm cm^{-3}$, the radiation efficiency is low; even for the most effective confinement of 10 PeV protons, $t_{\rm esc}\leq 250\,$yr  (?), thus $\kappa\leq 5\times 10^{-5}$. To explain the PeV $\gamma$-ray luminosity, $L_{\gamma}\approx 5\times 10^{31}\ \rm erg\ s^{-1}$, the acceleration power of $\sim 10$ PeV parent protons would need to be $\dot{W}_{\rm p}=W_{\rm p}/t_{\rm esc}=\kappa^{-1}L_{\gamma}\approx 10^{36}\ \rm erg\ s^{-1}$ or, in the case of a broad $E^{-2}$ type proton spectrum, an order of magnitude larger. These estimates are supported by numerical calculations presented in Fig.S5. They avoid exceeding the theoretical constraints on the proton fraction  (?) only if there is effective proton confinement in the nebula.

LHAASO  (?) is a complex of extensive air shower (EAS) detectors installed on Mt. Haizi (29^∘21’27.6” N, 100^∘08’19.6” E) at 4410 m above sea level, in Sichuan province, China.

In the local coordinate system, the 0.88 PeV event arrived at about 2 am on January 12th, 2020 (Beijing time), with a zenith angle of $(33.9\pm 0.2)^{\circ}$ according to the shower front reconstructed by the 395 registered scintillator counters, which are synchronized within 0.2 ns. Fig.S1 show the significance map around Crab Nebula region. A clear signal at 6.6 $\sigma$ is observed at energies above 400 TeV. The 0.88 PeV event is found just $0.21^{\circ}$ away from the celestial coordinates of the Crab Nebula, as illustrated in Fig.S1, consistent with the uncertainty in the reconstructed direction.

The event is only likely to be associated with the Crab Nebula if it can be identified as a $\gamma$-ray event. We used the muon detector (MD) array as a veto  (?). The 15 muons ($N_{\mu}$) was compared to the total number of particles measured by all counters ($N_{e}$), 4996, which indicates a very small muon content in this event. The probability of being a hadron induced shower, with such a small muon content, was evaluated by counting how many events have the ratio $N_{\mu}/N_{e}$ less than the ratio of 15/4996. It is found that only a few showers out of 12.3 million cosmic ray samples recorded by KM2A can have the reconstructed energy E${}_{rec}\geq$0.88 PeV. We estimate the chance probability of any single cosmic ray event to be recognized as a $\gamma$-ray event as $6.7\times 10^{-6}$, as illustrated in Fig.S2. To avoid the pollution of diffuse $\gamma$-rays produced within the Milky Way, only events coming from regions 12^∘ away from the Galactic plane are considered in the estimation. Taking into account that the Crab transits have specific occupancy of different zenith angles, this probability has been tested using events in 3 zenith angle intervals, i.e. 0^∘-20^∘, 20^∘-40^∘ and 20^∘-50^∘ and a negligible difference was found. In the direction of the Crab Nebula, KM2A recorded a total of 179 events with the angular distance from Crab Nebula less than 0.3^∘, which is slightly larger than the radius of the PSF of KM2A above 400 TeV. Using the probability given above, the number of expected background CR events is 0.0012. In other words, the chance probability of this event was not a photon is 0.1%.

The 395 scintillator counters measured the lateral distribution of the 4996 secondary particles in the shower. The core of this shower was reconstructed by fitting the distribution with the modified Nishimura-Kamata-Greisen (NKG) functional form, $N\Gamma(4.5-s)/\Gamma(s-0.5)/\Gamma(5-2s)(r/R_{M})^{s-2.5}(1+r/R_{M})^{s-4.5}$, where $N$ is the total number of secondary particles in the shower, $s$ is the shower age parameter, and $R_{M}$ is the Moliere unit of 136 m (?) as a measure of the width of the distribution at the altitude of LHAASO, and $r$ is the perpendicular distance from the detector to the shower axis. The uncertainty of the core location is about 3 m. Here, the particles density at 50 m from the core determined in the fitting has been used as the energy estimator and has resulted in the shower energy $E_{rec}=0.88\pm 0.11$ PeV by assuming the primary particle a $\gamma$-ray photon. As a pure electromagnetic cascade shower, the energy scale and the resolution established using existing simulation tools for air showers and detector responses. The resolution defined as the $(E_{rec}-E_{true})/E_{true}$, where $E_{true}$ is the input shower energy in the simulation, is found to have a Gaussian form with the systematical shift less than 1% and a resolution of 14% for showers above 100 TeV and arriving with angles less than 20^∘ from the zenith (?). The equivalent values for showers above 400 TeV and from a zenith angle less than $50^{\circ}$ are the systematic shift $\sim 1\%$ and the resolution of $\approx 18\%$, see Fig.S3.

To verify the estimation of the energy and its uncertainty of the 0.88 PeV photon, 10,000 events were simulated by using the geometric parameters of the measured event over a wide energy range. The distribution of input energies of events that have the total number of particles recorded by the registered scintillator counters and number of the counters are within $\pm$10% around the detected values, is used to estimate the photon energy and its error as ($0.866^{+0.103}_{-0.097}$) PeV, in a good agreement with the reconstructed energy.

One of the WFCTA telescopes, Telescope No.10, pointing at a zenith angle of 30^∘ and nearly toward west, saw the shower image in the lower 1/4 of its FoV starting from the zenith angle of 34^∘. Since the shower is quite far from the telescope, the image was not very much extended, only 11 pixels had signals above the threshold. However, the image is bright in terms of total number of photoelectrons, which was 8636. 30,000 $\gamma$-ray induced showers with the same geometric parameters were simulated in the energy range from 0.4 to 4 PeV. The generated showers reproduced the footprint on KM2A, WCDA and the image in the telescope. The most probable energy among the simulated showers with $8636\pm 100$ photoelectrons in their images is found to be $0.92^{+0.28}_{-0.20}$ PeV.

This same analysis was applied on all 5 events above 0.4 PeV. The results are summarized in Table S1.

The SED of Crab Nebula below 20 TeV was measured using data collected by WCDA-1 from September 2019 to October 2020. The total exposure of the Crab Nebula in this period of time was 343.5 transits of the Crab. More details about the performances of WCDA have been reported elsewhere (?). The independent data set collected in the operation of half KM2A for 314 days and three-quarter KM2A for 87 days from December 2019 to February 2021 has been used to measure the SED of the Crab Nebula above 10 TeV. More details about the measurement and detector performance have been reported elsewhere (?) for the half KM2A. Selected $>0.1\,$PeV $\gamma$-ray events according to the ratio $N_{\mu}/{N_{e}}<1/230$ were selected in a cone of 0.4^∘ centered at the Crab Nebula direction. The criteria ensure that 90% of signal is contained, giving the angular resolution better than 0.25^∘ of KM2A and the extension less than 0.07^∘ for the Crab Nebula. The background cosmic rays are suppressed to be in the ‘signal-dominated’ regime. In Table S2, the details about the number of on source photons and off source background events in the energy bins are listed.

The energy resolution of KM2A is found to be about 18% for all events with the zenith angle up to 50^∘ as shown in Figure S3 for photon showers above 0.4 PeV. A similar resolution function for lower energies has been reported  (?). Only a small non-Gaussian tails less than 1% is found in the range of 0.5 PeV and above. The systematic shift is less that 1.5%. Breaking into zenith angle ranges, the resolution is found better than 15% for events with the zenith angle less than 40^∘ and 30% for events in the range of $40^{\circ}<\theta<50^{\circ}$. Given the resolution, the energy is divided into 5 bins per decade with a bin width of $\Delta log_{10}E=0.2$ to ensure the bin-to-bin migration mainly happens in the adjacent bins according to the Gaussian-shaped resolution function. In the operation of KM2A, the numbers of the events in bins above 0.1 PeV are listed in Table S2. The number of background events N_b in each bin is estimated by using the direct-integration method. The flux is measured for each bin as a function of $E$, as presented in Fig.3. Many more technical details of the SED measurements can be found elsewhere (?).

Above 400 TeV, the Crab Nebula has been detected with the significance of 6.6 $\sigma$. The sky map around the Crab is shown in Fig. S1.

The standard one-zone scenario We mainly focus on the high-energy part of the electron spectrum which account for keV – MeV and TeV – PeV data, where a power-law type electron spectrum with a high-energy cutoff can match the SED. To avoid violating the radio – optical and GeV data, which should arise from low-energy electrons from a much larger volume of the nebula, a low-energy spectral cutoff or break is needed. Such a low-energy cutoff or break could naturally correspond to the bulk Lorentz factor of the cold pulsar wind just before the termination shock, which for the Crab Nebula is expected around $10^{6}$ (?). If we choose the break energy and the low-energy spectral index, we could also fit the low-energy emission phenomenologically. We here assume a broken-power-law function with a super-exponential high-energy cutoff for the electron spectrum, i.e,

where $\alpha$ is the spectral index after the spectral break, and $\Delta\alpha$ is the spectral difference before and after the break, $E_{b}$ is the break energy, $E_{0}$ is the cutoff energy. The parameter $\alpha$ should match simultaneously the X-ray and multi-TeV $\gamma$-ray spectra; it should be confined within a narrow range, close to $\alpha=3.4$, The parameters $\Delta\alpha$ and $E_{b}$ are important for fitting the low energy IC $\gamma$-ray and optical to radio synchrotron data, but are not directly linked to UHE energies. For any reasonable parameters set, $\Delta\alpha\leq 2$ and $E_{b}\leq 1$ TeV. $N_{0}$ is the normalization factor which can be determined once the total energy in the electron spectrum $W_{e}$ is given, via $W_{e}=\int E_{e}\frac{dN}{dE_{e}}dE_{e}$.

In the magnetic and radiation fields, electrons are radiatively cooled through the synchrotron radiation and the inverse Compton (IC) scattering. For the synchrotron radiation, we consider a turbulent magnetic field with Gaussian distribution of an average strength $B$  (?). For the IC radiation, we employ the 2.7 K CMBR, the interstellar radiation field  (?), the far-infrared radiation excess from the nebula with the temperature of 70 K and energy density of $0.5\rm\,eV~{}cm^{-3}$. The synchrotron radiation of the same population of electrons is also considered in order to calculate the synchrotron-self Compton process. The synchrotron radiation is assumed homogeneously generated within the nebula of a radius of 1.8 pc, and the volume-averaged photon density is enhanced by a factor of 2.24 with respect to the boundary of the nebula  (?). We then explore the parameter space evaluating the goodness of fit by the Pearson’s chi-squared test with the data above 1 keV. The best-fitting result (shown in Fig. 3) gives $\chi^{2}/dof=158.4/99$, among which the 11 data points of KM2A contributes a value of 36.9 corresponding to a $4\sigma$ deviation between the model prediction and the data.

The comparison between the broadband SED and the best-fitting model is shown in Fig. S6, where $\Delta\alpha=1.76$ and $E_{0}=0.76\,$TeV is employed. In the same figure, we show the contributions of electrons from different energy intervals to the SED at different energies. For the magnetic field $B=112$ G, the relation between the characteristic energies of the synchrotron and IC photons from the same electron, is shown in Fig. S7. For the IC scattering on an isotropic target photon field of blackbody/greybody radiation with temperature $T$, the relation between the electron energy and the mean energy of up-scattered photons $\bar{E}_{\gamma}$ is described by  (?):

where $b=E_{e}kT$, $k$ is the Boltzmann constant. All energies in this formula are expressed in unit of $m_{e}c^{2}$.

For $2.7$ K CMBR, a 2.3 PeV electron is required to produce 1.1 PeV photon as shown in Fig. S8. The characteristic energy of the synchrotron radiation is given by $\varepsilon_{\rm syn}=0.29\cdot 3E_{e}^{2}eB/4\pi m_{e}^{3}c^{5}$.

The synchrotron cooling limits the characteristics distance at which the highest energy electrons release all their energy: $l_{\rm rad}\leq 6\times 10^{17}(E_{e}/1\,{\rm PeV})^{-1}(B/100\mu{\rm G})^{-2}\,$cm. The multiwavelength modelling of the synchrotron and IC components of radiation does not allow a deviation of the average magnetic field from $100\,\mu$G implying that both the acceleration and radiation of PeV electrons take place in a compact region(s), attached, presumably, to the wind termination site. The short lifetime of electrons and their confinement in compact regions allow flux variability in both channels - IC $\gamma$-rays at PeV and synchrotron radiation at MeV energies on timescales of months. The synchrotron radiation could experience variability on shorter timescales caused, e.g. due to the synchrotron burning of electrons in locations with suddenly enhanced, up to 1 mG magnetic field, initiating the so-called Crab flares. Therefore, the MeV emission of the Crab Nebula measured previously  (?) may originate from a region which is different from the production region of PeV photons. It will relax the procedure of the spectral fitting given additional model parameters, and the result is shown in Fig. S4.

Another possible acceleration site is the wind termination, through the relativistic shock or shock-driven magnetic reconnection  (?, ?). The formed spectrum is expected to be power law with a slope between 2 and 3. Alternatively, hadronic emission could be produced by protons from the cold, ultrarelativistic pulsar wind with Lorentz factor $\Gamma_{\rm w}$  (?). In this scenario, the protons are monoenergetic, $E_{p}=m_{\rm p}c^{2}\Gamma_{\rm w}$. To produce 1 PeV $\gamma$-ray photon, the bulk Lorentz factor of the pulsar wind should exceed $10^{7}$. The latter is significantly larger than the “standard”, for the Crab pulsar, value of $10^{6}$ (?), although cannot be excluded for the small $e^{\pm}$ electron-positron multiplicity $k_{\pm}\sim 10^{3}$  (?).

In Fig. S5, we consider two types of proton spectra:

and

The normalization factor $N_{\rm 0,p}$ determines the total energy of protons. The target density for the hadronic interactions can be calculated from the the total mass of gas and dust in the nebula, $M=7.2\,M_{\odot}$ with $M_{\odot}$ being the solar mass  (?). For the radius of the Crab Nebula, $R=1.8$ pc, the mean number density of nucleons (independent of the gas composition) is estimated $10\,\rm cm^{-3}$. To calculate the spectra of $\gamma$-rays, we use a semi-analytical method  (?). The results are shown in Fig. S5.

The total energy in protons to explain the $\gamma$-ray flux at 1 PeV is $3\times 10^{47}\,$erg for the power-law distribution of protons, and $5\times 10^{46}\,$erg for monoenergetic 10 PeV protons. The particle escape time in the nebula determines the acceleration power of protons, $\dot{W}_{\rm p}=W_{\rm p}/t_{\rm esc}$. We consider two extreme cases: (i) protons escape ballistcally with the speed of light and (ii) propagate diffusively in the Bohm limit. The former case corresponds to the fastest escape and the timescale is comparable to the light crossing time of the entire nebula, i.e., $t_{\rm esc}=R_{\rm pwn}/c\simeq 6$ years. In the latter case, for the mean magnetic field in the nebula $B=300\,\mu$G, the Bohm diffusion coefficient is $D_{B}\simeq 10^{27}(E_{p}/10\,{\rm PeV})(B/300\,\mu\rm G)^{-1}\rm cm^{2}s^{-1}$. The diffusive escape time is estimated as $t_{\rm esc}=R_{\rm pwn}^{2}/4D_{B}\simeq 250(E_{p}/10\,{\rm PeV})^{-1}(B/300\,\mu\rm G)$ years. Correspondingly, the luminosity of 10 PeV protons estimated in the broad limits $\dot{W_{\rm p}}=W_{\rm p}/t_{\rm esc}\sim(7-260)\times 10^{36}(W_{p}/5\times 10^{46}\rm\,ergs)\,\rm erg~{}s^{-1}$. The required proton acceleration rate constitutes 1% of the present spin-down luminosity of the Crab pulsar only in the case of the most effective confinement, corresponding to the propagation in the extreme Bohm diffusion regime.