© The Authors 2023

1 Introduction

The Geminga pulsar (PSR J0633+1746) is a radio-quiet γ-ray source that was established as a pulsar in 1992, and at 250 pc is one of the closest pulsars to Earth (Bignami et al. 1983; Bertsch et al. 1992; Bignami & Caraveo 1992; Faherty et al. 2007). Searches for extended γ-ray emission around Geminga have been conducted ever since, proving unsuccessful for many years (Akerlof et al. 1993; Aharonian et al. 1999), until it was first detected at TeV energies by Milagro in 2007, with a diameter of 2.8° ± 0.8° (Abdo et al. 2007). This discovery was subsequently confirmed by the High Altitude Water Cherenkov (HAWC) collaboration, presenting evidence for significantly extended emission on angular scales of up to ~5.5° radius (Abeysekara et al. 2017b). However, a detection with imaging atmospheric Cherenkov telescopes (IACTs) remained elusive (Singh et al. 2009; Finnegan 2009; Ahnen et al. 2016). With a spin-down luminosity of Ė = 3.2 × 10³⁴ erg⁻¹, a spin period of P = 237 ms and a characteristic age of τ_c = 342 kyr, Geminga is one of the oldest pulsars around which extended very-high-energy γ-ray emission has been detected, providing evidence of energetic electron acceleration by middle-aged pulsars (τ_c ~ 0.1–1 Myr; Manchester et al. 2005)¹.

The morphology of the emission as seen with HAWC was found to indicate considerably slower diffusion than values typical for the interstellar medium (ISM), yielding diffusion coefficients a factor ~100 lower at 100 TeV (Abeysekara et al. 2017a). Geminga and the similar nearby pulsar PSRB0656+14 are considerably older pulsars, with longer periods and lower spin-down powers with respect to other pulsars associated with extended TeV γ-ray emission. It has been proposed that the γ-ray emission regions around these systems are in a different evolutionary stage to other TeV pulsar wind nebulae (PWNe) and they constitute a distinct source class of ‘TeV halos’ or ‘pulsar halos’ (Giacinti et al. 2020; Linden et al. 2017)². Within such halos, the γ-ray emission is due to inverse Compton (IC) scattering by electrons that have escaped from the PWN and the energy density of these electrons is lower than the ISM energy density, that is they do not dominate the surrounding region dynamically or energetically. As such, in contrast to the many PWNe detected with the High Energy Stereoscopic System (H.E.S.S.), with 12 firmly identified in the H.E.S.S. Galactic Plane Survey (H.E.S.S. Collaboration 2018), the detection of extended γ-ray emission around the Geminga pulsar with H.E.S.S. constitutes the first unambiguous detection of a pulsar halo at TeV energies by IACTs.

Extensive air showers (EAS) caused by cosmic rays are the primary background source of triggered events for a γ-ray analysis with IACTs. Despite state-of-the-art techniques in separating gamma-initiated from hadron-initiated EAS, there remains an irreducible background of gamma-like hadronic (mostly proton) showers – those in which a large fraction of the energy is transferred to neutral pions in the first interaction (Maier & Knapp 2007). As the rate of background events can vary based on the atmospheric conditions, observing direction, and hardware settings, methods to estimate and model the background level typically rely on counting gamma-like events within a region of the sky in the same dataset away from the target region (‘Off’), and subtracting this level from the region of interest (‘On’; Berge et al. 2007). In the case of the Geminga halo, the emission fills the field of view, such that there is no region free from γ-ray emission available for such a background estimation. This results in emission belonging to the halo being counted as background.

Previous IACT observations of the region have resulted in upper limits when probing angular scales in the range 0.1° to 0.3° (Finnegan 2009; Ahnen et al. 2016). Within the X-ray range a PWN has been identified with two lateral tails of length ~2–3′ and an axial tail of length 45″, motivating searches on these angular scales (Caraveo et al. 2003; Pavlov et al. 2010; Posselt et al. 2017). Analysis and detection of emission on larger scales is challenging for IACTs, for reasons outlined above (Berge et al. 2007). Dedicated analysis approaches developed for this sky region include matching observation conditions between observations of this region and of empty sky fields, in order to estimate the background across the full region of interest (Flinders 2015; Abeysekara 2019; Hona 2021).

As shown in a systematic study of differences between the analysis procedures of the water Cherenkov detector facility HAWC and the H.E.S.S. IACT array (Abdalla et al. 2021), the two experimental techniques are complementary approaches, with HAWC providing good sensitivity to high energies (E ≳ 10 TeV) and extended emission regions. In contrast, H.E.S.S. can provide detailed morphological and spectral studies, thanks to its few-arcmin angular and 10–15% energy resolution. Given the extent of the significant emission detected with HAWC around the Geminga pulsar (out to ~5° radius), it is not currently possible for IACTs to measure the full extent of the emission using standard background estimation techniques. However, as we also show in the analysis presented here, it is nevertheless possible for IACTs to significantly detect degree-scale extended emission whilst controlling the estimated background, and to start probing the morphological and spectral properties of the inner regions of such sources (Aharonian et al. 2022).

Searches for spectral variation, energy-dependent morphology or asymmetries to the emission are of particular interest to investigate properties of the large scale pulsar halo. Similarly, investigating the transport of energetic particles in the vicinity of the pulsar is important to establish whether the properties are consistent with suppressed diffusion as seen in previous analyses.

2 H.E.S.S. data and analysis

2.1 H.E.S.S. observations

H.E.S.S. is an array of five IACTs, located in the Khomas Highland of Namibia at 1800 m above sea level. Four of the telescopes (CT1–4) have mirror areas of 107 m² whilst the fifth (CT5) has a mirror area of 612 m² and correspondingly lower energy threshold (Holler et al. 2015). Due to its smaller field of view, of 3.2° (compared with the 5° field of view for the smaller telescopes), CT5 is not used in this analysis. The electronics of the cameras of the CT1–4 telescopes were upgraded in 2017 (Ashton et al. 2020). Observation data are collected with H.E.S.S. in ‘runs’ typically 28 min in length, during which telescopes are pointed towards a specific sky position.

H.E.S.S. observed the Geminga region using ‘wobble’ mode offsets of 0.5° and 0.7° around the Geminga pulsar (RA 06h33m54s, Dec + 17°46′13″) over two seasons in 2006-2008, for a total of 14.2 h of livetime (see Table 1; Aharonian et al. 2006). With these default wobble offsets, no γ-ray emission has been detected when using background estimation techniques probing radii of ≲0.3°. Further observations were taken during the first quarter of 2019, employing a pointing strategy with much wider pointing offsets of ±1.6° from the pulsar in RA and Dec, for a total livetime of 27.2 h. These offsets are among the widest ‘wobble’ mode pointing offsets yet used with IACTs, thereby pushing the limits of standard analysis techniques.

2.2 Analysis

Considerable advances in gamma-hadron separation were achieved by the use of template or model-based approaches. In this analysis, a sensitive likelihood-based template fitting analysis is used (Parsons & Hinton 2014) and the results are cross-checked using an independent calibration and analysis chain (de Naurois & Rolland 2009).

