# Spectrum and extension of the inverse-Compton emission of the Crab Nebula from a combined Fermi-LAT and H.E.S.S. analysis

F Aharonian, F. Ait Benkhali, J Aschersleben, H Ashkar, M Backes, A Baktash, V. Barbosa Martins, R Batzofin, Y Becherini, D Berge, et al.

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

## To cite this version:

F Aharonian, F. Ait Benkhali, J Aschersleben, H Ashkar, M Backes, et al. Spectrum and extension of the inverse-Compton emission of the Crab Nebula from a combined Fermi-LAT and H.E.S.S. analysis. Astronomy & Astrophysics - A&A, 2024, 686, pp.A308. 10.1051/0004-6361/202348651. hal-04529117

HAL Id: hal-04529117

https://hal.science/hal-04529117v1

Submitted on 4 Oct 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
Spectrum and extension of the inverse-Compton emission of the Crab Nebula from a combined Fermi-LAT and H.E.S.S. analysis

F. Aharonian<sup>1,2,3</sup>, F. Ait Benkhali<sup>4</sup>, J. Aschersleben<sup>5</sup>, H. Ashkar<sup>6</sup>, M. Backes<sup>7,8</sup>, A. Baktash<sup>9</sup>, V. Barbosa Martins<sup>10</sup>, R. Batzofin<sup>11</sup>, Y. Becherini<sup>12,13</sup>, D. Berge<sup>10,14</sup>, K. Bernlöhr<sup>2</sup>, B. Bi<sup>15</sup>, M. Böttcher<sup>8</sup>, C. Boisson<sup>16</sup>, J. Bolmont<sup>17</sup>, M. de Bony de Lavergne<sup>18</sup>, J. Borowska<sup>14</sup>, F. Bradascio<sup>18</sup>, M. Breuhaus<sup>2</sup>, R. Brose<sup>1</sup>, A. Brown<sup>19</sup>, F. Brun<sup>18</sup>, B. Bruno<sup>20</sup>, T. Bulik<sup>21</sup>, C. Burger-Scheidlin<sup>1</sup>, T. Bylund<sup>18</sup>, S. Caroff<sup>22</sup>, S. Casanova<sup>23</sup>, R. Cecil<sup>9</sup>, J. Celic<sup>20</sup>, M. Cerruti<sup>12</sup>, P. Chambery<sup>24</sup>, T. Chand<sup>8</sup>, S. Chandra<sup>8</sup>, A. Chen<sup>25</sup>, J. Chibueze<sup>8</sup>, O. Chibueze<sup>8</sup>, G. Cotter<sup>19</sup>, P. Cristofari<sup>16</sup>, J. Devin<sup>26</sup>, A. Djannati-Ataï<sup>12</sup>, J. Djuvsland<sup>2</sup>, A. Dmytriiev<sup>8</sup>, S. Einecke<sup>27</sup>, J.-P. Ernenwein<sup>28</sup>, S. Fegan<sup>6</sup>, K. Feijen<sup>12</sup>, M. Filipović<sup>29</sup>, G. Fontaine<sup>6</sup>, M. Füßling<sup>10</sup>, S. Funk<sup>20</sup>, S. Gabici<sup>12</sup>, Y. A. Gallant<sup>26</sup>, G. Giavitto<sup>10</sup>, D. Glawion<sup>20</sup>, J. F. Glicenstein<sup>18</sup>, J. Glombitza<sup>20</sup>, P. Goswami<sup>12</sup>, G. Grolleron<sup>17</sup>, M.-H. Grondin<sup>24</sup>, L. Haerer<sup>2</sup>, J. A. Hinton<sup>2</sup>, W. Hofmann<sup>2</sup>, T. L. Holch<sup>10</sup>, M. Holler<sup>30</sup>, D. Horns<sup>9</sup>, M. Jamrozy<sup>31</sup>, F. Jankowsky<sup>4</sup>, V. Joshi<sup>20</sup>, E. Kasai<sup>7</sup>, K. Katarzyński<sup>32</sup>, R. Khatoon<sup>8</sup>, B. Khélifi<sup>12</sup>, W. Klužniak<sup>33</sup>, Nu. Komin<sup>25</sup>, K. Kosack<sup>18</sup>, D. Kostunin<sup>10</sup>, A. Kundu<sup>8</sup>, R. G. Lang<sup>20</sup>, S. Le Stum<sup>28</sup>, F. Leitl<sup>20</sup>, A. Lemière<sup>12</sup>, M. Lemoine-Goumard<sup>24</sup>, J.-P. Lenain<sup>17</sup>, F. Leuschner<sup>15</sup>, A. Luashvili<sup>16</sup>, J. Mackey<sup>1</sup>, D. Malyshev<sup>15</sup>, D. Malyshev<sup>20</sup>, V. Marandon<sup>18</sup>, P. Marinos<sup>27</sup>, G. Martí-Devesa<sup>30</sup>, R. Marx<sup>4</sup>, A. Mehta<sup>10</sup>, M. Meyer<sup>9,*</sup>, A. Mitchell<sup>20</sup>, R. Moderski<sup>33</sup>, L. Mohrmann<sup>2,*</sup>, A. Montanari<sup>4</sup>, E. Moulin<sup>18</sup>, T. Murach<sup>10</sup>, M. de Naurois<sup>6</sup>, J. Niemiec<sup>23</sup>, P. O'Brien<sup>34</sup>, S. Ohm<sup>10</sup>, L. Olivera-Nieto<sup>2</sup>, E. de Ona Wilhelmi<sup>10</sup>, M. Ostrowski<sup>31</sup>, S. Panny<sup>30</sup>, M. Panter<sup>2</sup>, R. D. Parsons<sup>14</sup>, G. Peron<sup>12</sup>, D. A. Prokhorov<sup>35</sup>, G. Pühlhofer<sup>15</sup>, M. Punch<sup>12</sup>, A. Quirrenbach<sup>4</sup>, M. Regeard<sup>12</sup>, P. Reichherzer<sup>18</sup>, A. Reimer<sup>30</sup>, O. Reimer<sup>30</sup>, H. Ren<sup>2</sup>, M. Renaud<sup>26</sup>, B. Reville<sup>2</sup>, F. Rieger<sup>2</sup>, G. Roellinghoff<sup>20</sup>, B. Rudak<sup>33</sup>, V. Sahakian<sup>36</sup>, H. Salzmann<sup>15</sup>, M. Sasaki<sup>20</sup>, F. Schüssler<sup>18</sup>, H. M. Schutte<sup>8</sup>, J. N. S. Shapopi<sup>7</sup>, A. Specovius<sup>20</sup>, S. Spencer<sup>20</sup>, Ł. Stawarz<sup>31</sup>, R. Steenkamp<sup>7</sup>, S. Steinmassl<sup>2</sup>, C. Steppa<sup>11</sup>, K. Streil<sup>20</sup>, I. Sushch<sup>8</sup>, H. Suzuki<sup>37</sup>, T. Takahashi<sup>38</sup>, T. Tanaka<sup>37</sup>, R. Terrier<sup>12</sup>, M. Tluczykont<sup>9</sup>, N. Tsuji<sup>39</sup>, T. Unbehaun<sup>20,*</sup>, C. van Eldik<sup>20</sup>, M. Vecchi<sup>5</sup>, J. Veh<sup>20</sup>, C. Venter<sup>8</sup>, J. Vink<sup>35</sup>, T. Wach<sup>20</sup>, S. J. Wagner<sup>4</sup>, A. Wierzcholska<sup>23</sup>, M. Zacharias<sup>4,8</sup>, D. Zargaryan<sup>1</sup>, A. A. Zdziarski<sup>33</sup>, A. Zech<sup>16</sup>, S. Zouari<sup>12</sup>, N. Żywucka<sup>8</sup> (H.E.S.S. Collaboration), and A. Harding<sup>40</sup>

(Affiliations can be found after the references)

Received 17 November 2023 / Accepted 18 March 2024

ABSTRACT

The Crab Nebula is a unique laboratory for studying the acceleration of electrons and positrons through their non-thermal radiation. Observations of very-high-energy γ rays from the Crab Nebula have provided important constraints for modelling its broadband emission. We present the first fully self-consistent analysis of the Crab Nebula’s γ-ray emission between 1 GeV and ~100 TeV, that is, over five orders of magnitude in energy. Using the open-source software package GAMMAPY, we combined 11.4 yr of data from the *Fermi* Large Area Telescope and 80 h of High Energy Stereoscopic System (H.E.S.S.) data at the event level and provide a measurement of the spatial extension of the nebula and its energy spectrum. We find evidence for a shrinking of the nebula with increasing γ-ray energy. Furthermore, we fitted several phenomenological models to the measured data, finding that none of them can fully describe the spatial extension and the spectral energy distribution at the same time. Especially the extension measured at TeV energies appears too large when compared to the X-ray emission. Our measurements probe the structure of the magnetic field between the pulsar wind termination shock and the dust torus, and we conclude that the magnetic field strength decreases with increasing distance from the pulsar. We complement our study with a careful assessment of systematic uncertainties.

**Key words.** acceleration of particles – radiation mechanisms: non-thermal – gamma rays: general – ISM: individual objects: Crab Nebula

**1. Introduction**

The Crab Nebula, associated with the pulsar PSR B0531+21, represents the prototype of a pulsar wind nebula and has been extensively studied across all wavelengths of the electromagnetic spectrum (e.g. Hester 2008; Buehler & Blandford 2014; Amato & Olmi 2021). In the high-energy (HE; $100 \text{ MeV} < E <$

$100 \text{ GeV}$) and very-high-energy (VHE; $E > 100 \text{ GeV}$) wavebands, it stands out as one of the brightest objects in the sky. In fact, the Crab Nebula was the first γ-ray source to be unambiguously detected at very high energies from 1986–1988 with the Whipple telescope (Weekes et al. 1989). The radiation emitted by the nebula is mostly steady, although temporal variations have been observed in some wavebands (e.g. Wilson-Hodge et al. 2011; Abdo et al. 2011). In the VHE band, no flux variability has been observed so far, and the Crab Nebula has been used

<footnote>* Corresponding authors: H.E.S.S. Collaboration, e-mail: contact.hess@hess-experiment.eu.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** SED of the Crab Nebula. The synchrotron data points up to hard X-rays have been taken from Dirson & Horns (2023), and references therein. The HE and VHE data are from Buehler et al. (2012), Aharonian et al. (2004, 2006), Aleksić et al. (2015), Acciari et al. (2020), Meagher (2015), Amenomori et al. (2019), Abeysekara et al. (2019), and Cao et al. (2021), in the order as listed in the legend. A distance of 2 kpc has been assumed to convert to luminosity (see right-hand side axis).</figcaption></figure>

as a standard candle for detector calibration in this regime since its first detection (cf. Vacanti et al. 1991; Pühlhofer et al. 2003; Aharonian et al. 2006; Albert et al. 2008; Meyer et al. 2010; Abeysekara et al. 2017; Aharonian & An 2021; Abe et al. 2023).

We show in Fig. 1 the spectral energy distribution (SED) of the Crab Nebula, as measured by a large number of instruments. The distribution exhibits two broad peaks; both are generally accepted to be caused by a population of high-energy electrons<fnref n="1" /> that fills the nebula. The first peak is due to synchrotron emission that the electrons emit within the magnetic field of the nebula. The emission at energies close to the second peak is generated by high-energy electrons via the inverse Compton (IC) process. The general formation and structure of the nebula was already discussed by Rees & Gunn (1974), and a number of theoretical models describing the emission from high-energy electrons in the Crab Nebula have been put forward since then (e.g. Kennel & Coroniti 1984a; de Jager & Harding 1992; Atoyan & Aharonian 1996; Hillas et al. 1998; Aharonian et al. 2004; Meyer et al. 2010), including time-dependent models applicable to young pulsar wind nebulae (e.g. Torres et al. 2014; van Rensburg et al. 2018; for a review see also Amato & Olmi 2021, and references therein). Although usually considered sub-dominant over most of the energy range, γ-ray emission due to accelerated cosmic-ray nuclei has also been considered (e.g. Atoyan & Aharonian 1996; Arons 1998; Bednarek 2003). Recently, following the detection of PeV-energy photons from the Crab Nebula by the LHAASO detector (Cao et al. 2021), this possibility has received some renewed interest (Peng et al. 2022; Nie et al. 2022).

In recent years, it has furthermore become possible to spatially resolve the IC component of the Crab Nebula both in the HE and VHE bands, with extensions (68% containment radii) of $1.8' \pm 0.18'_{\text{stat}} \pm 0.42'_{\text{sys}}$ measured with the Fermi Large Area Telescope (LAT; Ackermann et al. 2018) and $1.31' \pm 0.07'_{\text{stat}} \pm 0.17'_{\text{sys}}$ measured with the High Energy Stereoscopic System (H.E.S.S.; Abdalla et al. 2020), respectively. These measurements provide new, crucial constraints for the theoretical modelling of the IC emission from the Crab Nebula. Generally, because higher-energy electrons cool more efficiently, the extension of the nebula is expected to decrease with increasing electron energy (Kennel & Coroniti 1984a; de Jager & Harding 1992; Atoyan & Aharonian 1996; Hillas et al. 1998; Aharonian et al. 2004). For example, in the X-ray domain – where the emission is due to synchrotron radiation of the highest-energy electrons in the system – the nebula exhibits a torus-like morphology that is only about $0.63' \times 0.3'$ in size (Weisskopf et al. 2000). In contrast, in the ultraviolet, the characteristic extension is approximately $1.6'$ (Abdalla et al. 2020, where we have converted from a containment fraction of 39–68% assuming a Gaussian morphology). However, besides the spatial distribution of the electrons, the apparent size of the nebula also depends on the extension of the nebula’s magnetic field (for synchrotron radiation) and its distribution of seed photons (for the IC component). At radio wavelengths the extension is indeed only ~1.9’, which is smaller than expected solely from the distribution of electrons (Yeung & Horns 2019). Recently, Dirson & Horns (2023) performed a combined spatial and spectral fit of the multi-wavelength (MWL) data of the Crab Nebula, using a model that features radially dependent electron, seed photon, and magnetic field densities, finding that they can describe the entire data set well.

In this work, we provide, for the first time, a simultaneous measurement of the energy spectrum and spatial extension of

<footnote n="1">¹ Throughout this paper, the term “electrons” refers to both electrons and positrons.</footnote>
the entire IC component of the Crab Nebula. To this end, we combined data from Fermi-LAT and H.E.S.S. at the level of detected γ-ray events in a joint likelihood fit, using the open-source software package GAMMAPY (Donath et al. 2023). In addition, we fitted phenomenological models based on those developed by Kennel & Coroniti (1984a), Meyer et al. (2010), and Dirson & Horns (2023) to the Fermi-LAT, H.E.S.S., and available synchrotron MWL data. We emphasise that in contrast to previous works, our analysis does not rely on fitting models to derived data products such as flux points, which depend on the assumed underlying spectral or spatial models used in their derivation.

