<figure><img src="image_1.png" /></figure>

HAL
open science

# Phenomenological modelling of the Crab Nebula's broad-band energy spectrum and its apparent extension

Ludmilla Dirson, Dieter Horns

► To cite this version:

Ludmilla Dirson, Dieter Horns. Phenomenological modelling of the Crab Nebula's broad-band energy spectrum and its apparent extension. Astronomy and Astrophysics - A&A, 2023, 671, pp.A67. 10.1051/0004-6361/202243578. [hal-03630649](https://hal.science/hal-03630649)

HAL Id: hal-03630649

https://hal.science/hal-03630649

Submitted on 10 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Phenomenological modelling of the Crab Nebula’s broadband energy spectrum and its apparent extension

L. Dirson<fnref n="1" />,<fnref n="2" /> and D. Horns<fnref n="1" />

Received 17 March 2022 / Accepted 21 December 2022

## ABSTRACT

**Context.** The Crab Nebula emits exceptionally bright non-thermal radiation across the entire wavelength range from the radio to the most energetic photons. So far, the underlying physical model of a relativistic wind from the pulsar terminating in a hydrodynamic standing shock has remained fairly unchanged since the early 1970s when it was first introduced. One of the predictions of this model is an increase in the toroidal magnetic field downstream from the shock where the flow velocity drops quickly with increasing distance until it reaches its asymptotic value, matching the expansion velocity of the nebula.

**Aims.** The magnetic field strength in the nebula is poorly known. Using the recent measurements of the spatial extension and improved spectroscopy of the gamma-ray nebula, it has become –for the first time – feasible to determine in a robust way both the strength as well as the radial dependence of the magnetic field in the downstream flow.

**Methods.** In this work, we introduce a detailed radiative model which was used to calculate the emission from non-thermal electrons (synchrotron and inverse Compton) as well as from thermal dust present in the Crab Nebula in a self-consistent way to compare it quantitatively with observational data. Special care was given to the radial dependence of the electron and seed field density.

**Results.** The radiative model was used to estimate the parameters related to the electron populations responsible for radio and optical/X-ray synchrotron emission. In this context, the mass of cold and warm dust was determined. A combined fit based upon a $\chi^2$ minimisation successfully reproduced the complete data set used. For the best-fitting model, the energy density of the magnetic field dominates over the particle energy density up to a distance of $\approx 1.3 r_s$ ($r_s$: distance of the termination shock from the pulsar). The very high energy (VHE: $E > 100$ GeV) and ultra-high energy (UHE: $E > 100$ TeV) gamma-ray spectra set the strongest constraints on the radial dependence of the magnetic field, favouring a model where $B(r) = (264 \pm 9) \mu G(r/r_s)^{-0.51 \pm 0.03}$. For a collection of VHE measurements during epochs of higher hard X-ray emission, a significantly different solution $B(r) = (167 \pm 5) \mu G(r/r_s)^{-0.29(+0.03,-0.06)}$ is found.

**Conclusions.** The high energy (HE: $E > 100$ MeV) and VHE gamma-ray observations of the Crab Nebula lift the degeneracy of the synchrotron emission between particle and magnetic field energy density. The reconstructed magnetic field and its radial dependence indicates a ratio of Poynting to kinetic energy flux $\sigma \approx 0.1$ at the termination shock, which is $\approx 30$ times larger than estimated up to now. Consequently, the confinement of the nebula would require additional mechanisms to slow the flow down through, for example, excitation of small-scale turbulence with possible dissipation of the magnetic field.

**Key words.** radiation mechanism: non-thermal – astroparticle physics – magnetohydrodynamics (MHD) – ISM: supernova remnants – gamma rays: general – dust, extinction

## 1. Introduction

The Crab Nebula and pulsar are the centrepieces of gamma-ray astronomy and high energy (HE) astrophysics in general. The discovery of its very high energy (VHE) emission by the pioneering observations with the Whipple air Cherenkov telescope (Weekes et al. 1989) marks a breakthrough in the field by establishing the Crab Nebula as the brightest steady gamma-ray source.

The current understanding of the nebula was already established several decades ago by Rees & Gunn (1974). The description of the dynamics of the nebula and its radiative model dominated by synchrotron emission of relativistic electrons (Pacini & Salvati 1973) has remained basically unchanged for the past 50 years. In this model, the nebula is powered by a continuous cold and ultra-relativistic pulsar wind. The wind terminates at a distance of about 0.13 pc (Weisskopf et al. 2012) to the pulsar in a standing shock. The position of the termination shock is commonly associated with the bright ring observed at soft

X-rays which has a major axis of $(13.3 \pm 0.2)''$ (Weisskopf et al. 2012)$^1$.

The slow expansion of the nebula requires the downstream flow to slow down. In an ideal magneto-hydrodynamic (MHD) setting as suggested by Kennel & Coroniti (1984), the ratio of magnetic to kinetic energy ($\sigma$ parameter) upstream of the shock is estimated to $\sigma = 0.003$ to match the boundary conditions of the slow expansion of the outer nebula of $\approx 1500 \text{ km s}^{-1}$. In the downstream flow, the toroidal (‘wound up’) magnetic field would increase due to magnetic flux conservation with increasing distance until equipartition is reached.

However, with the three-dimensional treatment of the MHD problem (Porth et al. 2013) it has been demonstrated that the magnetic field structure and magnetic field present in the plasma very likely deviate from the simplified model of Kennel & Coroniti (1984). Additionally, the conversion of ordered fields to turbulent fields which would be dissipated in

<footnote n="1">¹ Institute for Experimental Physics, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany
e-mail: dieter.horns@uni-hamburg.de</footnote>

<footnote n="2">² Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France
e-mail: ludmilla.dirson@astro.unistra.fr</footnote>

<footnote n="1">¹ The pulsar is slightly offset by 0.9" from the centre of the ellipse.</footnote>
non-ideal MHD plasma could change the magnetic field configuration substantially (Bucciantini et al. 2017; Tanaka et al. 2018).

The spatially resolved degree of polarisation and position angle in the radio (e.g., Bietenholz & Kronberg 1990) and optical (e.g., Hickson & van den Bergh 1990) indicate the presence of large-scale toroidal fields in the region of the torus, while at larger distances the field orientation becomes predominantly radial as traced out in the radio. The fraction of polarised X-ray emission from the Crab Nebula was measured to be (19.2 ± 1)% (Weisskopf et al. 1978; Feng et al. 2020; Liu & Wang 2021), which indicates that the magnetic field energy is roughly equally shared between an ordered and small-scale turbulent field (Bucciantini et al. 2017). A more detailed view of the spatial distribution of the magnetic fields will be possible with the upcoming observations of the Crab Nebula with the Imaging X-ray Polarimetry Explorer (IXPE; Weisskopf et al. 2016).

With the observation of (inverse Compton) generated gamma-ray emission, the ideal MHD flow-type solution (Kennel & Coroniti 1984) has been confirmed to provide a reasonable description of the nebula's gamma-ray spectrum (Atoyan & Aharonian 1996). Improved spectroscopy at VHEs demonstrates however, that the simple approach does not provide a reasonable description of the data anymore (Aleksić et al. 2015). Even though the ideal MHD flow solution and its siblings remain as benchmarks (see e.g., Abdalla et al. 2020), the observational data are of sufficient quality to test the underlying assumptions of the simple flow model already developed in the 1970s.

So far, it has not been possible to spatially resolve the magnetic field strength in the nebula or any other pulsar wind nebulae. High energy and VHE gamma-ray observations have been used to estimate an average magnetic field from the ratio of synchrotron to inverse Compton emission (Aharonian et al. 2004), which requires careful modelling of the seed photon fields in the nebula (Meyer et al. 2010). Even though this method is robust, observational details, which are required to match the precision of the observational data with sufficiently accurate model calculations, have been mostly omitted in past studies.

In this analysis, we take the next step and include the recent measurement of the spatial extension at gamma rays of the nebula (Yeung & Horns 2019; Abdalla et al. 2020) in combination with the most recent spatial and spectral multi-wavelength data to reconstruct, for the first time, the electron and magnetic field configuration in the nebula in a consistent phenomenological model without any underlying assumptions as to the flow properties.

The underlying parameters were estimated with a $\chi^2$ minimisation method. This is the first study of this type carried out on any gamma-ray emitting source. Thanks to the exceptional multi-wavelength observational data, the Crab Nebula remains a corner stone object for the study of non-thermal processes in relativistic plasma.

## 2. Observations

In the following, the heliocentric distance of the Crab Nebula is assumed to be $d_{\text{Crab}} = 2$ kpc (Trimble 1973). The basis of the compilation of multi-wavelength data used here is summarised by Aharonian et al. (2004), Meyer et al. (2010) and references therein. In the past decade, both new observations (see Table 1) as well as refined and relevant corrections on the extinction in the interstellar medium have become available.

<table><caption>**Table 1.** Summary of updated observations.</caption><thead><tr><th>Energy</th><th>Instrument</th><th>References</th></tr></thead><tbody><tr><td rowspan="2">Radio</td><td>WMAP</td><td>(1)</td></tr><tr><td>Planck for the HFI instrument</td><td>(2)</td></tr><tr><td rowspan="2">Sub-mm/IR</td><td>Herschel</td><td>(3)</td></tr><tr><td>WISE</td><td>(3)</td></tr><tr><td rowspan="3">X-ray to γ-ray</td><td>Crab</td><td>(4)</td></tr><tr><td>NuSTAR</td><td>(5)</td></tr><tr><td>Fermi-LAT</td><td>(6)</td></tr><tr><td rowspan="6">VHE</td><td>H.E.S.S.</td><td>(7)</td></tr><tr><td>VERITAS</td><td>(8)</td></tr><tr><td>MAGIC</td><td>(9)</td></tr><tr><td>Tibet ASγ</td><td>(10)</td></tr><tr><td>HAWC</td><td>(11)</td></tr><tr><td>LHAASO WCDA, KM2A</td><td>(12)</td></tr></tbody></table>

**Notes.** All other data shown e.g., in Fig. 7 are taken from Meyer et al. (2010) and references therein.

**References.** (1) Weiland et al. (2011), (2) Planck Collaboration XXVI (2016), Ritacco et al. (2018), (3) Gomez et al. (2012), De Looze et al. (2019), (4) Toor & Seward (1974), (5) Madsen et al. (2017), (6) Buehler et al. (2012), Yeung & Horns (2019), (7) Aharonian et al. (2006), Holler (2014), Holler et al. (2015), H. E. S. S. Collaboration (2014), (8) Wells (2019), Fig. 5.15, Aliu et al. (2014), (9) Aleksić et al. (2015), MAGIC Collaboration (2020), (10) Amenomori et al. (2019), (11) Abeysekara et al. (2019), (12) LHAASO Collaboration (2021).

### 2.1. Radio and infrared observations

The transition between the radio synchrotron component and the thermal dust emission is explored through observations with the Herschel (PACS and SPIRE), Spitzer (MIPS), while the transition from the dust emission to the optical/wind synchrotron emission is traced with WISE and Spitzer (IRAC) multi-band photometry. The integrated aperture-photometry (integrating over an ellipse centred on the pulsar position, with PA = 40° and semi-major and minor axes of 245" × 163") is taken from De Looze et al. (2019), Table 3. In addition to the continuum component, the emission lines contribute up to 18% of the continuum flux in the PACS $\lambda = 100$ µm pass band. The resulting data set is shown in Fig. 1 in combination with the best-fitting model. The new and additional data taken from the compilation of De Looze et al. (2019) are marked with white circles in Fig. 1.

### 2.2. Optical observations

In comparison to previous works (Meyer et al. 2010), the optical data used here have been updated to include a more recent estimate of the extinction along the line of sight towards the Crab Nebula. In the original compilation of Aharonian et al. (2006), the extinction was corrected using $A(V) = 1.5$ mag with a colour excess $E(B-V) = 0.47$ mag guided by the discussion of Véron-Cetty & Woltjer (1993). More recent 3D-reconstruction of the Galactic extinction with Gaia, Pan-STARRS 1, and 2MASS observations lead to a smaller value of $A(V) = 1.2$ mag (Green et al. 2019). We chose here the central value of $A(V) = (1.08 \pm 0.38)$ mag and $E(B-V) = 0.43$ mag determined by De Looze et al. (2019) from modelling the dust emission along the line of sight towards the Crab Nebula. The wavelength dependent extinction was calculated with the extinction curve from O'Donnell (1994). The resulting shape of the visual
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Observed spectral energy distribution (SED) from radio to optical photon energies. The additional data from the recent compilation of photometric data from De Looze et al. (2019) are marked with white circles. Optical/UV data are presented with and without extinction correction (see text for more details). The SED has been fit with the superposition of thermal dust emission (see Table 3 for a list of parameters and uncertainties) and the underlying synchrotron continuum (see Table 4).</figcaption></figure>

continuum emission follows a power law, such that $I_ν ∝ ν^{-s}$ with $s_O ≈ 0.8$. This is considerably harder than the X-ray spectrum with $s_X ≈ 1.3$. The far-UV spectroscopy of the entire nebula with the Hopkins UV telescope (Blair et al. 1992) confirms the extrapolation of the optical spectrum to the UV. The UV spectrum obtained in a slit aperture of 17" × 116" displayed in Fig. 1 was corrected by an ad hoc correction factor of 7.5 to match the overlapping photometric points of Hennessy et al. (1992).

Additionally, we extracted the spectral indices $s(9241–5364)$ obtained between $\lambda = 9241$ Å and $\lambda = 5364$ Å from Fig. 1 of Véron-Cetty & Woltjer (1993). The spectral indices were determined within bins of 10" × 10" covering the optical nebula. By averaging over annuli centred on the pulsar position, we calculated the mean and root mean square of $s_O$. The indices change from ≈0.6 in the centre to ≈1 at the periphery of the nebula (see Fig. 2). This is well-matched by the predicted behaviour of the synchrotron model as indicated by the superimposed blue line in the same figure.

## 2.3. X-ray observations

The Crab Nebula has a long history as a standard candle for X-ray astronomy (Toor & Seward 1974). Therefore, for many X-ray telescopes the instrumental response functions are tuned to reproduce a particular spectral shape that is deemed to be a canonical X-ray spectrum following a power law $dN/dE = N_0(E/\text{keV})^{-\Gamma}$, with $\Gamma = 2.11$ and $N_0 = 11$ in units of $\text{cm}^{-2} \text{s}^{-1} \text{keV}^{-1}$, readers can refer to Shaposhnikov et al. (2012), for example.

In this spirit, the X-ray spectrum of the Crab Nebula has been used to cross-calibrate different instruments (Kirsch et al. 2005; Weisskopf et al. 2010). Given the apparent brightness of the nebula, other weaker non-thermal X-ray sources are conveniently used for more sensitive instruments which suffer from saturation/pile-up for bright sources such as the Crab Nebula. In a recent study of the non-thermal X-ray source G21.5-0.9 (Tsujimoto et al. 2011), the findings of earlier comparisons were

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Predicted and observed optical spectra of the Crab Nebula. In concentric annuli, the mean and root mean square of the spectral indices $s_O (I_ν ∝ ν^{-s_O})$ have been determined from the data displayed in Fig. 1 of Véron-Cetty & Woltjer (1993). The superimposed curve is calculated from the intensity expected from the synchrotron model as described in Sect. 4.3 with the best-fitting parameters listed in Table 4. It is important to note that these data have not been used in the fit.</figcaption></figure>

confirmed that the energy spectra measured with different instruments in the energy range from 2 keV to 8 keV show relative differences in the normalisation of up to 20% and differences in the spectral index of $\Delta\Gamma = 0.1$.

Different to most other X-ray instruments in use, the pn-camera of the XMM-Newton X-ray observatory has been calibrated without relying on the Crab Nebula as a standard candle (Strüder et al. 2001). The normalisation $N_0 = 8.6$ (Kirsch et al. 2005) found in observations of the nebula with the pn-camera in burst mode is consistently smaller than the normalisation found with other instruments.

The challenging task to (cross-)calibrate the past, current, and future X-ray missions is carried out by a consortium that formed shortly after the publication of Kirsch et al. (2005). The International Astronomical Consortium for High-Energy Calibration (IACHEC) recognises these systematic differences (Madsen et al. 2021) but there is still no consensus which measurement can be considered to be more reliable.

In the previous works (Aharonian et al. 2006; Meyer et al. 2010), we chose to re-scale measurements at higher energies to match the XMM-Newton result. In view of the recent measurements of the nebula's X-ray spectrum with the NuSTAR telescope, we updated this approach. The absolute flux measurement reported by Madsen et al. (2017) uses stray light from the Crab Nebula in the focal plane of the NuSTAR telescope. Since the X-rays do not pass from the optical mirror assembly, the resulting flux determination does not suffer from the uncertainty of the energy dependent efficiency of the mirror assembly. Subsequently, the instrumental team of NuSTAR has changed their canonical Crab spectrum to $\Gamma = 2.103 \pm 0.001$ and $N_0 = 9.69 \pm 0.02$. This normalisation is 14% larger than the previously obtained value for the pn camera (Kirsch et al. 2005). Moreover, the NuSTAR measurement extends to hard X-rays and is consistent with, for example, the results from the coded mask imaging detector and spectrometer SPI on the INTEGRAL observatory (Jourdain & Roques 2009) as well as Earth occultation data taken with the BATSE instrument (Ling & Wheaton 2003). Both measurements indicate a break of the spectrum at an energy of approximately 100 keV.
Beyond the hard X-ray band, gamma-ray measurements rely on Compton scattering instead of photo-ionisation. The Comptel telescope on-board the Compton-Gamma-Ray-Observatory (CGRO) observed the Crab Nebula for 9 years (Kuiper et al. 2001). The resulting time-averaged energy spectrum between 0.75 MeV and 30 MeV follows a power law with photon index Γ = 2.227 ± 0.013 (Kuiper et al. 2001). While the power-law index is consistent with the one determined with INTEGRAL-SPI and BATSE data at lower energies, the flux normalisation of the Comptel measurement is systematically lower than the one obtained with BATSE and SPI data. From Fig. 31 of Kuiper et al. (2001), we estimated possible relative systematic deviations on the collection area from two different calibration methods to be 20%. In the data sample used to characterise the synchrotron component of the nebula for the fit, we introduced an ad hoc scaling of the Comptel flux by 1.3.