The Geminga region is observable with H.E.S.S. from November to April only at large zenith angles >40°, which raises the intrinsic energy threshold of the analysis. Additionally, depending on the method of background estimation and due to systematic uncertainties at the level of the statistical uncertainties or higher (caused for example by variation in the atmospheric conditions Hahn et al. 2014) we use a conservative set of selection cuts on reconstruction quality, known as hard selection cuts, such that only the best candidate γ-ray events are retained. This also has the effect of improving the signal-to-background ratio such that the uncertainty on the background normalisation is reduced. Correspondingly, the energy threshold is 1 TeV for the spectral analysis; whilst for the morphological analysis, we used an energy threshold of 0.5 TeV. This is due to the higher event reconstruction quality and an energy bias ofless than 10% required for the spectral analysis.

3 Detection of extended γ-ray emission

A systematic study of analysis differences and background estimation techniques between H.E.S.S. and HAWC was conducted and applied to the Galactic plane (Abdalla et al. 2021). By adapting the analysis of H.E.S.S. data to techniques suitable for angular scales comparable to those probed by HAWC, it was found that the detectability of large, extended sources could be improved. This enabled the measurement of γ-ray emission from a very extended region and, although the background analysis remains challenging, the H.E.S.S. data retains the advantage of being able to identify smaller scale structures, including potential point-like sources.

As the γ-ray emission is considerably more extended than most other single TeV γ-ray sources detected by IACTs to date, a variety of approaches are used to verify the detection and nature of the emission. Particular care needed to be taken with the background estimation, which is described in detail in Sects. 3.1 and 3.2. It is not yet possible to evaluate the true extent of the TeV γ-ray emission, as it extends beyond the sky region available with the current dataset.

Using an integration region of 1° radius, extended emission around the Geminga pulsar is detected at >6σ (evaluated using Li & Ma 1983) in the 2006–2008 dataset and similarly at ~9σ in the 2019 dataset. Due to the differences in wobble offset between the 2006–2008 and the 2019 observations, different background estimation techniques are used for the two datasets.

3.1 Background estimation: 2006–2008 dataset

Revisiting the Geminga observations taken with H.E.S.S. in 2006–2008, the detection of emission over a much larger angular scale is enabled by increasing the exclusion region³ radius from ~0.4° to 1.5°. The reflected region background method is well suited to spectral analysis, in which a region of the same size as the ON region yet reflected across the telescope pointing direction is used to estimate the background. However, this is not possible for cases where the ON region size (here of 1° radius) exceeds the pointing offset (here of ±0.5°−0.7°, see Table 1). Additionally, due to the large exclusion regions, a spectral analysis of the emission using the reflected background method is not possible on the Geminga region for any of the current H.E.S.S. datasets (Aharonian et al. 2001).

Employing the ring background method with a fixed ring thickness of 0.5° and a minimum radius in excess of 1.5°, for a 250 pc distance to Geminga, this corresponds to background estimation from radii >6.5 pc away from the pulsar (Berge et al. 2007). As the background counts are sampled from the same region of the sky in the ring background method, standard (std) selection cuts are used, whereas more conservative (hard) selection cuts are used for other methods of background estimation. With such a ring background subtraction, H.E.S.S. is sensitive to the difference in γ-ray emission between the ON region at small radii from the pulsar, and the emission at larger radii that is used to estimate the background. The limited field of view prevents using a background ring radius large enough such that the background region does not contain emission from the source. The measured flux is therefore a relative measurement and will underestimate the true flux (see also Mitchell et al. 2019).

Two different background methods, Ring and On-Off, are applied to evaluate the significance of γ-ray emission within a 1° radius region around the pulsar in this dataset, with a consistent detection obtained using both approaches (see Table 2). When using the On-Off background estimation, the entire field of view of the observations is considered as the On region, with Off data taken from observing runs matched for comparable conditions. The parameters used to match the Off data to On data include the zenith angle of the observations, the run duration, and the combination of telescopes participating. The presence of all four telescopes CT1–4 is required for this analysis in order to provide a smooth acceptance⁴ across the sky region. The Off runs used to estimate the background are extragalactic observations from 2004–2009 with no significant source detected in the field of view.

3.2 Background estimation: 2019 dataset

Following this detection, observations of the Geminga region were conducted in 2019 at large wobble offsets of 1.6° around the pulsar, out to which the camera acceptance remains at ≳60% of the on-axis acceptance (Aharonian et al. 2006). Once again, the emission from the region around Geminga is found to fill the available field of view, proving a considerable challenge for background estimation in the analysis. Therefore, two different methods are used; the so-called On-Off and field-of-view background approaches (Weekes et al. 1989; Berge et al. 2007).

In order to provide an estimate of the systematic uncertainties of the On-Off background estimator and cross-check our results, we used two independent matched lists of Off runs, with different tolerance levels in the matching criteria (see Table 3). Both lists are comprised of observations from 2017 to 2019. Unless otherwise specified, analysis results from the On-Off background approach are always presented using the better matched Off List 1, with Off List 2 contributing to the evaluation of the systematic uncertainty only (see Sect. 4.2).

The second background estimation approach used is the field-of-view (FoV) method, which uses an acceptance model to estimate the expected level of background counts throughout the field of view (Berge et al. 2007). Regions of significant γ-ray emission are excluded and the predicted number of background counts is normalised to the excess counts outside of these exclusion regions.

For both background estimation approaches, the background counts are normalised to the On counts in the region at radii θ > 3.2° away from the pulsar location (that is twice the offset of the observation positions), as shown in Fig. 1. Excess counts maps constructed using both background estimation approaches are shown in Fig. 2, using a 0.08° correlation radius. The observation positions used in 2006–2008 and in 2019 are indicated on the sky map shown in Fig. 3.

Some evidence for run-by-run variation in the Off count rate was investigated and found to be related to the atmospheric conditions, as indicated by the Cherenkov Transparency Coefficient, CTC (Hahn et al. 2014). The Pearson correlation coefficient for the Off count rate and CTC is ~0.7 for both Off Lists. A correction for the variation in atmospheric conditions between the On and Off runs is implemented via two methods; using the ratio of the CTC between On and Off runs; and using the ratio of the number of hadron-like events in On and Off runs. Both methods are found to remove the correlation between the off count rate and the CTC. This correction is applied prior to the normalisation of On counts to radii >3.2°. For the remainder of this paper only this dataset obtained in 2019 is used.

Fig. 1 Surface brightness profile of emission around the Geminga pulsar measured with HAWC (Abeysekara et al. 2017a). Shaded regions indicate the parts of the emission profile that are used as ON (radius θ ≲ 1.0°) and OFF normalisation (radius θ ≳ 3.2°) regions to estimate the background level and evaluate the significance in the analysis of the 2019 H.E.S.S. dataset. The region shown for normalisation of the OFF counts is only accessible with the 2019 dataset, due to the wider pointing strategy used.

3.3 Background estimation: Systematic tests of the On-Off and FoV methods

