# The starburst galaxy NGC 253 revisited by H.E.S.S. and Fermi-LAT

H. Abdalla, F. Aharonian, F. Ait Benkhali, E.O. Angüner, M. Arakawa, C.

Arcaro, C. Armand, M. Arrieta, M. Backes, M. Barnard, et al.

► To cite this version:

H. Abdalla, F. Aharonian, F. Ait Benkhali, E.O. Angüner, M. Arakawa, et al. The starburst galaxy NGC 253 revisited by H.E.S.S. and Fermi-LAT. Astronomy & Astrophysics - A&A, 2018, 617, pp.A73. 10.1051/0004-6361/201833202. hal-01833680

HAL Id: hal-01833680

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

Submitted on 22 Nov 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# The starburst galaxy NGC 253 revisited by H.E.S.S. and Fermi-LAT

H.E.S.S. Collaboration, H. Abdalla¹, F. Aharonian³,⁴,⁵, F. Ait Benkhali³, E. O. Angüner¹⁹, M. Arakawa³⁹, C. Arcaro¹, C. Armand²², M. Arrieta¹⁴, M. Backes¹,⁸, M. Barnard¹, Y. Becherini¹⁰, J. Becker Tjus¹¹, D Berge³⁴, S. Bernhard¹², K. Bernlöhr³, R. Blackwell¹³, M. Böttcher¹, C. Boisson¹⁴, J. Bolmont¹⁵, S. Bonnefoy³⁴, P. Bordas³, J. Bregeon¹⁶, F. Brun²⁴, P. Brun¹⁷, M. Bryan⁹, M. Büchele³³, T. Bulik¹⁸, T. Bylund¹⁰, M. Capasso²⁶, S. Caroff²⁷, A. Carosi²², S. Casanova²⁰,³, M. Cerruti¹⁵, N. Chakraborty³, S. Chandra¹, R. C. G. Chaves¹⁶,*, A. Chen²¹, S. Colafrancesco²¹, B. Condon²⁴, I. D. Davids⁸, C. Deil³, J. Devin¹⁶, P. deWilt¹³, L. Dirson², A. Djannati-Ataï²⁸, A. Dmytriiev¹⁴, A. Donath³, L. O’C. Drury⁴, J. Dyks³¹, K. Egberts³², G. Emery¹⁵, J.-P. Ernenwein¹⁹, S. Eschbach³³, S. Fegan²⁷, A. Fiasson²², G. Fontaine²⁷, S. Funk³³, M. Füßling³⁴, S. Gabici²⁸, Y. A. Gallant¹⁶, T. Garrigoux¹, F. Gaté²², G. Giavitto³⁴, D. Glawion²³, J. F. Glicenstein¹⁷, D. Gottschall²⁶, M.-H. Grondin²⁴, J. Hahn³, M. Haupt³⁴, G. Heinzelmann⁵, G. Henri²⁹, G. Hermann³, J. A. Hinton³, W. Hofmann³, C. Hoischen³²,**, T. L. Holch⁷, M. Holler¹², D. Horns², D. Huber¹², H. Iwasaki³⁹, A. Jacholkowska¹⁵,†, M. Jamrozy³⁵, D. Jankowsky³³, F. Jankowsky²³, L. Jouvin²⁸, I. Jung-Richardt³³, M. A. Kastendieck², K. Katarzyński³⁶, M. Katsuragawa⁴⁰, U. Katz³³, D. Kerszberg¹⁵, D. Khangulyan³⁹, B. Khélifi²⁸, J. King³, S. Klepser³⁴, W. Klužniak³¹, Nu. Komin²¹, K. Kosack¹⁷, S. Krakau¹¹, M. Kraus³³, P. P. Krüger¹, G. Lamanna²², J. Lau¹³, J. Lefaucheur¹⁷, A. Lemière²⁸, M. Lemoine-Goumard²⁴, J.-P. Lenain¹⁵, E. Leser³², T. Lohse⁷, M. Lorentz¹⁷, R. López-Coto³, I. Lypova³⁴, D. Malyshev²⁶, V. Marandon³, A. Marcowith¹⁶, C. Mariaud²⁷, G. Martí-Devesa¹², R. Marx³, G. Maurin²², P.J. Meintjes³⁷, A. M. W. Mitchell³, R. Moderski³¹, M. Mohamed²³, L. Mohrmann³³, E. Moulin¹⁷, T. Murach³⁴, S. Nakashima⁴⁰, M. de Naurois²⁷, H. Ndiyavala¹, F. Niederwanger¹², J. Niemiec²⁰, L. Oakes⁷, P. O’Brien³⁰, H. Odaka⁴⁰, S. Ohm³⁴,** M. Ostrowski³⁵, I. Oya³⁴, M. Padovani¹⁶, M. Panter³, R. D. Parsons³, C. Perennes¹⁵, P.-O. Petrucci²⁹, B. Peyaud¹⁷, Q. Piel²², S. Pita²⁸, V. Poireau²², A. Priyana Noel³⁵, D. A. Prokhorov²¹, H. Prokoph³⁴, G. Pühlhofer²⁶, M. Punch²⁸,¹⁰, A. Quirrenbach²³, S. Raab³³, R. Rauth¹², A. Reimer¹², O. Reimer¹², M. Renaud¹⁶, F. Rieger³,³⁸, L. Rinchiuso¹⁷, C. Romoli³, G. Rowell¹³, B. Rudak³¹, E. Ruiz-Velasco³, V. Sahakian⁶,⁵, S. Saito³⁹, D. A. Sanchez²², A. Santangelo²⁶, M. Sasaki³³, R. Schlickeiser¹¹, F. Schüssler¹⁷, A. Schulz³⁴, U. Schwanke⁷, S. Schwemmer²³, M. Seglar-Arroyo¹⁷, M. Senniappan¹⁰, A. S. Seyffert¹, N. Shafi²¹, I. Shilon³³, K. Shiningayamwe⁸, R. Simoni⁹, A. Sinha²⁸, H. Sol¹⁴, F. Spanier¹, A. Specovius³³, M. Spir-Jacob²⁸, Ł. Stawarz³⁵, R. Steenkamp⁸, C. Stegmann³²,³⁴, C. Steppa³², I. Sushk¹, T. Takahashi⁴⁰, J.-P. Tavernet¹⁵, T. Tavernier¹⁷, A. M. Taylor³⁴,**, R. Terrier²⁸, L. Tibaldo³, D. Tiziani³³, M. Tluczykont², C. Trichard²⁷, M. Tsirou¹⁶, N. Tsuji³⁹, R. Tuffs³, Y. Uchiyama³⁹, D. J. van der Walt¹, C. van Eldik³³, C. van Rensburg¹, B. van Soelen³⁷, G. Vasileiadis¹⁶, J. Veh³³, C. Venter¹, A. Viana³,⁴¹, P. Vincent¹⁵, J. Vink⁹, F. Voisin¹³, H. J. Völk³,**, T. Vuillaume²², Z. Wadiasingh¹, S. J. Wagner²³, P. Wagner⁷, R. M. Wagner²⁵, R. White³, A. Wierzcholska²⁰, A. Wörnlein³³, R. Yang³,**, D. Zaborov²⁷, M. Zacharias¹, R. Zanin³, A. A. Zdziarski³¹, A. Zech¹⁴, F. Zefi²⁷, A. Ziegler³³, J. Zorn³, N. Żywucka³⁵

