© The Authors 2024

1. Introduction

The Crab Nebula, associated with the pulsar PSR B0531+21, represents the archetype 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 MeV < E < 100 GeV) and very-high-energy (VHE; E > 100 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 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¹ 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).

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).

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}^{\prime }\pm {0.18}^{\prime }{}_{\text{stat}}\pm {0.42}^{\prime }{}_{\text{sys}}$ measured with the Fermi Large Area Telescope (LAT; Ackermann et al. 2018) and ${1.31}^{\prime }\pm {0.07}^{\prime }{}_{\text{stat}}\pm {0.17}^{\prime }{}_{\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′ × 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 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 ℒ is a measure of how well the model prediction agrees with the measurement and can be calculated as

$$\begin{array}{c}\hfill \mathcal{L}(\xi )={\displaystyle \prod _i}\mathcal{P}(n_i|{\nu }_i(\xi ))={\displaystyle \prod _i}\frac{{\nu }_i^{n_i}(\xi )}{n_i!}·exp(-{\nu }_i(\xi )),\end{array}$$(1)

with 𝒫(n_i|ν_i(ξ)) the Poisson probability to measure n_i events in cell i given the model prediction ν_i for a set of parameters ξ. Besides γ rays from the source under investigation, there are also contributions to n_i and ν_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

$$\begin{array}{c}\hfill \mathcal{C}\equiv -2ln\mathcal{L}=2{\displaystyle \sum _i}({\nu }_i-n_{i}ln({\nu }_i)+ln(n_i!))\end{array}$$(2)

(where the term ln(n_i!), being independent of the model parameters ξ, 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 γ 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 γ-ray induced particle showers in the atmosphere. Based on the recorded shower images, we reconstructed the direction and energy of the primary γ 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 culminates 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).

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² 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·ln(1-0.68)}\sigma \approx 1.5·\sigma$. From the stereo data set, we obtain a fitted extension (averaged over energy) of $1.62\prime \pm {0.05\prime }_{\mathrm{stat}}{}_{-0.24}^{+0.21}{\prime }_{\mathrm{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_{\text{68}}={1.30}^{\prime }\pm {0.07}^{\prime }{}_{\text{stat}}\pm {0.17}^{\prime }{}_{\text{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,

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}E}=N_0{\left(\frac{E}{E_0}\right)}^{-\alpha -\beta ln(E/E_0)},\end{array}$$(3)

with normalisation N₀, reference energy E₀, spectral index α, and curvature parameter β. From a joint fit to all H.E.S.S. data sets, we obtain N₀ = (4.06 ± 0.03)×10⁻¹¹ TeV⁻¹ cm⁻² s⁻¹ at E₀ = 1 TeV (fixed), α = 2.530 ± 0.006, and β = 0.086 ± 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 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%.

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σ 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).

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° ×10° 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 ϕ of all photons that have arrived within a 1° radius around the Crab Nebula, see Fig. 3. The off-pulse phase interval is determined with the “quiescent background” (QB) algorithm introduced by Meyer et al. (2019). The resulting interval is 0.51 ≤ ϕ ≤ 0.89, as indicated by the red-shaded region in Fig. 3.

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° of the Crab Nebula’s position. The shaded regions indicate the intervals for the Bridge emission (0.09 ≤ ϕ ≤ 0.22, not relevant here) and Off region (0.51 ≤ ϕ ≤ 0.89) derived with the QB algorithm. The grey shaded areas indicate regions found by the HOP algorithm, which identifies “watersheds” between local maxima (Meyer et al. 2019). There are shown here for illustration only.

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° 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⁵. 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°  pixel⁻¹ and eight energy bins per decade.

We initialised the ROI model with all point-like and extended γ-ray sources within 20° 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 templates for isotropic and Galactic diffuse emission are used⁶. 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₆₈ = 1.79′). As for the spectrum, the synchrotron part is modelled with a simple power law,

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}E}=N_0{\left(\frac{E}{E_0}\right)}^{-\mathrm{\Gamma }},\end{array}$$(4)