As the analysis of sources with highly extended γ-ray emission is inherently challenging for IACTs, we performed a range of tests to assess the robustness of the analysis techniques and to quantify systematic uncertainties. We performed an On-Off background analysis on two extragalactic regions with no expected γ-ray emission, using comparable parameters to the analysis of the Geminga region, that is with the same set of matching criteria for Off runs. For this purpose we used data obtained in 2017 and 2018 with the telescopes CT1–4 on two Dwarf Spheroidal galaxies, namely Reticulum II and Tucana III (Abdallah et al. 2020). The average zenith angles of these datasets are 42.4° and 39° respectively, comparable to the 43.5° average value of the Geminga dataset. No evidence for significant emission within the target region is found. The significance distributions for both dwarf spheroidals are found to be compatible with background, whereas the significance distribution for the Geminga region is asymmetric and exhibits a clear tail towards high significance. Also, the mean of the distribution is greater than zero for the Geminga region, reflecting the fact that the emission fills the field of view, whereas the distribution mean is compatible with zero for the test fields. This lack of signal from the test fields associated with dwarf spheroidal galaxies when using the same method provides confidence in the clear signal seen from the Geminga region. These empty fields are also analysed using the field-of-view background method; again, no significant emission is seen.

Figure 4 shows the ratio of On counts to background counts with distance from the FoV centre, with renormalisation to the background at a radial distance of >3.2° from the centre of the field of view. There are indications for a bias in the Ret II dataset with an On counts to background ratio lower than 1.0 out to beyond 3° (that is out to the radius which is used for normalisation to background). However, for the Tuc III dataset, the background normalisation results in approximately the expected level, averaging around 1.0. Correcting for a ~10% bias for the region within 1° radius of the pulsar would lead to a ~30% increase in significance for the values quoted in Table 2. Given the limited number of tests performed, we conservatively estimate a systematic error of the order of 30%. A similar bias is found using the FoV background method, averaging ~6% at radii <1°.

The field of view around the Geminga pulsar contains a bright star of magnitude 1.9, that leads to increased Night-Sky-Background (NSB) in the region. We investigated whether any dependence of the excess counts on the NSB in an On-Off analysis is seen, and found no correlation between the NSB rate and excess counts per run.

Fig. 2 Excess counts sky maps of the Geminga region above 0.5 TeV using the 2019 dataset with two different background estimation methods. Left: On-Off background estimation method. Right: field-of-view background estimation method. The maps are over-sampled with a 0.08º correlation radius, with contours at the level of 50 and 100 excess counts from maps over-sampled with a 0.5° correlation radius. The 68% containment PSF is shown for comparison and has a value of 0.06º, valid for the innermost θ < 1° region around the pulsar in which the significance is evaluated (indicated by a dashed circle). The Galactic plane is indicated by a green dashed line and the pulsar location with a green triangle, whilst an arrow indicates the proper motion direction of the pulsar (arbitrary length). Background normalisation is performed using regions at θ > 3.2°, indicated by a dotted circle.

Fig. 3 Observation positions corresponding to the 2006–2008 data set (magenta) and the 2019 dataset (green) at offsets of 0.7° and 1.6° from the pulsar respectively. The white dashed circle indicates the test On region with 1° radius around the pulsar, whilst the white dotted circle indicates the 3.2° radius region beyond which the background is normalised. Green and magenta circles indicate the radius around an observation position with the same offset as the pulsar. Counts map and contours are otherwise as in Fig. 2, left.

Fig. 4 Ratio of On counts to background counts using the On-Off background method as a function of radial distance from the centre of the field of view. Shown are two empty sky regions taken on the dwarf spheroidal galaxies Reticulum II and Tucana III, as well as Geminga data.

4 Analysis results

4.1 Morphology of the γ-ray emission

Figure 5 shows the On data and background counts as a function of radial distance from the pulsar. The estimated background from the FoV and On-Off background methods (using two different Off run lists) are shown, and normalised to the data at radii >3.2°, with a statistical uncertainty on the normalisation at the 1% level. An excess of On data above background counts can be seen for all three estimations of the background level, particularly towards the innermost radii. The levels of background counts from the different methods are broadly consistent, with the remaining mild discrepancies providing an indication of the level of background systematics, estimated as ≤5%. Figure 5 provides further confidence that the detected emission is significant above the level of Galactic diffuse emission, which would peak towards the Galactic plane (the Geminga pulsar being situated towards the Galactic anti-centre at a galactic latitude of 4.27°).

The extension of the γ-ray emission is of particular relevance to studies of energetic particle transport in the region. Therefore, we tested the extent of the region of significant emission, firstly by varying the integration region radius (the radius around the pulsar within which the significance of the γ-ray emission is computed). Once the significant γ-ray emission is fully contained, the significance curve will start to flatten with increasing radius.

Figure 6 shows the significance with radius for the two independent background estimates. Notably, neither of the curves flatten within a radius out to 3°, indicating that the significant γ-ray emission is not yet fully contained by the integration region. Secondly, the shape of the curves is consistent between the FoV and On-Off background methods; this is reassuring confirmation in the extent of significant emission provided by these fundamentally different approaches. Lastly, this analysis shows that the emission within a radius <0.25° around the pulsar is not significant. This further confirms that although concentrated towards the pulsar, the majority of the emission is spread out to larger radii.

Although the analysis is limited by systematic uncertainties, we attempt to quantify the significance of a possible asymmetry as follows. Firstly, an acceptance-corrected azimuthal profile of the emission within a 3° radius around the pulsar is constructed. The peak of the azimuth profile is found to be at 105° ± 3° measured anticlockwise from the north, and is used to divide the region into two halves. From each semicircular half, radial surface brightness profiles of the emission, corrected for acceptance, are constructed in the two independent regions. Indications for a higher level of emission in one region compared to the other are seen at radii >1° (for example Fig. 2), however this is found to be neither statistically significant nor independently verified in the cross-check analysis.

To quantify whether the emission is significantly offset from the pulsar, despite the inability to measure the true extent, we evaluated the barycentre and 68% containment radius of the excess counts, with acceptance correction applied. Indications for an offset of the emission centroid from the pulsar are found for all background methods, yet are nevertheless compatible with the pulsar location when systematic errors are taken into account. The emission centroid has an offset of 0.6°, at RA 99.1° ± 0.1° ± 0.5° and Dec 17.7° ± 0.1° ± 0.5° where errors are statistical and systematic (as estimated from the difference between the background methods). Evaluating the containment radius in three energy bands (0.5–2 TeV, 2–8 TeV, and 8–40 TeV) no evidence for significant energy-dependent morphology is found. These containment radii in energy bands are summarised in Table 4.

Although TeV-bright γ-ray emission may be expected on angular scales ≲0.1° corresponding to the size of the X-ray PWN, with this analysis we see no indications of a separate component at these scales. Identifying such a separate emission component is challenging given the evolved state of the system, the proper motion of the pulsar (which is fast compared to the cooling time of electrons at the nearby distance of Geminga) and the predominance of inverse Compton emission from electrons already escaped from the PWN region.