*(Affiliations can be found after the references)*

Received 10 April 2018 / Accepted 24 May 2018

## ABSTRACT

**Context.** NGC 253 is one of only two starburst galaxies found to emit γ-rays from hundreds of MeV to multi-TeV energies. Accurate measurements of the very-high-energy (VHE; $E > 100$ GeV) and high-energy (HE; $E > 60$ MeV) spectra are crucial to study the underlying particle accelerators, probe the dominant emission mechanism(s) and to study cosmic-ray interaction and transport.

**Aims.** The measurement of the VHE γ-ray emission of NGC 253 published in 2012 by H.E.S.S. was limited by large systematic uncertainties. Here, the most up to date measurement of the γ-ray spectrum of NGC 253 is investigated in both HE and VHE γ-rays. Assuming a hadronic origin of the γ-ray emission, the measurement uncertainties are propagated into the interpretation of the accelerated particle population.

**Methods.** The data of H.E.S.S. observations are reanalysed using an updated calibration and analysis chain. The improved Fermi–LAT analysis employs more than 8 yr of data processed using pass 8. The cosmic-ray particle population is evaluated from the combined HE–VHE γ-ray spectrum using NAIMA in the optically thin case.

**Results.** The VHE γ-ray energy spectrum is best fit by a power-law distribution with a flux normalisation of $(1.34 \pm 0.14^{\text{stat}} \pm 0.27^{\text{sys}}) \times 10^{-13} \text{ cm}^{-2} \text{ s}^{-1} \text{ TeV}^{-1}$ at 1 TeV – about 40% above, but compatible with the value obtained in Abramowski et al. (2012). The spectral index $\Gamma = 2.39 \pm 0.14^{\text{stat}} \pm 0.25^{\text{sys}}$ is slightly softer than but consistent with the previous measurement within systematic errors. In the Fermi energy range an integral flux of $F(E > 60 \text{ MeV}) = (1.56 \pm 0.28^{\text{stat}} \pm 0.15^{\text{sys}}) \times 10^{-8} \text{ cm}^{-2} \text{ s}^{-1}$ is obtained. At energies above ~3 GeV the HE spectrum is consistent with a power-law ranging into the VHE part of the spectrum measured by H.E.S.S. with an overall spectral index $\Gamma = 2.22 \pm 0.06^{\text{stat}}$.

**Conclusions.** Two scenarios for the starburst nucleus are tested, in which the gas in the starburst nucleus acts as either a thin or a thick target for hadronic cosmic rays accelerated by the individual sources in the nucleus. In these two models, the level to which NGC 253 acts as a calorimeter is estimated to a range of $f_{\text{cal}} = 0.1$ to 1 while accounting for the measurement uncertainties. The presented spectrum is likely to remain the most accurate measurements until the Cherenkov Telescope Array (CTA) has collected a substantial set of data towards NGC 253.

**Key words.** astroparticle physics – galaxies: starburst – gamma rays: galaxies

## 1. Introduction

Starburst galaxies are characterised by their high star-formation rate (SFR) and gas-consumption times of 1 Gyr or less. The

starburst phase typically lasts for a few hundred million years (see e.g. Kennicutt & Evans 2012; Krumholz 2014, and references therein). Supernova (SN) remnants are believed to be the main sources of the Galactic cosmic rays (CRs). Starburst galaxies with their enhanced SFR and SN rate provide a testbed to probe this paradigm. Furthermore, CRs are star-formation regulators and drive complex chemical reactions by penetrating deep into dense molecular cloud cores (e.g. Indriolo & McCall 2013). There is also increasing evidence that CRs play an

<footnote>* Funded by EU FP7 Marie Curie, grant agreement No. PIEF-GA-2012-332350.</footnote>

<footnote>** Corresponding authors: H.E.S.S. Collaboration, e-mail: contact.hess@hess-experiment.eu</footnote>

<footnote>† Deceased.</footnote>
important role in galaxy formation and evolution (Booth et al. 2013; Salem & Bryan 2014; Salem et al. 2016; Pakmor et al. 2016) by driving galactic winds along expanding magnetic loops (e.g. Breitschwerdt et al. 1991, 1993) that result from their excitation of a Parker instability in the disk (Parker 1966). In this process the CRs heat the outflowing gas through the non-linear Landau damping of the scattering Alfvén waves that are excited by the outward streaming CRs (e.g. Zirakashvili et al. 1996). This CR heating might even prevent low-mass star formation in regions of very high CR densities such as starburst galaxies (e.g. Papadopoulos & Thi 2013).

Observations of starburst galaxies at γ-ray energies provide a useful probe to test the physics of CRs in these objects: (i) They permit inference on the efficiency with which kinetic energy released in SN explosions is channelled via relativistic particles into γ-rays. (ii) γ-rays can be used to infer properties of the interstellar medium in starburst galaxies or probe energy partition between CRs, magnetic fields and radiation fields. (iii) Finally, γ-ray measurements can be used as an independent probe to test the paradigm of CR acceleration in SN remnant shocks.