with normalisation N₀, reference energy E₀, and spectral index Γ, 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 Γ = 2 and calculated its TS. No additional sources with $\sqrt{\mathrm{TS}}⩾5$ were found. The best-fit spectrum of the synchrotron component is shown in Fig. 4 and the parameters are provided in Table 2.

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.

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}^{\prime }\pm {{0.18}^{\prime }}_{\mathrm{stat}}{}_{-0.48}^{+0.54}{\prime }_{\mathrm{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_{\text{68}}={1.80}^{\prime }\pm {0.18}^{\prime }{}_{\text{stat}}\pm {0.42}^{\prime }{}_{\text{sys}}$, within the uncertainties.

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 (spectrum 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 𝒞 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 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₆₈ = (2.10 ± 0.18_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:

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}E}=N_0{\left(\frac{E}{E_0}\right)}^{-{\mathrm{\Gamma }}_1}{[1+{\left(\frac{E}{E_{\mathrm{break}}}\right)}^{\frac{{\mathrm{\Gamma }}_2-{\mathrm{\Gamma }}_1}{\beta }}]}^{-\beta },\end{array}$$(5)

with E₀ fixed at 1 TeV. A fit over the whole energy range yields the following parameters: N₀ = (3.35 ± 0.22)×10⁻¹⁰ TeV⁻¹ s⁻¹ cm⁻², Γ₁ = 1.61 ± 0.02, Γ₂ = 2.95 ± 0.02, E_break = 0.33 ± 0.02 TeV, and β = 1.73 ± 0.07. This model is preferred over the log-parabola model by a difference of 94 in the Akaike information criterion, AIC = 𝒞 + 2N_par (Akaike 1974), where the number of free model parameters is N_par = 3 for the log-parabola model and N_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.

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.

The best-fit position over the whole IC energy range is (184.5499 ± 0.0004_stat)° Galactic longitude and ( − 5.7825 ± 0.0004_stat)° Galactic latitude. This position is shifted to the north with respect to that of the Crab Pulsar by ≈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 of the extension. Potentially large systematic effects may otherwise 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.

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.

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 measurements 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₆₈ = (2.10 ± 0.18_stat)′. The extension measured with H.E.S.S. above 10 TeV, on the other hand, is r₆₈ = (0.82 ± 0.29_stat)′. In most models, as will be further detailed in the subsequent section, the optical radiation is produced by electrons that are also responsible for the IC emission in the Fermi-LAT energy range. In agreement with this, the position and extension of the emission measured with Fermi-LAT is roughly consistent with the outline of the nebula seen in the optical image. As the electrons radiating the X-ray synchroton 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% containment 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 modelling 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 incompleteness of the models.

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.

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 ± 0.2)″ from the pulsar (Weisskopf et al. 2012), which corresponds to r_s ≈ 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 γ_r, min and γ_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^−α. 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 σ and spin-down luminosity (see Appendix B for details).

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.

4.2. Fitting of the models

All three models give a prediction of the γ-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 γ-ray emission is slightly 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 γ-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 up-scattering, 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)⁷, to constrain the synchrotron flux predicted by our models (see also Fig. 1). This is done by computing a χ²-value between the predicted and observed synchrotron flux and extensions and adding this value as a penalty term to the total 𝒞-value of the fit of the IC spectrum,

$$\begin{array}{c}\hfill {\mathcal{C}}_{\mathrm{tot}}=-2ln{\mathcal{L}}_{\mathrm{tot}}={\mathcal{C}}_{\mathrm{IC}}+{\chi }_{\mathrm{SYN},\mathrm{flux}}^2+{\chi }_{\mathrm{SYN},\mathrm{ext}}^2,\end{array}$$(6)