Fig. 5 Radial profile without acceptance correction applied from uncorrected maps, showing counts per unit area with radial distance from the pulsar at energies >0.5 TeV, error bars are statistical. The different background methods show reasonable agreement.

Fig. 6 Significance with integration region radius around the pulsar for FoV and On-Off background estimation analyses. Emission within 1 degree radius of the pulsar is detected with a significance of ≳8σ with all analyses.

Fig. 7 Spectral energy distribution of the γ-ray emission within 1° radius of the Geminga pulsar. Spectra from two HAWC analyses of Geminga are shown for comparison (Abeysekara et al. 2017a,b), with flux normalisation scaled to match the sub-region from which the H.E.S.S. spectrum is extracted. Error bands include systematic uncertainties for H.E.S.S. but are statistical only for HAWC.

4.2 Spectral results

We perform a spectral analysis for a circular region with 1º radius centred on the pulsar using the On-Off background approach; the resulting spectrum is shown in Fig. 7 with parameter values for a best fit power law quoted in Table 5. Due to the quality cuts applied for a spectral analysis, the energy threshold of the analysis is 1 TeV. The spectrum obtained with HAWC for both the disc morphology of Abeysekara et al. (2017b) and the diffusion model of Abeysekara et al. (2017a) are also shown in Fig. 7, scaled for comparison to account for the different sizes of the regions probed. The scaling factor for the HAWC spectrum from Abeysekara et al. (2017b) corresponds to the ratio of the disc areas for a 2º radius and a 1º radius, the region probed by this H.E.S.S. analysis. For the HAWC spectrum from Abeysekara et al. (2017a), we rescale according to the ratio between the diffusion model integrated out to 1º radius and the full integral, assuming that the spectral index is not radially dependent.

Figure 7 shows that the scaled HAWC flux and the H.E.S.S. measurement are consistent at energies ≳5 TeV and compatible within the systematic errors of both measurements (Abeysekara et al. 2017a,b). From Abeysekara et al. (2017a), the best fit spectral measurement has a spectral index of 2.34 ± 0.07 and flux normalisation of TeV⁻¹ cm⁻² s⁻¹ at 20 TeV. We use this latter spectrum obtained from a dedicated analysis of HAWC data on the Geminga region for modelling purposes. The systematic errors of the H.E.S.S. spectral measurement are ~0.2 on the index and ~30% on the flux, obtained from the cross-check between the background models. The ~8% bias seen at radii <1.0° in Fig. 4 corresponds to a ~30% systematic on the flux, assuming these two contributions are uncorrelated, the total systematic error on the flux is of the order of 40%. The systematic errors of the HAWC spectral measurement are 0.2 on the index and 50% on the flux. We note that the measured HAWC flux is dependent on the morphological model assumed for the emission – here we show the disc model, the approach most similar to the H.E.S.S. spectral measurement, which is also approximately a disc with 1º radius.

Fig. 8 Radial profile constructed using three independent methods; from γ-ray maps with FoV and On-Off background methods (Off list 1), and from a spectral analysis in consecutive ring regions around the pulsar. The latter is limited to radii <2° from the pulsar and is shown without an additional normalisation to background at radii >3.2º applied.

4.3 Radial profiles

Acceptance-corrected radial profiles of the excess emission are constructed from the sky maps, using a circular region centred on the pulsar extending out to ~3.5° radius. The surface brightness profiles are constructed from acceptance-corrected excess maps with units of cm⁻² s⁻¹ deg⁻². For the conversion to surface brightness in units of TeV cm⁻² s⁻¹ deg⁻², a power law spectral model is used with index Γ = 2.76 (corresponding to the best-fit spectral parameters, see Table 5). These profiles are presented in Fig. 8.

Using the On-Off background method, spectral analyses are additionally performed in concentric rings of 0.5° width, out to a radius of 2°. The data is fit in each ring using a power law spectral model, from which radial profiles of the emission are constructed. The best fit model is integrated within the energy range 0.5–40 TeV to obtain the flux as a function of radius within a specific energy band. This radial profile is also shown in Fig. 8. The resulting radial surface brightness profiles are therefore found to be compatible between the normalised profiles from maps constructed with On-Off background and from maps constructed with the FoV background. The radial surface brightness profile from a spectral analysis is however limited to radii <2°. In order to normalise this radial profile constructed from spectra, we applied the same normalisation as that obtained for the radial profile for On-Off background constructed from sky maps.

These radial profiles also show the normalisation to background at radii >3.2°, which implies that the flux is likely underestimated over the region. The innermost radial bin corresponds to a region 0.2° around the pulsar, comparable in size to that of the region probed by X-ray instruments. A flux measurement is obtained from the surface brightness within this radial bin, assuming that the spectral index does not vary with radius. In Sect. 5, these radial profiles, together with the spectral energy distribution, are fit with a diffusion model.

5 Modelling

5.1 Diffusion model description

For the particle transport we consider energy-only dependent diffusion, energy losses, and a continuous point-like source of electrons: (1)

with N(E, r, t) the density of electrons per energy, D(E) the spatial diffusion coefficient, b(E) the energy loss rate, and Q(E, t)δ(r − r_s) the source term with Q(E, t) the rate of produced electrons per energy and r_s the source position. The assumption of a point-like source of injection is motivated by X-ray observations of the Geminga pulsar wind nebula (Pavlov et al. 2010) showing emission up to 0.06 pc whereas the observed TeV emission extends up to 14 pc. An advection term is not included in this expression for the particle transport, as this would generate asymmetries to the γ-ray emission that are not seen in this analysis.

5.1.1 Diffusion coefficient

The diffusion coefficient D(E) = D₀(E/E₀)^δ is assumed to follow a power law, with index δ in the range δ ∊ [0.3, 1], to cover the scale between Kolmogorov turbulence (Kolmogorov 1991), Kraichnan turbulence (Kraichnan 1965), and Bohmian diffusion (Bohm 1949). For this model, E₀ = 100 TeV as in Abeysekara et al. (2017a) is chosen to have a direct comparison of the diffusion coefficient with the HAWC results, whilst D₀ is kept as a free parameter of the fit.

5.1.2 Energy losses

Synchrotron and Inverse Compton energy losses are considered. We assume the presence of three different photon fields; the Cosmic Microwave Background, Infra-Red from dust, and Optical radiation from star light (see also Appendix A.1). Magnetic field values between 1 µG and 5 µG are tested, which correspond to inverse-Compton dominant and synchrotron dominant scenarios respectively.

The typical time of energy losses is defined as: (2)

where b(E′) represents the total energy-dependent energy losses. Eq. (2) can be defined as the maximum age of an electron given its current energy. This simple calculation enables us to assess the age of the particles that can be observed, as shown in Fig. 9. Given the H.E.S.S. energy threshold converted to electron energy, the maximal age of particles is between 140 kyr for the 1 µG case and 32 kyr for the 5 µG case.