The two archetypical starburst galaxies NGC 253 and M82 have been discovered to emit γ-rays with energies ranging from hundreds of MeV to several TeV (Acero et al. 2009; VERITAS Collaboration et al. 2009; Fermi-LAT Collaboration et al. 2010). Subsequently, detailed spectral studies of NGC 253 at TeV energies (Abramowski et al. 2012, HESS12 in the following) and the systematic search for GeV γ-ray emission from a sample of star-forming galaxies with Fermi-LAT, including NGC 253 and M82 (Ackermann et al. 2012), have been presented. Recently, NGC 253 has also been studied at hard X-rays with NuSTAR, soft X-rays with Chandra and at radio wavelengths with VLBA (Wik et al. 2014). Broadband spectral energy distribution (SED) modelling is performed with different approaches, ranging from semi-analytical one-zone models as described in e.g. Eichmann & Becker Tjus (2016), to three-dimensional (3D) steady state models (e.g. Persic et al. 2008; Rephaeli et al. 2010) and the treatment of time and space-dependent injectors (e.g. Torres et al. 2012). Starburst galaxies are also discussed as one of the possible source classes contributing to the astrophysical neutrino excess seen by the IceCube collaboration (Aartsen et al. 2014). NGC 253 remains one of the weakest detected TeV γ-ray sources. After three years of improvements to the calibration, reconstruction and analysis, a re-analysis of the γ-ray data, including a re-evaluation of systematic uncertainties, is presented in this work. In addition, 8 yr of Fermi-LAT data are analysed and the updated γ-ray spectrum from 60 MeV to more than 10 TeV is studied. As a result of both improvements, the discussion presented in HESS12 is revisited.

# 2. H.E.S.S. data analysis

## 2.1. H.E.S.S. data

H.E.S.S. is an array of imaging atmospheric Cherenkov telescopes located in the Khomas Highland of Namibia and detects cosmic γ-rays with energies ranging from ~0.1 to ~100 TeV. The data used in this work are identical to the data presented and analysed in HESS12. The observations carried out in 2005 and from 2007 to 2009 sum to 158 h of acceptance-corrected live time. NGC 253 has not been the target of new observations from H.E.S.S. since then. Significantly improved statistics would only be possible at the cost of unreasonably large amounts of observation time. For a detailed description of the instrument and data set, the reader is referred to the

original publication. The differences and improvements of the analysis methods compared to HESS12 are highlighted where applicable.

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

The results presented here and in HESS12 are based on a semi-analytical model of air showers for the event reconstruction and background suppression (de Naurois & Rolland 2009). This model analysis provides an improved angular resolution and a much better sensitivity compared to the classical Hillas parameter-based analysis. It is, however, susceptible to imperfections in the detailed modelling of the instrument response. Since the original publication, a small misalignment of the camera's position with respect to the telescope dish was found. This was not fully taken into account in the HESS12 analysis but has been accounted for in newer versions of the analysis software. The improvements in the pointing model have been verified using optical star positions and by application to known, strong γ-ray sources. The misaligned cameras resulted in a broadening of the point-spread function (PSF) and introduced a shift of the main discrimination variable in the model analysis. This shift led to misclassifications of γ-rays as background and resulted in an underestimation of the γ-ray flux of NGC 253. The same behaviour was uncovered and accounted for during the analysis of N157B in the Large Magellanic Cloud (H.E.S.S. Collaboration 2015, supplement). The cross-check analysis presented in HESS12 is not significantly affected by the imperfect pointing model. The resulting systematic difference in reconstructed γ-ray flux between the model analysis and the cross-check analysis was taken into account in the flux systematic uncertainty in HESS12. As the precision of the γ-ray spectrum presented in HESS12 was limited by the systematic flux uncertainty, and since the modelling of the camera positions were improved since then, a reanalysis of the NGC 253 H.E.S.S. data and revised VHE γ-ray spectrum using the same data set as used in HESS12 is presented.

### 2.2.1. VHE γ-ray spectrum

The data reanalysis was performed using the same analysis framework as in HESS12, namely the model analysis (de Naurois & Rolland 2009). An updated position, extension limit, light curve and γ-ray spectrum are derived. The source is detected with a slightly lower significance of $7.2\sigma$ compared to $8.4\sigma$ in HESS12. The updated source position is $RA = 0h47m32.54s \pm 0m11.2s$, $Dec = -25d17'25.4'' \pm 0'10.3''$ (J2000), which changed only marginally towards an even better agreement with the optical center of NGC 253 at $RA = 0h47m33.1s$, $Dec = -25d17'18''$. With the PSF being understood better, a possible extension of the γ-ray source is constrained to $\le 1.4'$ at the $3\sigma$ level compared to the previous limit of $\le 2.4'$. The new γ-ray spectrum, extracted at the best-fit position, is well described by a single power law, depicted in Fig. 1. The spectral points are given in Table 2. The flux normalisation $F_0(1\text{TeV}) = (1.34 \pm 0.14^{\text{stat}} \pm 0.27^{\text{sys}}) \times 10^{-13} \text{cm}^{-2} \text{s}^{-1} \text{TeV}^{-1}$ is ~40 higher, and the best-fit spectrum is with a spectral index of $\Gamma = 2.39 \pm 0.14^{\text{stat}} \pm 0.25^{\text{sys}}$ somewhat softer but consistent within errors compared to HESS12, where a spectral index of $\Gamma^{2012} = 2.14 \pm 0.18^{\text{stat}} \pm 0.30^{\text{sys}}$ and a normalisation at 1 TeV of $F_0^{2012} = (9.6 \pm 1.5^{\text{stat}} (+5.7, -2.9)^{\text{sys}}) \times 10^{-14} \text{TeV}^{-1} \text{cm}^{-2} \text{s}^{-1}$ were reported. Both spectral parameters are consistent within
<figure><img src="image_2.png" /></figure>

**Fig. 1.** H.E.S.S. and Fermi-LAT pass 8 γ-ray SEDs are shown in red and blue, respectively. All error bars represent $1\sigma$ statistical uncertainties. The upper limits are given at 95% confidence level. The red shaded area represents the $1\sigma$ confidence region of the H.E.S.S. fit with combined statistical and systematic uncertainties. The black box shows the $1\sigma$ confidence region of the H.E.S.S. fit from HESS12. The grey box shows the Fermi-LAT 3FGL best fit. We note that the Fermi-LAT measurement uncertainties are dominated by the low statistics. The systematic error of the Fermi-LAT points range from 5 to 20%. The blue area shows the best fit power-law to the Fermi-LAT pass 8 data.