The remainder of the paper is structured as follows. In Sect. 2, we first present separate analyses of the H.E.S.S. and Fermi-LAT data. The joint analysis detailed in Sect. 3 then provides a measurement of the energy spectrum and extension of the Crab Nebula across the energy range covered by both instruments. In Sect. 4, we describe the fit of several phenomenological models to our combined data set. The impact of systematic uncertainties is discussed in Sect. 5. Finally, we conclude in Sect. 6.

## 2. Data analysis

In parallel with combining the Fermi-LAT and H.E.S.S. data we performed separate analyses using the standard analysis tools of the respective instruments. These results are used to cross-check the joint analysis, which we carried out with the GAMMAPY package (v0.18.2; Donath et al. 2023; Deil et al. 2020).

### 2.1. Likelihood method

The analysis is performed as a three-dimensional (longitude, latitude, and energy), binned likelihood fit. The measured events are binned into a “counts” cube according to their reconstructed energy and sky direction. The spectral and spatial models are evaluated on a similar cube, which can have a different binning or larger dimensions for better accuracy. The model prediction is then forward-folded with the instrument response functions (IRFs: effective area, point spread function (PSF), and energy dispersion) in order to compute the number of predicted counts for each cell of the counts cube. The likelihood $\mathcal{L}$ is a measure of how well the model prediction agrees with the measurement and can be calculated as

$$ \mathcal{L}(\xi) = \prod_{i} \mathcal{P}(n_i|v_i(\xi)) = \prod_{i} \frac{\nu_i^{n_i}(\xi)}{n_i!} \cdot \exp(-v_i(\xi)), \quad (1) $$

with $\mathcal{P}(n_i|v_i(\xi))$ the Poisson probability to measure $n_i$ events in cell $i$ given the model prediction $v_i$ for a set of parameters $\xi$. Besides $\gamma$ rays from the source under investigation, there are also contributions to $n_i$ and $v_i$ from background processes, as will be described in more detail in the subsequent sections concerning the different instruments. By minimising the “Cash statistic” (Cash 1979), defined as

$$ C = -2 \ln \mathcal{L} = 2 \sum_{i} (v_i - n_i \ln(v_i) + \ln(n_i!)) \quad (2) $$

(where the term $\ln(n_i!)$, being independent of the model parameters $\xi$, is usually omitted), the fit maximises the likelihood and thus the agreement between model and data.

Because the synchrotron emission of the Crab Nebula measured in the X-ray regime is steady to within ±5% over the

period considered here (Jourdain & Roques 2020), we do not expect variations of the flux level or spatial appearance of the IC emission. We therefore do not take into account the time of data taking in our analysis, but perform a time-independent modelling.

### 2.2. H.E.S.S. data analysis

H.E.S.S. is an array of five imaging atmospheric Cherenkov telescopes (IACTs) located in Namibia at an altitude of 1800 m (Aharonian et al. 2006). It is sensitive to $\gamma$ rays in the energy range from several tens of GeV to ~100 TeV, where the exact achievable energy threshold depends on the observing conditions (in particular the zenith angle). The initial array, denoted H.E.S.S. Phase-I, was completed in 2003 and comprises four telescopes (CT1-4) with 108 m² mirror area each, arranged in a square grid. We refer to data taken with this configuration as “stereo” data. In 2012, a fifth, larger telescope (CT5) with 614 m² mirror area was added in the centre of the array, marking the beginning of H.E.S.S. Phase-II. The larger mirror area of CT5 can in favourable conditions lead to energy thresholds significantly below 100 GeV (e.g. Abdalla et al. 2017, 2018). When analysed separately, data taken with the CT5 telescope is referred to as “mono” data (Murach et al. 2015). In 2019, the camera of CT5 was replaced by a “FlashCam” prototype camera (Bi et al. 2021; Pühlhofer et al. 2021) that has been designed for the upcoming Cherenkov Telescope Array (CTA; Acharya et al. 2018). H.E.S.S. data taking is conducted in observation runs that typically last 28 min.

IACTs are built to measure the Cherenkov light emitted in $\gamma$-ray induced particle showers in the atmosphere. Based on the recorded shower images, we reconstructed the direction and energy of the primary $\gamma$ ray using the IMPACT reconstruction algorithm (Parsons & Hinton 2014). The dominant background consists of air showers initiated by charged cosmic rays, which we suppress by means of a multi-variate analysis (Ohm et al. 2009). In contrast to traditional methods that estimate the residual background from “OFF” regions in the same observation run, we used here a three-dimensional background model constructed from archival H.E.S.S. observations (this is a prerequisite for carrying out the three-dimensional likelihood analysis described at the beginning of this section; Mohrmann et al. 2019). The IRFs are obtained from extensive Monte-Carlo simulations, carried out with the SIM_TELARRAY package (Bernlöhr 2008). A correction for telescope efficiency, measured for each observation run using muons from atmospheric air showers, is applied as described in Mitchell (2016). Results have been cross-checked with independent calibration and reconstruction methods (de Naurois & Rolland 2009).

We used 114 observation runs (corresponding to an observation time of ≈50 h) of H.E.S.S. stereo (CT1-4) data taken between Nov. 2004 and Mar. 2015. Only observations with good data quality, excellent atmospheric conditions, zenith angles below 55°, and a maximum angle between the pointing position and the Crab Nebula of 1° have been included. Compared to the standard analysis configuration, we required a minimum of three telescopes for the reconstruction of each event, which leads to fewer detected events mostly at low energies but improves the angular resolution by ~20%, to 0.06° (68% PSF containment radius) at 1 TeV. At the same energy, the achieved energy resolution is 16% (standard deviation). Because the Crab Nebula dominates at a relatively large zenith angle of ~45° at the H.E.S.S. site, the energy threshold of the stereo data set is comparatively high (≈560 GeV).
<table><caption>**Table 1.** Overview of H.E.S.S. data.</caption><thead><tr><th>Data set</th><th>Year</th><th>Runs</th><th>Observation time (h)</th></tr></thead><tbody><tr><td rowspan="8">Stereo</td><td>2004</td><td>18</td><td>7.82</td></tr><tr><td>2005</td><td>6</td><td>2.81</td></tr><tr><td>2007</td><td>5</td><td>2.34</td></tr><tr><td>2008</td><td>4</td><td>1.43</td></tr><tr><td>2009</td><td>10</td><td>4.69</td></tr><tr><td>2013</td><td>58</td><td>25.19</td></tr><tr><td>2014</td><td>7</td><td>3.01</td></tr><tr><td>2015</td><td>6</td><td>2.80</td></tr><tr><td rowspan="3">Mono</td><td>2019</td><td>18</td><td>8.41</td></tr><tr><td>2020</td><td>42</td><td>19.44</td></tr><tr><td>2021</td><td>5</td><td>2.09</td></tr></tbody></table>

**Notes.** The separation into years is done only for illustrative purposes here; the data are not separated by year when they enter the analysis.

Additionally, we used 65 observation runs (≈30 h) of mono data taken with the CT5 telescope with its new camera, between Nov. 2019 and Oct. 2021. This data set has a lower energy threshold (≈240 GeV) compared to the stereo one, which increases the overlap with the energy range covered by the Fermi-LAT (see next section). Because the reconstruction of events uses information from only one telescope in this case, the performance in terms of angular resolution (0.13°) and energy resolution (29%) is worse compared to the stereo data set. A summary of the H.E.S.S. data sets can be found in Table 1.

The region of interest (ROI) for the H.E.S.S. data analysis is defined as a region of 5° × 5° around the position of the Crab Nebula, in Galactic coordinates. Because the extension of the Crab Nebula is at the limit of detectability for H.E.S.S., we chose a relatively fine spatial bin size of 0.01°, to achieve an accurate description of the PSF and to not smear the signal artificially. Furthermore, we binned the data in energy using 16 bins per decade. For better performance, we only considered events with a maximum offset of 1.5° from the pointing position in each observation run.

In a first step, we fitted the normalisation and spectral tilt<fnref n="2" /> of the background model template to each observation individually, where we excluded a circular region with radius 0.3° around the Crab Nebula. For the stereo data, we then binned the observations by zenith angle into three bins with edges [45, 47.5, 52.5, 55]°, so that the observations in each bin have nearly identical energy thresholds. The mono observations exhibit a smaller variation in zenith, so that this procedure is not necessary. For each bin, as well as for the mono data, we combined the respective observations to a stacked data set. In the stacking process, the counts, exposure, and predicted background are summed while the energy dispersion and point spread function are averaged. The stacking makes the analysis computationally less expensive at the cost of using slightly less accurate IRFs. The binning of observations in zenith angle prevents the averaging of IRFs across observation runs with too dissimilar observing conditions. We additionally smoothed the runs-averaged PSF with a 30'' one-dimensional Gaussian kernel, which represents the scale of pointing uncertainty of the telescopes (Braun 2007), where we assume that the pointing error varies from run to run. After this procedure, we are left with

three stereo data sets (for low, medium, and high zenith angles) as well as one mono data set.

Subsequently, we modelled the γ-ray emission of the Crab Nebula, that is, we performed a three-dimensional likelihood analysis as introduced in Sect. 2.1 to simultaneously fit the spatial morphology and the energy spectrum. We used a symmetric, two-dimensional Gaussian as spatial model for the intrinsic source distribution, which is then folded with the PSF. At this stage, we did not yet consider a possible energy dependence of the spatial model. All extension measurements reported in this work denote the 68% containment radius of the emission, which for a two-dimensional Gaussian with width σ is given by $r_{68} = \sqrt{-2 \cdot \ln(1 - 0.68)} \sigma \approx 1.5 \cdot \sigma$. From the stereo data set, we obtain a fitted extension (averaged over energy) of 1.62' ± 0.05'_{stat} ± 0.21'_{sys}. The systematic uncertainty is derived by repeating the analysis with a PSF that has been broadened (narrowed) by 5% with respect to the default one, independent of energy. The variation of 5% in width has been established as a realistic value in a study – with the same analysis configuration as employed here – of the bright blazar PKS 2155–304 (see e.g. Abramowski et al. 2010), which appears as a point-like source for the H.E.S.S. array. It amply covers the 30''-smoothing applied to the PSF earlier. Finally, we note that given the larger systematic uncertainty and considerably larger PSF of the mono data, we cannot constrain the extension with this data set.

The obtained value is larger than the one found in Abdalla et al. (2020), $r_{68} = 1.30' \pm 0.07'_{stat} \pm 0.17'_{sys}$, although we note that the deviation is far from significant when taking into account the systematic uncertainties. We attribute the difference partly to the fact that at the lower energy threshold achieved here (560 GeV instead of 700 GeV), the extension is indeed expected to be slightly larger. Nevertheless, we note that this comparison of results indicates that the systematic uncertainty may be slightly underestimated. The measured extension is compatible with the upper limits provided by the HEGRA telescopes (Aharonian et al. 2000, 2004).

For the spectral model, we used a log-parabola function,

$$ \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{-\alpha - \beta \ln(E/E_0)} \quad (3) $$

with normalisation $N_0$, reference energy $E_0$, spectral index $\alpha$, and curvature parameter $\beta$. From a joint fit to all H.E.S.S. data sets, we obtain $N_0 = (4.06 \pm 0.03) \times 10^{-11} \text{ TeV}^{-1} \text{ cm}^{-2} \text{ s}^{-1}$ at $E_0 = 1 \text{ TeV (fixed)}$, $\alpha = 2.530 \pm 0.006$, and $\beta = 0.086 \pm 0.005$; this result is shown as a black line in Fig. 2. For both the mono and stereo data, we subsequently derived separate sets of flux points by first performing the fit separately for each data set, and then re-adjusting the normalisation parameter in narrow energy ranges (where we have ensured that the best-fit models for the two sets are compatible with each other). Taking into account the relatively poor energy resolution of the mono data, we reduced the number of flux points (from 16 to 8 per decade in energy) for this data set in order to avoid correlations between the points. Both the stereo and mono flux points are also displayed in Fig. 2. The figure illustrates that the mono data set allows us to extend the measurement to lower energies, whereas the stereo data set extends up to almost ~100 TeV. A comparison with results from other instruments shows that our derived flux is largely compatible, with the noticeable deviations most likely explained by a systematic shift in the energy scale of the instruments. The spectral results have been confirmed in several one-dimensional analyses based on the more traditional reflected-background method (Fomin et al. 1994; Berge et al. 2007). One example of such a

<footnote n="2">² The spectral tilt parameter $\delta$ modifies the predicted background rate $R$ at energy $E$ as $R' = R \cdot (E/1 \text{ TeV})^{-\delta}$.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** H.E.S.S. flux points of the mono (CT5) analysis (red) and stereo (CT1-4) analysis (blue) in comparison to results from other instruments. The error bars show the 68% statistical uncertainty. Points less significant than $2\sigma$ are shown as 95% confidence level upper limits. The black line displays the result of fitting a log-parabola model to the H.E.S.S. data of this work. References (in the order of the legend): Aharonian et al. (2004, 2006), Aleksić et al. (2015), Meagher (2015), Amenomori et al. (2019), Abeysekara et al. (2019), and Cao et al. (2021).</figcaption></figure>

comparison is shown in Appendix A. Finally, we note that in contrast to the Fermi-LAT analysis (see next section), the γ-ray emission from the Crab Pulsar itself, despite extending into the TeV energy range (Ansoldi et al. 2016), can safely be neglected for our study. Adopting the measured spectrum of the pulsed emission from Ansoldi et al. (2016), even the contamination for the lowest-energy flux point derived with the H.E.S.S. mono data set is below 0.5%.

## 2.3. Fermi-LAT data analysis

The Fermi-LAT is a pair conversion telescope designed to measure γ rays with energies from 20 MeV to above 300 GeV (Atwood et al. 2009). We first performed a pulsar phase-resolved reference analysis of the Crab Nebula using the standard FERMITOOLS version 1.2.23³ and FERMIPY, version 0.19.0+14⁴ (Wood et al. 2017). This analysis serves as a comparison for the GAMMAPY analysis.

To this end, we selected γ rays that were observed between 4 August 2008 and 4 January 2020 in the energy range between 100 MeV and 3 TeV and within a $10^\circ \times 10^\circ$ ROI centred on the nominal Crab Nebula position provided in the fourth Fermi-LAT γ-ray source catalogue (4FGL; Abdollahi et al. 2020). As we are interested in emission of the nebula and not the pulsar, we proceeded to determine the off-pulse phase interval. To do so, we used an ephemeris prepared following the methods described in Ajello et al. (2022) and generated a histogram of the pulsar phases $\phi$ of all photons that have arrived within a $1^\circ$ radius around the Crab Nebula, see Fig. 3. The off-pulse phase inter-

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Distribution of pulsar phases of γ-ray candidate events detected with Fermi-LAT around the Crab Nebula. The distribution includes events with a reconstructed arrival direction that is within $1^\circ$ of the Crab Nebula's position. The shaded regions indicate the intervals for the Bridge emission ($0.09 \le \phi \le 0.22$, not relevant here) and Off region ($0.51 \le \phi \le 0.89$) derived with the QB algorithm. The grey shaded areas indicate regions found by the HOP algorithm, which identifies “water-sheds” between local maxima (Meyer et al. 2019). There are shown here for illustration only.</figcaption></figure>