where 𝒞_IC = −2lnℒ_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).

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{\mathrm{TS}}$, where TS ≡ −2ln(ℒ/ℒ_max) and ℒ_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).

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).

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 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{\mathrm{\Delta }{\mathcal{C}}_{\mathrm{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 = 𝒞_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₀(r/r_s)^−α, we find best-fit values of α = 0.47 and B₀ = 256 μG. The steepness of the radial profile of the magnetic field is thus compatible with the value of α = 0.51 ± 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 σ = 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 σ would lead to a high flow speed of ≳5000 km s⁻¹ 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\mathrm{pc}$, where r₀ ≈ 2 pc is the size of the nebula (Kennel & Coroniti 1984b). This is at odds with observed expansion and flow speeds of the order of ∼2000 km s⁻¹ and the shock radius at r_s = 0.13 pc. A value of σ compatible with these measurements is ruled out by our analysis with high significance, as it would lead to a much larger value of 𝒞_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 of distance to the centre. We find a value of σ ≈ 0.1 at the shock radius and ≈0.03 averaged across the whole nebula.

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 rs. For the variable B-field model B ∝ r−α, with α given in the legend. For the K&C model the magnetisation σ is specified in the legend.

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.

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.

The X-ray to γ-ray part of the synchrotron spectrum is not described well by any of the models. Specifically, the data between 10²⁰ − 10²² 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 γ_w, max, or constraints resulting from 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 parametrisation of the wind electron distribution. The first break in the spectrum, just beyond γ_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 (γ_w, min + γ)^s (compared to γ^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⁸. 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.

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.

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 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.

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.

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).

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.

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.

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 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 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 geo-magnetic cutoff was used to infer the absolute energy scale with an accuracy of ${}_{-5}^{+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 η between the energy scales is to evaluate the H.E.S.S. source model at a modified energy η ⋅ E, while the evaluation of the Fermi-LAT source model is unchanged (Meyer et al. 2010). Leaving η 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 η = 0.79 ± 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 η > 1. For the variable B-field model, we find η = 1.038, with only a very minor improvement of the fit (Δ𝒞_tot ≈ 2, implying that η 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 η = 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 of (1.9  ±  0.3)  ×  10⁻¹⁰ cm⁻² s⁻¹, which is consistent with the H.E.S.S. mono flux of (2.09  ±  02)  ×  10⁻¹⁰ cm⁻² s⁻¹ (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 χ² test, and determined the magnitude of a potential systematic error that would lead to an acceptable fit quality. For the spectral measurements in the γ-ray regime, we added a constant percentage of the flux in quadrature to the statistical uncertainty. To calculate the χ² values for the γ-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.48}^{+0.54}\prime$ (cf. Sect. 2.3) and at energies for which the H.E.S.S. data dominate we added an error of ${}_{-0.24}^{+0.21}\prime$ (cf. Sect. 2.2). Again, all of the systematic uncertainties have been added in quadrature to the statistical ones.

Table 4 summarises the χ² 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 χ² 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⁹. The χ² 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 χ² 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 synchroton 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 χ² ≈ 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.

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¹⁰. 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 σ 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.

Acknowledgments

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.

References

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.

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(γ, r)≡dN_el/dγdV, that is, the number of electrons per unit volume and Lorentz factor γ as a function of the distance to the nebula centre.

Given B(r) and n_el, we calculate the spectral volume emissivity $j_{\nu }^{\text{sync}}(\nu ,r)\equiv dE/dtdVd\nu d\mathrm{\Omega }$ for the synchrotron emission of a (pitch-angle averaged) chaotic magnetic field through (Blumenthal & Gould 1970; Aharonian et al. 2010)

$$\begin{array}{c}\hfill j_{\nu }^{\mathrm{sync}}(\nu ,r)=\frac{1}{4\pi }\frac{\sqrt{3}{e}^{3}B(r)}{m_e{c}^2}{\displaystyle \underset{1}{\overset{\infty }{\int }}}d\gamma n_{\mathrm{el}}(\gamma ,r)G\left(\frac{\nu }{{\nu }_c}\right),\end{array}$$(B.1)

where ν_c = 3eBγ²/4πm_ec is the critical frequency and G(x) is given by

$$\begin{array}{cc}\hfill G(x)& =\frac{x}{20}[(8+3x^2){\left(K_{1/3}(x/2)\right)}^2+\hfill \\ \hfill & x{K}_{2/3}(x/2)(2K_{1/3}(x/2)-3x{K}_{2/3}(x/2))],\hfill \end{array}$$(B.2)

with K_ξ the modified Bessel function of kind ξ. Throughout this paper, we use units of erg s⁻¹Hz⁻¹ cm⁻³ sr⁻¹ for j_ν.

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_{dust,  in} and r_{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

$$\begin{array}{c}\hfill j_{\nu }^{\mathrm{dust}}(\nu ,r)=\frac{3{\kappa }_{\mathrm{abs}}[\mathrm{\Theta }(r-r_{\mathrm{dust},\mathrm{in}})-\mathrm{\Theta }(r-r_{\mathrm{dust},\mathrm{out}})]}{4\pi (r_{\mathrm{dust},\mathrm{out}}^3-r_{\mathrm{dust},\mathrm{in}}^3)}{\displaystyle \sum _{i=1,2}}{M}_i{B}_{\nu }(T_i),\end{array}$$(B.3)

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

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{abs}}=2.15\times {10}^{-4}{\mathrm{cm}}^2{\mathrm{g}}^{-1}{\left(\frac{\lambda }{\mu \mathrm{m}}\right)}^{-1.3}.\end{array}$$(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 = (sync, dust), as a function of photon energy ϵ = hν is calculated through the integral (Atoyan & Aharonian 1996)

$$\begin{array}{c}\hfill n_{\mathrm{seed}}^t(ϵ,r)=\frac{r_0}{2hc}\frac{4\pi }{ϵ}{\displaystyle \underset{r_s/r_0}{\overset{1}{\int }}}dy\frac{y}{x}ln\frac{x+y}{|x-y|}{j}_{\nu }^t(\nu ,r_{0}y),\end{array}$$(B.5)

in units of photons eV⁻¹ cm⁻³ and where x = r/r₀ and y = r′/r₀, with r₀ the radius of the nebula and r_s = 0.13 pc (Weisskopf et al. 2012) the shock radius until which no emission is assumed.¹¹ The seed photon density for CMB photons is simply $n_{\mathrm{seed}}^{\mathrm{CMB}}=(4\pi /(hc))B_{\nu }(T_{\mathrm{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 synchtrotron and infrared dust emission, respectively.

The IC emissivity for up-scattered photons at frequency ν is then calculated using the integral over the IC kernel f_IC, which includes Klein-Nishina effects, and the seed photon density (Blumenthal & Gould 1970),

$$\begin{array}{cc}\hfill j_{\nu }^{\mathrm{IC}}(\nu ,r)& =\frac{3{\sigma }_{\mathrm{T}}hc}{4}\frac{h\nu }{4\pi }{\displaystyle \underset{1}{\overset{\infty }{\int }}}d\gamma \frac{n_{\mathrm{el}}(\gamma ,r)}{{\gamma }^2}\hfill \\ \hfill & \times {\displaystyle \underset{0}{\overset{\infty }{\int }}}dϵf_{\mathrm{IC}}(\nu ,ϵ,\gamma )\sum _t\frac{n_{\mathrm{seed}}^t(ϵ,r)}{ϵ},\hfill \end{array}$$(B.6)

with the usual expression for the full IC kernel,

$$\begin{array}{c}\hfill f_{\mathrm{IC}}(\nu ,ϵ,\gamma )=2qlnq+(1+2q)(1-q)+\frac{1}{2}\frac{{({\mathrm{\Gamma }}_{ϵ}q)}^2}{1+{\mathrm{\Gamma }}_{ϵ}q}(1-q),\end{array}$$(B.7)

where Γ_ϵ = 4ϵγ/(mc²) and q = hν/(Γ_ϵ(γmc² − hν)). The Thompson limit corresponds to Γ_ϵ ≪ 1 and IC scattering only occurs for 1/(4γ²)≤q ≤ 1.

The integrals are evaluated numerically and the code, fully written in python, is available online¹². The IC emissivity of the synchrotron emission involves four integrals for each frequency ν 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_{\nu }^{\text{sync}}$ and $n_{\mathrm{seed}}^{\mathrm{sync}}$.

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

$$\begin{array}{c}\hfill {\mathcal{L}}_{\nu }^t=\oint d\mathrm{\Omega }\int dV{j}_{\nu }^t(\nu ,r)={(4\pi )}^2{\displaystyle \underset{r_s}{\overset{r_0}{\int }}}dr{r}^2{j}_{\nu }^t(\nu ,r),\end{array}$$(B.8)

where the two factors of 4π 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_ν as a function of the angular distance θ 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_{\nu }^t$ is non-zero,

$$\begin{array}{c}\hfill I_{\nu }^t(\theta )=2{\displaystyle \underset{s_{\min }}{\overset{s_{\max }}{\int }}}ds{j}_{\nu }^t(\nu ,r(s)).\end{array}$$(B.9)

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

$$\begin{array}{c}\hfill r(s)=\sqrt{r_{\min }^2+s^2}.\end{array}$$(B.10)

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 rs (dark shaded region). We assume non-zero emission until the maximum radius of the nebula at r0. Here we illustrate two exemplary lines of sight s1 and s2 at angular offsets θ1 and θ2 with respect to the centre (dashed lines). s1 represents a line which intersects the shock whereas s2 does not. We integrate the emissivity jν 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.

In this expression, r_min = d sin θ is the smallest distance from the LOS to the centre of the Crab Nebula, which depends on the observing angle θ 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

$$\begin{array}{c}\hfill s_{\min }=\{\begin{array}{cc}\sqrt{r_s^2-r_{\min }^2}& r_{\min }\le r_s\hfill \\ 0& \text{otherwise}\hfill \end{array}.\end{array}$$(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 θ is large enough so that the LOS misses the shock radius. The upper integration bound s_max is set to the intersection with r₀ where the emissivities are negligible small.

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

$$\begin{array}{c}\hfill F_{\nu }^t({\theta }_0,{\theta }_1)={\displaystyle \int I_{\nu }}cos\theta d\mathrm{\Omega }=2\pi {\displaystyle \underset{{\theta }_0}{\overset{{\theta }_1}{\int }}}{I}_{\nu }cos\theta sin\theta d\theta .\end{array}$$(B.12)

The extension of our emission model is defined as the angle θ₆₈ which contains 68 % of the flux, that is,

$$\begin{array}{c}\hfill 0.68=\frac{F_{\nu }^t(0,{\theta }_{68})}{F_{\nu }^t(0,{\theta }_{\max })},\end{array}$$(B.13)

where the maximum angle θ_max is defined through tanθ_max = r₀/d.

For a full description of the emission model, what remains to be defined are the electron spectrum n_el(γ, 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,

$$\begin{array}{c}\hfill n_{\mathrm{el}}(\gamma ,r)=n_{\mathrm{radio}}(\gamma ,r)+n_{\mathrm{wind}}(\gamma ,r).\end{array}$$(B.14)

The spectrum of the radio electrons is a simple power law with index s_r between γ factors γ_r, min and γ_r, max with a super-exponential cutoff,

$$\begin{array}{cc}\hfill n_{\mathrm{radio}}(\gamma ,r)& =\frac{n_{r,0}}{{\rho }_r^3}{\gamma }^{-s_r}exp[-{\left(\frac{\gamma }{{\gamma }_{r,\max }}\right)}^{{\beta }_{\min }}]\hfill \\ \hfill & \times \mathrm{\Theta }(\gamma -{\gamma }_{r,\min })exp(-\frac{r^2}{2{\rho }_r^2})\hfill \\ \hfill & \times \mathrm{\Theta }(r-r_s),\hfill \end{array}$$(B.15)

where Θ(x) is the Heavyside step function. The spatial dependence of the radio electrons is given by a Gaussian function of constant width ρ_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

$$\begin{array}{cc}\hfill n_{\mathrm{wind}}(\gamma ,r)& =\frac{n_{w,0}}{{\rho }_w{(\gamma )}^3}exp[-{\left(\frac{{\gamma }_{w,\min }}{\gamma }\right)}^{{\beta }_{\min }}]exp[-{\left(\frac{\gamma }{{\gamma }_{w,\max }}\right)}^{{\beta }_{\max }}]\hfill \\ \hfill & \times [{\left(\frac{\gamma }{{\gamma }_{w,1}}\right)}^{-s_{w,1}}{\left(\frac{{\gamma }_{w,1}}{{\gamma }_{w,2}}\right)}^{-s_{w,2}}(1-\mathrm{\Theta }(\gamma -{\gamma }_{w,1}))\hfill \\ \hfill & +{\left(\frac{\gamma }{{\gamma }_{w,2}}\right)}^{-s_{w,2}}\mathcal{B}(\gamma ,{\gamma }_{w,1},{\gamma }_{w,2})\hfill \\ \hfill & +{\left(\frac{\gamma }{{\gamma }_{w,2}}\right)}^{-s_{w,3}}\mathrm{\Theta }(\gamma -{\gamma }_{w,2})]\hfill \\ \hfill & \times exp(-\frac{r^2}{2{\rho }_w{(\gamma )}^2})\hfill \\ \hfill & \times \mathrm{\Theta }(r-r_s),\hfill \end{array}$$(B.16)

where ℬ(x, a, b) is the boxcar function, ℬ(x, a, b) = Θ(x − a)−Θ(x − b). The energy dependence of the spatial extension is modelled through a simple power law

$$\begin{array}{c}\hfill {\rho }_w(\gamma )={\rho }_{w,0}{\left[{\left(\frac{\gamma }{9·{10}^5}\right)}^2\right]}^{-{\alpha }_w}.\end{array}$$(B.17)

For the free parameters γ_r, max, γ_w, min, and γ_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 β_min = 2.8 (this is the best-fit value obtained by Meyer et al. 2010) and β_max = 2.

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

$$\begin{array}{c}\hfill B(r)=B_0{\left(\frac{r}{r_s}\right)}^{-\alpha },\end{array}$$(B.18)

with α ≥ 0. The case of α = 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 σ at the shock, which is given by the ratio of the electromagnetic and particle energy flux,

$$\begin{array}{c}\hfill \sigma =\frac{B_s^2}{4\pi n{u}_s\gamma m_e{c}^2},\end{array}$$(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:

$$\begin{array}{cc}\hfill n_s(\gamma )& =n_{0,w}exp[-{\left(\frac{{\gamma }_{w,\min }}{\gamma }\right)}^{{\beta }_{\min }}]exp[-{\left(\frac{\gamma }{{\gamma }_{w,\max }}\right)}^{{\beta }_{\max }}]\hfill \\ \hfill & \times [{(1+\frac{\gamma }{{\gamma }_{w,\min }})}^{s_{w,1}}(1-\mathrm{\Theta }(\gamma -{\gamma }_{w,2}))\hfill \\ \hfill & +\frac{{(1+\gamma /{\gamma }_{w,\min })}^{s_{w,3}}}{{(1+{\gamma }_{w,2}/{\gamma }_{w,\min })}^{s_{w,3}-s_{w,1}}}\mathrm{\Theta }(\gamma -{\gamma }_{w,2}))].\hfill \end{array}$$(B.20)

Similar to the wind electron spectrum for the static models, β_min and β_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, γ′,

$$\begin{array}{c}\hfill \gamma (z)=\frac{{\gamma }^{\prime }}{{(v{z}^2)}^{1/3}+{\gamma }^{\prime }/{\gamma }_{\max }},\end{array}$$(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 σ, and γ_max is the maximum γ factor at distance z, see, for example, Eq. (6) in Atoyan & Aharonian (1996) (which also depends on vz²). The wind electron spectrum is then found to be

$$\begin{array}{c}\hfill n_{\mathrm{wind}}(\gamma ,z)={(v{z}^2)}^{-4/3}{\left(\frac{{\gamma }^{\prime }}{\gamma (z)}\right)}^2{n}_s({\gamma }^{\prime }).\end{array}$$(B.22)

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

$$\begin{array}{c}\hfill B_s=\sqrt{\frac{L_{\mathrm{spin}-\mathrm{down}}}{c{r}_s^2}\frac{\sigma }{1+\sigma }},\end{array}$$(B.23)

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

$$\begin{array}{c}\hfill B(z)=B_s\frac{3(1-4\sigma )z}{v{z}^2}.\end{array}$$(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 ‘simplex’ 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.

All Tables

All Figures

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).

In the text

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σ 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).

In the text

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° of the Crab Nebula’s position. The shaded regions indicate the intervals for the Bridge emission (0.09 ≤ ϕ ≤ 0.22, not relevant here) and Off region (0.51 ≤ ϕ ≤ 0.89) derived with the QB algorithm. The grey shaded areas indicate regions found by the HOP algorithm, which identifies “watersheds” between local maxima (Meyer et al. 2019). There are shown here for illustration only.

In the text

	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.
In the text

	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.
In the text

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.

In the text

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.

In the text

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.

In the text

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).

In the text

	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 rs. For the variable B-field model B ∝ r−α, with α given in the legend. For the K&C model the magnetisation σ is specified in the legend.
In the text

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.

In the text

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.

In the text

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.

In the text

	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.
In the text

	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.
In the text

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).

In the text

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 rs (dark shaded region). We assume non-zero emission until the maximum radius of the nebula at r0. Here we illustrate two exemplary lines of sight s1 and s2 at angular offsets θ1 and θ2 with respect to the centre (dashed lines). s1 represents a line which intersects the shock whereas s2 does not. We integrate the emissivity jν 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.

In the text