the previously estimated systematic uncertainties. The relative statistical errors are slightly reduced due to the higher reconstructed γ-ray flux. We note that the systematic uncertainties are now comparable to the statistical uncertainties.

### 2.2.2. Estimation of systematic uncertainties

In HESS12, systematic uncertainties were estimated using a cross-check analysis, which accounted for systematic differences of the calibration, reconstruction and background subtraction. As this cross-check analysis proved to be unaffected by the imperfections in the modelling of the camera positions described above, the updated analysis presented here is compared to the original cross-check analysis. In HESS12, the difference in the flux normalisation between the two analysis chains was found to be 50%, while the re-analysed flux normalisations agree within 2%. The difference of best-fit spectral indices is on the 10% level. As a 2% difference is likely not representative for the real systematic uncertainty caused by different calibration chains, event reconstruction and background subtraction procedures, additional tests for systematic effects within the primary analysis framework have been performed.

A test for systematic shifts in the reconstructed γ-ray flux caused by the run selection has been performed by applying the original data quality criteria (e.g. Aharonian et al. 2006) used in HESS12 in comparison to a revised data quality selection. We found the data selection has an impact on the reconstructed flux at a level of 10% and 3% in flux normalisation and spectral index, respectively, for this faint source.

The applied γ-ray selection cuts might also introduce systematic effects. To estimate the impact of the chosen cuts, the data set was analysed using two different cut configurations: one designed for a low-energy threshold and optimised for spectral studies (Standard), the other optimised for weak sources (like NGC 253) and morphological studies, with a higher energy threshold (Faint). The differences between the analyses with the two cut configurations are 5% in the spectral index and 13% in the normalisation, and represent an estimate of the systematic uncertainty associated to the specific choice of the cut configuration.

The atmosphere is an integral part of an imaging atmospheric Cherenkov telescope and varies over time. The assumed atmospheric density profile influences the amount of light predicted to be seen by each camera. The light yield is uncertain by ~10 (Aharonian et al. 2006). In order to estimate the effect of this uncertainty on the spectral parameters, the fit was repeated using response functions that were shifted by ±10% in energy. The

<table><caption>**Table 1.** Estimated systematic uncertainties of the H.E.S.S. observations towards NGC 253.</caption><thead><tr><td>Origin of uncertainty</td><td>Spectral index</td><td>Normalisation (%)</td></tr></thead><tbody><tr><td>Reconstruction, calibration & analysis</td><td>±0.19</td><td>2</td></tr><tr><td>Run selection</td><td>±0.07</td><td>10</td></tr><tr><td>Selection cuts</td><td>±0.11</td><td>13</td></tr><tr><td>Atmospheric modelling</td><td>±0.09</td><td>10</td></tr><tr><td>Total systematic uncertainty</td><td>±0.25</td><td>19</td></tr></tbody></table>

resulting uncertainties are 10% and ±0.09 for the flux normalisation and the spectral index, respectively.

All uncertainties obtained for the flux normalisation and spectral index are summarised in Table 1. The error bars for the H.E.S.S. flux points shown in Fig. 1 only represent the statistical uncertainties. In this figure, the red shaded area indicates the combined statistical and systematic uncertainties of the best-fit power-law model. The black contour depicts the error region as derived in HESS12.

## 3. Fermi-LAT data analysis

Since the publication of HESS12, the dataset accumulated by Fermi-LAT has increased by a factor of two. In addition, the release of the pass 8 data products (Atwood et al. 2013) allows for an additional gain in sensitivity and performance, especially at the lower end of the Fermi-LAT energy range of around 100 MeV. As this is the energy region where differences between hadronic and leptonic emission processes are best visible, a new analysis of more than 8 yr' Fermi-LAT data was performed.

### 3.1. Fermi-LAT data

Fermi-LAT observations towards NGC 253 were selected in the time period of MET 239557417-MET 507108410 (corresponding to 4 August 2008–26 January 2017), more than 8 yr in total and double the data that was used in HESS12. The standard Fermi Science Tools<fnref n="1" /> were employed to carry out the data analysis. A square region-of-interest (ROI) of 15° × 15° was defined around the optical center of NGC 253. In order to suppress albedo background events from the Earth's limb, events arriving during times in which the ROI was observed under unfavourable zenith angles were not included in the analysis. Specifically,

<footnote n="1">¹ v10r0p5, http://fermi.gsfc.nasa.gov/ssc</footnote>
<table><caption>**Table 2.** The γ-ray spectral data as displayed in Figs. 1 and 2.</caption><thead><tr><th>Min. – Mean – Max.<br>Energy (TeV)</th><th>Energy Flux<br>(erg s<sup>-1</sup> cm<sup>-2</sup>)</th></tr></thead><tbody><tr><td>(6.00 – 8.45 – 12.0) × 10<sup>-5</sup></td><td>&lt;1.52 × 10<sup>-12</sup></td></tr><tr><td>(1.20 – 1.70 – 2.39) × 10<sup>-4</sup></td><td>(1.06 ± 0.62) × 10<sup>-12</sup></td></tr><tr><td>(2.39 – 3.39 – 4.76) × 10<sup>-4</sup></td><td>(1.86 ± 0.33) × 10<sup>-12</sup></td></tr><tr><td>(4.76 – 6.79 – 9.50) × 10<sup>-4</sup></td><td>(1.48 ± 0.21) × 10<sup>-12</sup></td></tr><tr><td>(0.95 – 1.36 – 1.89) × 10<sup>-3</sup></td><td>(1.64 ± 0.20) × 10<sup>-12</sup></td></tr><tr><td>(1.89 – 2.72 – 3.78) × 10<sup>-3</sup></td><td>(1.30 ± 0.21) × 10<sup>-12</sup></td></tr><tr><td>(3.78 – 5.43 – 7.54) × 10<sup>-3</sup></td><td>(1.06 ± 0.25) × 10<sup>-12</sup></td></tr><tr><td>(0.75 – 1.09 – 1.50) × 10<sup>-2</sup></td><td>(3.79 ± 2.42) × 10<sup>-13</sup></td></tr><tr><td>(1.50 – 2.17 – 3.00) × 10<sup>-2</sup></td><td>(1.10 ± 0.51) × 10<sup>-12</sup></td></tr><tr><td>0.03 – 0.08 – 0.23</td><td>(2.65 ± 1.71) × 10<sup>-13</sup></td></tr><tr><td>0.26 – 0.35 – 0.43</td><td>(2.97 ± 1.25) × 10<sup>-13</sup></td></tr><tr><td>0.43 – 0.55 – 0.71</td><td>(2.26 ± 0.77) × 10<sup>-13</sup></td></tr><tr><td>0.71 – 0.88 – 1.17</td><td>(2.05 ± 0.64) × 10<sup>-13</sup></td></tr><tr><td>1.17 – 1.43 – 1.93</td><td>(2.88 ± 0.70) × 10<sup>-13</sup></td></tr><tr><td>1.93 – 2.33 – 3.17</td><td>(2.67 ± 0.78) × 10<sup>-13</sup></td></tr><tr><td>3.17 – 3.93 – 5.75</td><td>&lt;1.34 × 10<sup>-13</sup></td></tr><tr><td>5.75 – 7.07 – 10.4</td><td>&lt;1.99 × 10<sup>-13</sup></td></tr></tbody></table>