val is determined with the “quiescent background” (QB) algorithm introduced by Meyer et al. (2019). The resulting interval is $0.51 \le \phi \le 0.89$, as indicated by the red-shaded region in Fig. 3.

In the 4FGL, the total emission from the Crab Nebula is modelled with two super-imposed source models for the high-energy tail of the synchrotron emission and the emission due to IC scattering, respectively. For the combination with H.E.S.S. data, we are mostly interested in the IC part of the spectrum. As the synchrotron and IC components are of similar intensity around 1 GeV, we proceeded in an iterative fashion: we first optimised an ROI model over the entire energy range. We then fixed the synchrotron component to the best-fit model and performed a second analysis of the nebula above 1 GeV, which is dominated by IC radiation. For this second analysis, we also used a finer spatial binning (see below) in order to be able to determine the extension of the IC nebula.

For the initial ROI model, we selected photons within the off-pulse interval and we further excluded events that have arrived at a zenith angle above $90^\circ$ in order to mitigate contamination of γ rays originating from Earth's limb. We also excised periods of bright γ-ray bursts and solar flares that have been detected with a test statistic (TS) > $100^5$. We used the Pass 8 IRFs (Atwood et al. 2013) and selected γ rays passing the P8R3 SOURCE event selection. We chose a spatial binning of $0.04^\circ$ pixel$^{-1}$ and eight energy bins per decade.

We initialised the ROI model with all point-like and extended γ-ray sources within $20^\circ$ of the ROI centre as listed in the 4FGL (except the Crab Pulsar), from which we also took the initial spectral source parameters. Additionally, the standard tem-

<footnote n="3">³ http://fermi.gsfc.nasa.gov/ssc/data/analysis/software</footnote>

<footnote n="4">⁴ http://fermipy.readthedocs.io</footnote>

<footnote n="5">⁵ The TS is defined as $TS = -2 \ln(\mathcal{L}_1/\mathcal{L}_0)$, i.e. the log-likelihood ratio between the maximised likelihoods $\mathcal{L}_1$ and $\mathcal{L}_0$ for the hypotheses with and without an additional source, respectively (Mattox et al. 1996).</footnote>
plates for isotropic and Galactic diffuse emission are used<fnref n="6" />. In the Fermi-LAT energy range, the synchrotron component is expected to exhibit only a very small spatial extension that is not resolvable with Fermi-LAT, and therefore modelled as a point-like source. For the IC component, we used the 4FGL default source morphology, which is given by a radial Gaussian with a width of $σ = 0.0198°$ ($r_{68} = 1.79°$). As for the spectrum, the synchrotron part is modelled with a simple power law,

$$ \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{-\Gamma}, \qquad (4) $$

with normalisation $N_0$, reference energy $E_0$, and spectral index $\Gamma$, whereas the IC component is described with a log-parabola as defined in Eq. (3). After an initial optimisation, we freed the spectral normalisation and spectral shape parameters of sources within 10° from the ROI centre as well as all spectral parameters of the isotropic and diffuse emission templates. We froze all spectral parameters for sources with TS < 5. After successful convergence of the fit, we generated a TS map to search for additional point sources. For each pixel in the ROI, we added a putative point source with a power-law spectrum with index $\Gamma = 2$ and calculated its TS. No additional sources with $\sqrt{TS} \ge 5$ were found. The best-fit spectrum of the synchrotron component is shown in Fig. 4 and the parameters are provided in Table 2.

The best-fit ROI model was then used to perform an analysis above 1 GeV. In terms of data selection, we followed the official Fermi-LAT recommendations and relaxed the zenith angle cut from 90° to 100°. We shrank the ROI to 6° × 6° and following Ackermann et al. (2018), we chose a finer spatial binning of 0.025° pixel⁻¹ that we apply in the calculation of the livetime and exposure cubes. Furthermore, we made use of PSF event types (Atwood et al. 2009). The data are split into four sets according to the quality of the reconstruction of the γ-ray arrival directions. Each set is analysed with its own specific IRFs and the four sets are combined by multiplying their respective likelihood values. Accordingly, the isotropic background templates for these event types have been used. We also limited the maximum energy of the analysis to the edge of the energy bin that contains the highest energy photon, that is 1.78 TeV. The spectral parameters of the synchrotron component have been fixed. The resulting spectrum is presented in Fig. 4 and the best-fit spectral parameters are listed in Table 2. They agree well with the spectral parameters of the fit over the whole energy range.