Even though this scaling is larger than the systematic uncer-
tainty of the Comptel flux, it is well within the range of the com-
bined systematic uncertainty when cross-calibrating between the
very different types of detector. Clearly, the MeV energy range is
most challenging and even though heroic efforts are under way
to re-analyse the complete Comptel data set (Strong & Collmar
2019), a new mission is required to improve the situation in
this poorly explored energy range (see e.g., the COSI mission
Tomsick & COSI Collaboration 2022).

An additional concern in constructing the X-ray SED of the Crab Nebula is the variability of its X-ray (Greiveldinger & Aschenbach 1999; Ling & Wheaton 2003; Wilson-Hodge et al. 2011) and gamma-ray flux (Tavani et al. 2011; Abdo et al. 2011). The X-ray variability has been observed to occur mainly between 2001 and 2010 with an apparent relative amplitude of 2(4)% in the energy band <15(15–50) keV on timescales of ≈3 years. In the preceding years of monitoring the Crab Nebula with the RXTE mission, the flux appears to remain constant (Wilson-Hodge 2012). The X-ray spectrum softened between 2009 and 2011 with a photon index changing from 2.14 to 2.17 as measured with RXTE-PCA (Wilson-Hodge 2012). We show in Fig. A.1 the combined hard X-ray light-curve between 2000 and beginning of 2022.

When comparing the recent average flux with the histori-
cal measurements of about 50 years ago (Toor & Seward 1974),
possible long-term flux variations do not exceed the variations
observed during the 20 years covered by Wilson-Hodge (2012).
The origin of the variability is not clear but a systematic effect
can be firmly rejected as the sole reason given that the variability
has consistently been observed with independent instruments.

On a shorter timescale and at larger energies, the Crab Neb-
ula shows rapid variations of the flux between 100 MeV and
approximately 1 GeV. The flux can both increase by a factor
of ≈8 over timescales of days (Buehler et al. 2012) and it can
also drop below the average flux within a similar timescale
(Yeung & Horns 2020).

For the characterisation of the time-average synchrotron spectral energy distribution we assume that the X-ray variability is averaged out during the long exposure. It is to be expected that the variability measured at soft to hard X-rays should translate in a similar variability in the VHE gamma-ray flux at TeV energies. So far, multi-year observations of the nebula with ground-based gamma-ray instruments such as HEGRA have only constrained the relative variability to be smaller than 10% (Aharonian et al. 2006).

The more volatile situation at the end point of the syn-
chrotron spectrum is unlikely to change the inverse Compton
emission of the nebula. This has been confirmed by simulta-

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** SED of the Crab Nebula centred on the X-ray energy range. The best-fitting model (blue line, see details on the fitting procedure in Sect. 4 and Table 4) is shown along with measurements from different X-ray observatories. The normalisation of the spectrum measured with *XMM-Newton* (yellow band) has been scaled up by 16% to be consistent with the absolute flux measurement with *NuSTAR* (magenta band). The COMPTEL gamma-ray flux has been scaled up by 30% to match the extrapolation of the SPI spectrum.</figcaption></figure>

neous observations of flares with Fermi-LAT and ground-based
gamma-ray telescopes: during these observations no correlated
flux changes have been found (H. E. S. S. Collaboration 2014;
Aliu et al. 2014). There may be however a measurable change of
the PeV gamma-ray flux during high flux or flaring states (see
Sect. 4.4.4).

In Fig. 3, the data points used for the fit are shown as
blue markers together with measurements listed above. The blue
line indicates the best-fitting synchrotron model which includes
the Fermi-LAT data points between 70 MeV and 500 MeV (see
below).

With the scaling of the flux measurements to match the abso-
lute flux measurement of NuSTAR applied to the XMM-Newton
and COMPTEL measurement, we obtain a consistent spectral
measurement from 0.2 keV to 30 MeV. The spectrum in this
energy range can be well characterised with a Band-like model
(Band et al. 1993)