**Notes.** The Fermi and H.E.S.S. spectral points are separated by the horizontal line. The energy flux uncertainties correspond to the statistical uncertainties only. The flux limits are calculated for a confidence interval of 95%.

times in which the spacecraft was tilted more than 52° and in which the ROI was only observable at zenith angles >90° were neglected. To avoid biases due to energy-dispersion effects at low energies, a pre-selection of γ rays was performed that restricts the energy range to a minimum of 30 MeV.

The HE γ-ray light curve of NGC 253 was monitored with the high cadence long-term monitoring tool *FlaapLUC* (Lenain 2017) which did not reveal significant variability.

### 3.2. New HE γ-ray spectrum

The spectral analysis was performed based on the P8_R2_v6 version of the post-launch instrument response functions (IRFs). Both front and back converted photons were selected. A binned likelihood analysis using the *gtlike* tool in the energy range from 60 MeV to 500 GeV was performed. All known sources as described in the 3FGL, as well as the diffuse galactic background *iso_P8R2_SOURCE_V6* and the isotropic diffuse background γ-ray emission (*gll_iem_v06*) were included in the fit. The flux normalisations and spectral indices of all contributing sources in the ROI were left free during the fit. NGC 253 is detected with a TS value of 480, corresponding to roughly 22σ. The flux above 60 MeV $F(E > 60$ MeV) is $1.56 \pm 0.28^{\text{stat}} \pm 0.15^{\text{sys}} \times 10^{-8} \text{ cm}^{-2} \text{ s}^{-1}$. Assuming the spectrum follows a power-law, the best-fit spectral index is $2.09 \pm 0.07^{\text{stat}} \pm 0.05^{\text{sys}}$.

To derive the SED, the energy range from 60 MeV to 30 GeV was divided into nine equally spaced $log_{10}(E)$ energy bins, while the higher-energy part from 30 to 300 GeV was combined into a single bin due to the low count statistics at these energies. The fluxes obtained in these bins are shown in Fig. 1 and were obtained with a likelihood fit carried out in each bin accounting for the energy dispersion. The spectral points are given in Table 1. In the first energy bin ranging from 60 to 120 MeV NGC 253 is not detected significantly. Therefore an upper limit at a confidence level of 95% is derived. All higher-energy spectral

points have TS values larger than 4, which corresponds to a significance of more than 2σ. A photon with an energy of 214 GeV was detected within 0.1 degrees from NGC 253, which limits the highest-energy bin to 230 GeV. At energies above ~3 GeV, the Fermi-LAT SED is very well described by a power-law extending into the entire H.E.S.S. energy domain. A power-law fit to all data points at energies above 3 GeV is found to yield a spectral slope of $2.22 \pm 0.06^{\text{stat}}$.

## 4. Results and discussion of the combined HE- and VHE-gamma ray spectrum

### 4.1. Cosmic-ray luminosity and propagation in the starburst

From the combined Fermi-H.E.S.S. γ-ray observations in the energy range from 0.1 GeV to 3 TeV, the inferred integrated γ-ray luminosity is estimated to be $L_γ = 1.19 \pm 0.35^{\text{stat}} \times 10^{40} \text{ erg s}^{-1}$. Adopting a fiducial CR luminosity for the nucleus of the starburst galaxy region of $L_{CR} = 10^{41} \text{ erg s}^{-1}$, it is immediately apparent that of the order of ~10% of such a CR luminosity must be transferred to γ-rays. The motivation for this fiducial CR luminosity comes from the inferred Milky Way's CR luminosity, which is estimated to lie within the range $0.6-3 \times 10^{41} \text{ erg s}^{-1}$ (Drury 2012).

Further consideration of the reference CR luminosity value for the starburst region comes from an estimation of the power fed into the CR population by SNe in this system, $L_{CR} = \Theta E_{SN} v_{SN} \approx 1.6 \times 10^{41} \text{ erg s}^{-1}$, where a fixed fraction $\Theta \approx 0.1$ of the supernova remnant (SNR) kinetic power is fed into CRs, a total SNR kinetic energy $E_{SN} = 10^{51} \text{ erg}$ is assumed to be released in each event, and the SNR rate within the starburst region of NGC 253 is $v_{SN} \approx 0.05 \text{ yr}^{-1}$, motivated from radio, infrared (IR), and optical observations of NGC 253 and taken from Engelbracht et al. (1998), Van Buren & Greenhouse (1994), Ohm & Hinton (2013) for the distance of 3.5 Mpc (Dalcanton et al. 2009, see HESS12). This estimated value is within the estimated range of $(0.6-3) \times 10^{41} \text{ erg s}^{-1}$ for the CR luminosity of the Milky Way under the same assumptions (Drury 2012).