Finally, we obtained the extension of the nebula using a two-dimensional Gaussian model as before. In this procedure, the nebula's sky coordinates are additional free parameters, as are the spectral normalisation parameters of the diffuse background sources. The spectral parameters of the nebula itself are frozen to their best-fit values. We obtain a 68% containment radius of $r_{68} = 2.22' \pm 0.18'_{-0.48}^{+0.54'}_{\text{sys}}$, where the systematic uncertainties are derived by using a bracketing method for the PSF (Ackermann et al. 2018). We find a TS value for the extension of TS$_{ext}$ = 47.6 with the nominal PSF, which reduces to TS$_{ext}$ = 16.3 if the larger bracketing PSF is used instead. These results are compatible with the extension found by Ackermann et al. (2018), $r_{68} = 1.80' \pm 0.18'_{-0.42}^{+0.42'}_{\text{sys}}$, within the uncertainties.

<figure><img src="image_6.png" /><figcaption>Fig. 4. Fermi-LAT spectra of synchrotron and IC emission. The IC emission is derived with both FERMIPY and GAMMAPY and the spectra agree remarkably well. Details on the GAMMAPY analysis of Fermi-LAT data are given in Sect. 3.1. The error bars on the flux points show the 68% statistical uncertainties.</figcaption></figure>

### 3. Combined analysis

In contrast to previous analyses, in which final data products such as flux points or extension measurements from different instruments have been analysed together, here we combine the Fermi-LAT and H.E.S.S. data at the event level. The advantage of this approach is that we are able to exploit more information and ensure consistency of spectral and morphological models between the different instruments. In the combined analysis we simultaneously fitted one three-dimensional source model (specrum and morphology) to the (binned) event data of Fermi-LAT and H.E.S.S. The analysis is carried out with GAMMAPY, where we used the data sets which we generated and cross-checked for the individual analyses and summed their respective Cash statistic values in the combined analysis. The total C value is then minimised, which corresponds to a joint fit of the model to the Fermi-LAT and H.E.S.S. data sets.

#### 3.1. Preparation of Fermi-LAT data in GAMMAPY

During the FERMIPY analysis the events are binned into a counts cube, and the IRFs (exposure, energy dispersion, and PSF) are calculated. We used the data products generated during the Fermi-LAT analysis to set up the combined analysis within GAMMAPY, where we generated four separate Fermi-LAT data sets, corresponding to the four PSF classes. While the analyses with FERMIPY and GAMMAPY are in principle very similar, there are some differences which require adaptations in the data production steps. These are described in the following.

First, GAMMAPY only calculates model fluxes of sources within the ROI, whereas the standard Fermi-LAT analysis tools allow the inclusion of source models outside the ROI, which may contribute to the observed counts in the ROI especially at low energies, where the Fermi-LAT PSF is very broad. Consequently, the ROI size of the Fermi-LAT data sets has been

<footnote n="6">⁶ We used the file gll_iem_v07.fits or the Galactic diffuse emission and the file iso_P8R2_SOURCE_V6_v07.txt for the isotropic diffuse component; see: http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html</footnote>
<table><caption>**Table 2.** Best-fit parameters of synchrotron and IC emission models of the Crab Nebula derived from Fermi-LAT data.</caption><thead><tr><th></th><th>N<sub>0</sub><br>(10<sup>-7</sup> MeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>E<sub>0</sub><br>(MeV)</th><th>Γ or α</th><th>β</th></tr></thead><tbody><tr><td>Synchrotron</td><td>2.47 ± 0.06</td><td>50.6</td><td>3.638 ± 0.024</td><td>—</td></tr><tr><td>IC (full energy range)</td><td>(4.75 ± 0.10) × 10<sup>-6</sup></td><td>10<sup>4</sup></td><td>1.762 ± 0.012</td><td>0.053 ± 0.006</td></tr><tr><td>IC (E ≥ 1 GeV)</td><td>(4.93 ± 0.11) × 10<sup>-6</sup></td><td>10<sup>4</sup></td><td>1.768 ± 0.013</td><td>0.054 ± 0.008</td></tr><tr><td>IC (E ≥ 1 GeV, GAMMAPY)</td><td>(4.88 ± 0.11) × 10<sup>-6</sup></td><td>10<sup>4</sup></td><td>1.763 ± 0.013</td><td>0.056 ± 0.008</td></tr></tbody></table>

increased to 10° × 10°, such that the 4FGL sources outside the inner (6° × 6°) analysis region are also evaluated by GAMMAPY. Additionally, to correctly take into account the energy dispersion at the boundaries of the analysed “reconstructed” energy range with GAMMAPY, we need to be able to evaluate the IRFs at “true” γ-ray energies that extend beyond this range. As it is not possible to separately define the ranges to be considered for reconstructed and true energy in FERMIPY, we increased the energy range by two bins at the low-energy end and by one bin at the high-energy end. The added spatial and energy bins were then again masked in the GAMMAPY analysis, so that exactly the same regions and energy bins as in the FERMIPY analysis contribute to the likelihood calculation.

To verify that our analysis setup in GAMMAPY is correct, we repeated the analysis of the Fermi-LAT data previously carried out with FERMIPY (cf. Sect. 2.3), using exactly the same model components. In the fit, the model parameters of the Crab Nebula’s IC emission as well as normalisation of the isotropic diffuse and Galactic diffuse backgrounds and spectral index of the Galactic diffuse model are left free, while all parameters of the other source models in the ROI are fixed to the best-fit values from the FERMIPY analysis. With this procedure, we obtain an excellent agreement between the FERMIPY and the GAMMAPY analysis. For the spectrum, this is illustrated in Fig. 4 and summarised in Table 2. For the extension, the GAMMAPY fit yields $r_{68} = (2.10 \pm 0.18_{\text{stat}})$, which is consistent with the FERMIPY result given at the end of Sect. 2.3.

## 3.2. Measuring the Crab Nebula’s combined energy spectrum and extension

Through the combined fit, we can now consistently measure the extension of the Crab Nebula across the combined energy range of Fermi-LAT and H.E.S.S.. To this end, we fitted our models jointly to the Fermi-LAT and H.E.S.S. data sets as defined in the previous sections. As in the separate analyses, we used a two-dimensional Gaussian as the spatial model. Because the log-parabola model does not yield a satisfactory description of the IC component across the entire energy range covered by our combined fit, we used a smoothly broken power-law model for the spectral component:

$$ \frac{dN}{dE} = N_0 \left(\frac{E}{E_0}\right)^{-\Gamma_1} \left[ 1 + \left(\frac{E}{E_{\text{break}}}\right)^{\frac{\Gamma_2-\Gamma_1}{\beta}} \right]^{-\beta}, \quad (5) $$

with $E_0$ fixed at 1 TeV. A fit over the whole energy range yields the following parameters: $N_0 = (3.35 \pm 0.22) \times 10^{-10} \text{ TeV}^{-1} \text{ s}^{-1} \text{cm}^{-2}$, $\Gamma_1 = 1.61 \pm 0.02$, $\Gamma_2 = 2.95 \pm 0.02$, $E_{\text{break}} = 0.33 \pm 0.02$ TeV, and $\beta = 1.73 \pm 0.07$. This model is preferred over the log-parabola model by a difference of 94 in the Akaike information criterion, $AIC = C + 2N_{\text{par}}$ (Akaike 1974), where the number of free model parameters is $N_{\text{par}} = 3$ for the

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Comparison of log-parabola and smoothly broken power law spectral models. Both are fitted to the Fermi-LAT and H.E.S.S. data sets, while the LHAASO flux points are not included in the fit and only shown for comparison.</figcaption></figure>

log-parabola model and $N_{\text{par}} = 5$ for the smoothly broken power-law model. Both models are shown in Fig. 5, together with the Fermi-LAT and H.E.S.S. flux points derived in this work. A comparison of the models with the LHAASO flux points, also shown in the plot, indicates that an extension of the analysis beyond 100 TeV would require the addition of another break to the spectral model.

The best-fit position over the whole IC energy range is $(184.5499 \pm 0.0004_{\text{stat}})^{\circ}$ Galactic longitude and $(-5.7825 \pm 0.0004_{\text{stat}})^{\circ}$ Galactic latitude. This position is shifted to the north with respect to that of the Crab Pulsar by $\approx 28''$. In what follows, the only parameter of interest is the extension parameter of the Gaussian spatial model. However, in order to take into account possible correlations with other parameters, the following parameters are also left free: amplitude and position of the Crab Nebula source model, as well as the normalisation parameters of the Fermi-LAT isotropic and Galactic diffuse background models and the H.E.S.S. background models. The rest of the spectral parameters as well as the rest of the Fermi-LAT ROI model remain fixed. To obtain an energy-dependent measurement of the extension, we optimised the free parameters separately in different energy intervals, where we used a binning of two bins per decade in the observed energy range between 1 GeV and 100 TeV. As an important cross-check, we also compared the best-fit position of the Crab Nebula spatial model between all energy bands, finding that they are all compatible with each other within the uncertainties. Furthermore, the fitted extensions remain unchanged if the centre of the spatial model is fixed to a common position for all energy bands.

As the PSF of the H.E.S.S. mono data set is not understood well enough, we did not include this set in the measurement
<figure><img src="image_8.png" /></figure>

**Fig. 6.** 68% containment radius of the γ-ray emission from the Crab Nebula measured with the combined Fermi-LAT and H.E.S.S. data analysis. The error bars denote the 1σ statistical uncertainties and do not include systematic effects. For comparison, the red points show the extension measured with Fermi-LAT by Yeung & Horns (2019, 2021), whereas the purple point displays the extension previously measured with H.E.S.S. (Abdalla et al. 2020). In grey we show the systematic uncertainty of the H.E.S.S. and Fermi-LAT measurements, respectively (cf. Sects. 2.2 and 2.3). The bottom panel shows the square root of the TS value of the measured extension, which gives an indication of the (statistical) significance of the extension with respect to the null hypothesis of a point-like source.

of the extension. Potentially large systematic effects may other-
wise dominate the extension measurement between 240 GeV
and 560 GeV (where the stereo data set has its threshold), due to
the large number of events detected at these energies.

In Fig. 6 we show the measured extension as a function of the γ-ray energy. Below the threshold of the H.E.S.S. stereo data set, 560 GeV, only the Fermi-LAT data contribute. To obtain a more significant extension, we have merged the first and last two bins in this range. Above 560 GeV, the H.E.S.S. data immediately dominate the fit due to the much larger effective area and better angular resolution. The highest significance for an extension is reached in the energy bin ranging from 1–3 TeV; this range also dominates the energy-integrated extension measurement with H.E.S.S. presented in Sect. 2.2.

The extensions measured in the ranges dominated by the
Fermi-LAT and H.E.S.S. data, respectively, connect smoothly
and are furthermore consistent with published results. Overall,
we find strong evidence for an extension that decreases with
energy. As we subsequently see in Sect. 4, this result is in very
good agreement with theoretical expectations.

The decreasing extension is also illustrated in Fig. 7, where
we show a measurement of the Crab Nebula in the optical regime
and at X-ray energies, and compare this to our extension mea-
surements with Fermi-LAT and H.E.S.S.. As noted earlier, the fit
of the Fermi-LAT data above 1 GeV yields an extension of r<sub>68</sub> =
(2.10 ± 0.18<sub>stat</sub>)'. The extension measured with H.E.S.S. above
10 TeV, on the other hand, is r<sub>68</sub> = (0.82 ± 0.29<sub>stat</sub>)'. In most

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

**Fig. 7.** Optical image of the Crab Nebula in green (credit: NASA, ESA/Hubble), overlaid with an X-ray image in red, taken with Chandra (0.3–10 keV; credit: NASA/CXC/SAO, https://chandra.harvard.edu/photo/openFITS/crab.html). The circles show the 68% containment radii of the γ-ray emission measured with Fermi-LAT (>1 GeV; blue) and H.E.S.S. (>10 TeV; yellow), where the bands illustrate the statistical uncertainties. The crosses indicate the fitted position and its uncertainty for the respective data sets.

models, as will be further detailed in the subsequent section, the
optical radiation is produced by electrons that are also responsi-
ble for the IC emission in the *Fermi*-LAT energy range. In agree-
ment with this, the position and extension of the emission mea-
sured with *Fermi*-LAT is roughly consistent with the outline of
the nebula seen in the optical image. As the electrons radiating
the X-ray synchrotron photons also produce the highest-energy
IC emission, a similar agreement between the X-ray image and
the H.E.S.S. measurement above 10 TeV is expected. However,
while the centroid of the γ-ray emission coincides well with a
region of bright X-ray emission, we observe that the 68% con-
tainment radius measured with H.E.S.S. almost entirely encloses
the X-ray emission, which is an unexpected result. While the
argument here is rather qualitative, this point will be addressed
again in a more quantitative way in the discussion of the mod-
elling results in Sect. 4.3.3. There, we confirm that the models
investigated in this work fail to explain the difference between
the extension in X-rays and at the highest-energy γ rays. While
the difference could also be due to underestimated systematic
effects (which will be discussed in Sect. 5), it hints at an incom-
pleteness of the models.

4. Modelling of the Crab Nebula

In this section we introduce the phenomenological models that
we investigate in this work. We present the results of fitting
them to the Fermi-LAT, H.E.S.S., and MWL data and discuss
the implications.

4.1. Description of the tested models

In total, we tested three time-independent leptonic models that
differ in the modelling of the high-energy electrons and the radial
dependence of the magnetic field. All models are radially symmetric synchrotron self-Compton (SSC) models, in which the low-energy photons are produced by synchrotron emission and the high-energy photons by IC scattering, where the synchrotron photons themselves serve as one of the targets. We assume two electron distributions: relic “radio” electrons at lower energies and freshly injected “wind” electrons at higher energies, where the latter are believed to be accelerated at a standing shock. The shock is commonly associated with a bright X-ray ring feature located at a distance of $r_s = (13.3 \pm 0.2)''$ from the pulsar (Weisskopf et al. 2012), which corresponds to $r_s \approx 0.129(2)$ pc for a distance of $d = 2$ kpc between Earth and the pulsar, which we take as the centre of the nebula in our radially symmetric models. The electron density is both a function of energy and distance from the nebula centre. The spatial distribution of the radio electrons is assumed constant in energy since the synchrotron cooling times are larger than the age of the nebula. These electrons could be the result of an injection during the early phases of the pulsar spin down (Atoyan 1999). For the wind electrons, the spatial distribution is allowed to change with energy. All models assume three seed photon fields for the IC scattering: the cosmic microwave background (CMB, 2.7 K), far-infrared emission from nearby dust (39 K and 149 K; see Dirson & Horns 2023), and the synchrotron photons.

For the description of the radial dependence of the magnetic field, we follow (1) Meyer et al. (2010) for a model where the magnetic field is constant throughout the nebula (constant B-field model); (2) Dirson & Horns (2023), where the magnetic field decreases (with adjustable rate) with increasing distance from the centre of the nebula (variable B-field model); and (3) Kennel & Coroniti (1984b) for a magneto-hydrodynamic (MHD) flow model (K&C model). The models are described in full detail in Appendix B. In the following, we refer to the former two models as static models in contrast to the K&C model, which is a result of the MHD flow equations. All fit parameters of the models are summarised in Table C.1.

In all models, the radio electron density follows a simple power law in energy between the Lorentz factors $\gamma_{r,min}$ and $\gamma_{r,max}$, where a super-exponential cutoff occurs. This spectral distribution is multiplied with a Gaussian to model the spatial dependency from the centre of the Crab Nebula. For the wind electrons, the constant and variable B-field models assume an electron density parameterised for different energies and radii, as shown in Fig. 8. Specifically, the wind electron distribution is modelled with a double-broken power law with super-exponential cut-offs at low and high energies. These two models are only distinguished by the spatial shape of the magnetic field, which is either constant or, for the variable B-field model, proportional to $r^{-\alpha}$. In the K&C model, on the other hand, the wind-electron spectrum is the result of an injection of electrons at the pulsar wind termination shock, which then evolve (together with the B-field) according to the solution of the static MHD flow equations matching the boundary condition of flow speed and position of the shock. In addition to the parameters of the electron injection spectrum, the model depends on the magnetisation $\sigma$ and spin-down luminosity (see Appendix B for details).

## 4.2. Fitting of the models

All three models give a prediction of the $\gamma$-ray intensity as a function of energy and angular separation from the Crab Pulsar, which we forward-folded with the IRFs and fitted to the Fermi-LAT and H.E.S.S. binned event data sets as introduced in the previous sections. As the centroid of the $\gamma$-ray emission is slightly

<figure><img src="image_10.png" /><figcaption>Fig. 8. Electron densities of models, with parameters as given in Table C.1. The solid, dotted and dashed black lines show the volume-averaged density for the variable B-field, the constant B-field and the K&C model, respectively. The coloured lines show the electron density for the variable B-field model at different radii. The radio electrons (at energies below the “dip”) are modelled with an energy-independent extension, therefore the density drops with radius independent of energy. On the other hand, the density of the freshly injected wind electrons decreases more rapidly for higher energies.</figcaption></figure>

offset from the pulsar position (cf. Sect. 3.2), we shifted the model prediction spatially such that its centre coincides with the centroid of the $\gamma$-ray emission in order to ensure a proper convergence of the fit. This does not affect the predicted spectrum in a noticeable way. In the combined fit, all model components except the Crab Nebula SSC model (i.e. the background models for H.E.S.S. and the Galactic and isotropic diffuse emission models as well as the remaining sources in the ROI for Fermi-LAT) are frozen to pre-fitted values. Additionally, we fixed some of the model parameters to values employed by Dirson & Horns (2023), specifically all parameters relating to the dust model and the shock radius. The fixed parameters are marked as such in Table C.1.

As the synchrotron photon field is a seed field for the IC upscattering, we also need to ensure consistency with the observed synchrotron extensions which are measured with higher precision compared to the IC extensions. Therefore, we also included archival flux and extension measurements from radio frequencies to X-ray energies, taken from Dirson & Horns (2023)<fnref n="7" />, to constrain the synchrotron flux predicted by our models (see also Fig. 1). This is done by computing a $\chi^2$-value between the predicted and observed synchrotron flux and extensions and adding this value as a penalty term to the total C-value of the fit of the IC spectrum,

$$C_{tot} = -2 \ln \mathcal{L}_{tot} = C_{IC} + \chi_{SYN,flux}^2 + \chi_{SYN,ext}^2 \quad (6)$$

where $C_{IC} = -2 \ln \mathcal{L}_{IC}$ is the Cash statistic of the model fitted to the Fermi-LAT and H.E.S.S. data. Lacking a detailed understanding of the systematic uncertainties associated with each of the synchrotron flux measurements, we added in quadrature to the uncertainty specified for each data point a relative systematic uncertainty of 7% (Meyer et al. 2010, see also Sect. 5).

<footnote n="7">⁷ The data points are available at https://github.com/dieterhorns/crab_pheno, where we have omitted points that also do not appear in Figs. 1 and 7 of Dirson & Horns (2023).</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Overview showing the best-fit variable *B*-field model and the observed data in the IC regime, in five different energy bands. Each panel displays a 1° × 1° cutout around the position of the Crab Nebula. The first column gives the best-fit model prior to folding with the IRFs. The second column illustrates the predicted number of events based on source and background models (folded with the IRFs), smoothed with a top-hat kernel of 2.4' radius. The third column shows the observed number of events with the same smoothing. Finally, the fourth column displays residual significance maps (see main text for details).</figcaption></figure>

Compared to the analysis in Dirson & Horns (2023), we used data from only two instruments to constrain the IC emission. However, both the Fermi-LAT and H.E.S.S. data sets analysed here are much larger compared to those utilised by Dirson & Horns (2023). Furthermore, we are able to include the extension of the IC nebula self-consistently in the fit because we use the three-dimensional information of the Fermi-LAT and H.E.S.S. event data (i.e. a mis-match in extension would lead to strong residuals in the fit). In order to achieve the largest possible overlap of the energy ranges covered by the Fermi-LAT and H.E.S.S. data, we included the H.E.S.S. mono data set in the combined analysis. This is acceptable despite the systematic uncertainties associated with the PSF of the mono data set because this data set contributes to the fit mostly in a relatively narrow energy range (~240–560 GeV), and we do not attempt to measure the extension as a function of energy here (in contrast to the analysis presented in Sect. 3.2).

Figure 9 shows the best-fit spatial models and residuals in five different energy bands for the variable *B*-field model. The residual maps show the deviations between the total (i.e. background + source) model prediction and the observed number of events in terms of $\sqrt{\text{TS}}$, where $\text{TS} = -2 \ln(\mathcal{L}/\mathcal{L}_{\text{max}})$ and $\mathcal{L}_{\text{max}}$ is the likelihood corresponding to a model that perfectly predicts the observed number of events in each pixel. While not perfectly compatible with purely statistical fluctuations, the residual maps indicate that we achieve a satisfactory description of the extension of the nebula in the part of the IC spectrum probed with Fermi-LAT. However, in the part probed with H.E.S.S., especially in the 1–10 TeV range, we observe stronger residuals. The extension of the Crab Nebula is strongly constrained

by the synchrotron measurements and appears too small when compared with the H.E.S.S. data, leading to negative residuals at the position of the Crab Nebula that are surrounded by an unaccounted-for excess. For the other two models this trend is similar, with even slightly larger residuals. This indicates that our models cannot simultaneously describe both the synchrotron and IC extensions together with the broadband SED (however, see Sects. 4.3.3 and 5 for a discussion of systematic uncertainties).

### 4.3. Discussion of the models

We now turn to the discussion of the physical models and how well they fit our data. For each of the models, we give a quantitative overview of the fit results in Table 3, that is, we provide the Cash statistic minimised by the fit as well as the number of fitted parameters and the number of synchrotron data points. An overview of all best-fit model parameters is provided in Table C.1. The complexity of the models prevents a detailed assessment of the uncertainties on the model parameters. We therefore note that the best-fit values should not be regarded as measurements of the corresponding physical quantities, but simply as the set of values that yields the best agreement with the data.

Of the three models tested in this work, we find that the variable *B*-field model provides the best fit to the data. A comparison between the constant *B*-field model and the variable *B*-field model shows that the former is strongly disfavoured with respect to the latter, both in the synchrotron and the IC part of the fit. Since the models are nested (i.e. the variable *B*-field model can
<table><caption>**Table 3.** Overview of *C* and $\chi^2$ values obtained in the fits of our physical models.</caption><thead><tr><th>Model</th><th>ΔC<sub>tot</sub></th><th>ΔC<sub>IC</sub></th><th>χ<sup>2</sup><sub>SYN,flux</sub> (#)</th><th>χ<sup>2</sup><sub>SYN,ext</sub> (#)</th><th>N<sub>par</sub></th><th>ΔAIC</th></tr></thead><tbody><tr><td>Variable B-field model</td><td>0</td><td>0</td><td>152.7 (194)</td><td>60.3 (14)</td><td>16</td><td>0</td></tr><tr><td>Constant B-field model</td><td>374.0</td><td>264.1</td><td>246.0 (194)</td><td>76.9 (14)</td><td>15</td><td>372.0</td></tr><tr><td>Kennel & Coroniti</td><td>292.6</td><td>143.0</td><td>285.5 (194)</td><td>77.1 (14)</td><td>12</td><td>284.6</td></tr></tbody></table>

**Notes.** The $\Delta C_{tot}$, $\Delta C_{IC}$ (cf. Eq. (6)), and $\Delta AIC$ values are given with respect to the variable B-field model. The $\chi^2$-values are given together with the number of flux and extension points included in the fit. $N_{par}$ denotes the number of free parameters of the model. See main text for a definition of the different quantities.

be reduced to the constant B-field model for a particular choice of parameter values), we can perform a likelihood ratio test and invoke Wilks' theorem (Wilks 1938) to infer that the variable B-field model is preferred by $\sqrt{\Delta C_{tot}} \approx 19.3$ standard deviations. On the other hand, a comparison of the K&C model and the variable B-field model by means of a likelihood ratio test is not possible, since these models are not nested. Instead, for a statistical comparison of the K&C model and the variable B-field model, we use again the Akaike information criterion, now defined as $AIC = C_{tot} + 2N_{par}$. We conclude that the variable B-field model is also strongly favoured over the K&C model, with a difference in AIC of 284.6.

As a caveat to these comparisons, we note that the dust model has been optimised for the variable B-field model only, as a result of adapting the dust parameters from Dirson & Horns (2023). While the other models would potentially benefit from an additional optimisation of the dust parameters, we expect this effect to be minor and to have no impact on the conclusions drawn in this study. In the following, we discuss the models in terms of their predictions for the magnetic field, the broadband SED, and the nebula extension.

### 4.3.1. Magnetic field

As already stated in Dirson & Horns (2023), the VHE data provide the strongest constraint on the shape of the magnetic field. For the variable B-field model, in which the B-field is parameterised as $B(r) = B_0(r/r_s)^{-\alpha}$, we find best-fit values of $\alpha = 0.47$ and $B_0 = 256 \mu G$. The steepness of the radial profile of the magnetic field is thus compatible with the value of $\alpha = 0.51 \pm 0.03$ found by Dirson & Horns (2023), who have used other VHE data sets in their fit.

With the K&C model we find a magnetisation of $\sigma = 0.0214$, which is about four to seven times larger compared to previous estimates (Kennel & Coroniti 1984b; Meyer et al. 2010). Such high values of $\sigma$ would lead to a high flow speed of $\ge 5000 \text{ km s}^{-1}$ at the interface between the nebula and the supernova remnant and a predicted shock radius of $r_s \approx \sqrt{\sigma r_0} \approx 0.3 \text{ pc}$, where $r_0 \approx 2 \text{ pc}$ is the size of the nebula (Kennel & Coroniti 1984b). This is at odds with observed expansion and flow speeds of the order of $\sim 2000 \text{ km s}^{-1}$ and the shock radius at $r_s = 0.13 \text{ pc}$. A value of $\sigma$ compatible with these measurements is ruled out by our analysis with high significance, as it would lead to a much larger value of $C_{tot}$ and thus be disfavoured by the AIC (cf. Sect. 4.3).

Figure 10 shows that the B-field profiles of the K&C model and the variable B-field model agree quite well in the outer regions of the nebula ($r > 5r_s$). At the shock radius, the variable B-field model predicts more than twice the magnetic field strength compared to the K&C model. We can also compute a magnetisation for the variable B-field model by comparing the energy in the B-field against the energy in particles as a function

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Magnetic field strength as predicted by the different models, as function of the distance to the nebula centre in units of the shock radius $r_s$. For the variable B-field model $B \propto r^{-\alpha}$, with $\alpha$ given in the legend. For the K&C model the magnetisation $\sigma$ is specified in the legend.</figcaption></figure>

of distance to the centre. We find a value of $\sigma \approx 0.1$ at the shock radius and $\approx 0.03$ averaged across the whole nebula.

As discussed before, the constant B-field model performs significantly worse, implying that a B field that decreases with distance from the centre is clearly preferred. This confirms the findings of Dirson & Horns (2023), who have reached a similar conclusion.

### 4.3.2. Broadband SED

We compare the synchrotron flux prediction of the models in Fig. 11. Even though the parametrisation of the radio electrons is identical for all models, minor differences can be seen in the radio part of the SED (<3 THz). Compared to the variable B-field model, both the K&C model and the constant B-field model prefer a different transition between radio and wind electrons, that is, a larger cutoff energy for the radio electrons. This leads to a slight mis-match in the SED, but is required to fit the spatial extension in the IR and optical regime (see below). Since the parametrisation of the transition is not physically motivated, we will draw no further conclusions from the differences between the models in that regime.

The X-ray to $\gamma$-ray part of the synchrotron spectrum is not described well by any of the models. Specifically, the data between $10^{20}-10^{22} \text{ Hz}$ suggest a slightly harder spectrum followed by a stronger cutoff, which is not matched by the models. Possible reasons for this could be the parametrisation of the super-exponential cutoff at $\gamma_{w,max}$, or constraints resulting from
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** SED of the Crab Nebula, together with the predictions of the three models. The dotted lines show the emission produced by the radio electrons, while dashed lines are for the wind electrons (where synchrotron emission from all electrons is always included as IC seed photon field). The total emission is shown by the solid lines. The light blue line denotes infrared dust emission (identical for all models). The synchrotron data points are the same as in Fig. 1, with a 7% systematic uncertainty added in quadrature. The IC data are from this work (cf. Figs. 2 and 4). The bottom panel shows the ratio of the data points as well as the constant B-field and K&C model with respect to the variable B-field model.</figcaption></figure>

the measured extension in this part of the spectrum, but also in the high-energy part of the IC spectrum. Comparing the models among each other, we observe a very good agreement between the static models, while the K&C model behaves slightly differently. This is most likely due to the differences in the parametrization of the wind electron distribution. The first break in the spectrum, just beyond $\gamma_{w,min}$, helps the static models with the transition between the UV and X-ray flux measurements. The K&C model, on the other hand, has only one break in the wind spectrum as the attempt of adding a second break resulted in no significant improvement. Instead, the wind electrons in the K&C model follow the form $(\gamma_{w,min}+\gamma)^s$ (compared to $\gamma^s$ for the static models), which has a similar effect. It should be noted that the K&C model requires a harder electron index after the break to achieve agreement with the hard X-ray spectrum. Specifically, the spectral index of the electron spectrum hardens from 2.87 to 2.32, which is difficult to explain in the framework of the MHD flow model.

A zoom-in of the IC part of the SED is shown in Fig. 12.<fnref n="8" /> We stress again that the models have not been fitted to the Fermi-LAT and H.E.S.S. flux points shown in this figure (and Fig. 11), but to the binned event data of these instruments directly; the flux points are shown for comparison only here. As indicated by the error bars of the flux points, the data provide the strongest constraints at energies of a few GeV, and at around 1 TeV. As a result, the model predictions are very similar up to an energy of a few TeV. The main differences appear at energies above 10 TeV, where in particular the constant B-field model over-predicts the measured flux. The electrons responsible for this flux are also producing the highest-energy synchrotron flux,

where we find an excellent agreement between the static models. However, the variable B-field model predicts the same synchrotron flux with a softer electron spectrum and lower high-energy cutoff (cf. Fig. 8), because of the higher B-field strength near the shock radius. This leads to a stronger reduction of the IC flux with respect to the synchrotron flux. When extrapolating the IC flux above 100 TeV, the variable B-field model agrees with the LHAASO data (Cao et al. 2021) up to ~1 PeV, even though these data are not included in the fit.

The seed photon field contributions are shown in Fig. 12 for the variable B-field model, but are very similar for the other models. It can be seen that the SSC component follows the spectral shape across all energies. The dust component is similar in strength to the SSC component between ~10 GeV and ~10 TeV, while being suppressed at lower and higher energies. The CMB component, on the other hand, behaves in exactly the opposite way. This rather complex behaviour emerges as the seed photon fields do not only differ in their spectrum, but also their morphology. While the CMB is constant at all radii, the SSC component approximately follows a Gaussian profile and the dust component is modelled as a shell around the nebula. As an example, at several tens of TeV the extension of the electron population becomes smaller than the inner radius of the dust shell, leading to a reduced contribution to the IC emission from dust photons at these energies. Another way of illustrating the spatial distribution of the emission in close vicinity to the centre is displayed in Fig. 13, where we show the energy-integrated intensity profiles (emissivity integrated along the line of sight) for the thermal dust, synchrotron, and IC emission, as predicted by our best-fit models. It is evident that, compared to the constant B-field model, the two models with a variable B-field predict a narrower profile for the synchrotron emission. Despite this, the variable B-field model exhibits the flattest IC profile of all models. In

<footnote n="8">⁸ We note that the unsteady behaviour of the K&C model at the highest energies is due to numerical instabilities and therefore not physical.</footnote>
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** IC flux predictions of the three models, together with the Fermi-LAT and H.E.S.S. flux points derived in this work. The LHAASO flux points, taken from Cao et al. (2021), are only shown for comparison. The dashed, dotted and dash-dotted lines show the individual contributions of seed photon fields for the variable B-field model. The bottom panel shows the ratio of the data points as well as the constant B-field and K&C model with respect to the variable B-field model.</figcaption></figure>

other words, the rapidly decreasing B-field allows for a more extended electron distribution while maintaining a narrow synchrotron profile. This leads to a larger predicted extension in the IC regime, which matches the observations best.

### 4.3.3. Nebula extension

In Fig. 14, we show the extension of the nebula predicted by the models from the radio to the γ-ray domain, and compare this to the corresponding measurements. The extensions of the respective electron distributions are shown in Fig. 15. At energies for which the radio electrons dominate (≤10 meV for the synchrotron and ≤10 MeV for the IC part), the morphology of the model is independent of energy. When the more extended low-energy wind electrons start to dominate, the maximal extension is reached in the optical (around 100 MeV for the IC part). That the low-energy wind electrons exhibit a larger extension than the radio electrons can be physically motivated by energy-dependent diffusion, see for example Tang & Chevalier (2012).

First, we note that our measurement of a decreasing extension with increasing energy in the IC regime is generally well in line with the model predictions, and thus not a surprising result. Ultimately, this is a consequence of the already observed decrease of the size of the nebula in the synchrotron domain. Comparing the models in detail, we find that only the K&C model fully predicts the strong decrease in extension towards the X-ray energies, whereas the static models seem to be influenced more by the larger extensions measured at γ-ray energies. Accordingly, the K&C model under-predicts the extension across most of the IC regime, while the variable B-field model – albeit still predicting a too-small extension – provides a better match to the data here. We furthermore find strong correlations between the spectral and spatial parameters, as a narrower spatial distribution also decreases the total number of electrons

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Energy-integrated intensity predicted by the three models, as a function of distance from the centre of the Crab Nebula. The solid, dashed, and dotted line correspond to IC, synchrotron, and dust emission, respectively; the colours refer to the three different tested models. The intensity has been integrated above 1 meV for the synchrotron emission and above 1 GeV for the IC emission. The dark grey band up to $z=1$ indicates the region up to the shock radius, in which no emission is predicted by the models. The light grey band shows the dust shell, in which all the dust emission happens. The profiles show non-zero emission inside the shock region and dust emission outside the dust shell because the emissivity is integrated along the line of sight.</figcaption></figure>

and thus the γ-ray flux at the corresponding energies. It does not seem possible for any of the models – for the given parameters of the electron spectrum and the seed photon fields – to predict
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Comparison of the 68% flux containment predicted by the models as function of energy (dashed for the synchrotron part and solid for the IC part). The grey synchrotron extension points are taken from Dirson & Horns (2023), and references therein. The error bars of the data points in the IC domain denote statistical uncertainties only.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Comparison of the 68% containment radii of the electron energy densities of the three models. The radio electrons dominate the low-energy part of the spectrum where the extension is not energy-dependent. The extension of the wind electrons changes with energy, following, in case of the static models, a simple power law. The spectral transition of the two distributions is shown in Fig. 8.</figcaption></figure>

the measured IC extensions together with the broadband SED. We see two possible reasons for that.

First, the discrepancies could be due to simplifying assumptions made by the models. For instance, all models assume spherical symmetry, which is clearly an oversimplification. Additionally, three-dimensional and non-ideal MHD calculations suggest deviations from the predictions of the K&C model (Porth et al. 2014; Bucciantini et al. 2017; Tanaka et al. 2018; Lyutikov et al. 2019; Luo et al. 2020). In particular, non-ideal MHD models which predict turbulent magnetic fields to be present can resolve a number of problems of the K&C model (Luo et al. 2020). This would necessarily alter the magnetic field structure in the nebula and thereby the relation between synchrotron and IC photons. This could potentially lead to a better agreement with the data.

However, a spatially resolved calculation of such models is so far not available and the current non-ideal MHD models cannot reproduce the full SED in all its details. Another option would be the presence of an additional electron population that extends beyond the synchrotron nebula, for example electrons that have already escaped the nebula or electrons moving inwards from the so-far unobserved outer shock. The IC scattering of such a population with the CMB could potentially lead to a larger IC extension.

The second possibility is systematic uncertainties, for example a mis-modelling of the H.E.S.S. PSF beyond the estimated systematic uncertainty, which could cause the extension measured in the IC regime to appear larger than it actually is. Systematic uncertainties could also include an unaccounted factor between true and reconstructed energies of the IC photons. However, the required magnitude of such systematic effects appears to be unrealistically high to bring the measurements in agreement with the model predictions. We elaborate on this in the next section.

## 5. Systematic uncertainties and energy scale

Several systematic uncertainties could impact our analysis. For example, the absolute energy scales of Fermi-LAT and H.E.S.S. could differ, as the energy reconstruction of the two instruments is calibrated in different ways. Specifically, for Fermi-LAT a beam test in combination with a measurement of the geomagnetic cutoff was used to infer the absolute energy scale with an accuracy of ±2% (Ackermann et al. 2012). Reaching this level of accuracy is not possible with H.E.S.S., where the atmosphere effectively is part of the detector and variations in atmospheric conditions can lead to a change in the number of Cherenkov photons reaching the telescopes. While in our selection of observation runs we have excluded observations carried out under sub-optimal conditions, distributions of the “Cherenkov Transparency Coefficient” (Hahn 2014) for the selected observations indicate that variations of about 10% remain. A difference in energy scale between the instruments could thus lead to a systematic shift between the spectrum of the Crab Nebula measured with Fermi-LAT and H.E.S.S. (cf. also Nigro et al. 2019).

One method of incorporating a scale factor $\eta$ between the energy scales is to evaluate the H.E.S.S. source model at a modified energy $\eta \cdot E$, while the evaluation of the Fermi-LAT source model is unchanged (Meyer et al. 2010). Leaving $\eta$ as an additional free “nuisance” parameter during the fit determines the difference in energy scale between the two instruments. This method, however, depends on the assumed spectral model, and a model that fits the data poorly can falsely indicate a scale factor significantly different from unity if this improves the agreement with the data. For example, adopting a log-parabola spectrum leads to a factor $\eta = 0.79 \pm 0.03$ simply because the model is not suitable for the whole energy range. Other models with the possibility of a sharper break in the spectrum result in values $\eta > 1$. For the variable B-field model, we find $\eta = 1.038$, with only a very minor improvement of the fit ($\Delta C_{\text{tot}} \approx 2$, implying that $\eta$ is consistent with a value of 1). As we cannot be certain about the correctness of our models, we therefore chose the approach of fixing $\eta = 1$ in our final analysis. This is further supported by looking at the integrated flux in the overlapping energy range between the Fermi-LAT data and the H.E.S.S. mono data set (i.e. 0.24–1.80 TeV). In this energy range, we do not see an indication for the energy shift as the integrated fluxes agree within statistical uncertainties. Specifically, summing the flux points in this range we find a Fermi-LAT flux
<table><caption>**Table 4.** Fit quality with and without including systematic uncertainties.</caption><thead><tr><th rowspan="2">Data set</th><th>w/o systematics</th><th rowspan="2">N<sub>dof</sub></th><th>w/systematics</th></tr><tr><th>χ<sup>2</sup> (sys. error)</th><th>χ<sup>2</sup> (sys. error)</th></tr></thead><tbody><tr><td>H.E.S.S. stereo (flux points)</td><td>45.5 (none)</td><td>35</td><td>33.6 (3%)</td></tr><tr><td>H.E.S.S. mono (flux points)</td><td>36.7 (none)</td><td>13</td><td>12.1 (4%)</td></tr><tr><td>Fermi-LAT (flux points)</td><td>41.7 (none)</td><td>26</td><td>27.7 (10%)</td></tr><tr><td>Synchrotron (flux points)</td><td>235.0 (5%)</td><td>194</td><td>152.7 (7%)</td></tr><tr><td>Fermi-LAT + H.E.S.S. (extension)</td><td>121.5 (none)</td><td>6</td><td>8.9 (see text)</td></tr><tr><td>Synchrotron (extension)</td><td>60.3 (n.a.)</td><td>13</td><td>13.0 (8")</td></tr></tbody></table>

**Notes.** The $\chi^2$ values are calculated for the respective data set and the variable B-field model. For each data set, we conservatively estimate the number of degrees of freedom ($N_{dof}$) as the number of data points minus one. The magnitude of the systematic uncertainties (specified in parentheses) are chosen such that a good fit quality is obtained, that is, $\chi^2 \sim N_{dof}$, and thus do not represent estimates of the actual systematic uncertainty of each data set. The first section of the table summarises the flux measurements, where the systematic uncertainty is specified as a relative error on the flux level. The extension uncertainty, given in the second section of the table, is increased by a constant error for each data set (see main text for details). We note that we still add a systematic uncertainty of 5% to the synchrotron flux points for the 'w/o systematics' case, as the $\chi^2$-value would otherwise increase to >10$^6$, due to the small statistical uncertainties of some of the points. Similarly, the extension measurements in the synchrotron regime already contain systematic uncertainties, which are however not easily separable and thus not stated explicitly here (see Dirson & Horns 2023 for details).

of $(1.9 \pm 0.3) \times 10^{-10} \text{ cm}^{-2} \text{ s}^{-1}$, which is consistent with the H.E.S.S. mono flux of $(2.09 \pm 0.2) \times 10^{-10} \text{ cm}^{-2} \text{ s}^{-1}$ (statistical uncertainties only).

Nonetheless, we have aimed to estimate the systematic uncertainty associated with our flux measurements. To do so, we assessed the quality of the fit of the variable B-field model by means of a $\chi^2$ test, and determined the magnitude of a potential systematic error that would lead to an acceptable fit quality. For the spectral measurements in the $\gamma$-ray regime, we added a constant percentage of the flux in quadrature to the statistical uncertainty. To calculate the $\chi^2$ values for the $\gamma$-ray data sets, we used the asymmetric errors of each flux point, meaning that we used the positive error if a point falls below the model prediction and vice versa. In the synchrotron range, the systematic error is also added as a percentage of the measured flux in quadrature to the statistical errors (as in the fitting process, see Sect. 4). For the synchrotron extension measurements, the assumed systematic error is a constant angle. Finally, in the IC regime, we used the systematic errors of the extensions previously cited in the single-instrument analyses. That is, at energies for which the Fermi-LAT data dominate we added an error of $^{+0.54'}_{-0.48}$(cf. Sect. 2.3) and at energies for which the H.E.S.S. data dominate we added an error of $^{+0.21'}_{-0.24}$(cf. Sect. 2.2). Again, all of the systematic uncertainties have been added in quadrature to the statistical ones.

Table 4 summarises the $\chi^2$ values for the variable B-field model with and without the added systematic uncertainties, together with the estimated number of degrees of freedom ($N_{dof}$). For the case of a single data set to which a model has been fitted, $N_{dof}$ is equal to the number of data points minus the number of model parameters (when neglecting correlations between the parameters). We note that in the present case, with the models being fitted to multiple data sets at once, it is not straightforward to properly state the exact $N_{dof}$ for each of the individual sets. The reason for this is that the models have not been adjusted to any of the individual sets alone, but to their combination. The high degree of correlation between some of the model parameters presents an additional complication. In Table 4, we therefore provide as a very conservative estimate of $N_{dof}$ the number of points in the data set minus one (taking into account one overall scale factor that remains even if all model parameters are perfectly correlated). We find that without added systematic

uncertainties, the $\chi^2$ values are all well above $N_{dof}$, indicating an unsatisfactory model description or the presence of systematic uncertainties. For the H.E.S.S. stereo and mono flux points, an additional systematic flux uncertainty of 3% and 4%, respectively, leads to an acceptable fit quality. This flux uncertainty would correspond to an uncertainty on the energy scale of the instrument of 1.2% and 1.6%, respectively. For the Fermi-LAT data set, a systematic flux uncertainty of 10% is required, corresponding to a 5% uncertainty on the energy scale. We stress that these values are not generally valid estimates for the energy scale uncertainty of H.E.S.S. and Fermi-LAT, but only the level of systematic scaling required to achieve an acceptable goodness-of-fit in our analysis. The variations are, however, well within the typically adopted systematic uncertainties of the two instruments$^9$. The $\chi^2$ value of the synchrotron flux points depends sensitively on the assumed systematic error because of the relatively small statistical uncertainties on some of the points. We find that adding a 7% flux uncertainty leads to a $\chi^2$ value well below $N_{dof}$, while 5% is not sufficient. Compatible best-fit parameters are obtained with either assumed systematic error.

Regarding the extension of the nebula, it is evident that without consideration of systematic uncertainties, the predicted extension in both the IC and synchrotron regime is strongly at odds with the measurements. For the Fermi-LAT and H.E.S.S. data sets, taking into account the previously estimated systematic uncertainties leads to a marginally acceptable fit quality. The assumed error bars of the synchrotron data points, on the other hand, already contain systematic uncertainties (see Dirson & Horns 2023). Nevertheless, an additional uncertainty of 8", which is approximately twice as large as the uncertainty used in the fit (~4" on average), would be required to achieve $\chi^2 \approx N_{dof}$. This illustrates again the previously discussed inability of the model to simultaneously describe the SED and the extension measurements. The simplifying assumptions of the model already discussed in the preceding Sect. 4.3.3, as for example radial symmetry and a simple power law describing the electron distribution size evolution with energy, may not be justified.

<footnote>$^9$ See e.g. Aharonian et al. (2006) for H.E.S.S. and https://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html for Fermi-LAT.</footnote>
## 6. Conclusion

In this work, we present a joint Fermi-LAT and H.E.S.S. analysis of the Crab Nebula, carried out with the open-source analysis package GAMMAPY. Below, we first summarise our findings regarding the technical part of our work – related to the joint fit of the Fermi-LAT and H.E.S.S. data – before we discuss the implications of our analysis on emission models of the Crab Nebula.

First, we have demonstrated that applying the open-source analysis package GAMMAPY to Fermi-LAT and H.E.S.S. data, we are able to obtain results that are consistent with the traditionally used analysis chains. Specifically, for Fermi-LAT, we show that we can almost exactly reproduce results obtained with the FERMIPY package after some minor adjustments in the data reduction procedure. In the context of this work, GAMMAPY has the advantage that it enables a combination of the Fermi-LAT data with that of H.E.S.S., and it offers the possibility to include customised physical models in the analysis. For H.E.S.S., on the other hand, GAMMAPY allows us to apply the three-dimensional likelihood analysis method, that is, to model the spectrum and morphology of the Crab Nebula simultaneously under full consideration of the Poisson event statistics. A key ingredient for this is a three-dimensional model for the residual hadronic background. We find consistent results with more traditional methods that estimate this residual background from the analysed observations themselves. Second, by combining the Fermi-LAT and H.E.S.S. data in a joint-likelihood analysis, we provide a fully consistent spectro-morphological analysis of the Crab Nebula over five decades in energy. We confirm previously found indications for an extension of the nebula in the IC regime that decreases with increasing energy – a result that agrees well with theoretical expectations. We discuss the possibility of taking into account possible systematic effects (e.g. different energy scales between instruments) in the form of nuisance parameters in the likelihood fit. Finally, to be able to fit physical emission models to our data that are consistent with other MWL observations, we demonstrate that it is possible to take into account MWL constraints through the addition of penalty terms to the fit statistic.

This has enabled us to consistently model the synchrotron and IC part of the emission from the Crab Nebula. Specifically, we have confronted three phenomenological models for the non-thermal emission of the Crab Nebula with our measurements: a constant B-field model based on Meyer et al. (2010), a variable B-field model based on Dirson & Horns (2023), and the K&C model by Kennel & Coroniti (1984b). We note that these models are open-source and publicly available<fnref n="10" />. All models predict the γ-ray emission not only as a function of energy, but also as a function of the distance from the centre of the nebula. Hence, by projecting the symmetric model on the sky, we were able to fit a three-dimensional model of the Crab Nebula emission directly to the measured event data from Fermi-LAT and H.E.S.S.. We conclude that, within the range of the estimated systematic uncertainties, none of the models are able to describe the full SED together with the extension of the nebula at different energies. We firmly rule out the constant B-field model, as it completely fails to describe the decreasing extension with energy and the spectrum at the same time. The K&C model describes the SED and the extension of the synchrotron nebula well, but under-predicts the size of the IC nebula. Furthermore, it requires a spectral hardening in the electron spectrum to fit

the X-ray data, and the best-fit value of the magnetisation $\sigma$ is at odds with the position of the shock and the expansion velocity of the nebula. Finally, the variable B-field model achieves the best agreement with the data. This indicates that the magnetic field strength in the nebula decreases with increasing distance from the pulsar. However, compared to the model prediction, the new H.E.S.S. measurements of the extension presented here are in tension with the extension measurements in the X-ray domain. As a consequence, the variable B-field model does not achieve the same fit quality as in Dirson & Horns (2023). In summary, even though we find strong evidence that the extension of the IC emission decreases with energy, it still appears too large compared to the distribution of the synchrotron photons, which are supposedly produced by the same electrons. This could suggest more complicated spatial profiles of the electron distribution and the magnetic field as predicted, for example, in non-ideal MHD calculations (e.g. Porth et al. 2014; Luo et al. 2020). Additionally, a non-spherical geometry (such as a toroidal or flattened profile) will also change the seed photon density, which could lead to a better agreement between the model predictions and data.

As a final note, we remark that the IC spectrum above ~30 TeV is not well constrained by the data used in this work. Hence, differences between the models occur predominantly in this regime. In future, including data from instruments optimised for the highest energies, such as LHAASO, could help in placing further constraints on the models and in determining the shape of the magnetic field in more detail. Moreover, data from the upcoming Cherenkov Telescope Array (CTA; Acharya et al. 2018; Hofmann & Zanin 2023) with its unprecedented angular resolution will allow us to measure the extension of the Crab Nebula in the IC regime with higher precision.

**Acknowledgements.** We thank Matthew Kerr for providing the pulsar timing information for our Fermi-LAT data set. The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG), the Helmholtz Association, the Alexander von Humboldt Foundation, the French Ministry of Higher Education, Research and Innovation, the Centre National de la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU), the Commissariat à l'énergie atomique et aux énergies alternatives (CEA), the U.K. Science and Technology Facilities Council (STFC), the Irish Research Council (IRC) and the Science Foundation Ireland (SFI), the Knut and Alice Wallenberg Foundation, the Polish Ministry of Education and Science, agreement no. 2021/WK/06, the South African Department of Science and Technology and National Research Foundation, the University of Namibia, the National Commission on Research, Science & Technology of Namibia (NCRST), the Austrian Federal Ministry of Education, Science and Research and the Austrian Science Fund (FWF), the Australian Research Council (ARC), the Japan Society for the Promotion of Science, the University of Amsterdam and the Science Committee of Armenia grant 21AG-1C085. We appreciate the excellent work of the technical support staff in Berlin, Zeuthen, Heidelberg, Palaiseau, Paris, Saclay, Tübingen and in Namibia in the construction and operation of the equipment. This work benefited from services provided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation. The Fermi-LAT Collaboration acknowledges support for LAT development, operation and data analysis from NASA and DOE (United States), CEA/Irfu and IN2P3/CNRS (France), ASI and INFN (Italy), MEXT, KEK, and JAXA (Japan), and the K.A. Wallenberg Foundation, the Swedish Research Council and the National Space Board (Sweden). Science analysis support in the operations phase from INAF (Italy) and CNES (France) is also gratefully acknowledged. This work performed in part under DOE Contract DE-AC02-76SF00515. The authors gratefully acknowledge the compute resources and support provided by the Erlangen National High Performance Computing Center (NHR@FAU). This research made use of the ASTROPY (https://www.astropy.org; Astropy Collaboration 2013, 2018, 2022), MATPLOTLIB (https://matplotlib.org; Hunter 2007) and IMINUIT (https://iminuit.readthedocs.io; Dembinski et al. 2020) software packages.

<footnote n="10">¹⁰ See https://github.com/me-manu/crabmeyerpy for the implementation we have used.</footnote>
## References

Abdalla, H., Abramowski, A., Aharonian, F., et al. 2017, A&A, 600, A89

Abdalla, H., Aharonian, F., Ait Benkhali, F., et al. 2018, A&A, 620, A66

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

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

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33

Abe, H., Abe, K., Abe, S., et al. 2023, ApJ, 956, 80

Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017, ApJ, 843, 39

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

Abramowski, A., Acero, F., Aharonian, F., et al. 2010, A&A, 520, A83

Acciari, V. A., Ansoldi, S., Antonelli, L. A., et al. 2020, A&A, 635, A158

Acharya, B. A., Agudo, I., Al Samarai, I., et al. 2018, Science with the Cherenkov Telescope Array (World Scientific Publishing)

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

Ackermann, M., Ajello, M., Baldini, L., et al. 2018, ApJS, 237, 32

Aharonian, F. A., Akhperjanian, A. G., Barrio, J. A., et al. 2000, A&A, 361, 1073

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

Aharonian, F., An, Q., Axigeku, et al. 2021, Chinese Phys. C, 45, 025002

Ajello, M., Atwood, W. B., Baldini, L., et al. 2022, Science, 376, 521

Akaike, H. 1974, IEEE Trans. Automat. Contr., 19, 716

Albert, J., Aliu, E., Anderhub, H., et al. 2008, ApJ, 674, 1037

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

Amato, E., & Olmi, B. 2021, Universe, 7, 448

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

Ansoldi, S., Antonelli, L. A., Antoranz, P., et al. 2016, A&A, 585, A133

Arons, J. 1998, Mem. Soc. Ast. Ital., 69, 989

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167

Atoyan, A. M. 1999, A&A, 346, L49

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

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071

Atwood, W., Albert, A., Baldini, L., et al. 2013, arXiv e-prints [arXiv:1303.3514]

Bednarek, W. 2003, A&A, 407, 1

Berge, D., Funk, S., & Hinton, J. 2007, A&A, 466, 1219

Bernlöhr, K. 2008, Astropart. Phys., 30, 149

Bi, B., Barcelo, M., Bauer, C. W., et al. 2021, Proc. 37th Int. Cosmic Ray Conf. – PoS(ICRC2021), 395, 743, https://pos.sissa.it/395/743

Blumenthal, G. R., & Gould, R. J. 1970, Rev. Mod. Phys., 42, 237

Braun, I. 2007, Ph.D. Thesis, Ruprecht-Karls-Universität Heidelberg, Germany

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

Buehler, R., & Blandford, R. 2014, Rep. Prog. Phys., 77, 066901

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

Cao, Z., Aharonian, F., An, Q., et al. 2021, Science, 373, 425

Cash, W. 1979, ApJ, 228, 939

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

de Neurois, M., & Rolland, L. 2009, Astropart. Phys., 32, 231

Deil, C., Donath, A., Terrier, R., et al. 2020, https://doi.org/10.5281/zenodo.4701500

Dembinski, H., Ongmongkolkul, P., Deil, C., et al. 2020, https://doi.org/10.5281/zenodo.3949207

Dirson, L., & Horns, D. 2023, A&A, 671, A67

Donath, A., Terrier, R., Remy, Q., et al. 2023, A&A, 678, A157

Fomin, V. P., Stepanian, A. A., Lamb, R. C., et al. 1994, Astropart. Phys., 2, 137

Hahn, J., de los Reyes, R., Bernlöhr, K., et al. 2014, Astropart. Phys., 54, 25

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Hester, J. J. 2008, ARA&A, 46, 127

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

Hofmann, W., & Zanin, R. 2023, arXiv e-prints [arXiv:2305.12888]

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

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

Kennel, C. F., & Coroniti, F. V. 1984a, ApJ, 283, 710

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

Lam, S. K., Pitrou, A., & Seibert, S. 2015, Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, 1

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

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

Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396

Meagher, K. 2015, Proc. 34th Int. Cosmic Ray Conf., 236, 792

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

Meyer, M., Scargle, J. D., & Blandford, R. D. 2019, ApJ, 877, 39

Mitchell, A. M. W. 2016, Ph.D. Thesis, Ruperto-Carola-University Heidelberg Germany, https://doi.org/10.11588/heidok.00021768

Mohrmann, L., Specovius, A., Tiziani, D., et al. 2019, A&A, 632, A72

Murach, T., Gajdus, M., & Parsons, R. D. 2015, Proc. 34th Int. Cosmic Ray Conf. – PoS(ICRC2015), 236, 1022, https://pos.sissa.it/236/1022

Nie, L., Liu, Y., Jiang, Z., & Geng, X. 2022, ApJ, 924, 42

Nigro, C., Deil, C., Zanin, R., et al. 2019, A&A, 625, A10

Ohm, S., van Eldik, C., & Egberts, K. 2009, Astropart. Phys., 31, 383

Parsons, R., & Hinton, J. 2014, Astropart. Phys., 56, 26

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

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

Pühlhofer, G., Bolz, O., Götting, N., et al. 2003, Astropart. Phys., 20, 267

Pühlhofer, G., Bernlöhr, K., Bi, B., et al. 2021, Proc. 37th Int. Cosmic Ray Conf. – PoS(ICRC2021), 395, 764, https://pos.sissa.it/395/764

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

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

Tang, X., & Chevalier, R. A. 2012, ApJ, 752, 83

Torres, D. F., Cillis, A., Martín, J., & de Oña Wilhelm, E. 2014, J. High Energy Astrophys., 1, 31

Vacanti, G., Cawley, M. F., Colombo, E., et al. 1991, ApJ, 377, 467

van Rensburg, C., Krüger, P. P., & Venter, C. 2018, MNRAS, 477, 3853

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

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

Weisskopf, M. C., Hester, J. J., Tennant, A. F., et al. 2000, ApJ, 536, L81

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

Wilks, S. S. 1938, Ann. Math. Stat., 9, 60

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

Wood, M., Caputo, R., Charles, E., et al. 2017, arXiv e-prints [arXiv:1707.09551]

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

Yeung, P. K. H., & Horns, D. 2021, ApJ, 909, 221

1. Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland

2. Max-Planck-Institut für Kernphysik, PO Box 103980, 69029 Heidelberg, Germany

3. Yerevan State University, 1 Alek Manukyan St, Yerevan 0025, Armenia

4. Landessternwarte, Universität Heidelberg, Königstuhl, 69117 Heidelberg, Germany

5. Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands

6. Laboratoire Leprince-Ringuet, École Polytechnique, CNRS, Institut Polytechnique de Paris, 91128 Palaiseau, France

7. University of Namibia, Department of Physics, Private Bag 13301, Windhoek 10005, Namibia

8. Centre for Space Research, North-West University, Potchefstroom 2520, South Africa

9. Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, 22761 Hamburg, Germany

10. Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany

11. Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Strasse 24/25, 14476 Potsdam, Germany

12. Université de Paris, CNRS, Astroparticule et Cosmologie, 75013 Paris, France

13. Department of Physics and Electrical Engineering, Linnaeus University, 351 95 Växjö, Sweden

14. Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany

15. Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, 72076 Tübingen, Germany

16. Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université Paris Cité, 5 Pl. Jules Janssen, 92190 Meudon, France

17. Sorbonne Université, Université Paris-Diderot, Sorbonne Paris Cité, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Énergies, LPNHE, 4 Place Jussieu, 75252 Paris, France

18. IRFU, CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette, France
19 University of Oxford, Department of Physics, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK

20 Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Nikolaus-Fiebiger-Str. 2, 91058 Erlangen, Germany

21 Astronomical Observatory, The University of Warsaw, Al. Ujazdowskie 4, 00-478 Warsaw, Poland

22 Université Savoie Mont Blanc, CNRS, Laboratoire d'Annecy de Physique des Particules – IN2P3, 74000 Annecy, France

23 Instytut Fizyki Jądrowej PAN, ul. Radzikowskiego 152, 31-342 Kraków, Poland

24 Université Bordeaux, CNRS, LP2I Bordeaux, UMR 5797, 33170 Gradignan, France

25 School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, Johannesburg 2050, South Africa

26 Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, 34095 Montpellier Cedex 5, France

27 School of Physical Sciences, University of Adelaide, Adelaide 5005, Australia

28 Aix-Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

29 School of Science, Western Sydney University, Locked Bag 1797, Penrith South DC, NSW 2751, Australia

30 Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstraße 25, 6020 Innsbruck, Austria

31 Obserwatorium Astronomiczne, Uniwersytet Jagielloński, ul. Orla 171, 30-244 Kraków, Poland

32 Institute of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland

33 Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland

34 Department of Physics and Astronomy, The University of Leicester, University Road, Leicester LE1 7RH, UK

35 GRAPPA, Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

36 Yerevan Physics Institute, 2 Alikhanian Brothers St., 0036 Yerevan, Armenia

37 Department of Physics, Konan University, 8-9-1 Okamoto, Higashinada, Kobe, Hyogo 658-8501, Japan

38 Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, 5-1-5 Kashiwa-no-Ha, Kashiwa, Chiba 277-8583, Japan

39 RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan

40 Theoretical Division, Los Alamos National Laboratory Los Alamos, Alamos, NM 87545, USA
**Appendix A: Comparison with the reflected regions background method**

In this section we compare the H.E.S.S. flux points derived with the three-dimensional (3D) likelihood analysis as described in this paper to those derived with a more traditional one-dimensional (1D) analysis, in which only the energy dimension is considered. The main difference between the two methods lies in the estimation of the residual hadronic background. While the 3D analysis requires a background model, the background for the 1D analysis is estimated using ‘Off’ regions in the field of view (‘reflected regions method’, see Fomin et al. 1994; Berge et al. 2007). The 1D analysis has been carried out using the H.E.S.S.-internal ‘HAP’ pipeline. Because the Crab Nebula is slightly extended we choose a comparatively large ‘On’ region of 0.1° radius. For technical reasons, we also split the H.E.S.S. stereo data set into two subsets for this test, corresponding to the different phases of the instrument. The comparison of the derived spectra is shown in Fig. A.1. At low energies, the GAMMAPY points tend to be higher than those obtained with HAP. We attribute this partly to the fact that the 1D analysis corrects for leakage of γ rays outside of the On region under the assumption of a point-like source, whereas the Crab Nebula appears increasingly extended towards lower energies.

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

**Fig. A.1.** Comparison of the Crab Nebula γ-ray flux measured with H.E.S.S. using different analysis techniques. The blue points with circle markers show the result as derived in the three-dimensional likelihood analysis with GAMMAPY described in this work. The orange points with square markers have been computed with the H.E.S.S. analysis program (HAP), using a traditional one-dimensional analysis based on the ‘reflected background’ method. (a) Data from H.E.S.S. Phase-I (2004–2009, cf. Table 1). (b) Data from H.E.S.S. Phase-II (2013–2015).
## Appendix B: Physical emission models for the Crab Nebula

We model the broadband electromagnetic emission of the Crab Nebula through a radially symmetric steady-state model of the electron distribution, inspired by the constant B-field model introduced by Meyer et al. (2010), which has been extended to a spatially variable B-field model in Dirson & Horns (2023). An alternative description is the radially symmetric MHD flow solution obtained by Kennel & Coroniti (1984b). For comparison with the data, we need to calculate the predicted emission as a function of energy and angular separation between the line of sight and the nebular centre (this is equal to the two-dimensional sky position in the radially symmetric case). The basic ingredients to do so are the radial profile of the magnetic field, $B(r)$, and the electron distribution, $n_{el}(\gamma, r) \equiv dN_{el}/dydV$, that is, the number of electrons per unit volume and Lorentz factor $\gamma$ as a function of the distance to the nebula centre.

Given $B(r)$ and $n_{el}$, we calculate the spectral volume emissivity $j_v^{sync}(v, r) \equiv dE/dtdVdvd\Omega$ for the synchrotron emission of a (pitch-angle averaged) chaotic magnetic field through (Blumenthal & Gould 1970; Aharonian et al. 2010)

$$j_v^{\text{sync}}(v, r) = \frac{1}{4\pi} \frac{\sqrt{3}e^3 B(r)}{m_e c^2} \int_1^\infty dy n_{\text{el}}(\gamma, r) G\left(\frac{\nu}{v_c}\right), \quad (B.1)$$

where $v_c = 3eBy^2/4\pi m_e c$ is the critical frequency and $G(x)$ is given by

$$G(x) = \frac{x}{20} \left[ (8 + 3x^2)(K_{1/3}(x/2))^2 + xK_{2/3}(x/2)(2K_{1/3}(x/2) - 3xK_{2/3}(x/2)) \right], \quad (B.2)$$

with $K_\xi$ the modified Bessel function of kind $\xi$. Throughout this paper, we use units of erg s$^{-1}$Hz$^{-1}$ cm$^{-3}$ sr$^{-1}$ for $j_v$.

For the dust emission, we follow Dirson & Horns (2023), who model the emission as a mixture of two dust populations at different temperatures $T_i$ and total masses $M_i$, $i = 1, 2$. The two populations are contained in a shell with inner and outer radii $r_{\text{dust, in}}$ and $r_{\text{dust, out}}$, respectively, with respect to the pulsar position (which is equal to the centre of the nebula in our model). The density of dust within the shell is assumed to be constant. Assuming the two populations predominantly consist of amorphous carbon dust grains, the emissivity is given by

$$j_v^{\text{dust}}(v, r) = \frac{3\kappa_{\text{abs}}[\Theta(r - r_{\text{dust, in}}) - \Theta(r - r_{\text{dust, out}})]}{4\pi(r_{\text{dust, out}}^3 - r_{\text{dust, in}}^3)} \sum_{i=1,2} M_i B_v(T_i), \quad (B.3)$$

where $\Theta(r)$ is the Heavyside step function, and $B_v(T)$ is the intensity of the black body emission at temperature $T$. The wavelength-dependent absorption coefficient is given by

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

Both the synchrotron and dust emission together with the CMB act as seed photon fields for IC scattering of the electrons producing the synchrotron emission. The seed photon density, $n_{\text{seed}}^t$, $t = (\text{sync}, \text{dust})$, as a function of photon energy $\epsilon = h\nu$ is calculated through the integral (Atoyan & Aharonian 1996)

$$n_{\text{seed}}^{t}(\epsilon, r) = \frac{r_0}{2hc} \frac{4\pi}{\epsilon} \int_{r_s/r_0}^{1} dy \frac{y}{x} \ln \frac{x+y}{|x-y|} j_v^t(v, r_0 y), \quad (B.5)$$

in units of photons eV$^{-1}$ cm$^{-3}$ and where $x = r/r_0$ and $y = r'/r_0$, with $r_0$ the radius of the nebula and $r_s = 0.13$ pc (Weisskopf et al. 2012) the shock radius until which no emission is assumed.<fnref n="11" /> The seed photon density for CMB photons is simply $n_{\text{seed}}^{\text{CMB}} = (4\pi/(hc))B_v(T_{\text{CMB}})/(h\nu)$. We do not consider any additional photon fields such as optical line emission from the filaments (e.g. Meyer et al. 2010) or interstellar radiation fields, as these fields will give a subdominant contribution compared to the optical synchrotron and infrared dust emission, respectively.

The IC emissivity for up-scattered photons at frequency $\nu$ is then calculated using the integral over the IC kernel $f_{\text{IC}}$, which includes Klein-Nishina effects, and the seed photon density (Blumenthal & Gould 1970),

$$j_v^{\text{IC}}(v, r) = \frac{3\sigma_{\text{T}}hc}{4} \frac{h\nu}{4\pi} \int_1^{\infty} d\gamma \frac{n_{\text{el}}(\gamma, r)}{\gamma^2} \times \int_0^{\infty} d\epsilon f_{\text{IC}}(v, \epsilon, \gamma) \sum_i \frac{n_{\text{seed}}^t(\epsilon, r)}{\epsilon}, \quad (B.6)$$

with the usual expression for the full IC kernel,

$$f_{\text{IC}}(v, \epsilon, \gamma) = 2q \ln q + (1 + 2q)(1 - q) + \frac{1}{2} \frac{(\Gamma_{\epsilon} q)^2}{1 + \Gamma_{\epsilon} q} (1 - q), \quad (B.7)$$

where $\Gamma_\epsilon = 4\epsilon\gamma/(mc^2)$ and $q = h\nu/(\Gamma_\epsilon(ymc^2 - h\nu))$. The Thompson limit corresponds to $\Gamma_\epsilon \ll 1$ and IC scattering only occurs for $1/(4\gamma^2) \le q \le 1$.

The integrals are evaluated numerically and the code, fully written in python, is available online<fnref n="12" />. The IC emissivity of the synchrotron emission involves four integrals for each frequency $\nu$ and radius $r$, see Eqs. (B.1), (B.5), and (B.6). In order to speed up the calculations, we make heavy use of numerical routines of numpy, scipy, and numba (Harris et al. 2020; Virtanen et al. 2020; Lam et al. 2015) and use multi-dimensional spline interpolations for $j_v^{\text{sync}}$ and $n_{\text{seed}}^{\text{sync}}$.

The specific luminosity for component $t = (\text{sync}, \text{dust}, \text{IC})$ is found from the integration over radius,

$$\mathcal{L}_v^t = \oint d\Omega \int dV j_v^t(v, r) = (4\pi)^2 \int_{r_s}^{r_0} dr r^2 j_v^t(v, r), \quad (B.8)$$

where the two factors of $4\pi$ come from the volume and solid angle integrations, respectively. For the analysis of H.E.S.S. and Fermi-LAT data, we are also interested in the spatial profile of the emission for which we have to calculate the specific intensity $I_v$ as a function of the angular distance $\theta$ from the centre of the nebula. The intensity is given as an integral over the line of sight (LOS) $s$ ($s = 0$ is the position on the LOS which is closest to the centre of the nebula) through the nebula in which $j_v^t$ is non-zero,

$$I_v^t(\theta) = 2 \int_{s_{\text{min}}}^{s_{\text{max}}} ds j_v^t(v, r(s)). \quad (B.9)$$

<footnote n="11">¹¹ One should note that the actual size of the nebula is unknown as the outer shock is not observed (e.g. Hester 2008). The size of the observed synchrotron nebula is a free parameter in our fit (see below) and $r_0$ acts as an arbitrary upper integration bound that should be large enough to accommodate the different nebula sizes probed in the fit. Here, we set it to 3.6 pc, which is twice the value assumed by Atoyan & Aharonian (1996) for the visible nebula size.</footnote>

<footnote n="12">¹² https://github.com/me-manu/crabmeyerpy</footnote>
<figure><img src="image_19.png" /><figcaption>**Fig. B.1.** Illustration of the LOS integration. The pulsar in the centre is surrounded by a region with no emission, up to the shock radius $r_s$ (dark shaded region). We assume non-zero emission until the maximum radius of the nebula at $r_0$. Here we illustrate two exemplary lines of sight $s_1$ and $s_2$ at angular offsets $\theta_1$ and $\theta_2$ with respect to the centre (dashed lines). $s_1$ represents a line which intersects the shock whereas $s_2$ does not. We integrate the emissivity $j_v$ along the parts of the lines that are marked in solid orange. Because of the radial symmetry of our model, these parts contain exactly half of the emission expected from the full LOS.</figcaption></figure>

This is illustrated in Fig. B.1. The distance $r(s)$ from each point on the LOS to the centre of the Crab Nebula calculates as

$$r(s) = \sqrt{r_{\min}^2 + s^2}. \qquad (B.10)$$

In this expression, $r_{\min} = d \sin \theta$ is the smallest distance from the LOS to the centre of the Crab Nebula, which depends on the observing angle $\theta$ and the distance $d$ from the earth to the nebula. Consequently, the emissivity along the LOS is symmetric with respect to the $s = 0$ point which results in the factor two in front of the integral if we choose

$$s_{\min} = \begin{cases} \sqrt{r_s^2 - r_{\min}^2} & r_{\min} \le r_s \\ 0 & \text{otherwise} \end{cases}. \qquad (B.11)$$

The different cases correspond to a scenario where the LOS intersects with the shock radius and the integration starts at the point of the intersection, and another scenario where $\theta$ is large enough so that the LOS misses the shock radius. The upper integration bound $s_{\max}$ is set to the intersection with $r_0$ where the emissivities are negligible small.

The flux emitted from a point on the outer sphere of the nebula into a solid angle $\Omega$ is

$$F_v^t(\theta_0, \theta_1) = \int I_v \cos \theta d\Omega = 2\pi \int_{\theta_0}^{\theta_1} I_v \cos \theta \sin \theta d\theta. \qquad (B.12)$$

The extension of our emission model is defined as the angle $\theta_{68}$ which contains 68 % of the flux, that is,

$$0.68 = \frac{F_v^t(0, \theta_{68})}{F_v^t(0, \theta_{\max})}, \qquad (B.13)$$

where the maximum angle $\theta_{\max}$ is defined through $\tan \theta_{\max} = r_0/d$.

For a full description of the emission model, what remains to be defined are the electron spectrum $n_{el}(\gamma, r)$ and the magnetic field profile $B(r)$. The tested models are described in the following subsections.

## B.1. Variable B-field model

This is a phenomenological model developed in Meyer et al. (2010) and Dirson & Horns (2023), who expanded the previous works of de Jager & Harding (1992) and Hillas et al. (1998) to describe the Crab Nebula’s SED and extension. Two distinct electron populations are assumed: radio and wind electrons. The wind electrons are constantly injected by the pulsar wind and accelerated at the wind’s termination shock, whereas the origin of the radio electrons is still unclear (Atoyan & Aharonian 1996). As the name suggests, the radio electrons are responsible for the synchrotron emission from radio frequencies to sub-millimetre / optical wavelengths whereas the wind electrons produce synchrotron emission at higher frequencies. Both wind and radio electron densities are assumed to be zero for radii smaller than the shock radius. The total electron spectrum is given by the sum of the two populations,

$$n_{\text{el}}(\gamma, r) = n_{\text{radio}}(\gamma, r) + n_{\text{wind}}(\gamma, r). \qquad (B.14)$$

The spectrum of the radio electrons is a simple power law with index $s_r$ between $\gamma$ factors $\gamma_{r,\min}$ and $\gamma_{r,\max}$ with a super-exponential cutoff,

$$n_{\text{radio}}(\gamma, r) = \frac{n_{r,0}}{\rho_r^3} \gamma^{-s_r} \exp\left[-\left(\frac{\gamma}{\gamma_{r,\max}}\right)^{\beta_{\min}}\right] \times \Theta(\gamma - \gamma_{r,\min}) \exp\left(-\frac{r^2}{2\rho_r^2}\right) \times \Theta(r - r_s), \qquad (B.15)$$

where $\Theta(x)$ is the Heavyside step function. The spatial dependence of the radio electrons is given by a Gaussian function of constant width $\rho_r$. The wind electrons, on the other hand, are modelled with a double broken power law with super-exponential cut-offs at low and high energies. As the radio wind electrons, the spatial extension of the wind electrons is assumed to follow a Gaussian whose width, however, also depends on energy. Putting everything together one finds

$$n_{\text{wind}}(\gamma, r) = \frac{n_{w,0}}{\rho_w(\gamma)^3} \exp\left[-\left(\frac{\gamma_{w,\min}}{\gamma}\right)^{\beta_{\min}}\right] \exp\left[-\left(\frac{\gamma}{\gamma_{w,\max}}\right)^{\beta_{\max}}\right] \times \left[\left(\frac{\gamma}{\gamma_{w,1}}\right)^{-s_{w,1}} \left(\frac{\gamma_{w,1}}{\gamma_{w,2}}\right)^{-s_{w,2}} (1 - \Theta(\gamma - \gamma_{w,1})) + \left(\frac{\gamma}{\gamma_{w,2}}\right)^{-s_{w,2}} \mathcal{B}(\gamma, \gamma_{w,1}, \gamma_{w,2}) + \left(\frac{\gamma}{\gamma_{w,2}}\right)^{-s_{w,3}} \Theta(\gamma - \gamma_{w,2})\right] \times \exp\left(-\frac{r^2}{2\rho_w(\gamma)^2}\right) \times \Theta(r - r_s), \qquad (B.16)$$
where $\mathcal{B}(x, a, b)$ is the boxcar function, $\mathcal{B}(x, a, b) = \Theta(x - a) - \Theta(x - b)$. The energy dependence of the spatial extension is modelled through a simple power law

$$ \rho_w(\gamma) = \rho_{w,0} \left[ \left( \frac{\gamma}{9 \cdot 10^5} \right)^2 \right]^{-\alpha_w} . \qquad (B.17) $$

For the free parameters $\gamma_{r,\max}$, $\gamma_{w,\min}$, and $\gamma_{w,\max}$ we do not use abrupt cutoffs as those would lead to discontinuities on the numerical integration grid. Instead we use super-exponential cutoffs which are not motivated physically but help the smoothness of the likelihood landscape. For simplicity, the super-exponential indices are fixed to $\beta_{\min} = 2.8$ (this is the best-fit value obtained by Meyer et al. 2010) and $\beta_{\max} = 2$.

For the magnetic field profile, we follow Dirson & Horns (2023) and choose

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

with $\alpha \ge 0$. The case of $\alpha = 0$ corresponds to the constant B-field model studied, for example, in Meyer et al. (2010).

## B.2. MHD flow model

This model follows the solution to the steady-state MHD equations under the assumption of spherical symmetry and describes the flow of the electron plasma from the shock to the nebula's boundary. The MHD flow and the magnetic field are determined by the magnetisation $\sigma$ at the shock, which is given by the ratio of the electromagnetic and particle energy flux,

$$ \sigma = \frac{B_s^2}{4\pi n_s \gamma m_e c^2}, \qquad (B.19) $$

where $B_s$ is the magnetic field, $n$ the particle density, and $u_s = \sqrt{\gamma^2 - 1}$ the relativistic four speed (all quantities are at the shock radius).

The (injected) electron spectrum is another free parameter in the model. As in the previous model, we use two distinct electron populations, with the radio electrons modelled in the same way as in Section B.1. On the other hand, the spatial and spectral distribution for the wind electrons are the result of solving

the MHD equations for some electron population injected at the termination shock. While Kennel & Coroniti (1984b) considered a power law type injection, this was generalised by Atoyan & Aharonian (1996) to a more complex shape, which we also assume here with an additional break in the spectrum:

$$ n_s(\gamma) = n_{0,w} \exp\left[-\left(\frac{\gamma_{w,\min}}{\gamma}\right)^{\beta_{\min}}\right] \exp\left[-\left(\frac{\gamma}{\gamma_{w,\max}}\right)^{\beta_{\max}}\right] \\ \times \left[ \left(1 + \frac{\gamma}{\gamma_{w,\min}}\right)^{s_{w,1}} (1 - \Theta(\gamma - \gamma_{w,2})) \right. \\ \left. + \frac{(1 + \gamma/\gamma_{w,\min})^{s_{w,3}}}{(1 + \gamma_{w,2}/\gamma_{w,\min})^{s_{w,3}-s_{w,1}} \Theta(\gamma - \gamma_{w,2})} \right]. \qquad (B.20) $$

Similar to the wind electron spectrum for the static models, $\beta_{\min}$ and $\beta_{\max}$ are fixed to 2.8 and 2.0, respectively. From the MHD flow solution, the energy of electrons at some distance $z = r/r_s$ from the shock can be related to the electron Lorentz factor at injection, $\gamma'$,

$$ \gamma(z) = \frac{\gamma'}{(vz^2)^{1/3} + \gamma'/\gamma_{\max}}, \qquad (B.21) $$

where $v(z) = u(z)/u_s$ is the four-velocity of the electrons (see Eqs. (A7a) - (A7d) in Kennel & Coroniti 1984a) which depends on the magnetisation $\sigma$, and $\gamma_{\max}$ is the maximum $\gamma$ factor at distance $z$, see, for example, Eq. (6) in Atoyan & Aharonian (1996) (which also depends on $vz^2$). The wind electron spectrum is then found to be

$$ n_{\text{wind}}(\gamma, z) = (vz^2)^{-4/3} \left( \frac{\gamma'}{\gamma(z)} \right)^2 n_s(\gamma'). \qquad (B.22) $$

For small values of the magnetisation, the magnetic field in the MHD solution at the shock is given by

$$ B_s = \sqrt{\frac{L_{\text{spin-down}}}{cr_s^2} \frac{\sigma}{1 + \sigma}}, \qquad (B.23) $$

where $L_{\text{spin-down}}$ is the spin-down luminosity of the pulsar. Downstream of the shock, the magnetic field in the nebula evolves as

$$ B(z) = B_s \frac{3(1 - 4\sigma)z}{vz^2}. \qquad (B.24) $$
**Appendix C: Best-fit parameters of physical models**

We present the best-fit parameters for all considered models in Table C.1. The best-fit values were found using the 'simple' method of the SHERPA fitting backend of GAMMAPY. The Nelder-Mead algorithm converged into slightly lower minima

than the 'migrad' method from the MINUIT backend. As we found the estimation of uncertainties on the fit parameters to not work reliably – presumably due to the complexity of the model and the large number of parameters – we specify only the best-fit values here.

<table><caption>**Table C.1.** Best-fit parameter values of the physical models. Parameters that are not free in the fit are indicated by '(fixed)' after their value.</caption><thead><tr><th>Parameter</th><th>variable B-field model</th><th>constant B-field model</th><th>Kennel & Coroniti</th></tr></thead><tbody><tr><td>ln(n<sub>r,0</sub>)</td><td>117.170</td><td>117.69</td><td>118.766</td></tr><tr><td>ln(γ<sub>r,min</sub>)</td><td>3.09 (fixed)</td><td>3.09 (fixed)</td><td>3.09 (fixed)</td></tr><tr><td>ln(γ<sub>r,max</sub>)</td><td>11.599</td><td>12.35</td><td>12.625</td></tr><tr><td>s<sub>r</sub></td><td>-1.5439</td><td>-1.649</td><td>-1.7419</td></tr><tr><td>ρ<sub>r</sub> ["]</td><td>88.3</td><td>80.40</td><td>88.64</td></tr><tr><td>ln(n<sub>w,0</sub>)</td><td>76.822</td><td>76.8315</td><td>-27.625</td></tr><tr><td>ln(γ<sub>w,min</sub>)</td><td>12.841</td><td>12.69</td><td>12.8366</td></tr><tr><td>ln(γ<sub>w,1</sub>)</td><td>15.26</td><td>14.24</td><td>—</td></tr><tr><td>ln(γ<sub>w,2</sub>)</td><td>19.197</td><td>19.35379</td><td>17.96</td></tr><tr><td>ln(γ<sub>w,max</sub>)</td><td>22.115</td><td>22.371</td><td>22.251</td></tr><tr><td>β<sub>min</sub></td><td>2.8 (fixed)</td><td>2.8 (fixed)</td><td>2.8 (fixed)</td></tr><tr><td>β<sub>max</sub></td><td>2 (fixed)</td><td>2 (fixed)</td><td>2 (fixed)</td></tr><tr><td>s<sub>w,1</sub></td><td>-3.117</td><td>-2.75</td><td>—</td></tr><tr><td>s<sub>w,2</sub></td><td>-3.3928</td><td>-3.1764</td><td>-2.8695</td></tr><tr><td>s<sub>w,3</sub></td><td>-3.782</td><td>-3.5118</td><td>-2.316</td></tr><tr><td>ρ<sub>w,0</sub> ["]</td><td>98.14</td><td>78.94</td><td>—</td></tr><tr><td>α<sub>w</sub></td><td>0.12544</td><td>0.11973</td><td>—</td></tr><tr><td>B<sub>0</sub> [μG]</td><td>256.4</td><td>126.39</td><td>—</td></tr><tr><td>r<sub>s</sub> ["]</td><td>13.4 (fixed)</td><td>13.4 (fixed)</td><td>13.4 (fixed)</td></tr><tr><td>α</td><td>-0.4691</td><td>—</td><td>—</td></tr><tr><td>σ</td><td>—</td><td>—</td><td>0.021396</td></tr><tr><td>ln(L<sub>spin-down</sub>[erg/s])</td><td>—</td><td>—</td><td>88.716</td></tr><tr><td>r<sub>dust,in</sub> [pc]</td><td>0.55 (fixed)</td><td>0.55 (fixed)</td><td>0.55 (fixed)</td></tr><tr><td>r<sub>dust,out</sub> [pc]</td><td>1.53 (fixed)</td><td>1.53 (fixed)</td><td>1.53 (fixed)</td></tr><tr><td>log<sub>10</sub>(M<sub>1</sub>/M<sub>☉</sub>)</td><td>-4.4 (fixed)</td><td>-4.4 (fixed)</td><td>-4.4 (fixed)</td></tr><tr><td>log<sub>10</sub>(M<sub>2</sub>/M<sub>☉</sub>)</td><td>-1.2 (fixed)</td><td>-1.2 (fixed)</td><td>-1.2 (fixed)</td></tr><tr><td>T<sub>1</sub> [K]</td><td>149 (fixed)</td><td>149 (fixed)</td><td>149 (fixed)</td></tr><tr><td>T<sub>2</sub> [K]</td><td>39 (fixed)</td><td>39 (fixed)</td><td>39 (fixed)</td></tr></tbody></table>

**Notes.** For an explanation of the parameters, see Appendix B.