Fig. 9 Typical energy loss time versus electron energy. If the particle loss timescale is longer than the pulsar age (here chosen as 340 kyr), then the typical time tE(E), (see Eq. (2)) is taken as the pulsar age. The difference to the HAWC model is due to the approximation of energy losses made in the HAWC model.

5.1.3 Source term

The total energy released by the pulsar since its birth, W₀, is defined as: (3)

where T_* is the current pulsar age and L is the pulsar luminosity function (see Appendix A.2). We assume that the source term Q(E, t), (defined as the rate of electrons per energy produced by the pulsar) follows a power law with an exponential cut off, and its normalisation follows the same time-dependence as the pulsar luminosity: (4)

where τ₀ is the initial spin-down timescale, n is the braking index, and E_c the energy cut-off of electron spectra (see Appendix A.2 for the definitions of τ₀ and E_c). For simplicity, E_c is assumed to be constant during the lifetime of the pulsar. The normalisation of the source term can be related to the total energy released by the pulsar by: (5)

where the efficiency parameter η is the proportion of the pulsar spin-down energy converted to electron-positron pairs that escaped from the close source environment to the extended emission region.

5.2 Solution

The solution to the diffusion energy loss equation as described in Di Mauro et al. (2019) is: (6)

where E_s(E, t₀, t) is defined as the energy of an electron that was produced at time t₀ and observed at time t with energy E. This value is obtained by solving numerically: (7)

while the diffusion radius λ(t₀, t, E) is defined as: (8)

Since the data we are using to constrain this model are radial profiles, we compute the mean density over the angle between r and r_s of electrons observed today at a given energy and distance from the current pulsar position r can be written as: (9)

with T_* being the current pulsar age, where E_s(E, t₀, T_*) is defined as the energy of an electron that was produced at time t₀ and observed at time T_* with energy E. Since r_s(t) is varying as a function of time, the impact of the asymmetric morphology is taken into account in our radial profiles. In addition to the source parameters, the density of electrons depends on the diffusion parameters (D₀ and δ) and electron cooling parameters (B field and IC photon field energy density and temperature). The parameter r_s allows us to introduce the effect of the pulsar proper motion, and is defined as: r_s(t) = V_T(T_* − t) where V_T is the transverse speed of the pulsar.

As shown in Fig. 10, the diffusion radius peaks at the energy where the cooling time for electrons is equal to the pulsar age. For the Geminga pulsar with a characteristic age of T_c = 342 kyr, this occurs at ~ 1 TeV. At lower energies, the electrons are uncooled and the diffusion radius λ is limited by the diffusion parameters; whereas at higher energies, the diffusion radius is limited by the cooling time for electrons in the ambient magnetic field. Figure 10 demonstrates that the diffusion radius in the energy range of H.E.S.S. is much larger than the H.E.S.S. field of view, such that the full extent of the emission cannot be measured.

5.3 Model joint fit with H.E.S.S., HAWC and XMM-Newton data

The solution of the diffusion model provides us with the density of electrons per energy interval. Two more steps are needed to compare the diffusion model with the multi-wavelength data: Firstly, the electron density needs to be projected onto a 2D plane and converted to a density versus angle based on the pulsar distance. Secondly, the obtained electron density profile is then converted to the γ-ray flux through Inverse Compton (IC) scattering. This step makes use of the Naima package (Zabalza 2015), with a photon field for IC approximated as three black-body components defined in Appendix A.1.

The diffusion model described above is simultaneously fit to radial profiles and spectral energy distribution (SED) using a chi-square minimisation. The normalisation of the diffusion coefficient D₀ is used as a free parameter while the energy cut-off of the electron spectrum E_c is tested as a free or fixed parameter in order to evaluate its significance in the final result. We performed this analysis using the combination of H.E.S.S. (radial profile and SED within 1°) and HAWC 2017 (radial profile) data (see Sect. 4.2). The HAWC SED from Abeysekara et al. (2017a) is not used in the fit but is superimposed to the model SED obtained with an integration radius of 10°. The GeV part of the Geminga halo is assumed to be highly offset compared to the actual pulsar position due to the proper motion of the pulsar. A detection of such emission is claimed by Di Mauro et al. (2019), where a template fit, including the effect of pulsar proper motion, enabled the detection of emission between 10 GeV and 100 GeV with an offset position and extent consistent with that expected from the diffusion model fit to HAWC data.

We decided not to include these data because our single zone model is not valid on the full extent of the claimed Fermi-LAT emission. Indeed, a single zone of diffusion for the Fermi-LAT extent scale cannot take into account the positron excess of AMS-02 as explained in Martin et al. (2022). A two zone diffusion model is beyond the scope of this paper, since the emission measured with H.E.S.S. is not extended enough to constrain such a model. In order to constrain the synchrotron contribution to the SED, we compared keV upper limits to the model SED, derived using XMM-Newton data in a 10′ region around the Geminga pulsar as described in Liu et al. (2019).

In order to explore the phase space whilst keeping D₀ and E_c free in the fit, we conduct a parameter scan over five variables (n, η, α, B, δ), because a global minimisation procedure is found to be degenerate. These scanned parameters are listed in Table 6. All possible combinations of parameters are tested, leading to 243 different combinations to explore the full parameter space. Values for the parameters characterising the Geminga pulsar are taken from the ATNF pulsar catalogue Manchester et al. (2005). The initial period of the pulsar is assumed to be 15 ms since its value is not affecting results above the H.E.S.S. energy threshold. Two of the free parameters are further constrained by validity conditions: P₀ < P_* and η < 1.

The fits are performed for the results using both On-Off background lists separately, as well as for the main and crosscheck analyses. For each fit result, the p-value is computed and is used to define the goodness of fit. A p-value >0.003, corresponding to three standard deviations is taken, as a threshold to define a good fit. This can be considered as a fairly relaxed threshold since the systematic uncertainties of both H.E.S.S. and HAWC data are important.

The p-value distribution as a function of the electron efficiency from the fixed E_c scan is presented in Fig. 11, where only those good fit models with p-value >0.003 are retained. The combination of H.E.S.S. and HAWC data excludes an efficiency of injection and electron-positron pair conversion lower or equal to 1% at the 5σ level. However, no significant preference is found between the values scanned for the power law index of the diffusion coefficient δ, the ambient magnetic field B, injection index α and braking index n.

The fitted diffusion coefficients are comparable with the one obtained by HAWC (Abeysekara et al. 2017a) as shown in Fig. 12. The diffusion coefficients derived from the Boron-to-Carbon ratio (labelled B/C diffusion hereafter) are obtained under the assumption that the galactic magnetic halo size is 1 kpc, 4 kpc, and 16 kpc respectively for the three values indicated (Genolini et al. 2019). The H.E.S.S. result confirms a significantly lower diffusion coefficient than that obtained with the cosmic ray Boron-to-Carbon ratio. Moreover, the mean statistical uncertainty of the fits improved from 27% with HAWC only to 17% with the combination of H.E.S.S. and HAWC.