The rather hard overall differential γ-ray spectrum observed for this system up to the highest detected energies is an indication against diffusion-dominated transport of the CRs in the starburst region, a scenario which would be analogous to the conventional diffusion picture for diffuse CRs in the ISM of our Galaxy. Indeed, the high velocity of the starburst wind (see HESS12) rather suggests an advection-dominated transport, and therefore a spectrum of γ-rays emitted from the starburst region whose form is close to that of the source charged-particle spectrum, at least at energies above a few GeV, where the form of the hadronic γ-ray energy spectrum should roughly follow that of the generating charged particles. How closely the form of the γ-ray spectrum follows that of the CRs depends on the density of the gas in the starburst region. We consider two scenarios, which represent two extreme cases for this system.

If the gas in the system acts as a thick target, CRs will lose all their energy in the starburst region through pp interactions. In this regime, the rise of the pp interaction cross-section with energy softens the spectrum of CRs in the system relative to their source spectrum. This softening, however, is naturally compensated for by the γ-ray emission process itself, resulting in the photon index matching that of the source CR spectral slope. To ascertain the best-fit spectral slope for this case, a power-law fit to all data points from 3 GeV to 3 TeV was performed. This approach neglects the systematic uncertainties of both measurements, resulting in small statistical uncertainties on the obtained fit parameters. Energies below 3 GeV are not considered as the proton kinematics
<figure><img src="image_3.png" /><figcaption>Fig. 2. The γ-ray SED obtained with H.E.S.S. and Fermi-LAT are shown together. In addition, the best-fit γ-ray spectrum from pion decay of inelastically scattered protons is shown as a black solid line. The grey shaded areas highlight the 1 and 3σ confidence regions of the fit. The lower panel shows the residuals of the measurement with respect to the best fit, normalised to the measurements statistical uncertainty. Blue and red data points correspond again to the measurements obtained with Fermi-LAT and H.E.S.S., respectively.</figcaption></figure>

start to impact the results, leading to a departure from the power-law description at energies close to the pion production threshold. The best-fit spectral index was found to be 2.22 ± 0.06<sup>stat.</sup>

If the gas in the starburst region is considered to be a thin target for CRs, particles are able to escape the starburst region via advection before losing a significant fraction of their energy. In this regime, the spectral slope of CRs in the system is not altered relative to that in their source. Once again, however, the growth of the *pp* cross-section with energy results in higher energy CRs more efficiently losing their energy than lower energy CRs, hardening the γ-ray spectrum produced. In turn, correcting for this over-representation of high energy γ-rays yields a CR spectrum that is softer than the γ-ray spectrum. In order to estimate the CR spectral shape under these assumptions, a description of starburst nucleus of NGC 253 as well as the inelastic proton scattering cross-section is necessary. The cross-section and branching ratios can be forward folded with a proton test distribution. A tool that allows for all of this is NAIMA (Zabalza 2016), a tool which employs the Markov chain Monte Carlo methods from Goodman & Weare (2010) implemented in emcee (Foreman-Mackey et al. 2013). It fits the parameters of the proton test distributions based on the measured γ-ray spectrum, utilising the Kafexhiu et al. (2014) parameterisation of the *pp* interaction differential cross-sections. Assuming a distance of 3.5 Mpc (Dalcanton et al. 2009) and utilising the cross-section of the *pp* energy losses and pion production from PYTHIA 8, protons were simulated in a kinetic energy range from 0.1 GeV to 0.5 PeV according to a power-law in momentum with spectral index α and normalisation N₀ at the reference momentum p₀ of the form

$$N(E) = \frac{N(p_0)}{\beta c} \times \left(\frac{p}{p_0}\right)^{-\alpha}, \qquad (1)$$

where E is the total energy of the proton, p is the proton momentum, and β the proton velocity in units of c. The fit was done using both, the updated H.E.S.S. and Fermi-LAT spectrum, using statistical uncertainties only. Also, the upper limits

are taken into account in the fit. The resulting best-fit spectral index obtained is α = 2.46 ± 0.03<sup>stat.</sup>. The γ-ray spectrum that is produced via pion decay from this proton distribution is depicted in Fig. 2. The χ²/N<sub>d.o.f.</sub> of the best fit of 0.97 for 12 degrees of freedom, as well as the residuals indicate a good fit to the data. The total energy available in protons above a kinetic energy of 290 MeV (the pion production threshold) is (2.0 ± 0.2<sup>stat.</sup>) × 10<sup>53</sup> erg.

Since the systematic uncertainties on the measured flux normalisation in the H.E.S.S. energy range are as large as the statistical ones, the fit was performed again by shifting the H.E.S.S. data points. For this, all measured flux points were increased (or decreased) by 20%. The resulting best fit parameters for the proton spectrum varied by ~3%<sup>sys</sup><sub>H.E.S.S. norm</sub> and ±0.03<sup>sys</sup><sub>H.E.S.S. norm</sub> for the normalisation and the spectral index, respectively. Additionally, the impact of the exact choice of *pp*-interaction parameterisation was tested by repeating the fit using alternative descriptions of *pp*-interaction processes from QGSJET and SYBILL instead of PYTHIA8. We found that the choice of the parameterisation adds another uncertainty at the same level as the statistical ones, namely ~6%<sup>sys</sup><sub>interaction</sub> and ±0.03<sup>sys</sup><sub>interaction</sub> for the flux normalisation and the spectral index respectively. If the H.E.S.S. measurement is excluded from the fit, the best-fit spectral index is around 2.6 and undershoots the H.E.S.S. measurement.

A limiting factor in the effort to recover the CR source spectrum comes from the complications introduced into the in-situ CR picture, from a consideration of the potential competing particle propagation and energy loss channels. The source spectrum will have an index between those obtained in the *thick* and *thin* one-dimensional (1D) models as long as convection dominates the particle propagation in the starburst region and only pion-decay γ-ray production matters.

## 4.2. Cosmic-ray calorimetry in the starburst nucleus

The level at which the starburst system acts as a CR calorimeter, f<sub>cal</sub>, is defined by the ratio of power that is channelled into pion production relative to the total amount of potential power
available for pion production. A complete comparison of the calorimetric level therefore requires the further estimation of the fraction of CR energy in the population with energies above the threshold for pion production, $f_π$. We note that γ-ray observations are only sensitive to this high-energy component. As was discussed in HESS12, $f_π$ is reasonably estimated with $f_π \approx 3 - \Gamma$ (see Appendix in HESS12) which is based on the assumption of a simple power-law extrapolation of the CR momentum spectrum with index $\Gamma$ over the whole particle energy range.