$$
\frac{dN}{dE} = N_0 \left\{ \begin{array}{ll} \left( \frac{E}{100\text{keV}} \right)^{\Gamma_1} e^{-E/E_0} & \text{if } E \le E_0(\Gamma_1 - \Gamma_2), \\ \left[ \frac{(\Gamma_1-\Gamma_2)E_0}{100\text{keV}} \right]^{\Gamma_1-\Gamma_2} e^{\Gamma_2-\Gamma_1} \left( \frac{E}{100\text{keV}} \right)^{\Gamma_2} & \text{else} \end{array} \right. \quad (1)
$$

with parameters for the photon indices Γ₁ = -2.095(4), Γ₂ = -2.34(1), E₀ = 1.36(17) MeV and normalisation at 100 keV: N₀ = 6.40(8) × 10⁻⁴ cm⁻² s⁻¹ keV⁻¹. The resulting χ²(d.o.f.) = 97(158) indicates a good fit of the Band function. Similar parameters (Γ₁ = -1.99(1), Γ₂ = -2.31(2), E₀ = 531(30), N₀ = 7.5(2)×10⁻⁴ cm⁻² s⁻¹ keV⁻¹) have been obtained from analysing 17 years of INTEGRAL SPI data (Jourdain & Roques 2020), see also Fig. 3.

2.4. HE and VHE gamma-ray observations

Since the work of Meyer et al. (2010), roughly nine more
years of Fermi-LAT observations have been presented in recent
publications. Therefore, the statistical uncertainty on the flux
above 100 GeV has been reduced by a factor of three. The
updated instrumental response function of Fermi-LAT extends
the HE reach from formerly 300 GeV to 500 GeV. The system-
atic uncertainty on the energy scale has been reduced from pre-
viously ΔE/E = +5% (see e.g., Abdo et al. 2009) to ΔE/E = +2%
(Ackermann et al. 2012).
In a recent investigation of the gamma-ray extension of the Crab Nebula, Yeung & Horns (2019) analysed the LAT data accumulated over ≈9.1 years with a refined model for the Crab pulsar’s spectrum. They obtain a >5 GeV spectrum of FL8Y J0534.5+2201i (the PWN component) which essentially matches the off-pulse spectrum reported by Buehler et al. (2012) but with improved statistical uncertainties and wider energy reach.

The MAGIC collaboration has published an updated energy spectrum of the Crab nebula overlapping the energy range covered with Fermi-LAT and reduced uncertainties at energies below 100 GeV (Aleksić et al. 2015). The measurement of the flux from the Crab Nebula with the MAGIC instruments spans for the first time from 50 GeV to 30 TeV. In a subsequent analysis of data taken with the MAGIC telescope at ultra-large zenith angles, the energy reach has been extended up to 100 TeV (MAGIC Collaboration 2020).

The VHE end of the Nebula’s spectrum has recently become a focus of interest with observational data taken with non-imaging air shower detectors with improved analysis of data taken with the HAWC water Cherenkov detector (Abeysekara et al. 2019) or with long-term observations with the Tibet ASγ air shower array which extends the energy spectrum up to 450 TeV (Amenomori et al. 2019). The latest addition to results on the Crab nebula obtained with the LHAASO KM2A array extends the spectrum beyond 1 PeV (LHAASO Collaboration 2021).

The High Energy Stereoscopic System (H.E.S.S.) collaboration has presented their preliminary results of H.E.S.S. phase II observations of the Crab Nebula Holler et al. (2015) superseding the results obtained with the previous setup (phase I) which consisted of four imaging Cherenkov telescopes (IACTs) with a mirror surface area of about 100 m² each and an energy threshold of approximately 100 GeV. Due to the visibility of the Crab Nebula from the southern hemisphere, the energy spectrum obtained with the phase I instruments from the Crab Nebula covers the energy range from 440 GeV up to 40 TeV (Aharonian et al. 2006). The large telescope of the H.E.S.S. phase II with a mirror surface of 600 m² has lead to a lowered energy threshold of 230 GeV for observations of the Crab Nebula (Holler et al. 2015).

Finally, the VERITAS collaboration has reported on observations of the Crab Nebula (55 h before and 46 h after the camera update in 2012; Wells 2019). In a previous analysis, slightly more data were used (Meagher & VERITAS Collaboration 2016), but the investigation of Wells (2019) indicate that the calibration of the two different levels of amplification for each channel need to be adapted. More detailed information on the VHE data sets is provided in Appendix A.

## 3. Non-thermal emission model

The specific luminosity $L_ν$ from the Crab Nebula was calculated by integrating the synchrotron and IC emissivities ($j_ν^{Sy}$, $j_ν^{IC}$) over the volume of the nebula (assuming optically thin plasma and spherical symmetry):

$$L_ν = 4π \int_V d^3r (j_ν^{Sy}(r) + j_ν^{IC}(r)), \quad (2)$$

where the specific emissivity for inverse Compton emission ($j_ν^{IC}$) and synchrotron emission ($j_ν^{Sy}$) is given by

$$j_ν^{IC}(r) = \frac{1}{4π} \int_{γ} dγ n_{el}(r, γ) \frac{3}{4} \frac{\sigma_T c}{γ^2} h^2 v \int_{hv/4γ^2}^{hv} \frac{dε}{ε} f_{IC}(ε, v, γ) n_{seed}(r, ε), \quad (3)$$

$$j_ν^{Sy}(r) = \frac{1}{4π} \int_γ dγ n_{el}(r, γ) \frac{\sqrt{3}e^3 B(r)}{mc^2} G\left(\frac{ν}{ν_c}\right), \quad (4)$$

with $n_{seed}(r, ε)$ and $n_{el}(r, γ)$ being the differential number densities of seed photons and electrons respectively.

The resulting specific intensity $I_ν(θ)$ was calculated by integrating along a line of sight subtending an angle $θ$ with the centre of the nebula. The source region is optically thin. The external absorption in the interstellar medium was used to correct the apparent flux measurements in the optical/UV (de-reddening, see Sect. 2.2), X-ray (photo-ionisation, see Sect. 2.3), as well as at ultra-high energy (>100 TeV) gamma-rays. The latter effect was re-evaluated here taking into account the anisotropy of the local radiation field which leads to an increase of the optical depth in comparison to previous works (see Appendix B for more details).

The function $G(x)$ used for the synchrotron emissivity for an isotropic electron distribution in random magnetic fields (Crusius & Schlickeiser 1986) has been recently expressed in terms of modified Bessel functions as given in Eq. (D5) in Aharonian et al. (2010). For the numerical calculations, we used the convenient approximation of Eq. (D7) in Aharonian et al. (2010).

The function $f_{IC}$ is the general function for scattering of electrons in an isotropic photon gas derived by Jones (1968; Eq. (44)).

In order to calculate the emissivity of the electrons we need to determine the electron distribution, seed fields, and magnetic field structure to reproduce all available observations. While the observed synchrotron emission is mostly degenerate for the choice of magnetic field and particle distribution, the observed inverse Compton emission breaks the degeneracy to allow us to reconstruct both the strength and the spatial variation of the magnetic field in the nebula.

### 3.1. Distribution of electrons

The spectral and spatial distribution of electrons in the nebula is inferred by matching a phenomenological model of the electron distribution and its synchrotron emission to observations. The electron distribution was kept as simple as possible and as complex as necessary to match the observations.

We assumed two different populations of electrons, so-called radio and wind electrons which do not share a common origin. While the wind electrons are most likely accelerated at the wind termination shock, the origin of the radio emitting electrons is not well understood. This includes both the generation of sufficient number of electrons/positrons as well as their acceleration. Recent suggestions include stochastic acceleration (Tanaka & Asano 2017), for example, in reconnection layers in the turbulent nebula (Lyutikov et al. 2019).

#### 3.1.1. Spatial distribution of electrons

For photon energies below about 70 keV, the synchrotron emission is well-resolved by the available instruments. The observed (projected) structures are characterised by small-scale features (e.g., wisps) which evolve dynamically on timescales of weeks as well as large scale features including the X-ray bright ring, torus, and jets which remain stable for longer time-periods of years.

For the purpose of modelling the energy spectrum and extension of the inverse Compton emission of the nebula, we
followed an approach similar to de Jager & Harding (1992), Hillas et al. (1998), Meyer et al. (2010) where the electron distribution is found by matching the expected synchrotron emission to the observed broadband spectral energy distribution (SED). In a second step, the magnetic field was determined by adjusting the electron distribution (mainly the normalisation) to match the observed inverse Compton emission to the data. In order to keep the model simple, the source was assumed to be spherical. Given the coarse angular resolution of the instruments used to observe the nebula at the HE end, this simplification is acceptable.

An important benefit of the assumption on spherical symmetry is the simple characterisation of the extension of the nebula. Here, we used the apparent angular extension $\theta_{68}$ which is defined by the angular radius which encompasses 68% of the total flux. For a helio-centric distance $d_{Crab}$ of the Crab Nebula, the spatial radius $R_{68}$ relates to the apparent size $R_{68} \approx \theta_{68} d_{Crab}$. We assumed a distance $d_{Crab} = 2$ kpc.

We collected available observations in order to determine $\theta_{68}$ for different photon energies. The resulting values of $\theta_{68}$ are summarised in Table 2. The *Chandra* data points<fnref n="2" /> indicate the 68% containment angular radius $\theta_{68}$ of the Crab Nebula for different energy bands. For each energy band (0.5–1.2; 1.2–2; 2–7; 7–10 keV), the centroid of the background subtracted image was calculated. The value of $\theta_{68}$ for each energy band was determined using the method of encircled count fraction (ecf_CALC task of the CIAO package v4.13). The uncertainty on $\theta_{68}$ is dominated by systematic uncertainties which were estimated by checking against effects of vignetting and background contamination. The NuSTAR X-ray telescope has been used to image the X-ray nebula for the first time in a broad energy range from 3 keV–78 keV (Madsen et al. 2015). We estimated $\theta_{68}$ from their Fig. 13 in three energy ranges. In the lowest energy band (3 keV–5 keV), the extension derived from the NuSTAR observations is about 5% smaller than the values obtained with the *Chandra* image, well within the estimated uncertainties of the two measurements.

The differential electron number density $dn = n_{el}(\gamma, r)dydV$ was assumed to follow a multiple broken power law in Lorentz factor $\gamma$ (see below) with a spatial distribution approximated by a radial Gaussian with a scale length $\rho(\gamma)$ that depends on the electrons' energy:

$$n_{el}(\gamma, r) = n(\gamma)e^{-\frac{r^2}{2\rho(\gamma)^2}}, \quad (5)$$

such that the predicted synchrotron emission matches the observed synchrotron extension.

The observed synchrotron size is mostly decreasing with increasing frequency, reflecting the cooling of the electrons as they propagate outwards. Therefore, we assumed for the wind-electrons a power law for $\rho(\gamma)$:

$$\rho(\gamma) = \frac{\rho_0}{3600''} \cdot \frac{\pi}{180^\circ} \cdot \frac{d_{Crab}}{2\,\text{kpc}} \left( \frac{B_0}{264\,\mu\text{G}} \left( \frac{\gamma}{9 \times 10^5} \right)^2 \right)^{-\beta}, \quad (6)$$

with $\rho_0$ (in arc sec), $B_0$, and $\beta > 0$ free parameters in the fit to the data (see below).

Different to the wind electrons, for the low energy, radio-emitting electrons, the radial scale length $\rho_r$ is independent of the electrons' Lorentz factor. The value of $\rho_r$ was left to vary freely in the fit to match the measured synchrotron extension.

<table><caption>**Table 2.** Apparent extension of the Crab Nebula.</caption><thead><tr><th>ν [Hz]</th><th>θ<sub>68</sub>["]</th><th></th></tr></thead><tbody><tr><td>5 × 10<sup>9</sup></td><td>117 ± 1</td><td>VLA (1)</td></tr><tr><td>150 × 10<sup>9</sup></td><td>115 ± 2</td><td>NIKA (2)</td></tr><tr><td>1.8 × 10<sup>12</sup></td><td>112 ± 5</td><td>PACS 160 µm (3)</td></tr><tr><td>4.3 × 10<sup>12</sup></td><td>112 ± 5</td><td>PACS 60 µm (3)</td></tr><tr><td>1.4 × 10<sup>13</sup></td><td>137 ± 5</td><td>WISE 22 µm (4)</td></tr><tr><td>2.5 × 10<sup>13</sup></td><td>119 ± 5</td><td>WISE 12 µm (4)</td></tr><tr><td>6.5 × 10<sup>13</sup></td><td>140 ± 5</td><td>WISE 4.6 µm (4)</td></tr><tr><td>8.9 × 10<sup>13</sup></td><td>141 ± 5</td><td>WISE 3.3 µm (4)</td></tr><tr><td>1.03 × 10<sup>15</sup></td><td>99 ± 5</td><td>UVOT UVW1 (5)</td></tr><tr><td>2.66 × 10<sup>17</sup></td><td>46 ± 4</td><td>Chandra (0.5–1.5) keV (6)</td></tr><tr><td>3.7 × 10<sup>17</sup></td><td>43 ± 3</td><td>Chandra (1.5–3.0) keV (6)</td></tr><tr><td>8 × 10<sup>17</sup></td><td>39 ± 3</td><td>Chandra (3–6) keV (6)</td></tr><tr><td>1.2 × 10<sup>18</sup></td><td>37 ± 3</td><td>NuSTAR (5–6) keV (7)</td></tr><tr><td>2.4 × 10<sup>18</sup></td><td>33 ± 3</td><td>NuSTAR (8–12) keV (7)</td></tr><tr><td>1.2 × 10<sup>19</sup></td><td>27 ± 3</td><td>NuSTAR (35–80) keV (7)</td></tr><tr><th>E<sub>γ</sub> [GeV]</th><th>θ<sub>68</sub>["]</th><th></th></tr><tr><td>7</td><td>121 ± 47</td><td>Fermi-LAT (8)</td></tr><tr><td>14</td><td>209 ± 22</td><td>(8)</td></tr><tr><td>28</td><td>125 ± 26</td><td>(8)</td></tr><tr><td>56</td><td>114 ± 14</td><td>(8)</td></tr><tr><td>112</td><td>95 ± 29</td><td>(8)</td></tr><tr><td>214</td><td>83 ± 26</td><td>(8)</td></tr><tr><td>990</td><td>121 ± 38</td><td>(8)</td></tr><tr><td>1050</td><td>79 ± 10</td><td>H.E.S.S. (9)</td></tr></tbody></table>

**Notes.** The apparent angular size of the nebula was calculated by integrating the observed intensity over the solid angle up to $\theta_{68}$ which is the percentile that encloses 68% of the total flux. Whenever needed, point sources were excluded and background subtracted.

**References.** (1) Bietenholz et al. (2004), (2) Fig. 5 of Ritacco et al. (2018), (3) Herschel Science Archive (HSA), (4) Infrared Science Archive (IRSA), (5) XMM-Newton Science Archive (XSA), (6) *Chandra* Data archive, (7) Madsen et al. (2015), (8) Yeung & Horns (2019), (9) H. E. S. S. Collaboration (2020).

### 3.1.2. Spectral distribution of electrons

We considered two distinctly different populations dubbed radio electrons as they are mostly emitting in the radio band and so-called wind electrons which are producing the bulk of the emission from optical to gamma-rays via synchrotron processes:

$$n_{el}(r, \gamma) = n_{\text{radio}}(r, \gamma) + n_{\text{wind}}(r, \gamma). \quad (7)$$

The differential number density of each population was assumed to follow as a power law for the radio electrons

$$n_{\text{radio}}(r, \gamma) = N_{r,0} \rho_r^{-3} e^{-\frac{r^2}{2\rho_r^2}} \gamma^{-s_r} H(\gamma; \gamma_0, \gamma_1), \quad (8)$$

and a broken power law for the wind electrons

$$n_{\text{wind}}(r, \gamma) = N_{w,0} \rho(\gamma)^{-3} e^{-\frac{r^2}{2\rho(\gamma)^2}} \left[ \left( \frac{\gamma}{\gamma_{w1}} \right)^{-s_1} \left( \frac{\gamma_{w1}}{\gamma_{w2}} \right)^{-s_2} H(\gamma; \gamma_{w0}, \gamma_{w1}) + \left( \frac{\gamma}{\gamma_{w2}} \right)^{-s_2} H(\gamma; \gamma_{w1}, \gamma_{w2}) + \left( \frac{\gamma}{\gamma_{w2}} \right)^{-s_3} H(\gamma; \gamma_{w2}, \gamma_{w3}) \right] \quad (9)$$

respectively. The auxiliary function $H(\gamma; \gamma_0, \gamma_1) := \Theta(\gamma_1 - \gamma)\Theta(\gamma - \gamma_0)$ with the Heaviside function $\Theta(x)$ defines the intervals of the broken power law.

<footnote n="2">² *Chandra* ObsId 13151.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 4. Electron number density $n_{el}$, as a function of $\gamma$ for different values of the distance to the centre of the Crab Nebula $r$ and the volume average. The parameters used here are listed in Table 4.</figcaption></figure>

A representative electron number density $n_{el}$ at various distances and volume averaged is shown in Fig. 4. Note, the distribution of electrons follows a phenomenological Ansatz guided by the notion that electrons lose energy (radiative and adiabatic losses) while expanding in the nebula. However, without further knowledge about the properties and dynamics of the flow and the various and competing processes at work, the structure of the nebula can only be generically explained through the idea of a shock-accelerated electron distribution in a (ideal) MHD flow of a relativistic plasma as suggested for example by Rees & Gunn (1974) and following works (e.g., Kennel & Coroniti 1984).

In the absence of a specific model that matches available observations, we followed the approach of reconstructing the energy distribution of electrons and their spatial distribution from the observations. In this approach, we do not attempt to follow the temporal evolution of the nebula through its complex and largely unknown history.

## 3.2. Seed photon fields

The seed photon fields in the nebula were self-consistently calculated and matched to the observed intensity and brightness of the nebula. Note, in many previous attempts to model the nebula, the seed photon fields are assumed to be mostly homogeneous in the nebula and not specifically modelled to match the observations.

Assuming an isotropic specific volume emissivity $j_v$ (in $4\pi$ solid angle) in a spherical source, the spectral number density $dn = n_{seed}(r, \varepsilon)dVd\varepsilon$ of seed photons at a distance $r$ to the centre of the nebula is:

$$n_{\text{seed}}(r, \varepsilon) = \frac{4\pi}{\hbar\varepsilon} \frac{1}{2c} \int_{r_{\text{min}}}^{r_{\text{max}}} \frac{r_1}{r} j_v(r_1, \varepsilon) \ln\left(\frac{r+r_1}{|r-r_1|}\right) dr_1. \quad (10)$$

The three main seed photon fields contributing to the IC scattering of the relativistic electrons in the Crab Nebula are: (i) the synchrotron radiation; (ii) the FIR 'excess' radiation which is attributed to the dust emission; (iii) the 2.7K cosmic microwave background radiation (CMB).

## 3.3. Thermal emission from dusty plasma

Observations at sub-millimetre and far-infrared wavelengths show an excess flux in comparison to the extrapolated synchrotron continuum alone. The excess is commonly modelled with the thermal emission of spherical, radiatively heated dust grains in thermal equilibrium (Marsden et al. 1984; Temim et al.

2006; Gomez et al. 2012; Temim & Dwek 2013; De Looze et al. 2019). A carbon dominated composition of the dust has been favoured in past studies where the temperature of the dust is either self-consistently modelled in combination with a dust grain composition (e.g., Temim & Dwek 2013) or estimated by fitting a two-temperature population to the data (e.g., De Looze et al. 2019). The dust emitting region is concentrated in the vicinity of the optical filaments which are forming a shell with inner radius of $r = 0.55$ pc (Čadež et al. 2004).

The observed emission in the wavelength band from $\lambda \approx 6$ µm to $\lambda \approx 1000$ µm is the superposition of continuum synchrotron and thermal dust emission as well as line-emission from excited atoms in the plasma. In Fig. 1, the observed continuum emission is shown. In all previous works, a quantitative analysis of the dust emission has relied upon an ad hoc assumption on the underlying continuum emission.

Following the argument of Owen & Barlow (2015), we assumed the dust to be predominantly composed of amorphous carbon grains. The emissivity of dust grains in thermal equilibrium with a mass density $\rho_d$ in the emitting volume is given by

$$j_v = \rho_d \kappa_{\text{abs}} B_v(T), \quad (11)$$

with $B_v(T)$ denoting the black body spectral radiant emission at temperature $T$ and $\kappa_{\text{abs}}$ the mass absorption coefficient specific to the dust grain (Temim & Dwek 2013). In order to calculate the amount of seed photons present in the nebula which is relevant for the inverse Compton emissivity, we chose to fit a minimum set of free parameters to achieve an acceptable fit to the data. The minimal model found here is a mixture of dust at two different temperatures and with different values of mass ($M_1, T_1, M_2, T_2$). An additional parameter is the outer radius $r_{\text{out}}$ of the dust population. The value of $r_{\text{out}}$ needs to be adjusted to match the observed extension which is the result of the superposition of the underlying synchrotron continuum emission and the thermal emission of the dust. With the emissivity

$$j_v = \frac{3\kappa_{\text{abs}}}{4\pi(r_{\text{out}}^3 - r_{\text{in}}^3)} [M_1 B_v(T_1) + M_2 B_v(T_2)] H(r; r_{\text{in}}, r_{\text{out}}) \quad (12)$$

of the two temperature dust model for a fixed set of grain parameters (see next paragraph and Table 3) we achieved an acceptable description of the data (see Fig. 1).

The mass absorption coefficient $\kappa_{\text{abs}}$ depends upon the geometry of the dust grain and the optical constants for the bulk material. Here, we used the measured optical constants of the ACAR sample (Zubko et al. 1996) and calculated the resulting values of $\kappa_{\text{abs}}$ assuming a spherical symmetry and following the Mie-scattering theory in a series expansion. The underlying calculation was carried out with the numerical algorithm described in Bohren & Huffman (1998), Appendix A. A more thorough self-consistent modelling of the dust emission by radiatively heated dust grains is left for further work. The grains' diameter distribution was assumed to follow a power law between the minimum grain radius $a_{\text{min}} = 10^{-3}$ µm and $a_{\text{max}} = 0.1$ µm such that the probability density of finding a grain between $a$ and $a + da$ is given by $p(a)da = p_0 a^{-s}$ with $s = 4$ for $a$ in the interval $[a_{\text{min}}, a_{\text{max}}]$ (Temim & Dwek 2013). In the relevant range of wavelengths (from $\lambda = 1$ µm to $\lambda = 10^3$ µm), the resulting opacity is conveniently approximated with

$$\kappa_{\text{abs}} = 2.15 \times 10^4 \frac{\text{cm}^2}{\text{g}} \left(\frac{\lambda}{\mu\text{m}}\right)^{-1.3}, \quad (13)$$

with relative deviations between the model (Temim & Dwek 2013) and the parametrisation of less than 5% between $\lambda = 10$ µm and $\lambda = 200$ µm.
<table><caption>**Table 3.** Parameters of the dust model.</caption><thead><tr><th>Dust parameter</th><th>Value</th></tr></thead><tbody><tr><td>Grain material</td><td>Amorph. Carbon</td></tr><tr><td>log<sub>10</sub>(M<sub>1</sub>/M<sub>☉</sub>)</td><td>-4.4 ± 0.1</td></tr><tr><td>T<sub>1</sub> [K]</td><td>149 ± 8</td></tr><tr><td>log<sub>10</sub>(M<sub>2</sub>/M<sub>☉</sub>)</td><td>-1.2 ± 0.1</td></tr><tr><td>T<sub>2</sub> [K]</td><td>39 ± 2</td></tr><tr><td>r<sub>in</sub> [pc]</td><td>0.55</td></tr><tr><td>r<sub>out</sub> [pc]</td><td>1.53 ± 0.09</td></tr><tr><td>a<sub>min</sub>, a<sub>max</sub> [µm]</td><td>10<sup>-3</sup>, 10<sup>-1</sup></td></tr><tr><td>s</td><td>4</td></tr></tbody></table>

The thermal line emission from the excited plasma was not considered here. In the visual band, the relative contribution of the line emission to the overall flux is approximately 30% when comparing the sum of the equivalent width of the lines listed in Table 2 of Smith (2003) with the wavelength-band of the measured spectrum.

### 3.4. Magnetic field structure

The widely used model for the confinement of the Crab Nebula by Kennel & Coroniti (1984) is based upon a toroidal magnetic field without any turbulence forming. Three-dimensional relativistic ideal MHD simulations supports the notion of disordering of the magnetic field (Porth et al. 2013) which has been argued to account for the slowing of the expansion velocity to match the observations (Tanaka et al. 2018) and the comparably small polarisation observed at X-rays (Bucciantini et al. 2017; Weisskopf et al. 1978). In non-ideal MHD plasma flows, dissipation of magnetic fields would need to be considered as well (Tanaka et al. 2018; Luo et al. 2020). In this type of extension of the pulsar wind nebula scenario, several of the fundamental short comings of the toroidal field model of the Crab Nebula (Kennel & Coroniti 1984) can be overcome as argued by Luo et al. (2020).

In the following, we assumed a model where the magnetic field is dominated by small-scale turbulence. The root mean square value $B(r)$ follows a power law for its radial dependence:

$$B(r) = B_0 \left(\frac{r}{r_s}\right)^{-\alpha}, \qquad (14)$$

with $B_0$ the field strength at the distance of the termination shock $r_s = 0.13$ pc and $\alpha \ge 0$ free parameters of the model. The value $\alpha = 0$ corresponds to the case of a constant magnetic field as for example assumed by Meyer et al. (2010). The case $\alpha = -1$ corresponds to an MHD flow dominated by kinetic energy ($\sigma \ll 1$).

## 4. Parameter estimation

The complete model description of the non-thermal electron distribution (Eqs. (8) and (9)), the spatial distribution, temperature, and mass of the dust grains (Eq. (12)) as well as the strength and radial dependence of the magnetic field (Eq. (14)) requires one to choose values for the underlying parameters. The complete set of parameters were combined in a tuple $\Psi = (\Psi_1, \dots, \Psi_{21})$ (see Table 4 for the complete list of parameters and best fitting values found in Sect. 4.4.2).

The subset of parameters $\Psi_1, \dots, \Psi_{14}$ relates to the electron distribution and is closely related to the choice of the magnetic field with normalisation $\Psi_{20} = B_0$ and power-law index $\Psi_{21} = \alpha$ for the radial distribution (see Eq. (14)). The parameters $\Psi_{15}, \dots, \Psi_{19}$ characterise the dust distribution and the amount and temperature of the two dust populations (see Eq. (12)).

### 4.1. Data sets used

The method to estimate the best-fitting values $\Psi$ is based upon the minimisation of the sum of $\chi^2$ values calculated for the two compilations of observational data. $D_{\text{sync}}$ and $D_{\text{IC}}$ (see the following Sect. 4.2 on the definition of the test statistics).

The observational data set $D_{\text{sync}}$ includes 203 flux measurements $f_v$ and uncertainties $\sigma_f^3$ as well as 15 measurements of the angular extent (see Table 2) $\theta_{68}^{\text{obs}}(v)$ and corresponding uncertainties $\sigma_{\theta}(v)$. The resulting synchrotron data set

$$D_{\text{sync}} = D_{\text{sync,SED}} \cup D_{\text{sync,ext}}, \qquad (15)$$

is the union of the spectral measurements

$$D_{\text{sync,SED}} = \{(f_{v,i}, \sigma_{f,i}) | i = 1, \dots, 203\}, \qquad (16)$$

and the measurement of the angular (apparent) extension

$$D_{\text{sync,ext}} = \{(\theta_{68}(v_j), \sigma_{\theta}(v_j)) | j = 1, \dots, 15\}. \qquad (17)$$

The spectral data set $D_{\text{sync,SED}}$ comprises a variety of observations collected with different instruments. In many instances, the resulting measurements are mutually inconsistent indicating systematic uncertainties (see also Sect. 2 for more details). The quoted uncertainties on the flux measurements $\sigma_{f,i}$ are limited to the statistical uncertainties.

While below energies of 1 GeV, the dominant contribution to the Crab Nebula's luminosity is synchrotron emission of electrons and thermal emission from radiatively heated dust grains, at energies larger than 1 GeV, the dominant emission process is inverse Compton scattering of the same electron population on various seed photon fields present in the nebula.

For energies above 1 GeV, we compiled a data set $D_{\text{IC}}$

$$D_{\text{IC}} = D_{\text{IC,SED}} \cup D_{\text{IC,ext}}, \qquad (18)$$

which includes the measurements of the energy spectrum using Fermi-LAT data (25 data points: $D_{\text{IC,SED}}$) and the angular extension ($D_{\text{IC,ext}}$: 8 measurements as listed in Table 2) combining Fermi-LAT (Yeung & Horns 2019) and H.E.S.S. measurements (H. E. S. S. Collaboration 2020). The spectral measurements obtained with H.E.S.S. and other ground-based instruments were included separately in the fitting procedure (see Sect. 4.4.2).

### 4.2. Cost functions

The approach followed here used combinations of cost functions $\chi^2(D_a, \Psi)$ calculated for the data sets $D_a$, where $a$ indicates the underlying observations as given in Eqs. (15)–(18); for example, for the synchrotron SED:

$$\chi^2_{\text{sync,SED}} = \chi^2(D_{\text{sync,SED}}, \Psi) = \sum_{i=1}^{203} \left( \frac{L_{v,i}(\Psi) - 4\pi d_{\text{Crab}}^2 f_{v,i}}{4\pi d_{\text{Crab}}^2 (\sigma_{f,i}^2 + (0.05 f_{v,i})^2)} \right), \qquad (19)$$

<footnote n="3">³ See the supplementary table available under https://github.com/dieterhorns/crab_pheno</footnote>
where $L_{v,i}$ is the model synchrotron luminosity at frequency $v_i$. The statistical uncertainty on the flux measurement was adjusted by a systematic uncertainty, added in quadrature. Without further detailed knowledge about the contributions to the systematic uncertainty of individual measurements, we chose to assume that the flux measurements are subject to relative uncertainty of 5%<fnref n="4" />. The best-fitting parameters do not change within the accuracy of the minimising routine when varying the relative systematic uncertainty in the range from 0% to 7%.

The value of the model luminosity $L_{v,i}$ was calculated for a set of parameters $\Psi$ (see Table 4). Furthermore, we calculated the $\chi^2$-statistics for the angular extension:

$$ \chi^2_{\text{sync,ext.}} = \chi^2(\mathcal{D}_{\text{sync,ext}}, \mathbf{\Psi}) = \sum_{j=1}^{15} \left( \frac{\theta_{68}(v_j, \mathbf{\Psi}) - \theta_{68}(v_j)}{\sigma_{\theta}(v_j)} \right)^2, \quad (20) $$

where $\theta_{68}(v_j, \mathbf{\Psi})$ was calculated from the intensity $I_{v,j}(\theta)$, such that

$$ 0.68 = \frac{\int_{0}^{\theta_{68}(v_j)} d\theta \sin\theta I_{v,j}(\theta)}{\int_{0}^{\theta_{\max}} d\theta \sin\theta I_{v,j}(\theta)}. \quad (21) $$

The maximum value $\theta_{\max} = 250''$ was chosen such that the calculated angular extension $\theta_{68}$ does not change within the numerical accuracy for a choice of a larger value of $\theta_{\max}$. The measured extensions and uncertainties are listed in Table 2. The uncertainties on the measured angular extension were estimated from the data and do not require additional adjustments as the resulting $\chi^2_{\text{sync,ext}}$ (see Sect. 4.3) indicates an acceptable goodness of fit.

The value of

$$ \chi^2_{\text{sync}} = \chi^2_{\text{sync,SED}} + \chi^2_{\text{sync,ext}} \quad (22) $$

is the contribution to the cost function related to the synchrotron data set.

Similarly to the calculation of $\chi^2_{\text{sync}}$, the inverse Compton data set $\mathcal{D}_{\text{IC}}$ was used to calculate separately $\chi^2_{\text{IC,SED}}$ and $\chi^2_{\text{IC,ext}}$ and the sum of the two

$$ \chi^2_{\text{IC}} = \chi^2_{\text{IC,SED}} + \chi^2_{\text{IC,ext}}, \quad (23) $$

which is relevant to estimate the parameters related to the strength and radial distribution of magnetic fields in the nebula.

The main systematic uncertainty of the flux measurement with Fermi-LAT (see also Sect. 2) is related to the energy calibration of the calorimeter. This known uncertainty was taken into account by adding a contribution to $\chi^2_{\text{IC,SED}}$ which relates to a relative scaling of the energy (frequency):

$$ v \mapsto v' = v(1 + \zeta), $$

such that

$$ \chi^2_{\text{IC,SED}} = \sum_{i=1}^{25} \left( \frac{L_{v,i}(\mathbf{\Psi}) - 4\pi d_{\text{Crab}}^2 f_{v',i}}{4\pi d_{\text{Crab}}^2 \sigma_{f,i}} \right)^2 + \Theta(\zeta) \left( \frac{\zeta}{0.02} \right)^2 + \Theta(-\zeta) \left( \frac{\zeta}{-0.05} \right)^2. \quad (24) $$

The estimated systematic uncertainty allows for a relative shift in energy $\zeta$ to be a random deviate with an asymmetric normal distribution given by the quoted uncertainty on $\zeta = 0^{+0.02}_{-0.05}$

<table><caption>**Table 4.** Best-fitting parameters, and $\chi^2$-values using VHE data (MAGIC, VERITAS, Tibet Asy, LHAASO).</caption><thead><tr><th>Parameter</th><th>Best-fitting values<br>(68% c.l.)</th></tr></thead><tbody><tr><td colspan="2"><em>Radio electrons</em></td></tr><tr><td>Ψ<sub>1</sub> = s<sub>r</sub></td><td>1.54 ± 0.03</td></tr><tr><td>Ψ<sub>2</sub> = ln(N<sub>r,0</sub>)</td><td>114.7 ± 0.2</td></tr><tr><td>Ψ<sub>3</sub> = ln(γ<sub>1</sub>)</td><td>11.4 ± 0.1</td></tr><tr><td>Ψ<sub>4</sub> = ρ<sub>r</sub> ['']</td><td>89 ± 3</td></tr><tr><td colspan="2"><em>Wind electrons</em></td></tr><tr><td>Ψ<sub>5</sub> = s<sub>1</sub></td><td>3.1 ± 0.2</td></tr><tr><td>Ψ<sub>6</sub> = s<sub>2</sub></td><td>3.45 ± 0.01</td></tr><tr><td>Ψ<sub>7</sub> = s<sub>3</sub></td><td>3.77 ± 0.04</td></tr><tr><td>Ψ<sub>8</sub> = ln(γ<sub>w0</sub>)</td><td>12.7 ± 0.2</td></tr><tr><td>Ψ<sub>9</sub> = ln(γ<sub>w1</sub>)</td><td>15.6 ± 0.8</td></tr><tr><td>Ψ<sub>10</sub> = ln(γ<sub>w2</sub>)</td><td>19.2 ± 0.2</td></tr><tr><td>Ψ<sub>11</sub> = ln(γ<sub>w3</sub>)</td><td>22.3 ± 0.03</td></tr><tr><td>Ψ<sub>12</sub> = ln(N<sub>w,0</sub>)</td><td>73.8 ± 0.5</td></tr><tr><td>Ψ<sub>13</sub> = β</td><td>0.15 ± 0.01</td></tr><tr><td>Ψ<sub>14</sub> = ρ<sub>0</sub>['']</td><td>99 ± 4</td></tr><tr><td colspan="2"><em>Dust parameters</em></td></tr><tr><td>Ψ<sub>15</sub> = r<sub>out</sub> [pc]</td><td>1.53 ± 0.09</td></tr><tr><td>Ψ<sub>16</sub> = log<sub>10</sub>(M<sub>1</sub>/M<sub>☉</sub>)</td><td>-4.4 ± 0.1</td></tr><tr><td>Ψ<sub>17</sub> = log<sub>10</sub>(M<sub>2</sub>/M<sub>☉</sub>)</td><td>-1.2 ± 0.1</td></tr><tr><td>Ψ<sub>18</sub> = T<sub>1</sub> [K]</td><td>149 ± 8</td></tr><tr><td>Ψ<sub>19</sub> = T<sub>2</sub> [K]</td><td>39 ± 2</td></tr><tr><td colspan="2"><em>Magnetic field parameters</em></td></tr><tr><td>Ψ<sub>20</sub> = B<sub>0</sub> [μG]</td><td>264 ± 9</td></tr><tr><td>Ψ<sub>21</sub> = α</td><td>0.51 ± 0.03</td></tr><tr><td colspan="2"><em>Goodness of fit</em></td></tr><tr><td>χ<sup>2</sup><sub>sync,SED</sub> (d.o.f.)</td><td>182 (184)</td></tr><tr><td>χ<sup>2</sup><sub>sync,ext</sub> (d.o.f.)</td><td>16 (15)</td></tr><tr><td>χ<sup>2</sup><sub>IC,SED</sub> (d.o.f.)</td><td>22 (23)</td></tr><tr><td>χ<sup>2</sup><sub>IC,ext</sub> (d.o.f.)</td><td>24 (8)</td></tr><tr><td>χ<sup>2</sup><sub>VHE</sub> (d.o.f.)</td><td>41 (55)</td></tr><tr><td>χ<sup>2</sup><sub>tot</sub> (d.o.f.)</td><td>285 (285)</td></tr></tbody></table>

**Notes.** See Sect. 4.4.2 on details of the fit procedure. The value of $\gamma_0$ for the radio electrons has been kept fixed to be $\gamma_0 = 22$ (higher values are excluded by the low frequency radio measurements; smaller values are not constrained by observations).

(Ackermann et al. 2012). The value of $\zeta$ was found whenever evaluating $\chi^2_{\text{IC,SED}}$ (see Sect. 4.4.1).

Finally, we included the contribution of the VHE spectral measurements to the cost function with $\chi^2_{\text{VHE}}$ (see Sect. 4.4.2).

In the following sections, we follow three different minimisation approaches, where the most general fit is presented in Sect. 4.4.2. Firstly, in Sect. 4.3, we minimise the cost function $\chi^2_{\text{sync}}$ (see Eq. (22)) on a grid of values for $B_0$ and $\alpha$, demonstrating, that the synchrotron data set provides only a weak constraint on the parameter $\alpha$ and is degenerate for $B_0$. Secondly, we proceed in Sect. 4.4 with the inverse Compton component which allows us to lift the degeneracy with respect to $B_0$. In Sect. 4.4.1, the sum of $\chi^2_{\text{IC}}$ (Eq. (23)) and $\chi^2_{\text{sync}}$ is evaluated on the same grid of $B_0$ and $\alpha$ as in the previous approach. Finally, in the most general (time-consuming) fit in Sect. 4.4.2, the cost function $\chi^2_{\text{sync}} + \chi^2_{\text{IC}}$ is used for marginalising over 20 parameters ($\Psi_1...20$). The resulting model is used to calculate $\chi^2_{\text{VHE}}$ and

<footnote n="4">⁴ In the previous works of Meyer et al. (2010), the relative systematic uncertainty was assumed to be 7%.</footnote>
consequently $\chi^2_{\text{tot}}(\alpha) = \chi^2_{\text{sync}} + \chi^2_{\text{IC}} + \chi^2_{\text{VHE}}$ is used to find the best estimate for $\alpha$.

### 4.3. Fitting the synchrotron spectrum and extension

Since we are mainly interested in the best-fitting values of $\Psi_{20} = B_0$ and $\Psi_{21} = \alpha$, we marginalised over $\Psi_{1...19}$ such that we used the parameters $\Psi_{1...19}$ which minimise $\chi^2_{\text{sync}}$ for a given pair of $\alpha, B_0$. The resulting

$$ \chi^2_{\text{sync}}(\alpha, B_0) = \chi^2(\mathcal{D}_{\text{sync}}, \alpha, B_0, \Psi_{1...19}), \quad (25) $$

depends on $\alpha$ and $B_0$ only.

As a demonstration of the approach, we show in Fig. 5a the increment $\Delta\chi^2_{\text{sync}}(\alpha, B_0)$ relative to the value of the global minimum $\chi^2_{\text{sync}}(\tilde{\alpha}, \tilde{B}_0)$

$$ \Delta\chi^2(\alpha, B_0) = \chi^2_{\text{sync}}(\alpha, B_0) - \chi^2_{\text{sync}}(\tilde{\alpha}, \tilde{B}_0). \quad (26) $$

For a constant value of $\alpha$, the resulting $\chi^2_{\text{sync}}(\alpha, B_0)$ does not change substantially over the range of $B_0 \in [50 \, \mu\text{G}, 1025 \, \mu\text{G}]$ used in the scan with a step-size of $\Delta B_0 = 25 \, \mu\text{G}$.

We find a global minimum for $\chi^2_{\text{sync}}$ at $\tilde{\alpha} = 0.9$ with an uncertainty of $\sigma_{\alpha}(95\%\text{c.l.}) \approx 0.8$ (see bottom left panel, Fig. 5c).

The parameter $B_0$ is the normalisation of the magnetic field at the position of the shock $r_s$. The resulting field $\langle B \rangle$ averaged over the volume from $r_s$ to $r_{\text{PWN}} = 3 \, \text{pc}$ was calculated to be

$$ \langle B \rangle = \frac{1}{V} \int dV B(r) = \frac{3 B_0 r_s^\alpha r_{\text{PWN}}^{3-\alpha} - r_s^{3-\alpha}}{3-\alpha r_{\text{PWN}}^3 - r_s^3}. \quad (27) $$

In Fig. 5a we plot contours for $\langle B \rangle / \mu\text{G} \in \{30, 60, 90\}$ for orientation.

For a fixed value of $\alpha$, $\chi^2_{\text{sync}}$ is degenerate for any choice of the normalisation $B_0$. However, the energy in magnetic field

$$ W_B = \int dV \frac{B(r)^2}{8\pi} \quad (28) $$

increases for increasing $B_0$ while the normalisation of the electron spectrum decreases and correspondingly, the energy in particles $W_e$

$$ W_e = \int dV \int d\gamma m_e c^2 \gamma n_{\text{el}}(r, \gamma) \quad (29) $$

decreases monotonically with increasing $B_0$.

The minimum energy configuration for the total energy

$$ W_{\text{tot}} = W_e + W_B \quad (30) $$

is reached for $B_0 = 514 \, \mu\text{G}$, which is slightly lower than the equipartition ($W_e = W_B$) field $B_{\text{eq}} = 558 \, \mu\text{G}$.

The best-fit result which minimises $\chi^2_{\text{sync}} = 195.2(197 \, \text{d.o.f.})$ is reached for $\tilde{\alpha} = 0.9$ and large uncertainties which go well beyond the range of parameters scanned here.

The fitting of the model to the synchrotron data $\mathcal{D}_{\text{sync}}$ provides an overall satisfying goodness of fit. The parameters related to the electron distribution and dust properties can be estimated accurately.

We identified the combinations of $B_0$ and $\alpha$ which lead to a minimum energy configuration for each scanned value of $\alpha$. The resulting contour $B_0 = B_{\text{min}}$ is indicated as a cyan dashed line in Fig. 5.

Even though the configuration of the system is expected to be close to the minimum or equipartition value of $B_0$, without additional information the normalisation $B_0$ is degenerate with the number density of the electrons and $\alpha$ is only weakly constrained. The latter parameter is mostly degenerate with the radial distribution of the electron density ($\Psi_{13} = \beta$).

Additional information such as the inverse Compton emission needs to be incorporated in the fit to break these degeneracies and determine the parameters $B_0$ and $\alpha$.

### 4.4. Fitting the inverse Compton spectrum and extension

The inverse Compton emissivity of the electrons is proportional to the product of number density of electrons and the energy density of the seed photon field (in the Thomson limit of inverse Compton scattering).

The observation of inverse Compton emission and its spatial extension breaks the degeneracy of the model parameters. Different to a fit to the synchrotron component alone (see previous section), no additional constraints as for instance equipartition or minimum energy are required to determine $\alpha$ and $B_0$.

#### 4.4.1. HE SED and extension

In order to make the fitting numerically tractable and to get an overview on the best-fitting parameters, we carried out a minimisation following these steps$^5$:

1. Choose $\alpha \in [0, 1.2]$, starting at 0, incremental by 0.05

2. For each value of $\alpha$, marginalise over $B_0$: $\chi^2(\alpha) = \chi^2(\alpha, \tilde{B}_0)$

3. Repeat steps (1) and (2) and locate the global minimum $\chi^2(\tilde{\alpha})$.

The minimisation carried out in step (2) relates to $\chi^2 = \chi^2_{\text{sync}} + \chi^2_{\text{IC}}$ with $\chi^2_{\text{IC}} = \chi^2_{\text{IC,SED}} + \chi^2_{\text{IC,ext}}$ given by

$$ \chi^2_{\text{IC,SED}} = \chi^2(\mathcal{D}_{\text{IC,SED}}, \alpha, B_0, \Psi_{1...19}, \zeta), \quad (31) $$

and

$$ \chi^2_{\text{IC,ext}} = \chi^2(\mathcal{D}_{\text{IC,ext}}, \alpha, B_0, \Psi_{1...19}). \quad (32) $$

Note, the values of the parameter estimates related to the electron and dust population $\tilde{\Psi}_{1...19}$ were found by minimising $\chi^2_{\text{sync}}$ separately for each choice of parameters $\alpha, B_0$. This leads to an acceptable computation time for a parameter scan.

The resulting $\chi^2_{\text{IC}}(\alpha)$ shows a minimum for $\tilde{\alpha} = 0.8$ shown in Fig. 5d (bottom right panel). While $\chi^2_{\text{sync}}$ remains essentially constant over the entire scan range, $\chi^2_{\text{IC}}$ rises quickly towards the lower and upper boundary of the scan range. In order to accommodate the dynamical range of $\Delta\chi^2_{\text{IC}}$, the right hand panels of Fig. 5 show $\log_{10}(\Delta\chi^2_{\text{IC}} + 1)$.

The best estimates for $\tilde{\alpha} = 0.8^{+0.16}_{-0.13}$ and $\tilde{B}_0 = (497 \pm 43) \, \mu\text{G}$ are close to the minimum energy condition ($B_{0,\text{min}} = 514 \, \mu\text{G}$, see also cyan lines in Figs. 5c,d and 6).

For an overview of the best-fitting value $\tilde{B}_0$, we also carried out a scan in $\alpha, B_0$ shown in Fig. 5b. The solid yellow contour in the same Figure indicates the 95% confidence level parameter region favoured here. The shape of the favoured parameter region follows the shape of the cyan line which traces the minimum energy configuration.

The fit to the synchrotron component is in comparison to the previous fit to the synchrotron component only slightly worse with $\chi^2_{\text{sync}} = 196.4$ instead of $\chi^2_{\text{sync}} = 195.2$ found in the previous

<footnote>$^5$ In the subsequent section, a more general approach is chosen.</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Results of the fitting procedure. Left column: (a) $\Delta\chi^2_{\text{sync}}(\alpha, B_0)$ and (c) $\Delta\chi^2_{\text{sync}}(\alpha)$. Right column: (b) $\Delta\log_{10}(\chi^2_{\text{IC}}(\alpha, B_0) + 1)$ and (d) $\log_{10}(\Delta\chi^2_{\text{IC}}(\alpha) + 1)$. In the top panels, for each combination of $\alpha, B_0$, the other parameters of the model are marginalised; in the bottom panels, this includes $B_0$. The dotted white lines in the top panels indicate the contours for constant values of volume averaged $\langle B \rangle = 30, 60, 90$ µG. The cyan dashed line indicates the combinations $B_0$ and $\alpha$ with the minimum total energy $W_{\text{tot}}$ in electrons and magnetic field. The yellow solid contour marks the 95% confidence level error box for $B_0, \alpha$. The grey boxes in the bottom panel indicate the 95% confidence interval for $\alpha$. In the bottom right diagram, the cyan dashed line indicates the minimal total energy $W_{\text{tot}}$ for a given $\alpha$. For a discussion on the relevance of the co-location of the global minimum $\chi^2_{\text{IC}}$ with the corresponding global minimum of $W_{\text{tot}}$, see the text.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Volume-integrated total energy $W_{\text{tot}}$ equals the sum of the energy in relativistic electrons $W_e$ and magnetic field $W_B$. For any best-fitting set of parameters, the energy partition is a function of the normalisation $B_0$ of the magnetic field at the shock distance $r_s$. The best-fitting value is found to be $\tilde{B}_0 = (497 \pm 43)$ µG (90% c.l. uncertainty, including systematic uncertainties). This value is slightly smaller than the minimal energy configuration at $B_0 = 514$ µG.</figcaption></figure>

section. The inverse Compton component is fit with $\chi^2_{\text{IC,SED}} = 32.0$ and $\chi^2_{\text{IC,ext}} = 17.2$. The total $\chi^2(\text{d.o.f.}) = 245.6(228)$ indicates still an acceptable fit with $p(\chi^2 > 245.6, \text{d.o.f.} = 228) =$

0.2. However, when considering the fit to the inverse Compton component only, the resulting $\chi^2_{\text{IC}}(\text{d.o.f.}) = 49.2(29)$ indicates that the fit is not acceptable with $p(\chi^2 > 49.2, \text{d.o.f.} = 29) = 0.01$.

Note, the best fit value for $\tilde{\alpha}$ is consistent for the two different estimates from fitting the synchrotron and inverse Compton data.

### 4.4.2. VHE SED

The fitting procedure outlined above has been primarily optimised for the synchrotron emission of the Crab Nebula (SED, extension). The inverse Compton component of the SED was mainly used to constrain independently from the electron spectral shape the magnetic field $B_0$ and the related parameter $\alpha$.

Both, the synchrotron as well as the inverse Compton SED at energies below the peak in the SED are fairly degenerate for different choices of the radial dependence of the magnetic field (see also the following Sect. 4.4.3). The resulting uncertainty of the power-law index $\alpha$ used for $B(r) \propto r^{-\alpha}$ is comparably large with $\sigma_{\alpha} \approx 0.15$, even after including the Fermi-LAT spectral information. Furthermore, when constraining the electron distribution by matching the synchrotron emission to the SED and determining in a separate step the best-fitting values of $\alpha$ and $B_0$ (see previous section), the resulting goodness of fit to the Fermi-LAT data is not acceptable.
<table><caption>**Table 5.** Energy-scaling factors $\zeta_i$ and resulting $\chi^2$-values before and after changing the energy-scale.</caption><thead><tr><th>Data set<br>(α̃ = 0.51)</th><th>Scaling factor ζ<sub>i</sub></th><th>Statistical error σ<sub>ζ</sub><br>68% c.l.</th><th>χ<sup>2</sup><sub>before</sub>(d.o.f.)</th><th>χ<sup>2</sup><sub>after</sub>(d.o.f.)</th></tr></thead><tbody><tr><td>Fermi-LAT</td><td>1.00</td><td>Fixed</td><td>21.7 (25)</td><td>21.7 (25)</td></tr><tr><td>MAGIC</td><td>0.891</td><td>0.006</td><td>221.0 (14)</td><td>11.7 (13)</td></tr><tr><td>VERITAS</td><td>1.050</td><td>0.003</td><td>469.9 (12)</td><td>9.1 (11)</td></tr><tr><td>MAGIC (VLZA)</td><td>1.141</td><td>0.025</td><td>55.6 (7)</td><td>2.6 (6)</td></tr><tr><td>Tibet ASγ</td><td>0.918</td><td>0.017</td><td>22.8 (10)</td><td>6.4 (9)</td></tr><tr><td>LHAASO</td><td>0.964</td><td>0.006</td><td>41.0 (17)</td><td>11.6 (16)</td></tr><tr><td>Combined</td><td></td><td></td><td>832.0 (85)</td><td>63.1 (80)</td></tr></tbody></table>

**Notes.** Energy-scaling factors $\zeta_i$ for the VHE spectra are derived for the best-fitting model with $\tilde{\alpha} = 0.51$ (the corresponding best fit parameters $\Psi$ are given in Table 4).

The accurate measurement of the IC component at energies covered with ground based VHE gamma-ray telescopes bears the potential of determining $\alpha$ with greater precision. This is evident from the wide range of flux values and slopes of gamma-ray spectra at energies beyond several 100 GeV predicted for the interval of $\alpha$ considered here (see Sect. 4.4.3). Therefore, the VHE measurement of the SED reduces the uncertainty on the value of $\alpha$.

We modified the method used in the previous section by including the LAT spectral information to estimate the parameters $\Psi_{1...20}$. This numerically expensive approach was necessary to improve the goodness-of-fit of the inverse Compton method and to include the VHE measurement.

The fit is carried out in the following steps:

1. Choose $\alpha \in [0.4, 0.8]$, increment by 0.005

2. Minimise $\chi^2(\alpha) = \chi^2_{\text{sync}} + \chi^2_{\text{IC}}$ to estimate $\tilde{\Psi}_{1...20}$

3. Calculate for each $\alpha$ the predicted VHE spectrum and the resulting $\chi^2_{\text{VHE}}$.

4. Repeat until the minimum $\chi^2(\tilde{\alpha})$ is found to estimate the parameter $\alpha$.

The third step involves both individual VHE data sets $D_{\text{VHE},i}$ as well as their combinations.

$$D_{\text{VHE}} = \bigcup_{i} D_{\text{VHE},i}. \qquad (33)$$

A complete list of the data sets and individual fit results are presented in Appendix A and Table A.1. In the following, we considered a combination of five data sets obtained with MAGIC (normal zenith angles as well as very large zenith angles), VERITAS, Tibet ASγ, and LHAASO. With this combination of observations, a consistent model was found. Another consistent solution was found for the combination of the data from the HEGRA and HAWC observations (see Appendix C), while the spectra obtained with the H.E.S.S. I and II telescope arrays cannot be fit well by the model presented here (see Table A.1 for the individual fit results).

For each data set $D_{\text{VHE},i}$, we introduced a relative energy scaling $\zeta_i$ as an additional nuisance parameter. This adjustment of the relative energy scale compensates for systematic uncertainties related to the absolute energy calibration of ground based air Cherenkov telescopes.

The scaling factors $\zeta_i$ determined for the VHE data are listed in Tab. 5. The values are close to unity with the largest correction applied to the MAGIC data set taken at very large zenith angles (MAGIC VLZA: $\zeta = 1.14 \pm 0.03$). The value found here is well within the error budget of the systematic uncertainty related to

the energy scale which is estimated to be 17% for this data set (MAGIC Collaboration 2020). Note, the dramatic reduction of the $\chi^2$ between $\chi^2_{\text{before}} = 832$ before the scaling and $\chi^2_{\text{after}} = 63.1$ after the scaling (see Tab. 5 for more details). The introduction of the five additional parameters reduces the number of degrees of freedom from 85 to 80 for the VHE data sets used here.

The best-fitting parameters were used to calculate the broadband SED for $\tilde{\alpha} = 0.51 \pm 0.03$ and $\tilde{B}_0 = (264 \pm 9) \mu\text{G}$ (for a complete list of parameters used, see Table 4). The result is shown in Fig. 7. The overall fit from radio to PeV energies traces the observational data with a resulting $\chi^2_{\text{SED}} = \chi^2_{\text{sync,SED}} + \chi^2_{\text{IC,SED}} + \chi^2_{\text{VHE}} = 182+24+41 = 247$ with $184+23+55 = 262$ degrees of freedom (see Table 4 for the resulting parameter estimates as well as the goodness of fit summary). In the residual shown in the bottom panel of Fig. 7, deviations are visible in the MeV spectrum as well as at PeV energies. We return to possible explanations in Sect. 4.4.4.

In more detail, the SED at HE and VHE gamma-rays is shown in Fig. 8. The VHE data sets selected are consistent after the relative energy scales have been slightly adjusted to a common energy-scale provided by the LAT data. The five different VHE data sets overlap in the energy range covered both with the LAT data as well as among each other. The resulting goodness of fit is quantified by $\chi^2_{\text{VHE}} = 41$ for 55 degrees of freedom. Even without additional systematic uncertainties related to the flux measurement, the statistical uncertainties are not sufficiently small to further constrain or even rule out the model used here.

The spatial extension of the model emission is matching the measured extension quite well. This is shown in Fig. 9. The resulting $\chi^2_{\text{ext}} = \chi^2_{\text{sync,ext}} + \chi^2_{\text{IC,ext}} = 16 + 24 = 40$ for $15 + 8 = 23$ data points (the number of degrees of freedom is reduced by the number of parameters of $\ge 2$) indicates at first glance a poor fit: $p(\chi^2 > 40, \text{d.o.f.} = 21) \approx 1.2 \times 10^{-5}$. However, the large value of $\chi^2_{\text{IC,ext}}$ is mainly driven by one (outlying) measurement at $E = 14$ GeV. When removing that individual point, the resulting $\chi^2_{\text{ext}}(\text{d.o.f.}) = 21(20)$.

After including $D_{\text{VHE}}$, the constraint on the magnetic field distribution improved considerably. The best-fitting estimate for $\tilde{\alpha} = 0.51 \pm 0.03$ has a remarkably small statistical uncertainty. The value of $\alpha$ is mainly constrained by the accurate measurement of the SED at VHE energies. The measurement of the spatial extension is less constraining (see also Sect. 4.4.3). The difference in the best-fitting values for $\alpha$ when comparing the results obtained in Sect. 4.4.1 and the ones found here is not significant, but will be discussed in the context of an extension of the model in Sect. 5.5.
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Best-fitting model in comparison to the data. The upper plot shows the synchrotron and inverse Compton (IC) radiation by radio (R) and wind (W) electrons produced in the spatially varying magnetic field and seed photon field (blue and magenta coloured broken lines indicate the contributions from IC scattering on individual seed photon fields). The blue data points are used in the fit to determine the parameters $\Psi_{1...20}$ listed in Table 4. The VHE data have been scaled using the values of $\zeta_i$ determined for each data-set (see text for further details) to adjust the energy scale for prediction given by the global minimum $\chi^2(\hat{\alpha} = 0.51)$. In the lower plot, the relative residuals are shown (statistical uncertainties only).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Comparison of the best-fitting model with the HE and VHE data. The detailed view confirms the consistency of the selected data sets (after scaling) and the good agreement of the model to the data. The selected VHE data sets overlap with the LAT data and each data set overlaps with at least one other VHE data set.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Extension of the model emission is characterised by the radial size $\theta_{68}$. The model curve is the calculated extension which matches well the observational data.</figcaption></figure>

### 4.4.3. Varying $\alpha$: Effect on the IC SED and extension

As discussed in the previous sections, there is a unavoidable degeneracy of $\chi^2_{sync}$ for variations of parameters related to the magnetic field ($\alpha, B_0$) and the electron spectrum.

As an illustration of this degeneracy with respect to the choice of $\alpha$, we show in Fig. 10 the best-fitting SEDs, colour-coded for different values of $\alpha$. Each SED is the result of the fitting procedure outlined in Sect. 4.4.1 above. For clarity, the data points are not included.

The resulting SEDs indicate that the synchrotron component does not change in a noticeable way when varying the value of $\alpha$ while the inverse Compton emission varies in a systematic way: For a constant magnetic field ($\alpha = 0$), the energy spec-

trum is harder than for models where the magnetic field drops with increasing radius ($\alpha > 0$).

The model with the largest value of $\alpha = 1.2$ covered in the scan has the softest spectrum. The differences between the models with varying $\alpha$ is most pronounced at the energies beyond the peak energy of approximately 100 GeV.

Similar to the SED, the angular extension of the synchrotron nebula is highly degenerate with respect to variations of the parameter $\alpha$. This is demonstrated in Fig. 11.

Different to the SED which shows the largest variations of the SED for the VHE part when changing the value of $\alpha$, the angular extension varies the strongest for the part of the inverse Compton emission covered with HE gamma-ray observations.

At several GeV of gamma-ray energy, for a constant magnetic field, the extension is the lowest (about 100′′) and it
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Overall spectral energy distribution for the estimated best-fitting parameters $\tilde{\Psi}_{1\dots20}$ (after minimising $\chi^2_{\text{IC}}$ for fixed values of $\alpha$) is shown. The value of $\alpha$ characterises the radial decrease of the magnetic field $B \propto r^{-\alpha}$ (see Eq. (14)). While the synchrotron part of the spectral energy distribution is almost degenerate with respect to the choice of $\alpha$, the inverse Compton dominated part (right panel) shows noticeable differences up to the peak at approximately $10^{11}$ eV. Note, in the left panel the sum of synchrotron and inverse Compton emission is shown, while in the right panel, the SED of the inverse Compton emission is additionally superimposed.</figcaption></figure>

systematically increases with increasing $\alpha$, such that for $\alpha = 1.2$, the extension reaches a value of about 150$''$. At the highest energy part, the variation of the extension for different values of $\alpha$ is effectively constant.

The underlying reason for the observed changes in the SED is related to the magnetic field relevant for electrons of different energies. The most energetic electrons are located close to the shock which features the largest magnetic field. Consequently, a reduced number density of these electrons is required to match the synchrotron emission which in turn leads to a reduced inverse Compton emission, softening the VHE spectrum with increasing value of $\alpha$. The energy range covered with Fermi-LAT is less affected by changes of $\alpha$ as the underlying electrons populate the entire volume and therefore produce synchrotron emission in an averaged magnetic field which does not change for the best-fitting model when varying $\alpha$.

#### 4.4.4. An additional ultra-high energy component

The energy spectrum measured at ultra-high energy (UHE) gamma-rays above several 100 TeV with the LHAASO array (LHAASO Collaboration 2021) apparently deviates from the predicted shape of the inverse Compton emission at the highest energy bin centred on $E_\gamma = 1.26$ PeV$^6$. It has been speculated that the PeV emission could be due to an additional hadronic component, similar to the one possibly present in the Vela-X region (Horns et al. 2006). Previous upper limits on the presence of a hadronic wind have constrained the fraction of spin-down power diverted to the acceleration of relativistic hadrons to be less than ≈20% up to PeV energies (Aharonian et al. 2004). The PeV excess could be produced by as little as 1% of the spin down power assuming efficient confinement of the protons in the Crab Nebula (LHAASO Collaboration 2021).

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Angular extension from the fitting procedure described in the text. After marginalising $\tilde{\Psi}_{1\dots20}$, the predicted angular extension of the synchrotron nebula matches the data well for values of $\alpha = 0, \dots, 1.1$ with very small variations visible in the transition between thermal and non-thermal emission at near- to far infrared. At higher energies, the predicted size differs mostly at the low-energy end of the inverse Compton emission.</figcaption></figure>

Here, we invoked a second electron distribution similar to previous attempts of explaining the flaring component of the Crab Nebula (see e.g., Lobanov et al. 2011) with the energetic electrons located in the vicinity of the pulsar or close to the shock. The finding of a fast (several days) dimming in the Nebula's emission (Yeung & Horns 2020) indicates that most of the highest energy synchrotron emission (>75%) is produced in a compact region.

Motivated by these findings and the particular shape of the spectrum measured from the Crab Nebula with the COMPTEL instrument, we introduced a secondary PeV electron component which (on average) matches the highest energy part of the synchrotron spectrum. In the framework of the model developed here, we reduced the maximum energy of the wind electrons from $\gamma_{w3} = 4.8 \times 10^9$ (E$_{w3}$ = 2.5 PeV) to $\gamma_{w3} = 2.5 \times 10^9$ (E$_{w3}$ = 1.3 PeV). This way, the presumably steady synchrotron emission cuts off below 10 MeV (see Fig. 12).

<footnote n="6">⁶ When considering the actual photon statistics, the Poissonian probability to observe 2 or more photons in the energy bins around 1 PeV is predicted to be 15.6% when calculating the expected number of events with the best fitting model and the collection area given by LHAASO Collaboration (2021), see Appendix D for more details.</footnote>
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** A suggested model for the SED including a secondary electron contribution located at the shock position. The low and high state spectra as measured with Fermi-LAT are shown for comparison (Yeung & Horns 2020). The secondary electrons' spectrum during average (high) state is a power-law with $dn/d\gamma \propto \gamma^{-2.3}$ between $\gamma_1 = 3 \times 10^9$ and $\gamma_2 = 7(18) \times 10^9$. The magnetic field in the region is chosen to be 85 µG. The total energy content of the electrons is $2(3) \times 10^{42}$ ergs. The dotted lines indicate the effect of absorption on the Galactic foreground emission (see also Appendix B).</figcaption></figure>

We added a secondary component which is located at a fixed distance where due to radiative cooling, the energy is quickly dissipated. In the model considered here, the magnetic field varies as $B(r) = 264 \mu G(r/r_s)^{-0.51}$. In order to fit the synchrotron and UHE part, the consistent combination of seed photon field and magnetic field strength would favour a position in the nebula at $r \approx 20r_s$. Since it seems unlikely that particle acceleration to PeV energies is ongoing in the nebula's periphery, we explored here the possibility of electron acceleration close to the shock at $r \approx r_s$, possibly in a region with weaker magnetic field strength (here, we chose a value of 85 µG). We considered the seed photon field at the shock self-consistently calculated in our model. We did not attempt a quantitative fit here but instead chose representative values in order to match the SED. The average additional electron spectrum was assumed to follow a power law with $dn/d\gamma \propto \gamma^{-2.3}$ between $\gamma_1 = 3 \times 10^9$ and $\gamma_2 = 7 \times 10^9$.

The resulting spectra for the three flux states: wind-only with a reduced cut-off energy, average, and flaring PeV component are shown in Fig. 12 in comparison to the same data as shown before.

In addition, we include the low and high-state energy spectra provided by Yeung & Horns (2020) to compare with the respective flux state from the model. The total energy required in the secondary electron distribution is moderate with $W_e = 2 \times 10^{42}$ ergs.

The computed SED matches the data quite well. At the MeV energy range, the COMPTEL measurements are described better by the superposition of the steady wind electron emission and the average secondary component. During times of lower activity,

the reduced cut-off in the wind electron spectrum matches the observed low-state spectrum Yeung & Horns (2020; dark green squares in Fig. 12).

The inverse Compton emission at energies beyond ≈400 TeV is dominated by the emission of the secondary component. The emission of the secondary component is produced by up-scattering of synchrotron (radio) seed photons and to a lesser extent by CMB photons.

We also included a model for the high flux state (light green curve in Fig. 12). The minimal modifications were an increase of maximum energy such that $\gamma_2 = 1.8 \times 10^{10}$ and an increase of the energy of the electrons to $3 \times 10^{42}$ ergs. The resulting synchrotron emission matches the observation during the high flux state with Fermi-LAT quite well (light green circles). The resulting inverse Compton component during the high flux state is also shown in Fig. 12. In both cases, the absorption of gamma-rays due to the anisotropic Galactic radiation field is included (see Appendix B for details about the calculation carried out here).

# 5. Discussion of the parameters

## 5.1. Radio electrons

The four parameters $\Psi_{1...4}$ related to the radio electrons are well constrained by the observational data. The total number of radio emitting electrons has been found to be $N_r \approx 10^{51}$ up to a maximum Lorentz factor of $\gamma_1 \approx 90000$ for $\gamma_0 = 22$. The lower bound $\gamma_0 = 22$ has been fixed to a value which is not in conflict with the low frequency radio measurements starting at 111 MHz in our sample. A possible smaller value of $\gamma_0$ would increase the
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Energy density of the seed photons as a function of the distance from the pulsar. The dust emissivity was modelled to be confined in a shell surrounding the optical nebula. As a comparison, a composite of X-ray (Chandra: blue, white)/optical(HST: purple)/ infrared (Spitzer: pink) composite image from the nebula is shown with approximately matching spatial scale. Background image credit: X-ray: NASA/CXC/SAO; Optical: NASA/STScI; Infrared: NASA-JPL-Caltech.</figcaption></figure>

total number of radio emitting electrons. The value found here for $N_r$ should be considered a lower limit.

The resulting total energy in radio electrons is $3.4 \times 10^{47}$ ergs (this value is determined by the value of $\gamma_1$). The spatial distribution of radio electrons is found to extend up to $\rho_r = 89''$ for the best-fitting model of the magnetic field $B(r) \propto r^{-0.51}$. The resulting extension is therefore constant for the radio synchrotron emission which is in good agreement with the radial extension measured to be unchanged between 5 GHz and 150 GHz (see also Table 2).

The upper bound of the radio electron spectrum $\gamma_1$ is anti-correlated with the amount of cold dust ($\Psi_{17}$). The comparably low value for $\gamma_1$ found here is favoured by the observational data between $\lambda = 100$ µm and 1000 µm (see Fig. 1) where the dust emission starts to dominate over the synchrotron emission.

## 5.2. Dust emission

As a consequence of the comparably low value of the upper end of the radio electrons' spectrum ($\gamma_1$), the amount of cold dust required here is larger than in previous estimates which were based upon a power-law extrapolation between the radio and optical. Using the same dust emissivities as De Looze et al. (2019), we therefore find for the cold dust mass a larger value ($\Psi_{17} = \log_{10}(M/M_\odot) = -1.2 \pm 0.1$) than previously ($\approx -1.7$). In a standard (explosive) nucleosynthesis model (Woosley & Weaver 1995), the progenitor star would have to be sufficiently massive ($M_{prog} \ge 12(18) M_\odot$) to produce the dust with 75(25)% condensation efficiency. The spatial distribution of the dust is confined to a shell with fixed inner radius 0.55 pc (see Sect. 3.3) and an

outer radius estimated to be $\Psi_{15} = 1.53 \pm 0.09$ pc. This corresponds to an angle of $\approx 150''$. This parameter is constrained by the angular size measured in the infra-red (see Table 2 and Fig. 9). The dust emissivity and resulting spatial distribution of dust generated seed photons (see Eq. (12)) was in turn used for the calculation of the IC emissivity (see Eq. (10)). The IC emission of wind electrons up-scattering the dust emission as seed photons dominates at the peak energy of a few 100 GeV. This is distinctly different from previous models where the dust seed photon field is assumed to maintain a constant energy density of 0.5...1 eV cm$^{-3}$ throughout the nebula. Since the dust is formed in a shell surrounding the pulsar wind nebula, the energy density of the dust seed photon field reaches a maximum of $\approx 1.5$ eV cm$^{-3}$ at a distance of $\approx 1$ pc (see Fig. 13).

## 5.3. Wind electrons

The spectral and spatial model of the wind electrons includes ten parameters ($\Psi_{5...14}$). Most of the electrons' energy is carried by the wind electrons at the lower boundary $\gamma_{w0} = 320 000$ (corresponding to $E_{w0} \approx 200$ GeV). The total wind electrons' energy in the nebula is very close to the energy in radio electrons ($3.4 \times 10^{47}$ ergs).

There is a considerable gap in energy between the two electron components ($\gamma_1 m_e c^2 \approx 50$ GeV vs. $\gamma_{w0} m_e c^2 \approx 200$ GeV). The de-reddening of the optical nebula emission has revealed the presence of a further break in the spectrum between the optical and X-ray part (spectral index changes from 0.7 to 1.1). The break in the X-ray spectrum at about 120 keV requires a second break in the electron spectrum, such that we fitted as particle indices $3.1(2) \rightarrow 3.45(1) \rightarrow 3.77(4)$. The breaks are located at energies of $\gamma_{w1} m_e c^2 \approx 3$ TeV and $\gamma_{w2} m_e c^2 \approx 110$ TeV. The upper bound of the energy of wind electrons is fixed by the cut-off in the average synchrotron spectrum at $\gamma_{w3} m_e c^2 \approx 2.5$ PeV. The power-law indices and the positions of the breaks are most likely intrinsic to the injected particle spectrum and the underlying acceleration mechanism. Similar results related on the position of the lower-energy break have been found by, for example Martín et al. (2012). Radiative cooling would not be sufficient to explain however the appearance of the higher energy break at 110 TeV.

The spatial distribution of the wind electrons was assumed to be a radial Gaussian with an extension that follows a power law $\rho(\gamma) = \rho_0 \gamma^{-\beta}$ with index $\beta = 0.15$ (see Sect. 3.1.1). The pre-factor $\rho_0$ for the wind electrons' distribution was found to be $(99 \pm 4)''$ which is larger than $\rho_r = (89 \pm 3)''$ characteristic for the radio electrons. This fit result is a consequence of the comparably large size ($\approx 145''$) found at near infrared in the WISE images (see Table 2). The large value of $\rho_0$ is in line with the extension measured at high energies with Fermi-LAT ($\approx 120''$).

The observation that the lowest energy wind electrons with energies of approximately 200 GeV are distributed within a larger volume than the radio electrons indicates that the electrons are not only moving convectively with the wind.

## 5.4. PeV electrons

The most energetic electrons at PeV energies are traced by their synchrotron emission at energies of about 100 MeV observed with Fermi-LAT. The light curve of the 80...800 MeV synchrotron flux has been found to include at least three different spectral states (high, medium, and low state) that can be identified from Fermi-LAT observations (Yeung & Horns 2020): The high state is evident from the observation of flaring activity which occur roughly once per year (Mayer et al. 2013;
Yeung & Horns 2020). In a flaring state, the flux increases during several days dramatically with a spectral hardening and a peak in the SED reaching up to several 100 MeV. The average high flux state spectrum as defined by Yeung & Horns (2020) is shown in Fig. 12. Despite large efforts (Weisskopf et al. 2013), no obvious counter-part to the flares has been identified, even though the optical variability (position, size, and brightness) of the inner knot does show some correlation with the changes of the gamma-ray flux (Rudy et al. 2015). This correlation is consistent with the prediction of Komissarov & Lyubarsky (2003) that the majority of the MeV emission, including the flaring activity, is produced in the inner knot, where supposedly most of the PeV electrons are accelerated.

The medium state has the highest duty cycle (>95%) and is therefore identical to the average flux state shown in the previous sections. It is characterised by a very soft energy spectrum which smoothly connects to the inverse Compton emission at approximately 1 GeV.

The low state has just recently been found to be equally short lived and as rare as the flaring state. The low state is characterised by a hard spectrum that follows the extrapolation of the IC spectrum down to an energy of about 200 MeV where the spectrum seems to become soft again (Yeung & Horns 2020) and see Fig. 12.

Motivated by the apparent excess emission at the UHE end of the spectrum, we explore a possible connection of the UHE gamma-ray excess found with LHAASO and the flux states found in the MeV emission. We therefore introduced (minimal) modifications to the previous best-fitting model by (a) reducing the maximum energy of the wind electrons from 2.5 PeV to 1.3 PeV and (b) adding a PeV electron component that would radiate on average a synchrotron spectrum between 1 MeV and about 800 MeV. We assumed the PeV electrons to scatter off the seed photon field close to the termination shock. With a magnetic field of 85 µG the resulting number of electrons radiate sufficient inverse Compton emission to match highest energy emission observed. The resulting average state is shown in Fig. 12. The model matches better the data specifically in the Comptel energy range between 1 and 20 MeV than in the previous model. The IC emission at PeV energies (after absorption) matches the LHAASO data points well, even though the predicted flux is at the low side. The IC emission is dominated by the scattering of seed photons produced by the radio electrons and to a lesser degree by CMB emission. This is a consequence of the distribution of seed photons in the nebula. Close to the termination shock, the radio synchrotron emission dominates<fnref n="7" />.

In this scenario, the low, medium, and high flux states can be successfully described by the dynamics of a single emission region close to the termination shock. The resulting variability at PeV energies when comparing the average and high state is modest (about a factor of 2) and not resolvable with the current influence sensitivity of for example LHAASO.

At first glance, it seems like a plausible scenario to include an additional PeV electron population that simultaneously explains the hard synchrotron spectrum between 1 and 20 MeV and the PeV excess found with LHAASO<fnref n="8" />. Similar scenarios have been suggested in several studies (see e.g., Khangulyan et al. 2020; LHAASO Collaboration 2021; Liu & Wang 2021; Peng et al. 2022).

However, the common origin of the high, average, and low state emission as suggested by Yeung & Horns (2020) leads to considerable tension in the magnetic field strength required for confinement and the value estimated above. The size of the variable emission region is limited to the light crossing horizon for the variability time-scale $t_{var} \approx 100$ ksec (the minimum timescales for the apparent disappearing of the nebula and for the flares are similar). This translates into a minimum magnetic field to confine the electrons radiating at $E_{sync}$ (Lobanov et al. 2011):

$$B > 500 \mu G \left( \frac{E_{sync}}{100 \text{ MeV}} \right)^{-1/3} \left( \frac{t_{var}}{100 \text{ ksec}} \right)^{-2/3} . \quad (34)$$

With a magnetic field sufficiently large to confine the energetic electrons, the expected inverse Compton emission would be too small to explain the PeV emission from the Crab Nebula. This makes it unlikely that the variable MeV emission is directly linked to the PeV emission. Other, more complicated scenarios would require additional particle populations (either leptonic or hadronic; Liu & Wang 2021; Peng et al. 2022). However, with the currently limited statistics of 2 photons detected at PeV energies and an expected number of 0.3 photons in our model, the tension between data and model at this point is at the level of one standard deviation.

### 5.5. Magnetic field

The radial dependence of the magnetic field downstream from the termination shock was assumed to follow a power law with a power-law index, such that $B(r) = B_0(r/r_s)^{-\alpha}$. The VHE gamma-ray observations are essential to determine the value of $\alpha$ as can be seen from Fig. 10. A more accurate measurement of the HE extension bears the potential to measure $\alpha$ independently of the spectral shape (see Fig. 11). The current precision of the gamma-ray extension measurement is however not sufficient to constrain $\alpha$ strongly. Taking the synchrotron and HE gamma-ray measurements alone, a value of $\tilde{\alpha} = 0.8^{+0.16}_{-0.13}$ and $\tilde{B}_0 = (497 \pm 43) \mu G$ is found. This fit is however constrained in the way, that the HE data are only used after the electron spectra are determined from the synchrotron SED. In the more general fit including the VHE gamma-ray data, the parameter $\alpha$ is constrained more accurately. The benchmark fit to the combined VHE spectra taken with MAGIC, VERITAS, Tibet ASY, and LHAASO favours a value of $\tilde{\alpha} = 0.51 \pm 0.03$ and $\tilde{B}_0 = (264 \pm 9) \mu G$. The resulting total energy in magnetic field in the nebula amounts to $W_B = 7.1 \times 10^{47}$ ergs which is the same energy as in particles present<fnref n="9" />.

The measured value of $\alpha = 0.51$ excludes the assumption of a wound up toroidal magnetic field (this corresponds to $\alpha = -1$). In order to avoid the winding up of the toroidal field in the slowing wind, the wind has to remain fast or the magnetic field energy is dissipated (e.g., through reconnection through tear instabilities) in the plasma.

The velocity of the downstream flow has been estimated using the ratio of X-ray brightness of the advancing and receding part of the torus. The resulting estimate of the flow velocity by Shibata et al. (2003), Mori et al. (2004) requires a faster flow than predicted in the model of Kennel & Coroniti (1984). Conversely, the role of dissipation of magnetic fields in the downstream flow has been discussed in the context of confinement of the nebula (see e.g., Tanaka et al. 2018) as well as particle acceleration (see e.g., Tanaka & Asano 2017; Lyutikov et al. 2019; Luo et al. 2020).

<footnote n="7">⁷ Due to the Klein-Nishina effect, the scattering with optical synchrotron emission is suppressed.</footnote>

<footnote n="8">⁸ Note, the observed 2 events at around 1 PeV have a probability of 15.6% to be observed, see Appendix D for details.</footnote>

<footnote n="9">⁹ The value depends upon the volume of integration, here we choose to integrate up to $25r_s$.</footnote>
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Comparison of the energy density $w_e(\theta)$ of the electrons and of the magnetic field $w_B(\theta)$. The two additional (green) lines indicate the extension $r_{68}$ for the synchrotron and for the inverse Compton component: the y-axes indicates the photon energy $h\nu$ in units of eV (synchrotron, solid line) or GeV (IC, dash-dotted line).</figcaption></figure>

The fit of the radiative model to the complete data set provides an estimate of the radial dependence of the energy density present in particles and in magnetic field. In Fig. 14, the energy density $w_B$ of the magnetic field dominates over the particle energy density in the downstream region until a distance of about $r \approx 1.3r_s$. The reconstructed magnetic field at the termination shock of $\dot{B}_0 = (264 \pm 9) \mu\text{G}$ corresponds to a value of $\sigma \approx 0.1$, substantially larger than the previous estimates of $\sigma \approx 0.003$. The shock is therefore not sufficiently strong (and mainly hydrodynamic) to slow down the relativistic wind in the downstream region.

Even though the nominal uncertainty on $\alpha$ is small ($0.51 \pm 0.0$), the VHE spectral measurements show significant differences in the shape. The estimates of $\alpha$ therefore differ when using different sets of VHE measurements (see also Appendix A and Table A.1). An alternative combination of measurements that are consistent with each other (see also Appendix C) favours a significantly different value of $\dot{\alpha} = 0.29^{+0.03}_{-0.06}$. Without further data, it is impossible to determine whether these measurements exclude each other or whether the underlying systematic uncertainty on the shape of the energy spectrum is larger than estimated by the instrumental teams. The contemporaneous data taken during a flaring episode of the nebula are consistent with each other (see Table A.1). However, when looking at the temporal sequence of the observations (see Fig. A.1) it appears conceivable that the dynamical changes visible in the relevant part of the pulsar wind between the termination shock and the torus could be caused or accompanied by a change in the dissipation rate of magnetic field energy in the flow.

In addition to the disparity between the best-fitting values for $\alpha$ obtained with different VHE measurements (0.29...0.51), there is also a noticeable (but less significant) difference to the best-fitting value found for $\dot{\alpha} = 0.8^{+0.16}_{-0.13}$ in Sect. 4.4.1. For that fit, the spectral data from Fermi-LAT were used only—without including VHE data which covers the SED beyond the peak of the inverse-Compton component (see Sect. 4.4.2). In the context of the discussion related to the potential variability of the magnetic field configuration (see above) and the possible correlation with the hard X-ray flux (see Fig. A.1) it is clear that the value of 0.8 is too far off the range of values found with the VHE data to be explained by variability. Putting aside the possibility of systematic uncertainties, there remains as explanation that the model is missing an additional ingredient. Without attempting

an additional fit here, a possible modification would be a more complicated model for the magnetic field strength and its variation in the nebula. An obvious modification could be applied to the radial dependence. The electrons contributing to the inverse Compton emission detected with Fermi-LAT are spatially more extended than the electrons radiating at VHE energies. In this sense, the value of $\alpha$ could smoothly change, steepening with increasing distance from a value of 0.5 to 1. In this way, a more accurate spectral measurement in the energy range covered with Fermi-LAT could establish a more detailed picture of the magnetic field distribution in the Crab Nebula.

## 6. Summary and conclusion

The observed broadband emission and spatial extension of the Crab Nebula have been used to reconstruct the contemporaneous properties of the relativistic electrons, magnetic field, and dust content in the pulsar wind nebula. The main findings can be summarised in the following points:

1. The assembled observational data for the Crab Nebula are unique as it covers an extremely broad energy range with an accuracy that is generally not appreciated when attempting to model the data.

2. The recent measurement of the spatial extension of the inverse Compton nebula provides insights into the magnetic field strength in the downstream region of the termination shock.

3. A detailed radiative model has been developed to estimate a minimal set of parameters related to the spatial and spectral distribution of electrons, dust, and the magnetic field under the assumption of radial symmetry.

4. The best-fitting model provides a good fit to the data. The $\chi^2$ minimisation approach provides an accurate estimate of all the underlying parameters without degeneracies.

5. The VHE gamma-ray data provide the strongest constraint on the magnetic field strength beyond the termination shock at $r_s$: $B(r) = B_0(r/r_s)^{-\alpha}$ with $\dot{\alpha} = 0.51 \pm 0.03$ and $\dot{B}_0 = (264 \pm 9) \mu\text{G}$.

6. The long-standing assumption introduced about 50 years ago that the downstream magnetic field increases in a kinetically dominated pulsar wind ($\sigma = 0.003$) is ruled out by the data.

7. The downstream flow is dominated by magnetic energy density up to $r \approx 1.3r_s$. For larger distances, the particle energy density dominates until equipartition is recovered in the periphery of the PWN.

8. VHE data taken during an increased hard X-ray flux state favour a value of $\dot{\alpha} = 0.29^{+0.03}_{-0.05}$, with an overall slightly worse goodness of fit.

9. The simultaneous fit of the synchrotron and dust emission provides an estimate of the total dust mass, temperature, and spatial distribution in the nebula. The resulting dust mass is approximately five times larger than previous estimates (De Looze et al. 2019). The larger dust mass increases the required efficiency for dust formation.

10. The production process of the two PeV photons detected from the Crab Nebula remains of uncertain origin. The simplest model would link the photon production to inverse Compton emission of PeV electrons which radiate 80–800 MeV synchrotron emission. Since the majority (>75%) of this energetic synchrotron emission has been found to be produced in a compact region (Yeung & Horns 2020) with fast variability, the resulting constraint on the minimum (confining) magnetic field (>500 µG) is in conflict with the comparably low value of magnetic field (≈85 µG) required to fit the inverse Compton and synchrotron emission simultaneously.
The results obtained here mainly rely on the assumptions of a spherical and static system. In principle, the model can be extended to axial symmetry and to include the effects of the flow in the radiative treatment. On the observational side, multi-instrument simultaneous VHE observations of the Crab Nebula are crucial to verify the effect of spectral variations. The upcoming spatially resolved observations of the X-ray polarisation of the Crab Nebula with the IXPE satellite (Weisskopf et al. 2016) will provide a measurement of the degree of small-scale turbulence in the flow (Bucciantini et al. 2017) which lead to a reduced fraction of polarised synchrotron emission.

**Acknowledgements.** Part of this work was supported by the German Deutsche Forschungsgemeinschaft, DFG project number 284334853: HO 3305/2 and by ANR-20-CE31-0010. This research has made use of NASA's Astrophysics Data System Bibliographic Services. The authors thank the developers of matplotlib (Hunter 2007) for making their code available on a free and open-source basis. DH acknowledges support by the DFG under Germany's Excellence Strategy - EXC 2121 "Quantum Universe" - 390833306. We thank the referee for the very useful reports and for encouraging us to investigate the differences found between the Fermi-LAT and VHE fit results.

## References

Abdalla, H., Aharonian, F., Ait Benkhali, F., et al. 2020, Nat. Astron., 4, 167

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, Phys. Rev. Lett., 102, 181101

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, Science, 331, 739

Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2019, ApJ, 881, 134

Ackermann, M., Ajello, M., Allafort, A., et al. 2012, Astropart. Phys., 35, 346

Aharonian, F., Akhperjanian, A., Beilicke, M., et al. 2004, ApJ, 614, 897

Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, A&A, 457, 899

Aharonian, F. A., Kelner, S. R., & Prosekin, A. Y. 2010, Phys. Rev. D., 82, 043002

Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2015, J. High Energy Astrophys., 5, 30

Aliu, E., Archambault, S., Aune, T., et al. 2014, ApJ, 781, L11

Amenomori, M., Bao, Y. W., Bi, X. J., et al. 2019, Phys. Rev. Lett., 123, 051101

Atoyan, A. M., & Aharonian, F. A. 1996, MNRAS, 278, 525

Band, D., Matteson, J., Ford, L., et al. 1993, ApJ, 413, 281

Bietenholz, M. F., & Kronberg, P. P. 1990, ApJ, 357, L13

Bietenholz, M. F., Hester, J. J., Frail, D. A., & Bartel, N. 2004, ApJ, 615, 794

Blair, W. P., Long, K. S., Vancura, O., et al. 1992, ApJ, 399, 611

Bohren, C., & Huffman, D. R. 1998, in Absorption and Scattering of Light by Small Particles, eds. C. Bohren, & D. R. Huffman (Wiley Science Paperback Series)

Bucciantini, N., Bandiera, R., Olmi, B., & Del Zanna, L. 2017, MNRAS, 470, 4066

Buehler, R., Scargle, J. D., Blandford, R. D., et al. 2012, ApJ, 749, 26

Čadež, A., Carramiñana, A., & Vidrih, S. 2004, ApJ, 609, 797

Crusius, A., & Schlickeiser, R. 1986, A&A, 164, L16

de Jager, O. C., & Harding, A. K. 1992, ApJ, 396, 161

De Looze, I., Barlow, M. J., Bandiera, R., et al. 2019, MNRAS, 488, 164

Feng, H., Li, H., Long, X., et al. 2020, Nat. Astron., 4, 511

Gomez, H. L., Krause, O., Barlow, M. J., et al. 2012, ApJ, 760, 96

Gould, R. J., & Schréder, G. P. 1967, Phys. Rev., 155, 1404

Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93

Greiveldinger, C., & Aschenbach, B. 1999, ApJ, 510, 305

Hennessy, G. S., O'Connell, R. W., Cheng, K. P., et al. 1992, ApJ, 395, L13

H. E. S. S. Collaboration (Abdalla, H., et al.) 2020, Nat. Astron., 4, 167

H. E. S. S. Collaboration (Abramowski, A., et al.) 2014, A&A, 562, L4

Hickson, P., & van den Bergh, S. 1990, ApJ, 365, 224

Hillas, A. M., Akerlof, C. W., Biller, S. D., et al. 1998, ApJ, 503, 744

Holler, M. 2014, PhD Thesis, Universität Potsdam, Germany

Holler, M., Berge, D., van Eldik, C., et al. 2015, ArXiv e-prints [arXiv:1509.02902]

Horns, D., Aharonian, F., Santangelo, A., Hoffmann, A. I. D., & Masterson, C. 2006, A&A, 451, L51

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Jones, F. C. 1968, Phys. Rev., 167, 1159

Jourdain, E., & Roques, J. P. 2009, ApJ, 704, 17

Jourdain, E., & Roques, J. P. 2020, ApJ, 899, 131

Kennel, C. F., & Coroniti, F. V. 1984, ApJ, 283, 694

Khangulyan, D., Arakawa, M., & Aharonian, F. 2020, MNRAS, 491, 3217

Kirsch, M. G., Briel, U. G., Burrows, D., et al. 2005, SPIE Conf. Ser., 5898, 22

Komissarov, S. S., & Lyubarsky, Y. E. 2003, MNRAS, 344, L93

Krimm, H. A., Holland, S. T., Corbet, R. H. D., et al. 2013, ApJS, 209, 14

Kuiper, L., Hermsen, W., Cusumano, G., et al. 2001, A&A, 378, 918

LHAASO Collaboration (Cao, Z., et al.) 2021, Science, 373, 425

Ling, J. C., & Wheaton, W. A. 2003, ApJ, 598, 334

Liu, R.-Y., & Wang, X.-Y. 2021, ApJ, 922, 221

Lobanov, A. P., Horns, D., & Muxlow, T. W. B. 2011, A&A, 533, A10

Luo, Y., Lyutikov, M., Temim, T., & Comisso, L. 2020, ApJ, 896, 147

Lyutikov, M., Temim, T., Komissarov, S., et al. 2019, MNRAS, 489, 2403

Madsen, K. K., Reynolds, S., Harrison, F., et al. 2015, ApJ, 801, 66

Madsen, K. K., Forster, K., Grefenstette, B. W., Harrison, F. A., & Stern, D. 2017, ApJ, 841, 56

Madsen, K. K., Burwitz, V., Forster, K., et al. 2021, ArXiv e-prints [arXiv:2111.01613]

MAGIC Collaboration (Acciari, V. A., et al.) 2020, A&A, 635, A158

Marsden, P. L., Gillett, F. C., Jennings, R. E., et al. 1984, ApJ, 278, L29

Martín, J., Torres, D. F., & Rea, N. 2012, MNRAS, 427, 415

Mayer, M., Buehler, R., Hays, E., et al. 2013, ApJ, 775, L37

Meagher, K. & VERITAS Collaboration 2016, Proc. 34th International Cosmic Ray Conference (ICRC2015), 30 July - 6 August, 2015, The Hague, The Netherlands, online at http://pos.sissa.it/cgi-bin/reader/conf.cgi?confid=236, 792

Meyer, M., Horns, D., & Zechlin, H.-S. 2010, A&A, 523, A2

Mori, K., Burrows, D. N., Hester, J. J., et al. 2004, ApJ, 609, 186

O'Donnell, J. E. 1994, ApJ, 422, 158

Owen, P. J., & Barlow, M. J. 2015, ApJ, 801, 141

Pacini, F., & Salvati, M. 1973, ApJ, 186, 249

Peng, Q.-Y., Bao, B.-W., Lu, F.-W., & Zhang, L. 2022, ApJ, 926, 7

Planck Collaboration XXVI. 2016, A&A, 594, A26

Popescu, C. C., Yang, R., Tuffs, R. J., et al. 2017, MNRAS, 470, 2539

Porter, T. A., Jóhannesson, G., & Moskalenko, I. V. 2017, ApJ, 846, 67

Porth, O., Komissarov, S. S., & Keppens, R. 2013, MNRAS, 438, 278

Rees, M. J., & Gunn, J. E. 1974, MNRAS, 167, 1

Ritacco, A., Macías-Pérez, J. F., Ponthieu, N., et al. 2018, A&A, 616, A35

Rudy, A., Horns, D., DeLuca, A., et al. 2015, ApJ, 811, 24

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

Shaposhnikov, N., Jahoda, K., Markwardt, C., Swank, J., & Strohmayer, T. 2012, ApJ, 757, 159

Shibata, S., Tomatsuri, H., Shimanuki, M., Saito, K., & Mori, K. 2003, MNRAS, 346, 841

Smith, N. 2003, MNRAS, 346, 885

Strong, A., & Collmar, W. 2019, Mem. Soc. Astron. Ital., 90, 297

Strüder, L., Briel, U., Dennerl, K., et al. 2001, A&A, 365, L18

Tanaka, S. J., & Asano, K. 2017, ApJ, 841, 78

Tanaka, S. J., Toma, K., & Tominaga, N. 2018, MNRAS, 478, 4622

Tavani, M., Bulgarelli, A., Vittorini, V., et al. 2011, Science, 331, 736

Temim, T., & Dwek, E. 2013, ApJ, 774, 8

Temim, T., Gehrz, R. D., Woodward, C. E., et al. 2006, AJ, 132, 1610

Tomsick, J., & COSI Collaboration 2022, 37th International Cosmic Ray Conference, 12-23 July 2021, Berlin, Germany, online at https://pos.sissa.it/cgi-bin/reader/conf.cgi?confid=395, 652

Toor, A., & Seward, F. D. 1974, AJ, 79, 995

Trimble, V. 1973, PASP, 85, 579

Tsujiimoto, M., Guainazzi, M., Plucinsky, P. P., et al. 2011, A&A, 525, A25

Véron-Cetty, M. P., & Woltjer, L. 1993, A&A, 270, 370

Weekes, T. C., Cawley, M. F., Fegan, D. J., et al. 1989, ApJ, 342, 379

Weiland, J. L., Odegard, N., Hill, R. S., et al. 2011, ApJS, 192, 19

Weisskopf, M. C., Silver, E. H., Kestenbaum, H. L., Long, K. S., & Novick, R. 1978, ApJ, 220, L117

Weisskopf, M. C., Guainazzi, M., Jahoda, K., et al. 2010, ApJ, 713, 912

Weisskopf, M. C., Elsner, R. F., Kolodziejczak, J. J., O'Dell, S. L., & Tennant, A. F. 2012, ApJ, 746, 41

Weisskopf, M. C., Tennant, A. F., Arons, J., et al. 2013, ApJ, 765, 56

Weisskopf, M. C., Ramsey, B., O'Dell, S., et al. 2016, SPIE Conf. Ser., 9905, 990517

Wells, R. M. 2019, PhD Thesis, Iowa State University, Graduate Theses and Dissertations, 1706

Wilson-Hodge, C. A. 2012, When a Standard Candle Flickers: Crab Nebula Variations in Hard X-rays, Tech. rep., NASA Marshall Space Flight Center Huntsville, USA

Wilson-Hodge, C. A., Cherry, M. L., Case, G. L., et al. 2011, ApJ, 727, L40

Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181

Yeung, P. K. H., & Horns, D. 2019, ApJ, 875, 123

Yeung, P. K. H., & Horns, D. 2020, A&A, 638, A147

Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321
# Appendix A: Individual VHE data sets

<figure><img src="image_17.png" /><figcaption>**Fig. A.1.** Observed hard X-ray light curve (RXTE-PCA data from (Wilson-Hodge et al. 2011), Swift-BAT (Krimm et al. 2013) available through the BAT transient monitor page) from the Crab Nebula. The observed flux varies on time-scales of several years. The observations at VHE (listed in Table A.1) cover different parts of the X-ray flux states.</figcaption></figure>

<table><caption>**Table A.1.** Summary of VHE gamma-ray data sets used here.</caption><thead><tr><th>Data set</th><th>Start</th><th>End</th><th>χ<sup>2</sup>(dof)</th><th>α</th><th>ζ</th><th>Ref.</th></tr></thead><tbody><tr><td>HEGRA (IACT)</td><td>1997</td><td>2002</td><td>18(15)</td><td>0.26<sup>+0.09</sup><sub>-0.04</sub></td><td>1.197 ± 0.005</td><td>(1)</td></tr><tr><td>H.E.S.S. (IACT)</td><td>2003</td><td>2005</td><td>24(10)</td><td>0.38<sup>+0.02</sup><sub>-0.09</sub></td><td>1.065 ± 0.005</td><td>(2)</td></tr><tr><td>MAGIC (IACT)</td><td>2009</td><td>2011</td><td>11(13)</td><td>0.495<sup>+0.05</sup><sub>-0.07</sub></td><td>0.896 ± 0.006</td><td>(3)</td></tr><tr><td>VERITAS (IACT)</td><td>2009</td><td>2013</td><td>9(11)</td><td>0.535<sup>+0.02</sup><sub>-0.06</sub></td><td>1.038 ± 0.002</td><td>(4)</td></tr><tr><td>HAWC (WCD)</td><td>6/2015</td><td>12/2017</td><td>7(8)</td><td>0.27<sup>+0.03</sup><sub>-0.08</sub></td><td>1.093 ± 0.003</td><td>(5)</td></tr><tr><td>MAGIC (VLZA)</td><td>2014</td><td>2018</td><td>2(6)</td><td>0.59 ± 0.25</td><td>1.09 ± 0.02</td><td>(6)</td></tr><tr><td>Tibet ASy (EAS, WCD)</td><td>2014</td><td>2017</td><td>4(9)</td><td>0.32 ± 0.15</td><td>1.06 ± 0.02</td><td>(7)</td></tr><tr><td>VERITAS (flare)</td><td>03/2013</td><td>03/2013</td><td>13(13)</td><td>0.42 ± 0.17</td><td>1.09 ± 0.01</td><td>(8)</td></tr><tr><td>H.E.S.S. II (flare)</td><td>03/2013</td><td>03/2013</td><td>17(31)</td><td>0.68 ± 0.3</td><td>0.96 ± 0.01</td><td>(9)</td></tr><tr><td>H.E.S.S. II (IACT)</td><td>11/2013</td><td>1/2014</td><td>93(37)</td><td>0.3<sup>+0.1</sup><sub>-0.05</sub></td><td>1.100 ± 0.007</td><td>(10)</td></tr><tr><td>LHAASO (EAS, WCD)</td><td>2019</td><td>2020</td><td>11(16)</td><td>0.56 ± 0.06</td><td>0.932 ± 0.006</td><td>(11)</td></tr></tbody></table>

**Notes.** Listed data sets include imaging air Cherenkov telescopes' (IACTs) observations as well as data taken with extended air shower front sampling (EAS) arrays with scintillators or water Cherenkov detectors (WCD). By observing sources with IACTs at very large zenith angles (VLZA) between 70° and 80°, the collection area at large energies beyond 10 TeV is substantially increased. The uncertainty on α is given for a 90 % confidence level and has been estimated from the discrete sampled values of α. For some data sets the χ² curve is undersampled while for others, the maximum range of α is not sufficient to cover the confidence interval.

**Reference.** (1) Aharonian et al. (2004), (2) Aharonian et al. (2006), (3) Aleksić et al. (2015), (4) Wells (2019), (5) Abeysekara et al. (2019), (6) MAGIC Collaboration (2020), (7) Amenomori et al. (2019), (8) Aliu et al. (2014), (i) H. E. S. S. Collaboration (2014), (9) Holler et al. (2015), Holler (2014), (10) LHAASO Collaboration (2021).

The observation of VHE ($E > 100$ GeV) gamma-rays is the domain of ground based air shower detection. In this energy range, the Crab Nebula is the brightest steady source in the sky and therefore, the detection and observation of the Crab Nebula is the litmus test for VHE gamma-ray detectors. The fact that the very first TeV gamma-ray source clearly detected was the Crab Nebula (Weekes et al. 1989) underlines the importance of this particular source for ground based gamma-ray astronomy. The routine monitoring of this standard candle has lead to a large number of published data sets in this energy range with

widely different detection techniques. We have selected observations which cover both a wide range of observation times (from 1997 to 2020) and different techniques (imaging air Cherenkov telescopes as well as shower front sampling with extended air shower arrays). This is not a complete collection but it is meant to provide a good cross section of the spectroscopic data from the past two decades.

The individual data sets are listed in Table A.1. The observational times as well as ranges of best-fitting values of α are included with the hard X-ray light curve from the nebula in Fig. A.1. The contemporaneous data sets are consistent within the uncertainties of α. When looking at the evolution of the hard X-ray flux and the changes of α, there is at least a trend visible: during higher flux states of hard X-rays, the VHE spectra tend to follow the shape given for a slower radial drop of the magnetic field. During lower flux states, the magnetic field seems to go down quicker with increasing distance.

# Appendix B: Pair-production opacity towards the Crab Nebula

The soft photon fields present in the Galactic disk include the cosmic microwave background (CMB) as well as dust and stellar emission. Energetic gamma-rays will be absorbed by pair-production processes when propagating in the soft photon background Gould & Schréder (1967). Using a recent radiation model of the Galaxy (Popescu et al. 2017), the optical depth for the line of sight towards the Crab Nebula is calculated to correct its apparent gamma-ray brightness for absorption (LHAASO Collaboration 2021). The CMB leads to a reduction by ≈ 24 % at $E_γ = 2$ PeV (assuming a distance of 2 kpc). At lower energies, the CMB photons are not sufficiently energetic to create pairs such that at $E_γ = 100$ TeV pairs are mainly produced with thermal dust emission. Given that the photon density of the dust emission is substantially smaller at the relevant wavelength of ≈ 100 µm, the resulting absorption is less pronounced than the CMB-related absorption at higher energies ($1 - \exp(-\tau_{γγ}(E_γ = 100 \text{ TeV})) \approx 1.5 \%$). However, the result on the opacity used by LHAASO Collaboration (2021) is not taking into account the anisotropy of the soft photon field.

Here, we provide an estimate of the previously neglected effect of anisotropy by using the following simplifications: (i) the relevant specific energy density $\lambda u_\lambda$ of the photon field is dominated by the dust component with a peak at $\lambda \approx 100$ µm, (ii) the angular distribution of the anisotropic background field follows the measured local intensity, (iii) both, the intensity as well angular distribution remain constant along the line of sight towards the Crab Nebula.

The resulting optical depth $\tau_{γγ}(E_γ)$ for the Crab Nebula located at a distance $d$ and in the direction along the normal vector $e_C$ is given by:

$$ \tau(E_γ) = \int_0^d dx \int_{\varepsilon_{min}}^{\varepsilon_{max}} d\varepsilon \, n(\varepsilon) \int d\Omega \, f(\mathbf{e}_{\Omega}) \sigma_{γγ}(s) (1 - \mathbf{e}_{\Omega} \cdot \mathbf{e}_C), \quad (B.1) $$

where $\sigma_{γγ}$ is the pair production cross section (see e.g. Eq. 1 of Gould & Schréder (1967)), $n(\varepsilon)$ is the total (sky integrated) differential number density of photons $dn = n(\varepsilon)d\varepsilon dV$. The angular distribution $f(\mathbf{e}_{\Omega})$ was normalised over the entire sphere, such that:

$$ \int d\Omega f(\mathbf{e}_{\Omega}) = 1. \qquad (B.2) $$
<figure><img src="image_18.png" /><figcaption>**Fig. B.1.** Energy-dependent attenuation of gamma rays for the line of sight towards the Crab Nebula: The solid and dashed lines are the attenuation calculated here (see text for details) for the anisotropic radiation field (solid) and for the assumption of an isotropic radiation field with identical energy density (dashed line). The dash-dotted line is the attenuation as calculated in the framework of a detailed Galactic dust emission model, however assuming an isotropic radiation field (LHAASO Collaboration 2021)</figcaption></figure>

The unit vector $\mathbf{e}_{\Omega}$ points towards a direction given in Galactic coordinates. This way, $\mu := \mathbf{e}_{\Omega} \cdot \mathbf{e}_{C}$ corresponds to the cosine of the angular distance between the photon directions of the gamma ray from the Crab Nebula and the soft photon from the Galactic background emission. The centre of momentum energy $s$ depends on the photon energy $E_{\gamma}$, the soft photon's energy $\varepsilon$ as well as $\mu$:

$$s = 2\varepsilon E_{\gamma}(1 - \mu). \qquad (B.3)$$

Following the approximation (i), the dust photon field was taken to be a single function

$$n(\varepsilon) = n_0 \left(\frac{\varepsilon}{\varepsilon_0}\right)^{3.5} \frac{1}{\exp(\varepsilon/\varepsilon_0) - 1}, \qquad (B.4)$$

with $n_0 = 1223.78 \text{ cm}^{-3} \text{ eV}^{-1}$ and $\varepsilon_0 = 1.7 \times 10^{-3} \text{ eV}$. This is a reasonably close approximation of the model proposed by Popescu et al. (2017) around the peak at $\lambda \approx 100 \text{ µm}$.

The angular distribution at $\lambda = 100 \text{ µm}$ was taken from the all-sky map from DIRBE measurements (Schlegel et al. 1998) with the zodiacal foreground and point sources removed. In order to highlight the effect of the anisotropy, we calculated the optical depth for both cases (anisotropic and isotropic photon field) with the same energy density $\lambda u_{\lambda}$ and compare the results in Fig. B.1. While the simplified model reproduces well the more detailed calculation, taking additional photon fields as well as the spatial dependence of the photon density (Popescu et al. 2017) into account, the anisotropic case increases the attenuation as expected. At $E_{\gamma} = 100 \text{ TeV}$, the flux is attenuated by about 3 % (as compared to about 1.7 % for the isotropic case). It is important to note, that this approach is particular to sources close to the Sun's position, such that the approximations (ii) and (iii) allow the simplified treatment of the anisotropic case. The general anisotropic case has been treated in Porter et al. (2017) (and references therein).

## Appendix C: Low $\alpha$ VHE solution

The combination of VHE data sets shown in Sect. 4.4.2 favours a value of $\alpha = 0.51 \pm 0.03$. However, the VHE data sets show individually quite different fitting results as demonstrated by the values in Table A.1. While there is a large number of possible combinations conceivable, there seems to be at least one additional combination of data sets that would favour a smaller value of $\alpha$ while still maintaining a good fit to the data. Excluding therefore the two data sets obtained with the H.E.S.S. telescopes which are not fit well by the model (see Table A.1), we combined the spectra measured with HEGRA, HAWC, and the Tibet ASγ experiments. The resulting estimate of $\alpha = 0.29^{+0.03}_{-0.06}$ is significantly smaller than the value found for the data sets shown in Sect. 4.4.2.

The resulting fit (shown in Fig. C.1) is still acceptable with $\chi^2_{VHE}(dof) = 30(32)$. The Tibet ASγ data requires an up-scaling of the energy scale ($\zeta = 1.07(2)$) to match the model while in the case of the $\alpha = 0.51$ fit, $\zeta = 0.92(2)$ (see Table C.1 for the results of this fit). The best-fitting parameters listed in Table C.2 are systematically different to the parameters listed for the solution with $\alpha = 0.51$ (see Table 4). For the solution shown here, the magnetic field at the shock is smaller and drops off slower. To compensate, the spatial particle distribution does not extend as far as in the case of $\alpha = 0.51$ (see parameters $\Psi_3, \Psi_{14}$). This in turn reduces the size of the IC nebula (see also Fig. 11). The resulting $\chi^2_{IC,ext}$ is therefore worse than in the solution with $\alpha = 0.51$. Combined with the worse VHE fit, the probability to obtain this large value of $\chi^2(dof) = 289(262)$ is only $p(> \chi^2, dof) = 0.12$. Even though the goodness of fit is noticeably worse than in the $\alpha = 0.51$ case, it is still acceptable.

<figure><img src="image_19.png" /><figcaption>**Fig. C.1.** When combining the energy spectra obtained with HEGRA, HAWC, and Tibet ASγ, the data-sets are consistent with each other (after scaling the energy) and can be described with the model reasonably well: $\chi^2_{VHE}(dof) = 30(32)$.</figcaption></figure>

<table><caption>**Table C.1.** Energy-scaling for the solution with $\tilde{\alpha} = 0.29$.</caption><thead><tr><th>Data set<br>(α̃ = 0.29)</th><th>ζ</th><th>σ<sub>ζ</sub></th><th>χ<sup>2</sup><sub>before</sub>(dof)</th><th>χ<sup>2</sup><sub>after</sub>(dof)</th></tr><tr><td></td><td colspan="4" style="text-align: center;">68 % c.l.</td></tr></thead><tbody><tr><td>Fermi-LAT</td><td>1.00</td><td>fixed</td><td>26.6(25)</td><td>26.6(25)</td></tr><tr><td>HEGRA</td><td>1.216</td><td>0.006</td><td>3050.7(16)</td><td>18.1(15)</td></tr><tr><td>HAWC</td><td>1.083</td><td>0.004</td><td>556.8(9)</td><td>7.7(8)</td></tr><tr><td>Tibet ASγ</td><td>1.073</td><td>0.020</td><td>21.4(10)</td><td>4.2(9)</td></tr><tr><td>Combined</td><td></td><td></td><td>3655.5(60)</td><td>56.6(57)</td></tr></tbody></table>
<table><caption>**Table C.2.** Best-fitting parameters for the solution with $\alpha = 0.29$.</caption><thead><tr><th>Parameter</th><th>best-fitting values<br>(68 % c.l.)</th></tr></thead><tbody><tr><td colspan="2"><em>Radio electrons</em></td></tr><tr><td>Ψ<sub>1</sub> = s<sub>r</sub></td><td>1.56 ± 0.03</td></tr><tr><td>Ψ<sub>2</sub> = ln(N<sub>r,0</sub>)</td><td>114.8 ± 0.2</td></tr><tr><td>Ψ<sub>3</sub> = ln(γ<sub>1</sub>)</td><td>11.5 ± 0.1</td></tr><tr><td>Ψ<sub>4</sub> = ρ<sub>r</sub> ["]</td><td>83 ± 3</td></tr><tr><td colspan="2"><em>Wind electrons</em></td></tr><tr><td>Ψ<sub>5</sub> = s<sub>1</sub></td><td>2.9 ± 0.1</td></tr><tr><td>Ψ<sub>6</sub> = s<sub>2</sub></td><td>3.3 ± 0.01</td></tr><tr><td>Ψ<sub>7</sub> = s<sub>3</sub></td><td>3.62 ± 0.04</td></tr><tr><td>Ψ<sub>8</sub> = ln(γ<sub>w0</sub>)</td><td>12.7 ± 0.2</td></tr><tr><td>Ψ<sub>9</sub> = ln(γ<sub>w1</sub>)</td><td>15.8 ± 0.8</td></tr><tr><td>Ψ<sub>10</sub> = ln(γ<sub>w2</sub>)</td><td>19.3 ± 0.2</td></tr><tr><td>Ψ<sub>11</sub> = ln(γ<sub>w3</sub>)</td><td>22.5 ± 0.04</td></tr><tr><td>Ψ<sub>12</sub> = ln(N<sub>w,0</sub>)</td><td>74.3 ± 0.5</td></tr><tr><td>Ψ<sub>13</sub> = β</td><td>0.14 ± 0.01</td></tr><tr><td>Ψ<sub>14</sub> = ρ<sub>0</sub> ["]</td><td>83 ± 3</td></tr><tr><td colspan="2"><em>Dust parameters</em></td></tr><tr><td>Ψ<sub>15</sub> = r<sub>out</sub> [pc]</td><td>1.57 ± 0.12</td></tr><tr><td>Ψ<sub>16</sub> = log<sub>10</sub>(M<sub>1</sub>/M<sub>☉</sub>)</td><td>-4.4 ± 0.1</td></tr><tr><td>Ψ<sub>17</sub> = log<sub>10</sub>(M<sub>2</sub>/M<sub>☉</sub>)</td><td>-1.2 ± 0.1</td></tr><tr><td>Ψ<sub>18</sub> = T<sub>1</sub> [K]</td><td>148 ± 8</td></tr><tr><td>Ψ<sub>19</sub> = T<sub>2</sub> [K]</td><td>39 ± 3</td></tr><tr><td colspan="2"><em>Magnetic field parameters</em></td></tr><tr><td>Ψ<sub>20</sub> = B<sub>0</sub> [μG]</td><td>167 ± 5</td></tr><tr><td>Ψ<sub>21</sub> = α</td><td>0.29<sup>+0.03</sup><sub>-0.06</sub></td></tr><tr><td colspan="2"><em>Goodness of fit</em></td></tr><tr><td>χ<sup>2</sup><sub>sync,SED</sub> (dof)</td><td>185 (184)</td></tr><tr><td>χ<sup>2</sup><sub>sync,ext</sub> (dof)</td><td>17 (15)</td></tr><tr><td>χ<sup>2</sup><sub>IC,SED</sub> (dof)</td><td>27 (23)</td></tr><tr><td>χ<sup>2</sup><sub>IC,ext</sub> (dof)</td><td>30 (8)</td></tr><tr><td>χ<sup>2</sup><sub>VHE</sub> (dof)</td><td>30 (32)</td></tr><tr><td>χ<sup>2</sup><sub>tot</sub> (dof)</td><td>289 (262)</td></tr></tbody></table>

**Appendix D. UHE gamma-ray event statistic**

In the main section, the expected number of events for the LHAASO observation (LHAASO Collaboration 2021) were used to compare with the observed number of events detected at PeV energies. The underlying calculation is performed using the exposure time ($T_{\text{obs}} = 10^7$ s) and the collection areas listed in the published data<fnref n="10" />. The resulting event numbers $N_i$ in bin $i$ are calculated by integrating the differential photon flux predicted by the model $\Phi_{\text{model}}(E)$ within the bin interval $E_{i,\text{low}}$ and $E_{i,\text{high}}$:

$$N_i = T_{\text{obs}} A_{\text{eff},i} \int_{E_{i,\text{low}}}^{E_{i,\text{high}}} dE \, \Phi_{\text{model}}(E/\zeta) e^{-\tau_{\gamma\gamma}(E/\zeta)}. \qquad (\text{D.1})$$

In Tab. D.1, the resulting event numbers are compared with the background subtracted excess events from the data. For the two bins just below and above 1 PeV, the expected number of events is $\lambda = 0.7$ while the number of observed events (background free) is $k = 2$. The probability to observe 2 or more events is calculated with the Poissonian distribution

$$p(k \ge 2) = 1-p_{\lambda}(0)-p_{\lambda}(1) = 1-\exp(-\lambda)(1+\lambda) = 15.6 \%. \quad (\text{D.2})$$

<table><caption>**Table D.1.** Event statistics for the LHAASO data set: comparison of model and data.</caption><thead><tr><th>log<sub>10</sub>(E<sub>i,low</sub>/TeV)</th><th>log<sub>10</sub>(E<sub>i,high</sub>/TeV)</th><th>N<sub>i</sub></th><th>N<sub>i,obs</sub></th></tr></thead><tbody><tr><td>1</td><td>1.2</td><td>3305</td><td>4011</td></tr><tr><td>1.2</td><td>1.4</td><td>1718</td><td>1979</td></tr><tr><td>1.4</td><td>1.6</td><td>601</td><td>649</td></tr><tr><td>1.6</td><td>1.8</td><td>255</td><td>278</td></tr><tr><td>1.8</td><td>2</td><td>145</td><td>151</td></tr><tr><td>2</td><td>2.2</td><td>53.2</td><td>53.6</td></tr><tr><td>2.2</td><td>2.4</td><td>19.1</td><td>22.5</td></tr><tr><td>2.4</td><td>2.6</td><td>6.3</td><td>5.8</td></tr><tr><td>2.6</td><td>2.8</td><td>1.9</td><td>2.9</td></tr><tr><td>2.8</td><td>3.0</td><td>0.5</td><td>1.0</td></tr><tr><td>3.0</td><td>3.2</td><td>0.2</td><td>1.0</td></tr><tr><td>3.2</td><td>4.0</td><td>0.02</td><td>0.0</td></tr></tbody></table>

<footnote n="10">¹⁰ http://english.ihep.cas.cn/lhaaso/pdl/202110/t20211026_286779.html</footnote>