A comparison is performed between models both with an energy cut-off fixed to 1 PeV and when leaving the energy cut-off free to vary. This value is chosen as the minimal cut-off value that does not significantly modify the spectral shape in the energy ranges of H.E.S.S. and HAWC. The cut-off is defined as significant if the model fit leaving E_c free is able to provide a good fit (p-value >0.003) and if the chi-square difference with the E_c fixed model is found to be higher than 25 (which corresponds to 5 standard deviations).

For the H.E.S.S. and HAWC combined data, 53 models passed this criterion showing preference for a sub-PeV cut-off. Figure 12 shows the fitted energy cut-off compared to the power law index of the injection spectrum and the fitted diffusion coefficient for good fits obtained with an energy cut-off either fixed to 1 PeV or left free to vary. A correlation can be observed between the spectral index and the energy cut-off in the top panel of Fig. 12, with a Pearson correlation coefficient of 0.64. A hard injection spectrum, with index lower than 2, implies the presence of a sub-PeV energy cut-off to reproduce the combined H.E.S.S. and HAWC dataset.

The expected synchrotron flux in the 10′ region around the pulsar from the good fit models is compared with the XMM-Newton upper limits with both a free cut-off and a fixed 1 PeV cut-off. For the models defined as a good fit of the H.E.S.S. and HAWC data, those with a fixed E_c of 1 PeV result in a synchrotron flux exceeding the upper limit. This leads to their exclusion, while for the free cut-off, the models with a magnetic field of 1 µG and a high-energy cut-off lower than 75 TeV are able to fulfil the XMM-Newton constraint. The impact of these two parameters on the SED is presented in the Fig. 13. A value higher than 1 µG overshoots the X-ray upper limits, as well as a higher injection cut-off energy than 75 TeV. A flux point for H.E.S.S. derived from the innermost radial bin of Fig. 8, a comparable angular scale of ~10′ is also shown for comparison. We note that due to the background normalisation applied, the absolute flux value is likely underestimated.

As an example, the best fit from all the scanned parameters with HAWC and H.E.S.S. data in terms of p-value is shown in the Fig. 14 compared with data. The list of scanned models that fulfiled the good fit and XXMM-Newton constraints is given in the Appendix B.

Fig. 10 Comparison of the predicted typical diffusion scale λ as a function of electron energy for different diffusion models and magnetic fields. The HAWC model is taken from Abeysekara et al. (2017a). It is clear that the emission is predicted to exceed the H.E.S.S. field of view over most of the H.E.S.S. energy range (here converted to electron energy).

Fig. 11 Distribution of p-value for the scanned parameter combinations presented in Table 6, with different values of η parameter, using the H.E.S.S. and HAWC datasets in combination with Ec fixed to 1 PeV. This distribution is presented as violin plot, illustrating kernel probability density. The width of the shaded area represents the proportion of the number of fits located at these position. The upper and lower bar represent the maximum and minimum values. Models with p-value <0.003 are excluded.

6 Discussion

This H.E.S.S. detection of extended TeV emission around the Geminga pulsar confirms the results reported by the water Cherenkov detectors Milagro, HAWC, LHAASO, and Tibet-ASγ with IACTs for the first time. Demonstrating this detection with IACT analysis techniques, potentially opens up the source class of highly extended ‘halos’ of TeV γ-ray emission around middle-aged pulsars to detailed studies by IACTs (Linden et al. 2017; Giacinti et al. 2020). In contrast to the known population of TeV pulsar wind nebulae, for which the TeV structures are already known to be typically larger than the X-ray synchrotron components (Aharonian et al. 1997), in pulsar halos, it is thought that the emission originates from escaped particles diffusing away from the canonical PWN.

Geminga is unambiguously in the halo stage of PWN environmental evolution and as such is the first clear example of extended halo-like TeV emission to be detected with H.E.S.S. or by IACT-based facilities in general (Giacinti et al. 2020). Previously, the limited field of view, pointed observations, and restrictive background estimation approaches of IACTs have prohibited such a measurement. However, this paper demonstrates that these aspects can be partly overcome, although a full characterisation of the γ-ray emission around Geminga remains challenging. Given the systematic uncertainties, no significant asymmetry in the morphology of the γ-ray emission is found. If an asymmetry were to be confirmed at a significant level, it could provide valuable information and constraints on the level of magnetic turbulence in the region (López-Coto & Giacinti 2018).

Given the observed extent of γ-ray emission around Geminga by the survey instruments HAWC and Fermi-LAT of 20–30 pc and ~ 100 pc respectively, we can anticipate that the true extent of TeV γ-ray emission within the H.E.S.S. energy range is ~50pc, or roughly double the current H.E.S.S. field of view, following the estimate shown in Fig. 10. In contrast, the observed extent in other wavebands is just 2′ in X-ray, 10″ in radio or a mere few arcseconds in optical (de Luca et al. 2006; Pellizzoni et al. 2011; Shibanov et al. 2006). It is hence apparent that the emission detected in these wavelength bands are produced by different regions within the PWN-halo system. Nevertheless, given the γ-ray evidence for a halo of escaped electrons far larger than the region indicated for energetic electrons in X-ray and radio observations to date, one may expect that a larger region of much weaker low surface brightness synchrotron emission exists, originating from this electron population. Faint X-ray emission associated with the wider halo of escaped electrons may be identifiable with instruments sensitive to low surface brightness emission, such as eROSITA (Predehl et al. 2021).

With the revised energy loss description and analytical solution from Di Mauro et al. (2019) as described in Sect. 5, a value for D₀ compatible with that of the HAWC analysis (Abeysekara et al. 2017a) is obtained. Three turbulence regimes are probed, Kolmogorov and Kraichnan turbulence, as well as Bohm diffusion, employing the diffusion index δ ∊ [0.3, 0.6, 1.] respectively (Kolmogorov 1991; Kraichnan 1965; Bohm 1949). No preferred value is found with the joint fit of H.E.S.S. and HAWC data, as well as no preference between the ambient magnetic field, injection index, and braking index over the scanned values. We find that the efficiency with which energy is converted to electrons has to be higher than 1% giving a direct lower limit to the proportion of electrons that escape into the pulsar wind nebula over the history of the pulsar. Even if the spectral index is not directly constrained, a spectral index lower than 2 implies the presence of a sub-PeV energy cut-off in order to explain the combined H.E.S.S. and HAWC data.

Considerable systematic uncertainties are also inherent in this analysis, as demonstrated by the different background estimation methods used. It was not possible to measure the true extent of the γ-ray emission around the Geminga pulsar, from the currently available H.E.S.S. data. Caution must therefore be urged in the interpretation of the diffusion model fit in particular: although the fit is performed simultaneously to both the radial profile and spectrum, the radial extent of detected emission is limited by the detector field of view. Given the uncertainty in where the level of background emission is reached (and necessary assumptions made for background normalisation), this does not enable us to make an absolute measurement of the flux or an assumption-free measurement of the diffusion coefficient.