An estimation of the calorimetric level $f_{cal} = L_π/L_{CR}(E_{CR} > E_π^{th})$ of the system, using the reference values and a value $f_π \sim 0.66$, intermediate between the two extreme cases $\Gamma = 2.22$ and 2.46, gives

$$f_{cal} \approx 0.34 \left( \frac{0.66}{f_π} \right) \left( \frac{L_γ}{1.2 \times 10^{40} \text{ erg s}^{-1}} \right) \left( \frac{1.6 \times 10^{41} \text{ erg s}^{-1}}{L_{CR}} \right), (2)$$

where $E_π^{th}$ is the threshold energy in pp collisions of CRs with the target gas, and the relation $L_π = 3L_γ$ has been assumed. It is important to note, however, that on top of the uncertainties in the determination of $L_γ$, considerable contributions to the uncertainty of this value exist in the determination of the luminosity of CRs with energies above the pion production threshold ($f_π L_{CR}$).

The adoption of the extrema in the estimated range of the values $L_γ$, $L_{CR}$, and $f_π$ in Eq. (2) can be used to estimate the subsequent uncertainty range in the calorimetric value estimation, $f_{cal}$. These uncertainty contributions in its derivation collectively broaden the overall uncertainty, with the corresponding range of values obtained for NGC 253 being $f_{cal} \approx 0.1-1$. As part of this derivation, we impose a ceiling limit of 1, since values beyond this level are considered unphysical. This result highlights that only a crude order of magnitude estimation of the calorimetric value is presently possible. As such, this result is compatible with the estimate of Wang & Fields (2017) found using older HE and VHE γ-ray spectra in a more detailed calculation of the thick target scenario.

Furthermore, it should be noted that this level of uncertainty prevents a true estimation of the underlying uncertainty in $L_{CR}$, which is difficult to estimate in the absence of direct CR (and their secondaries) observations in external galaxy systems. Our estimation of the level of uncertainty in this result should therefore be considered as a lower limit on this range, since these additional contributions would be expected to further broaden it. Such considerations highlight the difficulties faced in constraining the calorimetric value for starburst galaxies, and the essential role played by high-quality spectral measurements.

## 5. Conclusions

The observational and analysis results presented here strengthen the interpretation of the γ-ray emission in a hadronic scenario, as previously considered in Abramowski et al. (2012). Key supporting findings are based on the improved H.E.S.S. and Fermi-LAT analysis. The deeper understanding of the systematics at energies above 100 GeV help to better constrain the spectral shape. The Fermi-LAT pass 8 analysis and a factor two more statistics provide a more accurate measurement of the γ-ray emission in the energy range below 100 GeV.

The assumption that a population of protons is giving rise to the measured γ-ray spectrum through hadronic collisions provides an excellent fit to the data. Based on the presented analyses, the CR luminosity in the starburst nucleus was evaluated in two extreme scenarios assuming the gas to act as either a thin or thick target for the CRs. These two scenarios allowed us to

bracket the CR luminosity in NGC 253 with uncertainties of one order of magnitude. The calorimetric level of NGC 253 was calculated to lie in the transition between the two scenarios with allowed values ranging from $f_{cal} \approx 0.1$ to 1.

The presented spectra will remain the most precise measurements for the coming years. The VHE γ-ray spectrum will only be updated once the Cherenkov Telescope Array (CTA) has started operations and collected a sizeable set of observations towards NGC 253. The CTA measurement will provide more precise and detailed γ-ray data, potentially yielding an estimate of the extension of the region emitting VHE γ-rays. As demonstrated by CTA Consortium (2017; their Fig. 11.4), the current gap in high-quality γ-ray data from 50 to 200 GeV will be filled down to 100 GeV within 100 h of observations. At higher energies, CTA will be able to probe the presence or absence of a cut-off of the γ-ray spectrum and therefore able to probe the acceleration limit of the astrophysical particle accelerators in the starburst nucleus or the onset of γ-γ absorption in the dense radiation fields. Further significant improvements in the HE γ-ray domain will only be possible with missions like e-ASTROGAM, which could provide an accurate measurement in the energy range below 1 GeV (De Angelis et al. 2017) and probe the existence of a “pion-bump” as predicted by hadronic emission models.

**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 Alexander von Humboldt Foundation, the Deutsche Forschungsgemeinschaft, the French Ministry for Research, the CNRS-IN2P3, the U.K. Science and Technology Facilities Council (STFC), the Knut and Alice Wallenberg Foundation the National Science Centre, Poland grant no. 2016/22/M/ST9/00382, 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 Innsbruck University, the Austrian Science Fund (FWF), and the Austrian Federal Ministry for Science, Research and Economy, the University of Adelaide and the Australian Research Council, the Japan Society for the Promotion of Science and by the University of Amsterdam. We appreciate the excellent work of the technical support staff in Berlin, Heidelberg, Palaiseau, Paris, Saclay, 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.

## References

* Aartsen, M. G., Ackermann, M., Adams, J., et al. 2014, Phys. Rev. Lett., 113, 101101

* Abramowski, A., Acero, F., Aharonian, F., et al. 2012, ApJ, 757, 158

* Acero, F., Aharonian, F., Akhperjanian, A., et al. 2009, Science, 326, 1080

* Ackermann, M., Ajello, M., Allafort, A., et al. 2012, ApJ, 755, 164

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

* Atwood, W., Albert, A., Baldini, L., et al. 2013, in 2012 Fermi Symposium proceedings – eConf C121028 [arXiv: 1303.3514]

* Booth, C. M., Agertz, O., Kravtsov, A. V., & Gnedin, N. Y. 2013, ApJ, 777, L16

* Breitschwerdt, D., McKenzie, J., & Völk, H. 1991, A&A, 245, 79

* Breitschwerdt, D., McKenzie, J., & Völk, H. 1993, A&A, 269, 54

* CTA Consortium. 2017, Science with the Cherenkov Telescope Array [arXiv: 1709.07997]

* Dalcanton, J. J., Williams, B. F., Seth, et al. 2009, ApJS, 183, 67

* De Angelis, A., Tatischeff, V., Tavani, M., et al. 2017, Exp. Astron., 44, 25

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

* Drury, L. O. 2012, Astropart. Phys., 39, 52

* Eichmann, B., & Becker Tjus, J. 2016, ApJ, 821, 87