Our study of the X-ray upper limits in a 10′ region around the pulsar allows us to conclude that the magnetic field, in a one diffusion zone scenario, and assuming a constant magnetic field over the X-ray to γ-ray scale, has to be lower than 1 µG in the absence of a sub-PeV energy cut-off. An energy cut-off lower than 75 TeV is needed to account for a magnetic field of 1 µG; the absence of such a cut-off would imply a sub 1 µG magnetic field leading to a tension with the typical value for the normalisation of the magnetic field of the interstellar medium (2.3 µG ≲ B ≲ 6.1 µG Delahaye et al. 2010). This can be seen as another argument favouring the presence of a sub-PeV cut-off in the pulsar injection. The scenarios where the magnetic field is ~1 µG and a strong cut-off is present cannot be disentangled from a sub-µG magnetic field, but an extension towards higher energies of the current γ-ray observations would be able to distinguish these two scenarios.

Several models have considered diffusion as the dominant transport mechanism for trapped particles within PWNe (Tang & Chevalier 2012; Porth et al. 2016). However in older, more evolved systems the halo of energetic particles probes properties of the ambient medium. Initially, it was surprising that diffusion at larger distances away from the pulsar continued to be considerably below the Galactic average value. If this value is representative of properties of the intervening ISM, it would question the paradigm of pulsars as the origin of the cosmic-ray positron excess. Several theories have been proposed to reconcile these aspects, the most popular thus far being that the region of slow diffusion is constrained to the region around the pulsar in which accelerated electrons continue to have significant influence on the ISM (Evoli et al. 2018; Fang et al. 2018, 2019; Profumo et al. 2018). At larger distances, the lower particle energy densities reduce the influence, such that the halo emission gradually decreases to join smoothly with the ISM. Alternatively, it has been noted that although the implied diffusion coefficient is lower than the Galactic average value, it remains consistent with predictions of a spiral arm model with non-uniform diffusion coefficient throughout the Galaxy (Tang & Piran 2019).

In the future, observing these types of large extended sources will become more feasible with the Cherenkov Telescope Array (CTA, the next-generation IACT γ-ray observatory) due to the anticipated larger telescope field of view, of up to ~8–9° towards the highest energies (Acharya et al. 2013). Additionally, divergent pointing strategies are foreseen to cover a larger region of the sky at a single observation instant, using the same number of telescopes but trading an increased field of view for reduced sensitivity (de Naurois 2021; Donini et al. 2019). Whilst this is envisaged in particular for observations of transient phenomena where the location may be poorly constrained, there are obvious advantages of employing this mode to observe extended ‘halo’ phenomena. With these approaches and a much larger number of available telescopes, several of the challenges limiting the capabilities of current IACT arrays such as H.E.S.S. to observe pulsar halos and comparable phenomena at TeV energies will be overcome with CTA. Further survey instruments are also fore-seen, that may increase the number of confirmed pulsar halos dramatically in the near future (López-Coto et al. 2022).

Fig. 12 Distribution of models with p-value >0.003, presented as a violin plot, illustrating the kernel probability density. The width of the shaded area represents the proportion of number of fits located at these positions. The upper and lower bars represent the maximum and minimum values. Top: Fitted energy cut-off distribution as a function of injection index for the combined H.E.S.S. and HAWC dataset. The red horizontal line shows the value of energy cut-off chosen for the fixed Ec fit. Bottom: distribution of the fitted diffusion coefficient in the case of Ec = 1 PeV and Ec left free to vary, for a combined fit to H.E.S.S. and HAWC data. Solid lines indicate the value ranges for the B/C diffusion coefficient and the statistical uncertainty of HAWC diffusion coefficient obtained in Abeysekara et al. (2017a).

Fig. 13 Comparison of the model SED with the XMM-Newton upper limit for different parameters. A flux point for H.E.S.S. corresponding to the innermost radii ≲10′ is shown for comparison.

Fig. 14 Comparison of the fitted SED and profiles with data for the models with p-value >0.003. The best fit is obtained for n = 4.5, η = 0.1, α = 1.8, δ = 1.0, B = 1 µG and is highlighted compared to other models, with fitted parameters of cm2 s−1 and TeV. The uncertainty bands for H.E.S.S. indicates the systematic uncertainty with statistical uncertainties indicated by the error bars, while for HAWC the uncertainty band represents the systematic and statistical uncertainty.

7 Conclusions

With this study, H.E.S.S. confirms the presence of extended very-high-energy γ-ray emission around the Geminga pulsar, on angular scales reaching at least 3.2° radius. Two different methods of background estimation were used to evaluate the systematic uncertainties of the measurement; the field-of-view background method and the On-Off background method, with two independent lists of Off runs for the latter. No evidence for statistically significant asymmetries to the emission or energy-dependent morphology is found. Within a 1º radius region around the pulsar, a spectral analysis obtained a flux normalisation at 1 TeV of (2.8 ± 0.7) × 10⁻¹² cm⁻² s⁻¹ TeV⁻¹ with a best-fit power law spectral index of 2.76 ± 0.22. We fitted an electron diffusion model jointly to the H.E.S.S. data combined with HAWC and compare to XMM-Newton, taking the different spectral extraction regions into account. Thanks to the few-arcmin angular resolution of H.E.S.S., it was possible to extract a flux point for the innermost 10′ radius region around the pulsar, corresponding to the XMM-Newton upper limit. Due to the large number of free variables in the fit, a parameter scan was conducted to cycle over specific values for several quantities, whilst leaving the diffusion coefficient normalisation free. For the best-fit model, a normalisation of the diffusion coefficient of cm² s⁻¹ is preferred at an electron energy of 100 TeV, as well as a cut off energy of the electron spectrum TeV. This value is considerably lower than the Galactic average, yet consistent with results obtained by the HAWC collaboration (Abeysekara et al. 2017a). We find that the mean statistical uncertainty on the diffusion coefficient obtained from our model fit is 17% for the joint fit of H.E.S.S. and HAWC data, whereas it is 27% for the fit to HAWC data only.

This is a challenging analysis: due to the γ-ray emission extending beyond the available sky region, limitations apply such as a lower limit rather than absolute measurement of the size, and a relative rather than absolute flux measurement, likely underestimating the true γ-ray flux. Despite these caveats, we anticipate that this detection and study of extended γ-ray emission around the Geminga pulsar paves the way for further detailed studies of pulsar halos with current and future generation IACT facilities.

Acknowledgements

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.

Appendix A Diffusion model

Appendix A.1 Energy losses

The energy losses are assumed to be dominated by synchrotron losses and Inverse Compton losses. For the latter, we use the model developed in Delahaye et al. (2010, equation 32), which approximates the intermediate regime between Thomson and Klein-Nishina regimes. For the Inverse Compton losses, the three thermal radiation fields considered have the following properties (Moskalenko & Strong 1998): Cosmic Microwave Background (T = 2.7K, U_rad = 0.26 eV cm⁻³), Infra-Red (T = 20K, U_rad = 0.3 eV cm⁻³) and Optical (T = 5000K, U_rad = 0.3 eV cm⁻³).

The total energy losses are computed as: (A.1)

including losses due to synchrotron b_sync and IC scattering on the three aforementioned radiation fields. Due to the different temperatures and energy densities, the transition between the Thomson and Klein-Nishina IC scattering regimes occurs at different energies for the three radiation fields.

With increasing magnetic field strength, the energy range for which Synchrotron losses dominate over inverse Compton losses also increases, as shown in Fig. A.1.

Appendix A.2 Pulsar properties

The pulsar luminosity as a function of time is defined as: (A.2)

where L₀ is the initial luminosity of the pulsar, τ₀ the initial spin-down timescale, and n is the braking index. The pulsar period as a function of time is: (A.3)

and the pulsar age T_* is expressed as: (A.4)

with T_c the characteristic age and P_* the actual pulsar period. It is important to note here that T_* ≃ T_c for n = 3 and P_* ≫ P₀ if the pulsar is old enough to have had a significant increase of its rotation period. Given these previous equations, we can retrieve τ₀ solving the equations (A.3) and (A.4) with P(T_*) = P_*: (A.5)

Fig. A.1 Energy losses from IC and synchrotron processes for different magnetic fields used in the diffusion model. IC is the total inverse Compton losses over all three radiation fields considered.

Appendix B Results for the SED and profile fit

References

All Tables

All Figures

Fig. 1 Surface brightness profile of emission around the Geminga pulsar measured with HAWC (Abeysekara et al. 2017a). Shaded regions indicate the parts of the emission profile that are used as ON (radius θ ≲ 1.0°) and OFF normalisation (radius θ ≳ 3.2°) regions to estimate the background level and evaluate the significance in the analysis of the 2019 H.E.S.S. dataset. The region shown for normalisation of the OFF counts is only accessible with the 2019 dataset, due to the wider pointing strategy used.

In the text

Fig. 2 Excess counts sky maps of the Geminga region above 0.5 TeV using the 2019 dataset with two different background estimation methods. Left: On-Off background estimation method. Right: field-of-view background estimation method. The maps are over-sampled with a 0.08º correlation radius, with contours at the level of 50 and 100 excess counts from maps over-sampled with a 0.5° correlation radius. The 68% containment PSF is shown for comparison and has a value of 0.06º, valid for the innermost θ < 1° region around the pulsar in which the significance is evaluated (indicated by a dashed circle). The Galactic plane is indicated by a green dashed line and the pulsar location with a green triangle, whilst an arrow indicates the proper motion direction of the pulsar (arbitrary length). Background normalisation is performed using regions at θ > 3.2°, indicated by a dotted circle.

In the text

Fig. 3 Observation positions corresponding to the 2006–2008 data set (magenta) and the 2019 dataset (green) at offsets of 0.7° and 1.6° from the pulsar respectively. The white dashed circle indicates the test On region with 1° radius around the pulsar, whilst the white dotted circle indicates the 3.2° radius region beyond which the background is normalised. Green and magenta circles indicate the radius around an observation position with the same offset as the pulsar. Counts map and contours are otherwise as in Fig. 2, left.

In the text

	Fig. 4 Ratio of On counts to background counts using the On-Off background method as a function of radial distance from the centre of the field of view. Shown are two empty sky regions taken on the dwarf spheroidal galaxies Reticulum II and Tucana III, as well as Geminga data.
In the text

	Fig. 5 Radial profile without acceptance correction applied from uncorrected maps, showing counts per unit area with radial distance from the pulsar at energies >0.5 TeV, error bars are statistical. The different background methods show reasonable agreement.
In the text

	Fig. 6 Significance with integration region radius around the pulsar for FoV and On-Off background estimation analyses. Emission within 1 degree radius of the pulsar is detected with a significance of ≳8σ with all analyses.
In the text

	Fig. 7 Spectral energy distribution of the γ-ray emission within 1° radius of the Geminga pulsar. Spectra from two HAWC analyses of Geminga are shown for comparison (Abeysekara et al. 2017a,b), with flux normalisation scaled to match the sub-region from which the H.E.S.S. spectrum is extracted. Error bands include systematic uncertainties for H.E.S.S. but are statistical only for HAWC.
In the text

	Fig. 8 Radial profile constructed using three independent methods; from γ-ray maps with FoV and On-Off background methods (Off list 1), and from a spectral analysis in consecutive ring regions around the pulsar. The latter is limited to radii <2° from the pulsar and is shown without an additional normalisation to background at radii >3.2º applied.
In the text

	Fig. 9 Typical energy loss time versus electron energy. If the particle loss timescale is longer than the pulsar age (here chosen as 340 kyr), then the typical time tE(E), (see Eq. (2)) is taken as the pulsar age. The difference to the HAWC model is due to the approximation of energy losses made in the HAWC model.
In the text

	Fig. 10 Comparison of the predicted typical diffusion scale λ as a function of electron energy for different diffusion models and magnetic fields. The HAWC model is taken from Abeysekara et al. (2017a). It is clear that the emission is predicted to exceed the H.E.S.S. field of view over most of the H.E.S.S. energy range (here converted to electron energy).
In the text

Fig. 11 Distribution of p-value for the scanned parameter combinations presented in Table 6, with different values of η parameter, using the H.E.S.S. and HAWC datasets in combination with Ec fixed to 1 PeV. This distribution is presented as violin plot, illustrating kernel probability density. The width of the shaded area represents the proportion of the number of fits located at these position. The upper and lower bar represent the maximum and minimum values. Models with p-value <0.003 are excluded.

In the text

Fig. 12 Distribution of models with p-value >0.003, presented as a violin plot, illustrating the kernel probability density. The width of the shaded area represents the proportion of number of fits located at these positions. The upper and lower bars represent the maximum and minimum values. Top: Fitted energy cut-off distribution as a function of injection index for the combined H.E.S.S. and HAWC dataset. The red horizontal line shows the value of energy cut-off chosen for the fixed Ec fit. Bottom: distribution of the fitted diffusion coefficient in the case of Ec = 1 PeV and Ec left free to vary, for a combined fit to H.E.S.S. and HAWC data. Solid lines indicate the value ranges for the B/C diffusion coefficient and the statistical uncertainty of HAWC diffusion coefficient obtained in Abeysekara et al. (2017a).

In the text

	Fig. 13 Comparison of the model SED with the XMM-Newton upper limit for different parameters. A flux point for H.E.S.S. corresponding to the innermost radii ≲10′ is shown for comparison.
In the text

Fig. 14 Comparison of the fitted SED and profiles with data for the models with p-value >0.003. The best fit is obtained for n = 4.5, η = 0.1, α = 1.8, δ = 1.0, B = 1 µG and is highlighted compared to other models, with fitted parameters of cm2 s−1 and TeV. The uncertainty bands for H.E.S.S. indicates the systematic uncertainty with statistical uncertainties indicated by the error bars, while for HAWC the uncertainty band represents the systematic and statistical uncertainty.

In the text

	Fig. A.1 Energy losses from IC and synchrotron processes for different magnetic fields used in the diffusion model. IC is the total inverse Compton losses over all three radiation fields considered.
In the text