* Engelbracht, C., Rieke, M. J., Rieke, G. H., Kelly, D., & Achtermann, J. 1998, ApJ, 505, 639

* Fermi-LAT Collaboration, Abdo, A. A., Ackermann, M., Ajello, M., & et al. 2010, ApJ, 709, L152
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65

H.E.S.S. Collaboration. 2015, Science, 347, 406

Indriolo, N., & McCall, B. J. 2013, Chem. Soc. Rev., 42, 7763

Kafexhiu, E., Aharonian, F., Taylor, A. M., & Vila, G. S. 2014, Phys. Rev., D90, 123014

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Krumholz, M. R. 2014, Phys. Rep., 539, 49

Lenain, J.-P. 2017, Astrophys. Source Code Lib.

Ohm, S., & Hinton, J. 2013, MNRAS, 429, L70

Pakmor, R., Pfrommer, C., Simpson, C. M., & Springel, V. 2016, ApJ, 824, L30

Papadopoulos, P. P., & Thi, W.-F. 2013, Ap&SS Proc., 34, 41

Parker, E. 1966, ApJ, 145, 811

Persic, M., Rephaeli, Y., & Arieli, Y. 2008, A&A, 486, 143

Rephaeli, Y., Arieli, Y., & Persic, M. 2010, MNRAS, 401, 473

Salem, M., & Bryan, G. L. 2014, MNRAS, 437, 3312

Salem, M., Bryan, G. L., & Corlies, L. 2016, MNRAS, 456, 582

Torres, D. F., Cillis, A., Lacki, B., & Rephaeli, Y. 2012, MNRAS, 423, 822

Van Buren, D., & Greenhouse, M. A. 1994, ApJ, 431, 640

VERITAS Collaboration, Acciari, V. A., Aliu, E., Arlen, T., & et al. 2009, Nature, 462, 770

Wang, X., & Fields, B. D. 2017, MNRAS, 474, 4073

Wik, D. R., Lehmer, B. D., Hornschemeier, A. E., et al. 2014, ApJ, 797, 79

Zabalza, V. 2016, in Proc. 34th International Cosmic Ray Conference (ICRC2015), PoS 236, id.922

Zirakashvili, V., Breitschwerdt, D., Ptuskin, V., & Völk, H. 1996, A&A, 311, 113

15 Sorbonne Université, Université Paris Diderot, Sorbonne Paris Cité, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE, 4 Place Jussieu, 75252 Paris, France

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

17 IRFU, CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette, France

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

19 Aix Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

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

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

22 Laboratoire d'Annecy de Physique des Particules, Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LAPP, 74000 Annecy, France

23 Landessternwarte, Universität Heidelberg, Königstuhl, 69117 Heidelberg, Germany

24 Université Bordeaux, CNRS/IN2P3, Centre d'Études Nucléaires de Bordeaux Gradignan, 33175 Gradignan, France

25 Oskar Klein Centre, Department of Physics, Stockholm University, Albanova University Center, 10691 Stockholm, Sweden

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

27 Laboratoire Leprince-Ringuet, Ecole Polytechnique, CNRS/ IN2P3, 91128 Palaiseau, France

28 APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/Irfu, Observatoire de Paris, Sorbonne Paris Cité, 10, rue Alice Domon et Léonie Duquet, 75205 Paris Cedex 13, France

29 Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

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

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

32 Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Strasse 24/25, 14476 Potsdam, Germany
e-mail: choische@uni-potsdam.de

33 Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Erwin-Rommel-Str. 1, 91058 Erlangen, Germany

34 DESY, 15738 Zeuthen, Germany

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

36 Centre for Astronomy, Faculty of Physics, Astronomy and Informationics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland

37 Department of Physics, University of the Free State, PO Box 339, Bloemfontein 9300, South Africa

38 Heisenberg Fellow (DFG), ITA Universität Heidelberg, Germany

39 Department of Physics, Rikkyo University, 3-34-1 Nishi-Ikebukuro, Toshima-ku, Tokyo 171-8501, Japan

40 Japan Aerospace Exploration Agency (JAXA), Institute of Space and Astronautical Science (ISAS), 3-1-1 Yoshinodai, Chuo-ku, Sagamihara, Kanagawa 229-8510, Japan

41 Now at Instituto de Física de São Carlos, Universidade de São Paulo, Av. Trabalhador São-carlense, 400 – CEP 13566-590 São Carlos, SP, Brazil

<footnote n="1">¹ Centre for Space Research, North-West University, Potchefstroom 2520, South Africa</footnote>

<footnote n="2">² Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, 22761 Hamburg, Germany</footnote>

<footnote n="3">³ Max-Planck-Institut für Kernphysik, PO Box 103980, 69029 Heidelberg, Germany</footnote>

<footnote n="4">⁴ Dublin Institute for Advanced Studies, 31 Fitzwilliam Place Dublin 2, Ireland</footnote>

<footnote n="5">⁵ National Academy of Sciences of the Republic of Armenia, Marshall Baghramian Avenue, 24, 0019 Yerevan, Republic of Armenia</footnote>

<footnote n="6">⁶ Yerevan Physics Institute, 2 Alikhanian Brothers St., 375036 Yerevan, Armenia</footnote>

<footnote n="7">⁷ Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany</footnote>

<footnote n="8">⁸ University of Namibia, Department of Physics, Private Bag 13301 Windhoek, Namibia</footnote>

<footnote n="9">⁹ GRAPPA, Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands</footnote>

<footnote n="10">¹⁰ Department of Physics and Electrical Engineering, Linnaeus University, 351 95 Växjö, Sweden</footnote>

<footnote n="11">¹¹ Institut für Theoretische Physik, Lehrstuhl IV: Weltraum und Astrophysik, Ruhr-Universität Bochum, 44780 Bochum, Germany</footnote>

<footnote n="12">¹² Institut für Astro- und Teilchenphysik, Leopold-Franzens-Universität Innsbruck, 6020 Innsbruck, Austria</footnote>

<footnote n="13">¹³ School of Physical Sciences, University of Adelaide, Adelaide 5005, Australia</footnote>

<footnote n="14">¹⁴ LUTH, Observatoire de Paris, PSL Research University, CNRS, Université Paris Diderot, 5 Place Jules Janssen, 92190 Meudon, France</footnote>