OPEN ACCESS

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

CrossMark

Very-high-energy γ-Ray Emission from Young Massive Star Clusters in the Large Magellanic Cloud

F. Aharonian<sup>1,2,3</sup>, F. Ait Benkhali<sup>4</sup>, J. Aschersleben<sup>5</sup>, H. Ashkar<sup>6</sup>, M. Backes<sup>7,8</sup>, V. Barbosa Martins<sup>9</sup>, R. Batzofin<sup>10</sup>, Y. Becherini<sup>11,12</sup>, D. Berge<sup>9,13</sup>, K. Bernlöhr<sup>2</sup>, M. Böttcher<sup>8</sup>, J. Bolmont<sup>14</sup>, M. de Bony de Lavergne<sup>15</sup>, J. Borowska<sup>13</sup>, R. Brose<sup>1</sup>, A. Brown<sup>16</sup>, F. Brun<sup>15</sup>, B. Bruno<sup>17</sup>, C. Burger-Scheidlin<sup>1</sup>, S. Casanova<sup>18</sup>, J. Celic<sup>17</sup>, M. Cerruti<sup>11</sup>, T. Chand<sup>8</sup>, S. Chandra<sup>8</sup>, A. Chen<sup>19</sup>, J. Chibueze<sup>8</sup>, O. Chibueze<sup>8</sup>, G. Cotter<sup>16</sup>, P. Cristofari<sup>20</sup>, J. Devin<sup>21</sup>, A. Djannati-Atai<sup>11</sup>, J. Djuvsland<sup>2</sup>, A. Dmytriiev<sup>8</sup>, K. Egberts<sup>10</sup>, S. Einecke<sup>22</sup>, K. Feijen<sup>11</sup>, M. Filipovic<sup>23</sup>, G. Fontaine<sup>6</sup>, S. Funk<sup>17</sup>, S. Gabici<sup>11</sup>, Y. A. Gallant<sup>21</sup>, J. F. Glicenstein<sup>15</sup>, J. Glombitza<sup>17</sup>, G. Grolleron<sup>14</sup>, L. Haerer<sup>2</sup>, B. Heb<sup>24</sup>, J. A. Hinton<sup>2</sup>, W. Hofmann<sup>2</sup>, T. L. Holch<sup>9</sup>, D. Horns<sup>25</sup>, Zhiqiu Huang<sup>2</sup>, M. Jamrozy<sup>26</sup>, F. Jankowsky<sup>4</sup>, I. Jung-Richardt<sup>17</sup>, E. Kasai<sup>7</sup>, K. Katarzyński<sup>27</sup>, R. Khatoon<sup>8</sup>, B. Khélifi<sup>11</sup>, W. Kluźniak<sup>28</sup>, Nu. Komin<sup>19</sup>, K. Kosack<sup>15</sup>, D. Kostunin<sup>9</sup>, A. Kundu<sup>8</sup>, R. G. Lang<sup>17</sup>, S. Le Stum<sup>29</sup>, A. Lemière<sup>11</sup>, M. Lemoine-Goumard<sup>30</sup>, J.-P. Lenain<sup>14</sup>, F. Leuschner<sup>24</sup>, J. Mackey<sup>1</sup>, V. Marandon<sup>15</sup>, G. Martí-Devesa<sup>31</sup>, R. Marx<sup>4</sup>, A. Mehta<sup>9</sup>, A. Mitchell<sup>17</sup>, R. Moderski<sup>28</sup>, M. O. Moghadam<sup>10</sup>, L. Mohrmann<sup>2</sup>, A. Montanari<sup>4</sup>, E. Moulin<sup>15</sup>, M. de Naurois<sup>6</sup>, J. Niemiec<sup>18</sup>, S. Ohm<sup>9</sup>, L. Olivera-Nieto<sup>2</sup>, E. de Ona Wilhelmi<sup>9</sup>, M. Ostrowski<sup>26</sup>, S. Panny<sup>31</sup>, U. Pensec<sup>14</sup>, G. Peron<sup>11</sup>, G. Pühlhofer<sup>24</sup>, A. Quirrenbach<sup>4</sup>, S. Ravikularaman<sup>2,11</sup>, M. Regeard<sup>11</sup>, A. Reimer<sup>31</sup>, O. Reimer<sup>31</sup>, H. Ren<sup>2</sup>, M. Renaud<sup>21</sup>, B. Reville<sup>2</sup>, F. Rieger<sup>2</sup>, G. Rowell<sup>22</sup>, B. Rudak<sup>28</sup>, E. Ruiz-Velasco<sup>2</sup>, K. Sabri<sup>21</sup>, V. Sahakian<sup>32</sup>, H. Salzmann<sup>24</sup>, A. Santangelo<sup>24</sup>, M. Sasaki<sup>17</sup>, J. Schäfer<sup>17</sup>, F. Schüssler<sup>15</sup>, H. M. Schutte<sup>8</sup>, H. Sol<sup>20</sup>, S. Spencer<sup>17</sup>, Ł. Stawarz<sup>26</sup>, S. Steinmassl<sup>2</sup>, C. Steppa<sup>10</sup>, K. Streil<sup>17</sup>, I. Sushch<sup>8</sup>, A. M. Taylor<sup>9</sup>, R. Terrier<sup>11</sup>, M. Tsirou<sup>9</sup>, N. Tsuji<sup>33</sup>, C. van Eldik<sup>17</sup>, M. Vecchi<sup>5</sup>, C. Venter<sup>8</sup>, J. Vink<sup>34</sup>, S. J. Wagner<sup>4</sup>, R. White<sup>2</sup>, A. Wierzcholska<sup>18</sup>, M. Zacharias<sup>4,8</sup>, A. A. Zdziarski<sup>28</sup>, A. Zech<sup>20</sup>, and N. Żywucka<sup>8</sup>

(H.E.S.S. Collaboration)

<sup>1</sup> Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland

<sup>2</sup> Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany; contact.hess@hess-experiment.eu

<sup>3</sup> Yerevan State University, 1 Alek Manukyan St, Yerevan 0025, Armenia

<sup>4</sup> Landesfernwarte, Universität Heidelberg, Königstuhl, D 69117 Heidelberg, Germany

<sup>5</sup> Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands

<sup>6</sup> Laboratoire Leprince-Ringuet, École Polytechnique, CNRS, Institut Polytechnique de Paris, F-91128 Palaiseau, France

<sup>7</sup> University of Namibia, Department of Physics, Private Bag 13301, Windhoek 10005, Namibia

<sup>8</sup> Centre for Space Research, North-West University, Potchefstroom 2520, South Africa

<sup>9</sup> Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany

<sup>10</sup> Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Strasse 24/25, D 14476 Potsdam, Germany

<sup>11</sup> Université de Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France

<sup>12</sup> Department of Physics and Electrical Engineering, Linnaeus University, 351 95 Växjö, Sweden

<sup>13</sup> Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, D 12489 Berlin, Germany

<sup>14</sup> Sorbonne Université, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE, 4 place Jussieu, 75005 Paris, France

<sup>15</sup> IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

<sup>16</sup> University of Oxford, Department of Physics, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK

<sup>17</sup> Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Nikolaus-Fiebiger-Str. 2, 91058 Erlangen, Germany

<sup>18</sup> Institut Fizyki Jądrowej PAN, ul. Radzikowskiego 152, 31-342 Kraków, Poland

<sup>19</sup> School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, Johannesburg, 2050, South Africa; contact.hess@hess-experiment.eu

<sup>20</sup> Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université Paris Cité, 5 Pl. Jules Janssen, 92190 Meudon, France

<sup>21</sup> Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, F-34095 Montpellier Cedex 5, France

<sup>22</sup> School of Physical Sciences, University of Adelaide, Adelaide 5005, Australia

<sup>23</sup> School of Science, Western Sydney University, Locked Bag 1797, Penrith South DC, NSW 2751, Australia

<sup>24</sup> Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, D 72076 Tübingen, Germany

<sup>25</sup> Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, D 22761 Hamburg, Germany

<sup>26</sup> Obserwatorium Astronomiczne, Uniwersytet Jagielloński, ul. Orla 171, 30-244 Kraków, Poland

<sup>27</sup> Institute of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland

<sup>28</sup> Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland

<sup>29</sup> Aix Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

<sup>30</sup> Université Bordeaux, CNRS, LP2I Bordeaux, UMR 5797, F-33170 Gradignan, France

<sup>31</sup> Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstraße 25, 6020 Innsbruck, Austria

<sup>32</sup> Yerevan Physics Institute, 2 Alikhanian Brothers St., 0036 Yerevan, Armenia

<sup>33</sup> RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan

<sup>34</sup> GRAPPA, Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

Received 2024 May 13; revised 2024 July 2; accepted 2024 July 2; published 2024 July 19

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

The Tarantula Nebula in the Large Magellanic Cloud is known for its high star formation activity. At its center lies the young massive star cluster R136, providing a significant amount of the energy that makes the nebula shine so brightly at many wavelengths. Recently, young massive star clusters have been suggested to also efficiently produce very high-energy cosmic rays, potentially beyond PeV energies. Here, we report the detection of very-high-energy γ-ray emission from the direction of R136 with the High Energy Stereoscopic System, achieved through a multicomponent, likelihood-based modeling of the data. This supports the hypothesis that R136 is indeed a very powerful cosmic-ray accelerator. Moreover, from the same analysis, we provide an updated measurement of the γ-ray emission from 30 Dor C, the only superbubble detected at TeV energies presently. The γ-ray luminosity above 0.5 TeV of both sources is $(2-3) \times 10^{35} \text{ erg s}^{-1}$. This exceeds by more than a factor of 2 the luminosity of HESS J1646-458, which is associated with the most massive young star cluster in the Milky Way, Westerlund 1. Furthermore, the γ-ray emission from each source is extended with a significance of >3σ and a Gaussian width of about 30 pc. For 30 Dor C, a connection between the γ-ray emission and the nonthermal X-ray emission appears likely. Different interpretations of the γ-ray signal from R136 are discussed.

Unified Astronomy Thesaurus concepts: Young star clusters (1833); Massive stars (732); Large Magellanic Cloud (903); Gamma-ray astronomy (628)

1. Introduction

It has been known for many decades that cosmic rays (CRs) with extremely high energies reach us on Earth (Particle Data Group 2022). In recent years, observations of γ-rays with PeV energies from throughout the Galaxy have confirmed the long-standing hypothesis that CRs with multi-PeV energies are produced within the Milky Way (Tibet ASγ Collaboration 2021; LHAASO Collaboration 2023). Despite decades of searches, their precise origins are, however, still unresolved. While shock fronts of young supernova remnants (SNRs) have long been considered as the main acceleration sites of CR nuclei ("hadronic CRs"; e.g., Ginzburg & Syrovatskii 1964; Berezinskii et al. 1990), the potential of stellar winds to accelerate CRs was also realized early on (e.g., Cesarsky & Montmerle 1983). In the last few years, young massive star clusters (YMCs) have increasingly been discussed as potentially predominant sources of the highest-energy Galactic CRs (e.g., Aharonian et al. 2019; Morlino et al. 2021; Vieu & Reville 2023). If YMCs generate high-energy hadronic CRs, they are expected to also be sources of γ-rays, which are created predominantly in the decay of neutral pions that emerge when the CRs interact with ambient gas. This is referred to as the "hadronic scenario" for the generation of high-energy γ-ray emission. The hypothesis that YMCs are effective CR accelerators can therefore be tested through observations in the very-high-energy (VHE; $E > 0.1$ TeV) γ-ray domain.

In fact, H.E.S.S. Collaboration et al. (2022) have recently been able to associate the VHE γ-ray source HESS J1646-458 with Westerlund 1, the most massive young star cluster in our Galaxy, thus revealing it as a powerful particle accelerator. This does not yet constitute unequivocal evidence for the acceleration of hadronic CRs by the cluster, however, since the nature of the emitting particles remains ambiguous. In fact, for the case of Westerlund 1, Härer et al. (2023) have demonstrated that the morphology is inconsistent with the standard hadronic scenario, and that a model that explains the γ-ray emission as being due to inverse-Compton (IC) scattering of CR electrons (the "leptonic scenario") provides a more natural explanation of the High Energy Stereoscopic System (H.E.S.S.) measurements. Moreover, the exact acceleration site remains unidentified; proposals in the literature include shocks forming at the interaction of winds of massive stars inside the cluster (e.g., Bykov et al. 2013), the termination shock of the collective

cluster wind (Gupta et al. 2020; Morlino et al. 2021), and magnetic turbulences within the entire superbubble (SB) blown by the cluster wind (e.g., Vieu et al. 2022, and references therein). A definitive observational confirmation of any of these predictions is still lacking. Unfortunately, only about a handful of YMCs in the Milky Way have been detected in the VHE domain so far (see, e.g., H.E.S.S. Collaboration et al. 2022), and the association of the γ-ray emission with the star cluster is not always firm. The detection and observation of further YMCs with VHE γ-rays could therefore help to shed more light on their role as CR accelerators.

The Large Magellanic Cloud (LMC), containing many massive star clusters, is a promising target to search for γ-ray emission from YMCs. Indeed, it is host to 30 Dor C, an SB inflated by the LH 90 association of star clusters (Lucke & Hodge 1970; Lortet & Testor 1984) that is visible not only in the radio and optical domains (Mathewson et al. 1985) but also in nonthermal X-rays (Bamba et al. 2004; Smith & Wang 2004; Yamaguchi et al. 2009; Kavanagh et al. 2015, 2019; Lopez et al. 2020), indicating the presence of high-energy electrons. 30 Dor C is the only confirmed SB that has been detected in VHE γ-rays so far<fnref n="35" /> (H.E.S.S. Collaboration 2015). Located nearby, at the heart of the Tarantula Nebula and its central open cluster NGC 2070, lies the "super star cluster" R136, which is exceptionally rich in massive stars (Crowther et al. 2010). With an estimated age between 1 and 2 Myr (e.g., Massey & Hunter 1998; Brands et al. 2022), R136 is also relatively young, implying that only a few supernovae (SNe) are expected to have occurred since its birth (although some older massive stars have also been found in NGC 2070; see Sabbi et al. 2012; Schneider et al. 2018; Bestenlehner et al. 2020).

In this Letter, we report the discovery of VHE γ-ray emission from the direction of R136 with H.E.S.S., achieved through a multicomponent, likelihood-based modeling of the spatial and spectral distribution of γ-ray-like events. Moreover, from the same analysis, we provide updated results on the SB 30 Dor C. Throughout the Letter, we assume a distance to the LMC of 50 kpc (Pietrzyński et al. 2013).

<footnote n="35">³⁵ The γ-ray emission from Westerlund 1 is very likely also associated with the SB of that cluster. In that case, however, a firm detection of the SB at other wavelengths is lacking.</footnote>
## 2. Data Analysis

H.E.S.S. is a γ-ray observatory located in the Khomas highland of Namibia, at 1800 m above sea level. In its initial configuration that began operating in 2003, it consisted of four identical Cherenkov telescopes with 107 m² mirror area each, arranged in a square with 120 m side length (H.E.S.S. Collaboration 2006). In this configuration, H.E.S.S. is sensitive to γ-rays above a few hundred GeV. In 2012, the array was augmented with a fifth, larger telescope with 600 m² mirror area, extending the sensitivity range to energies below 100 GeV.

H.E.S.S. has collected an extensive data set on the Tarantula Nebula region. Observations are taken in individual “runs” of usually 28 minutes duration. Here, we have analyzed a total of 794 runs taken with all four of the initial telescopes, corresponding to an observation time of ≈360 hr. Of these, 301 runs (≈138 hr) have been taken between 2004 December 30 and 2013 February 7, while the remaining 493 runs (≈222 hr) date to between 2017 October 14 and 2022 February 2.<fnref n="36" /> Because a substantial fraction of the observations have been carried out prior to the installation of the fifth telescope, only uniform data from the four initial telescopes are being considered.

The detection of γ-rays proceeds via measuring the Cherenkov radiation that is emitted by secondary particles in the extensive air shower that is launched when the primary γ-ray impinges on the atmosphere of the Earth. We employ the IMPACT algorithm (Parsons & Hinton 2014) to reconstruct the incoming direction and energy of the primary particles (“events”) from the telescope camera images. Instead of the customary two telescopes, we require that every event is detected by at least three telescopes, which enhances the angular resolution of the instrument at the expense of a slight reduction in effective area. For optimal performance, we furthermore reject events with a reconstructed direction that deviates from the pointing direction of the telescopes by more than 1°.⁵. Owing to the relatively strict selection cuts, the achieved energy threshold for the data set is 0.5 TeV. Instrument response functions (IRFs) have been generated from extensive Monte Carlo simulations (Bernlöhr 2008). The suppression of “hadronic” background events, which consist primarily of air showers initiated by CR nuclei, is performed using the method described in Ohm et al. (2009). We perform high-level analysis of the data with GAMMAP Y0.18.2 (Deil et al. 2020; Donath et al. 2023), employing a three-dimensional binned likelihood fit (Mattox et al. 1996) in a 5° × 5° region of interest (ROI), with spatial pixels of size 0.02° × 0.02° (see Table 2 in Appendix A for details). In the analysis, the residual level of background in the final event sample is described using a model constructed from archival H.E.S.S. observations, following the procedure outlined in Mohrmann et al. (2019). We demonstrate in Appendix A that after an adjustment of the background model to each observation run, we obtain a good description of the residual hadronic background for the full data set.

After the adjustment of the background model, we proceed with modeling of the γ-ray emission, following an iterative procedure in which we continue to add sources to the ROI model until no significant γ-ray emission remains. This procedure is described in more detail in Appendix B. As a result, we obtain a best-fit spatial and spectral model for each source included in the model. As spatial models, we use two-dimensional, radially symmetric Gaussians with variable width $\sigma_{\text{Gauss}}$, while the energy spectra are modeled with either a power law or a log-parabola function. Furthermore, we employ the NAIMA package (Zabalza 2015) to fit physical spectral models of primary CR particles to our data. All modeling results have been cross-checked using a second analysis pipeline that is based on an independent calibration, event reconstruction, and event selection (de Naurois & Rolland 2009).

## 3. Results

We show a γ-ray flux map of the Tarantula Nebula region in Figure 1(a). The applied smoothing is representative of the H. E.S.S. angular resolution for the employed analysis configuration. The nebula is outlined by the H$\alpha$ emission contours from the Southern H-Alpha Sky Survey Atlas (SHASSA; Gaustad et al. 2001). As can be seen on the map, by far the brightest γ-ray source in this region is the pulsar wind nebula (PWN) N157B, also known as HESS J0537−691, which is associated with the pulsar PSR J0537−6910 (H.E.S.S. Collaboration 2012, 2015). To test for the presence of additional sources, we performed a spectromorphological modeling of the γ-ray emission, adding source components—beginning with a model for N157B—to the ROI model until no significant residual emission remains. The modeling, described in detail in Appendix B, reveals the presence of two additional sources: HESS J0535−691, previously associated with 30 Dor C (H.E. S.S. Collaboration 2015), and a new source, HESS J0538−691, whose association with R136 we argue for in this Letter. This is illustrated in Figure 1(b), where the emission from N157B as predicted by our best-fit ROI model has been subtracted, thus making evident the remaining emission from the two other sources. Best-fit parameter values for all modeled sources can be found in Appendix B, as can an explanation of the method we use to derive systematic uncertainties on these values. In what follows, we highlight the most relevant results.

### 3.1. Description of Source Models

N157B—We do not investigate N157B in detail in this Letter, but note that our analysis yields a nonzero extension of $\sigma_{\text{Gauss}}^{\text{N 157B}} = (0.82 \pm 0.20_{\text{stat}} \pm 0.18_{\text{sys}})'$ for this source. We caution, however, that the extended source model is preferred over a pointlike model by only ≈1.3$\sigma$—these seemingly contradictory results are due to the presence of the other two nearby sources, whose model components can absorb part of the emission for the case that N157B is less extended than suggested by our best fit. We therefore do not claim an extension and provide an upper limit of $\sigma_{\text{Gauss}}^{\text{N 157B}} < 1.14$ (95% confidence level, statistical uncertainties only). Whether extended or not, as we demonstrate in Appendix B, the results obtained for 30 Dor C and R136 do not depend strongly on the model assumed for N157B. The obtained spectrum for N157B, shown in Figure 7(a) in Appendix B, is compatible with our previously published result.

<footnote n="36">³⁶ The reduced exposure time of the pre-2014 data set with respect to H.E.S.S. Collaboration (2015) is due to stricter selection criteria, in particular due to the requirement that all four 107 m² telescopes participate in each run. This requirement also explains the omission of observations taken between 2013 and 2017, most of which have been taken with an incomplete array, largely due to an upgrade of the cameras of the initial telescopes that took place during this time (Ashton et al. 2020).</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** γ-ray flux maps of the target region of the analysis. The maps show the γ-ray flux $F$ in units of $10^{-8} \text{ cm}^{-2} \text{ s}^{-1} \text{ sr}^{-1}$, integrated above an energy of 0.5 TeV, assuming a power-law spectrum with index -2.5. Smoothing with a top-hat kernel of 0.07 radius has been applied. (a) Entire emission. (b) Residual emission after subtraction of the emission from N157B predicted by the best-fit ROI model. Dashed blue lines show flux contours at (2.5/7.5/12.5) × $10^{-8} \text{ cm}^{-2} \text{ s}^{-1} \text{ sr}^{-1}$ in panel (a) and (1.5/3) × $10^{-8} \text{ cm}^{-2} \text{ s}^{-1} \text{ sr}^{-1}$ in panel (b). Pixels with a negative excess after background subtraction are clipped at zero. The light green contour lines denote Hα emission as inferred by SHASSA (Gaustad et al. 2001).</figcaption></figure>

**30 Dor C**—For the first time, we find HESS J0535−691, associated with 30 Dor C, to be an extended γ-ray source. The best-fit Gaussian width is $\sigma_{\text{Gauss}}^{30 \text{ Dor C}} = (1.91 \pm 0.40_{\text{stat}} \pm 0.20_{\text{sys}})'$, which corresponds to $(27.8 \pm 5.8)_{\text{stat}} \pm 2.9_{\text{sys}})$ pc at the distance to the LMC. This model is preferred over one in which 30 Dor C is described as a pointlike source by $3.3\sigma$. The measured extension is of the same order as the observed size of the X-ray SB, as can be seen from Figure 2. The best-fit position deviates by 1.'1 from that previously obtained in H.E.S.S. Collaboration (2015); this is most likely due to the different analysis method used there (a two-dimensional, i.e., energy-integrated likelihood fit). We note that the new position is in better agreement with the center of the X-ray SB and the compact star clusters located there. The energy spectrum follows a power law with spectral index $\Gamma_{30 \text{ Dor C}} = -2.57 \pm 0.09_{\text{stat}}$.

**R136**—Lastly, HESS J0538−691 is detected as a new γ-ray source with a significance of $6.3\sigma$. The separation between the best-fit position and the location of the YMC R136 is only $\approx 20''$ (see Figure 3). Because there is no other plausible counterpart, we associate HESS J0538−691 with R136.$^{37}$ Similarly to 30 Dor C, we find a preference ($3.1\sigma$) for an extended source, with a Gaussian width of $\sigma_{\text{Gauss}}^{\text{R136}} = (2.30 \pm 0.54_{\text{stat}}^{+0.26}_{-0.21})'$, or $(33.5 \pm 7.9_{\text{stat}}^{+3.8}_{-3.1})_{\text{sys}}$ pc. For the energy spectrum, we find a power-law spectral index $\Gamma_{\text{R136}} = -2.54 \pm 0.15_{\text{stat}}$.

## 3.2. Spectral Results and Energy Requirements

In Figure 4, we show the energy spectra of 30 Dor C and R136 in terms of their γ-ray luminosity. Above the threshold energy of 0.5 TeV, the integrated luminosities are $L_{\gamma}^{30 \text{ Dor C}} \approx 2.9 \times 10^{35} \text{ erg s}^{-1}$ and $L_{\gamma}^{\text{R136}} \approx 2.2 \times 10^{35} \text{ erg s}^{-1}$. Remarkably, these values exceed the luminosity of Westerlund 1 above the same energy, $8 \times 10^{34} \text{ erg s}^{-1}$ (H.E.S.S. Collaboration et al. 2022), by a factor of 2–3.

Figure 4 also shows the γ-ray luminosities predicted by a leptonic and a hadronic model (see Appendix C for details). For the hadronic case, we find spectral indices of the primary proton spectrum of $\Gamma_{\text{p}}^{30 \text{ Dor C}} = -2.64 \pm 0.08_{\text{stat}}$ and $\Gamma_{\text{p}}^{\text{R136}} = -2.59 \pm 0.13_{\text{stat}}$. Assuming an extrapolation of the particle spectrum to 1 GeV with the same spectral index, this would imply a total energy requirement for protons of $W_{\text{p}}^{30 \text{ Dor C}} \approx 2.1 \times 10^{53} (n/1 \text{ cm}^{-3})^{-1} \text{ erg}$ and $W_{\text{p}}^{\text{R136}} \approx 1.1 \times 10^{53} (n/1 \text{ cm}^{-3})^{-1} \text{ erg}$, respectively, where n is the average target gas density. However, at least for 30 Dor C, such an extrapolation would violate the upper limit on the γ-ray flux in the 1–10 GeV range provided by the Fermi-LAT instrument (Fermi-LAT Collaboration 2016). To respect the limit, the primary proton spectrum would, for example, need to transition to a harder spectral index of −2 below ~1 TeV. In this case, one obtains a lower requirement of $W_{\text{p}}^{30 \text{ Dor C}} \approx 1.4 \times 10^{52} (n/1 \text{ cm}^{-3})^{-1} \text{ erg}$ (or, for comparison, although no limit from Fermi-LAT is available, $W_{\text{p}}^{\text{R136}} \approx 9.7 \times 10^{51} (n/1 \text{ cm}^{-3})^{-1} \text{ erg}$). These values can be treated as lower limits for the required energy in protons.

<footnote>$^{37}$ The position of HESS J0538−691 is also compatible with a nearby “clump” of stars (separated by a few parsecs; Sabbi et al. 2012) that do not belong to R136 itself but are part of the encompassing cluster NGC 2070. See also Section 4.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Optical image (credit: ESO; https://www.eso.org/public/images/eso1816a) showing the LH 90 association of star clusters. The yellow circles mark the positions of the compact clusters α-ζ found by Lortet & Testor (1984). The pink contour line outlines the SB 30 Dor C as visible in X-rays (1–2 keV) with XMM-Newton; the red circle denotes MCSNR J0536−6913, a putative SNR (Kavanagh et al. 2015). The white circle marker, solid line, and dotted line indicate the best-fit position, $1\sigma$ Gaussian radius, and 68% containment radius of the best-fit γ-ray model for 30 Dor C, respectively. The transparent band shows the statistical uncertainty on the $1\sigma$ radius, while the arrows denote the systematic uncertainty. The turquoise diamond marks the position reported in H.E.S.S. Collaboration (2015). The green contour lines denote bright spots of $^{12}$CO(1–0) emission measured with the Atacama Large Millimeter/submillimeter Array (Yamane et al. 2021).</figcaption></figure>

For the leptonic model, we include as IC target radiation fields the cosmic microwave background, infrared-to-optical radiation from dust and stars, and ultraviolet radiation specifically from the massive stars in the clusters themselves (see Appendix D for more details). We obtain primary electron spectral indices of $\Gamma_{e}^{30\,\text{Dor C}} = -3.27 \pm 0.11_{\text{stat}}$ and $\Gamma_{e}^{R136} = -3.19 \pm 0.17_{\text{stat}}$. Given the age of the star clusters and the H.E.S.S. energy range, this represents a population of cooled electrons, with spectra steepened compared to injection spectra. Assuming that the primary spectra extend down to at least 0.1 TeV, and taking into account the energy-dependent cooling due to IC scattering, we find a minimum required injection power for electrons of $L_{e}^{30\,\text{Dor C}} \approx 6.4 \times 10^{36} \text{ erg s}^{-1}$ for 30 Dor C and $L_{e}^{R136} \approx 3.9 \times 10^{36} \text{ erg s}^{-1}$ for R136. In the presence of a magnetic field of $B = 5 \, \mu\text{G}$, a value roughly representative for the LMC as a whole (Gaensler et al. 2005), additional synchrotron losses lead to larger requirements of $L_{e}^{30\,\text{Dor C}} \approx 8.5 \times 10^{36} \text{ erg s}^{-1}$ and $L_{e}^{R136} \approx 5.3 \times 10^{36} \text{ erg s}^{-1}$, respectively. For $B = 15 \, \mu\text{G}$, as derived by H.E.S.S. Collaboration (2015) within a leptonic model for 30 Dor C, the corresponding values are $L_{e}^{30\,\text{Dor C}} \approx 1.5 \times 10^{37} \text{ erg s}^{-1}$ and $L_{e}^{R136} \approx 9.5 \times 10^{36} \text{ erg s}^{-1}$, that is, starting to be dominated by synchrotron losses.

## 4. Discussion

For most of our discussion, we will assume a connection between the γ-ray emission and the YMCs LH 90³⁸ and R136.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Composite infrared image from Spitzer (credit: NASA/JPL-Caltech; https://www.spitzer.caltech.edu/image/ssc2020-06b-tarantula-nebula-spitzer-3-color-image) of the region around the star cluster R136. The nominal position of R136 (Høg et al. 2000) is marked with a gray cross. The white circle marker, solid line, and dotted line indicate the best-fit position, $1\sigma$ Gaussian radius, and 68% containment radius of the best-fit γ-ray model for R136, respectively. The transparent band shows the statistical uncertainty on the $1\sigma$ radius, while the arrows denote the systematic uncertainty. The green contour lines denote $^{12}$CO(1–0) emission (Johansson et al. 1998), while the light red contours indicate the positions of dense “knots” of molecular gas identified through $^{12}$CO(2–1) emission (Kalari et al. 2018).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** γ-ray luminosity of 30 Dor C and R136. Shown on the vertical axis is the power output per logarithmic energy interval, that is, $E^2 \cdot (dN/dE) \cdot 4\pi d^2$, with *d* the distance to the source. For comparison, the γ-ray luminosity of the YMC Westlund 1 is also shown (H.E.S.S. Collaboration et al. 2022). The solid and dashed lines display a leptonic and a hadronic model, respectively, which we have fitted to both sources using the NAIMA package (Zabalza 2015).</figcaption></figure>

This does not exclude scenarios in which an SN explodes inside the SB formed by the clusters; the inevitable interaction of the SN shock with the YMC and its environment distinguishes this case from a “normal,” isolated SNR (Badmaev et al. 2024). Alternative origins of the γ-ray emission from the newly detected source HESS J0538−691 will be discussed at the end of this section.

<footnote n="38">³⁸ While LH 90 is, strictly speaking, classified as an OB association (Testor et al. 1993), we will treat it as a YMC in our discussion, as it fulfills, e.g., the definition of Vieu & Reville (2023).</footnote>
<table><caption>**Table 1**

YMC Properties</caption><thead><tr><th>Property</th><th>LH 90</th><th>R136</th><th>References</th></tr></thead><tbody><tr><td>Cluster age [Myr]</td><td>4</td><td>1.5</td><td>1, 2</td></tr><tr><td>Wind power<sup>a</sup> [10<sup>38</sup> erg s<sup>-1</sup>]</td><td>1.5</td><td>10</td><td>3, 4</td></tr><tr><td>Wind velocity [km s<sup>-1</sup>]</td><td>3000</td><td>3000</td><td>5, 6</td></tr><tr><td>Average ISM density [cm<sup>-3</sup>]</td><td>100</td><td>100</td><td>7–11</td></tr><tr><td>Magnetic field [μG]</td><td>15</td><td>15<sup>b</sup></td><td>12</td></tr><tr><td>SB radius [pc]</td><td>74</td><td>56</td><td>...</td></tr><tr><td>Termination shock radius [pc]</td><td>7.9</td><td>8.7</td><td>...</td></tr><tr><td>2D Gaussian width [pc]</td><td>27.8</td><td>33.5</td><td>...</td></tr><tr><td>68% containment radius [pc]</td><td>42.0</td><td>50.5</td><td>...</td></tr><tr><td>Spectral index</td><td>-2.57</td><td>-2.54</td><td>...</td></tr><tr><td>Flux<sup>c</sup> [10<sup>-13</sup> cm<sup>-2</sup> s<sup>-1</sup>]</td><td>4.8</td><td>3.6</td><td>...</td></tr><tr><td>Luminosity<sup>c</sup> [10<sup>35</sup> erg s<sup>-1</sup>]</td><td>2.9</td><td>2.2</td><td>...</td></tr><tr><td>Req. power (pp) [10<sup>36</sup> erg s<sup>-1</sup>]</td><td>1.1</td><td>2.0</td><td>...</td></tr><tr><td>Req. power (IC) [10<sup>37</sup> erg s<sup>-1</sup>]</td><td>1.5</td><td>0.95</td><td>...</td></tr></tbody></table>

**Notes.** Entries in the first section of the table are assumed properties, based on information available in the literature. The expected size of the SB and its termination shock—assuming a spherical expansion—are derived from these properties, following Weaver et al. (1977) and Koo & McKee (1992). Entries in the third section summarize the properties of the γ-ray emission measured with H.E.S.S. (see Section 3; see Table 3 for fit results including uncertainties). The last two rows give the required power in primary CR protons in the hadronic (pp) scenario and the required power in primary CR electrons in the leptonic (IC) scenario, respectively, as determined from the H.E.S.S. measurements. For the hadronic case, we used the primary spectra with a break at 1 TeV (see Section 3), considered no escape or radiation losses, and assumed a continuous injection over the lifetime of the respective cluster.

<sup>a</sup> Averaged over the cluster lifetime.

<sup>b</sup> No literature estimate available.

<sup>c</sup> The integrated flux and luminosity are given above an energy of 0.5 TeV, which is the energy threshold of the H.E.S.S. analysis.

**References.** (1) Testor et al. (1993); (2) Crowther et al. (2016); (3) Kavanagh et al. (2015); (4) Crowther et al. (2010); (5) Kavanagh et al. (2019); (6) Brands et al. (2022); (7) Sano et al. (2017); (8) Yamane et al. (2021); (9) Johansson et al. (1998); (10) Indebetouw et al. (2013); (11) Kalari et al. (2018); (12) H.E. S.S. Collaboration (2015).

We provide in Table 1 an overview of the relevant properties of the YMCs and the surrounding medium, as well as results derived from the H.E.S.S. γ-ray observations. We emphasize that most of the listed properties of the clusters are highly uncertain. The assumed values should thus be viewed as estimates only.

For 30 Dor C, surrounding LH 90, we note in addition that the presence of a shell of nonthermal X-ray emission implies the existence of a fast-moving (≥3000 km s⁻¹) shock, which cannot be the forward shock of the expanding SB, as Hα observations indicate a shell expanding at only ≤100 km s⁻¹. Kavanagh et al. (2019) have concluded from this that the X-ray emission is due to an SNR shock wave that expands fast in the low-density (~10⁻³ cm⁻³) interior of the SB. An alternative explanation could be the termination shock of the collective wind from the LH 90 association of clusters.

Regarding R136, we are not aware of an estimate for the average magnetic field strength in the surroundings of the cluster. For ease of comparison, we assume here the same value adopted for 30 Dor C (15 μG; H.E.S.S. Collaboration 2015). Based on its total mass (~2.2 × 10⁴ M☉; Cignoni et al. 2015) and compactness, R136 is expected to exhibit a collective wind

and inflate an SB (see, e.g., Vieu & Reville 2023), similar to the case of 30 Dor C. However, no SB around R136 could yet be unambiguously identified—although a swept-up shell of gas (Chu & Kennicutt 1994) as well as diffuse thermal X-ray emission (e.g., Townsley et al. 2006) have been detected, which may be attributed to the working of a wind emanating from the cluster. Wang & Helfand (1991), on the other hand, have identified several Hα shells that seem to intersect at the position of R136. The lack of a spherical shell may be attributed to the inhomogeneity of the interstellar medium (ISM) around R136 (see Johansson et al. 1998; Kalari et al. 2018; and Figure 3). For the purpose of our discussion, we nevertheless compute and list in Table 1 the expected size of a putative SB around R136 and its termination shock.

It is interesting that in terms of their γ-ray emission, 30 Dor C and R136 look similar: both appear extended with a Gaussian width of around 30 pc and exhibit a relatively soft spectrum with a spectral index of around -2.6. This is despite the YMCs at their centers being rather unequal: R136 is younger but allegedly exerts a more powerful wind than the LH 90 association. However, it appears that these differences compensate in the sense that the expected size of the SB and position of the termination shock are comparable between the two cases. We thus cannot conclude that the γ-ray emission must originate from different processes or be created at different sites.

In both cases, the γ-ray emission extends well beyond the expected location of the termination shock of the collective cluster wind. This may point to a scenario that is different from the case of Westerlund 1, where the γ-ray emission was found to exhibit a ringlike structure with a radius similar to that of the termination shock (H.E.S.S. Collaboration et al. 2022). We note, however, that compared to Westerlund 1, R136 and 30 Dor C are located in a region with an ISM density that is, on average, 1 order of magnitude larger. Therefore, at least in a hadronic scenario, it is still conceivable that CR nuclei accelerated at the wind termination shock and interacting with gas clouds further away from the cluster are responsible for the γ-ray emission. In that case, on the other hand, one would expect the centroid of the γ-ray emission to coincide with the positions of the densest gas clouds (see Figures 2 and 3), whereas we find it to lie very close to the position of the YMC for both 30 Dor C and R136, somewhat disfavoring a hadronic origin of the emission.

The feasibility of the leptonic and hadronic emission scenario can also be scrutinized by comparing the power in primary CRs required to sustain the γ-ray emission (see Table 1) with the power provided by, for example, the cluster wind, noting that the large uncertainties associated with either prevent a detailed discussion. For 30 Dor C, we obtain in the leptonic scenario a ratio between these two quantities of ~10%. While this would be a surprisingly large efficiency for a leptonic accelerator, we stress that the wind power listed in Table 1 (1.5 × 10³⁸ erg s⁻¹) is an estimate for the *average* power over the cluster lifetime; Kavanagh et al. (2015) have estimated that the *current* power of just the known Wolf-Rayet stars in LH 90 is ~5 × 10³⁸ erg s⁻¹, which alleviates the requirements slightly. Nevertheless, it appears challenging in this scenario to entirely explain the γ-ray emission as resulting from the collective wind of the clusters. As already proposed previously, a recent SN in the LH 90 association may be providing the additional energy that is required to explain the
γ-ray signal (H.E.S.S. Collaboration 2015; Kavanagh et al. 2019). For R136, with its more powerful wind, we find a less demanding but still considerable efficiency of ~1% in the leptonic case.

Considering, on the other hand, a hadronic origin of the γ-ray emission, we deduce a minimum ratio between required power in protons and power provided by the cluster wind of ~0.7% for 30 Dor C and of ~0.2% for R136. We assume in this case that the γ-ray emission originates from interactions of CRs in the dense gas clouds that surround the respective SBs (see Figures 2 and 3); this is in line with the extension of the γ-ray emission approximately matching the expected size of the SB. The derived acceleration efficiencies are considerably lower than those typically derived in the framework of diffusive shock acceleration (DSA), which are O(10%) (e.g., Eichler 1979). In terms of energy requirements, the hadronic scenario thus appears viable, even for somewhat lower gas densities than assumed here. However, a hadronic scenario for 30 Dor C is disfavored, as it requires relatively large magnetic field strengths, which is in disagreement with the magnetic field estimate by Kavanagh et al. (2019). We also note that for both sources, a mixed leptonic–hadronic scenario is possible.

For the leptonic models, the inferred electron distributions ($\Gamma_e \approx -3.2$) are consistent with a synchrotron-cooled injection spectrum of $\Gamma_{e,inj} \approx -2.2$, close to the standard prediction of DSA. In the hadronic scenario, the derived spectral indices for the proton distributions ($\Gamma_p \approx -2.6$), while in line with those inferred from other massive star clusters, are steeper than the DSA prediction. This could be explained by energy-dependent escape from the emitting region, though it has been argued that steep spectra can also result from the acceleration process at wind termination shocks; see, for example, Webb et al. (1985). Future γ-ray observations may disentangle these different possibilities (see, e.g., CTA Consortium 2023).

Regarding the multiwavelength picture, one striking difference between 30 Dor C and R136 is the presence of a nonthermal X-ray shell in the former and the lack thereof in the latter (although a diffuse X-ray source coincident with R136 has been detected with eROSITA; Sasaki et al. 2022). This may be seen as a hint for a different origin of the γ-ray emission: R136 is too young for many stars to have exploded yet, and so its emission may be from stellar winds alone, whereas 30 Dor C may be powered by a combination of winds and a recent SN. There are, however, older massive stars in the encompassing cluster NGC 2070 (Sabbi et al. 2012; Bestenlehner et al. 2020), still within the putative SB of R136, and so SNe should have occurred within this region during the past ~$10^6$ yr, which calls for a different explanation for the dissimilar appearance of 30 Dor C and R136 in the X-ray domain.

Finally, we comment on the possibility of the γ-ray emission from HESS J0538−691 not being connected to R136 (i.e., neither to the collective cluster wind nor to SNe occurring inside the SB around the cluster). As PWNe constitute a large fraction of extended TeV γ-ray sources, it is natural to consider an association of HESS J0538−691 with a PWN. While this generally appears realistic in terms of the source extension and energy spectrum, there is no evidence for the presence of a PWN that could plausibly be associated with R136. It seems unlikely that a pulsar/PWN with a power output of ~$10^{37}$ erg s⁻¹ (see Section 3) leaves no trace at any other wavelength, despite the Tarantula Nebula being one of the most

deeply observed regions in any wave band. Another theoretically viable explanation would be a counterpart to HESS J0538−691 that is located in front of or behind R136 along the line of sight. However, as there is again no evidence of such an object, the association with R136 appears more likely.

## 5. Conclusion

We present new measurements of the VHE γ-ray emission from the Tarantula Nebula region in the LMC with the H.E.S.S. array of Cherenkov telescopes. Utilizing improved analysis techniques, we are able to resolve the emission into three distinct γ-ray sources. The brightest one, HESS J0537−691, is associated with the PWN N157B. The other two sources can be associated with YMCs and/or their surrounding SBs. Both sources are extremely luminous in γ-rays, exceeding even the most massive Galactic young star cluster Westerlund 1 in this regard.

We provide updated results on HESS J0535−691, associated with the SB 30 Dor C around the LH 90 OB association. The extension of the γ-ray emission—measured for the first time—is comparable to the size of the nonthermal X-ray shell around LH 90, suggestive of a common origin. A consideration of the energy requirements suggests that the emission is powered by a recent SN in this case. A lack of correlation between the γ-ray emission and the distribution of molecular gas disfavors a hadronic origin, although we cannot rule out hadronic contributions.

Furthermore, we report the discovery of a new γ-ray source associated with the YMC R136, labeled HESS J0538−691. This source is similar in terms of spatial extension and γ-ray energy spectrum to the source associated with 30 Dor C. However, the lack of an identified SB around R136 complicates the interpretation in this case. Given that R136 is likely to exhibit a strong collective cluster wind, both a leptonic and a hadronic origin of the γ-ray emission appear viable.

The detection of γ-ray emission from the direction of R136 adds to the growing list of YMCs associated with TeV emission. Despite still being small, the population shows quite some variety, in terms of both the γ-ray emission and its interpretation. Our analysis thus provides crucial information to understand the ability of YMCs to accelerate CRs better.

## Acknowledgments

We thank Hidetoshi Sano and Yumiko Yamane for providing gas maps of the 30 Dor C region.

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 Helmholtz Association, 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 Polish Ministry of Education and Science, agreement No. 2021/WK/06, the South African Department of Science and Innovation 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 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.

*Software:* Gammapy (Deil et al. 2020; Donath et al. 2023), Astropy (Robitaille et al. 2013; Price-Whelan et al. 2018, 2022), matplotlib (Hunter 2007), iminuit (Dembinski et al. 2020).

# Appendix A
## Analysis Procedure and Fit of Hadronic Background Model

In this appendix, we provide technical details about the H.E. S.S. data analysis. Furthermore, we present the results of fitting the model for the residual hadronic background to the observation runs.

The likelihood analysis is carried out in a three-dimensional geometry with specifications as listed in Table 2. For every bin $i$ in this “cube,” we calculate a number of expected events $\mu_i$, where we take into account contributions from the hadronic background model as well as from $\gamma$-ray source models, if present. The prediction for the latter is obtained by folding the spatial and spectral source model with the IRFs, that is, with the exposure, energy dispersion matrix, and point-spread function (PSF). The fit then proceeds by comparing $\mu_i$ to the number of actually observed events, $n_i$, simultaneously across all bins. More specifically, the best-fit models are obtained by adjusting the model parameters such that the quantity $-2\log(\mathcal{L})$ is minimized, with the likelihood $\mathcal{L} = \prod_i P(n_i|\mu_i)$ and $P(n_i|\mu_i)$ denoting the Poisson probability to observe $n_i$ events, given an expectation of $\mu_i$. The minimization is done numerically, where we have used the default fitting backend in GAMMAPY, IMINUIT (Dembinski et al. 2020). Two different models with optimized likelihoods $\mathcal{L}_0$ and $\mathcal{L}_1$ can be compared by means of a likelihood ratio test. In the limit of sufficient statistics and in case the parameter values are far from boundaries, the test statistic $\text{TS} = -2\log(\mathcal{L}_0/\mathcal{L}_1)$ follows a $\chi^2$ distribution with $k$ degrees of freedom if the two models are nested (i.e., model 1 can always be reduced to model 0 for a particular choice of parameter values), where $k$ is the difference in the number of model parameters between the models (Wilks 1938).

Although care is taken in the construction of the hadronic background model to predict the background rate in every given observation run as accurately as possible, it is typically necessary to adjust the model slightly to the actually observed

<table><caption>**Table 2**
Summary of Analysis Settings</caption><thead><tr><th>Setting</th><th>Value</th></tr></thead><tbody><tr><td>ROI center (J2000)</td><td>R.A. 5h35m28s25<br>decl. -69°16′13″ 08</td></tr><tr><td>ROI size</td><td>5° × 5°</td></tr><tr><td>Spatial pixel size</td><td>0.02 × 0.02</td></tr><tr><td>Maximum offset angle</td><td>1°5</td></tr><tr><td>Energy binning</td><td>16 bins decade<sup>-1</sup></td></tr><tr><td>Energy range</td><td>0.5–100 TeV</td></tr></tbody></table>

**Note.** The maximum offset angle denotes the maximum allowed angle between the reconstructed direction of every event and the pointing direction of the telescopes.

level of background in each run (see Mohrmann et al. 2019). To avoid a bias due to actual $\gamma$-ray emission in this procedure, we mask regions in the ROI that contain known $\gamma$-ray sources in this step (i.e., the corresponding spatial pixels do not enter the likelihood computation). Besides the sources discussed in detail in this work, this includes the SNR N132D (H.E.S.S. Collaboration 2021) and the $\gamma$-ray binary LMC P3 (H.E.S.S. Collaboration 2018). For every observation run, we fit two parameters of the background model: a global normalization parameter $\phi$ that scales the total background rate and a “spectral tilt” parameter $\delta$ that modifies the predicted rate $R$ at energy $E$ as $R' = R \cdot (E/1 \text{ TeV})^{-\delta}$. We show in Figure 5(a) an energy-integrated residual significance map for the full ROI obtained after the background model adjustment, where we have summed the predicted background and observed events of all observations and use the “Cash” statistic (Cash 1979) to determine the residual significance in each spatial pixel. In the absence of systematic biases, the distribution of significance values in all pixels outside the masked regions should follow a Gaussian distribution centered at zero and with a width of unity, reflecting statistical fluctuations around the expected rate. The inset in Figure 5(a) shows that this is very nearly the case for this analysis, indicating that an excellent description of the hadronic background has been achieved. In Figure 5(b), we provide a zoom-in of the target region of this work, with the location of three detected $\gamma$-ray sources indicated (see Section 3 and Appendix B).

For the subsequent modeling of the $\gamma$-ray emission in the target region, we divide (for technical reasons) the observations into six groups, where the first group contains all observations taken between 2004 and 2013 and the remaining observations are grouped into yearly observation periods. For each group, a “stacked” data set is created by summing the observed number of events, exposure, and predicted background events and averaging the energy dispersion matrix and PSF. The source modeling is then performed as a joint likelihood analysis across these six stacked data sets.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) Significance map for the entire ROI, above an energy of 0.5 TeV, after adjustment of the hadronic background model. The map has been smoothed with a top-hat kernel of 0.07 radius, as indicated in the bottom right corner. The dark gray lines enclose regions that have been excluded in the fit of the background, and the white dashed line indicates the target region of the analysis, shown in more detail on the right. The inset shows the distribution of significance values outside the exclusion regions (purple histogram) as compared to a normal distribution with unity width (green line). (b) Cutout of the significance map, showing the target region. The white rectangle denotes a slice along which γ-ray excess profiles have been computed (see Figure 9).</figcaption></figure>

## Appendix B
Source Fitting

In this appendix, we detail the modeling of the γ-ray sources in the Tarantula Nebula region as displayed in Figure 1. All best-fit parameter values of the final ROI model are summarized in Table 3. In what follows, we first describe the modeling procedure and its main results before we introduce our method of estimating systematic uncertainties for the fit parameters.

### B.1. Modeling Procedure and Results

The iterative modeling procedure is illustrated in Figure 6. For all components, we use as a spatial model a two-dimensional, symmetric Gaussian with width parameter $\sigma_{\text{Gauss}}$. The positions of the models are free to vary, i.e., not fixed to the nominal positions of the respective counterparts.

The first source we include in our model is the PWN N157B (H.E.S.S. Collaboration 2012, 2015). Its energy spectrum, shown in Figure 7(a), is modeled with a log-parabola function,

$$ \frac{dN}{dE} = \phi_0 \left(\frac{E}{E_0}\right)^{-\alpha - \beta \ln(E/E_0)}, \quad (\text{B1}) $$

where the reference energy $E_0 = 1$ TeV remains fixed in the fit. The log-parabola model is preferred over a simple power-law model with a significance of $7.5\sigma$. As shown in Figure 8, the best-fit position of the model is in excellent agreement with the center of the PWN as seen in X-rays.

After the inclusion of N157B, the largest excess emission is visible around 30 Dor C (see Figure 6(b)), which we add as a second source to our model. This leads to an improvement in the fit statistic of TS ≈ 138.9, which—assuming 5 additional degrees of freedom for the model—translates to a detection significance of $11\sigma$. The source appears significantly extended; a pointlike source is excluded with $3.3\sigma$ significance. We

<table><caption>**Table 3**
Best-fit Parameters of the γ-Ray Source Models</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td colspan="3" style="text-align: center;">N157B/HESS J0537-691</td></tr><tr><td>R.A.</td><td>deg</td><td>84.4394 ± 0.0048<sub>stat</sub><br>(5<sup>h</sup>37<sup>m</sup>45.5<sup>s</sup> ± 1.1<sup>s</sup><sub>stat</sub>)</td></tr><tr><td>decl.</td><td>deg</td><td>-69.1713 ± 0.0016<sub>stat</sub><br>(-69°10'17" ± 6<sup>stat</sup>")</td></tr><tr><td>σ<sub>Gauss</sub></td><td>deg</td><td>0.0137 ± 0.0033<sub>stat</sub> ± 0.0030<sub>sys</sub></td></tr><tr><td>φ<sub>0</sub></td><td>10<sup>-13</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup></td><td>8.69 ± 0.56<sub>stat</sub> ± 0.85<sub>sys</sub></td></tr><tr><td>α</td><td>...</td><td>2.03 ± 0.07<sub>stat</sub> ± 0.08<sub>sys</sub></td></tr><tr><td>β</td><td>...</td><td>0.311 ± 0.037<sub>stat</sub></td></tr><tr><td colspan="3" style="text-align: center;">30 Dor C/HESS J0535-691</td></tr><tr><td>R.A.</td><td>deg</td><td>84.021 ± 0.018<sub>stat</sub><br>(5<sup>h</sup>36<sup>m</sup>5.0<sup>s</sup> ± 4.3<sup>s</sup><sub>stat</sub>)</td></tr><tr><td>decl.</td><td>deg</td><td>-69.197 ± 0.006<sub>stat</sub><br>(-69°11'49" ± 22<sup>stat</sup>")</td></tr><tr><td>σ<sub>Gauss</sub></td><td>deg</td><td>0.0319 ± 0.0066<sub>stat</sub> ± 0.0034<sub>sys</sub></td></tr><tr><td>φ<sub>0</sub></td><td>10<sup>-13</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup></td><td>2.54 ± 0.37<sub>stat</sub> ± 0.44<sub>sys</sub></td></tr><tr><td>Γ</td><td>...</td><td>2.57 ± 0.09<sub>stat</sub></td></tr><tr><td colspan="3" style="text-align: center;">R136/HESS J0538-691</td></tr><tr><td>R.A.</td><td>deg</td><td>84.692 ± 0.038<sub>stat</sub> (5<sup>h</sup>38<sup>m</sup>46<sup>s</sup> ± 9<sup>s</sup><sub>stat</sub>)</td></tr><tr><td>decl.</td><td>deg</td><td>-69.103 ± 0.013<sub>stat</sub><br>(-69°06'11" ± 47<sup>stat</sup>")</td></tr><tr><td>σ<sub>Gauss</sub></td><td>deg</td><td>0.0384 ± 0.0090<sub>stat</sub> ± 0.0045<sub>sys</sub></td></tr><tr><td>φ<sub>0</sub></td><td>10<sup>-13</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup></td><td>1.90 ± 0.58<sub>stat</sub> ± 0.45<sub>sys</sub></td></tr><tr><td>Γ</td><td>...</td><td>2.54 ± 0.15<sub>stat</sub></td></tr></tbody></table>

**Note.** Coordinates are given for the epoch J2000. The flux normalization $\phi_0$ is specified at a reference energy of $E_0 = 1$ TeV. None of the sources exhibit significant emission above 30 TeV (see Figures 4 and 7). Statistical errors are at a 68% confidence level. Systematic errors have been derived as explained in Appendix B.2; no systematic errors are quoted in case they are found to be negligible with respect to the statistical error. Systematic errors do not include a systematic uncertainty in the pointing of the telescopes, which can be of the order of 10"–20" (Gillessen 2004).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Maps of residual significance of the analysis target region. (a) After the background model adjustment (see Appendix A). (b) After adding N157B to the ROI model. (c) After adding 30 Dor C to the ROI model. (d) After adding R136 to the ROI model. The blue crosses denote the best-fit position of the respective source models, while the blue circles indicate the best-fit Gaussian width.</figcaption></figure>

model the energy spectrum with a power-law function,

$$ \frac{dN}{dE} = \phi_0 \left( \frac{E}{E_0} \right)^{-\Gamma} \qquad (B2) $$

(with $E_0 = 1$ TeV again fixed), and find it to be in good agreement with our earlier measurement (H.E.S.S. Collaboration 2015), as shown in Figure 7(b). A log-parabola model for 30 Dor C improves the fit only marginally ($TS \approx 4.2$) and is thus not selected as the default model.

Figure 6(c) shows that residual emission is still visible after the addition of 30 Dor C, close to the location of the star cluster R136. Adding a model component for this source improves the fit by $TS \approx 53.6$, implying a detection significance of $6.3\sigma$. The extended source spatial model is preferred with respect to a pointlike model by $3.1\sigma$. As for 30 Dor C, the spectral model is a power-law function, displayed in Figure 7(c).

With N157B, 30 Dor C, and R136 included in the ROI model, the residual significance map (see Figure 6(d)) shows no further significant excess emission. The iterative modeling procedure hence stops at this point. As a demonstration of the good agreement between the observed data and the final model, we show in Figure 9 one-dimensional profiles along the slice displayed in Figure 5(b). For our final model, a comparison between the model prediction and the observed data by means of a $\chi^2$ test yields $\chi^2 = 28.9$ for 24 degrees of freedom.

We have modified our ROI model in various ways to test whether a significantly better-fitting model can be found. For example, as already alluded to in the main text, we have replaced the Gaussian spatial model for N157B with a pointlike one. Surprisingly, given the relatively small statistical error of $0.0033$ (at a 68% confidence level) on the measured extension of $\sigma_{\text{Gauss}}^{\text{N 157B}} = 0.0137$, we find the pointlike model to be disfavored by only $\approx 1.3\sigma$. We
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Spectral energy distribution for N157B (a), 30 Dor C (b), and R136 (c). The blue line and band show the best-fit spectral model with its statistical uncertainty. The flux points have been obtained by refitting the flux normalization $\phi_0$ in the corresponding energy range, keeping the other parameters of the spectral model of the respective source fixed. The published spectra are taken from H.E.S.S. Collaboration (2015).</figcaption></figure>

attribute this to the presence of the two other γ-ray sources, 30 Dor C and R136, which could be more extended than our best-fit model suggests, thus causing the likelihood profile for the extension parameter of N157B to flatten toward smaller values. Despite the marginal preference for the Gaussian spatial model for N157B, we use it in our final model in order to not obtain biased results for the extension of 30 Dor C or R136. For reference, if N157B is modeled as a pointlike source, we obtain $\sigma_{\text{Gauss}}^{30 \text{ Dor C}} = (0.°0334 \pm 0.°0066_{\text{stat}})$ and $\sigma_{\text{Gauss}}^{\text{R136}} = (0.°0431 \pm 0.°0089_{\text{stat}})$, which is compatible with our default result (see Table 3) within the uncertainties. Interpreted as an upper limit, our analysis yields $\sigma_{\text{Gauss}}^{\text{N157B}} < 0.°0190$ at a 95% confidence level (statistical uncertainties only).

We have also tested a model in which the Gaussian spatial model for N157B is allowed to be elongated. While this improves the fit quality considerably when only N157B is included as a γ-ray source (TS ≈ 9.1), this is no longer the case when 30 Dor C and R136 are included as well (TS ≈ 0.5). Conversely, the model with all three sources included is strongly preferred over that with only N157B, even when modeled as elongated (TS ≈ 183.4; note, however, that the two models are not nested in this case). This is still true when adding a component for 30 Dor C (in this case TS ≈ 50),

implying that also just the emission around R136 cannot be explained by allowing the model for N157B to be elongated.

## B.2. Estimation of Systematic Uncertainties

In this section, we explain our procedure for estimating the systematic uncertainties on the γ-ray source model parameters specified in Table 3. We study three different sources of systematic errors: a mismodeling of the instrument PSF, an incorrect estimation of the residual hadronic background, and a wrongly calibrated energy scale. For each effect, we determine the associated systematic uncertainty following a “bracketing” approach. That is, we vary the IRFs according to the assumed magnitude of the effect, repeat the source modeling, and quote the difference to the best-fit parameter values obtained with the nominal IRFs as systematic error. The total uncertainties stated in Table 3 have been computed by summing quadratically the errors obtained for each of the three effects. Where no systematic uncertainty is specified, the error derived with our method is negligible compared to the statistical uncertainty. This does not necessarily mean that there is no systematic uncertainty on the respective parameter in general. For instance, the fitted source positions are subject to possible pointing inaccuracies of the telescopes, which can be of the
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Chandra X-ray image of N157B (Chen et al. 2006), with properties of the best-fit γ-ray source model overlaid. The green circle marker, thin solid line, thick solid line, and thick dotted line indicate the best-fit position, position uncertainty (95% confidence level), $1\sigma$ Gaussian radius, and 68% containment radius, respectively. In addition, the transparent band shows the statistical uncertainty on the $1\sigma$ radius, while the arrows denote the systematic uncertainty.</figcaption></figure>

order of 10''–20'' (Gillessen 2004) and are not covered by the systematic effects we studied in detail here.

### B.2.1. PSF

We validate the accuracy of the PSF for the analysis configuration employed in this study using the bright γ-ray source PKS 2155–304, an active galactic nucleus that appears pointlike for H.E.S.S. (e.g., H.E.S.S. Collaboration 2005, 2010). Using the same configuration, we analyze 576 observation runs (worth ≈274.2 hr of observation time) on PKS 2155–304 taken between 2004 July 14 and 2012 November 10 above an energy threshold of 0.27 TeV. Observations taken during the very strong outburst of PKS 2155–304 in 2006 July (H.E.S.S. Collaboration 2007) have been excluded from the list of analyzed runs, such that the outcome is not dominated by observations taken over a few days only. We apply the same procedure of adjusting the hadronic background model to each observation run, as explained in Appendix A. Subsequently, all observations are combined to obtain a stacked data set.

We first model PKS 2155–304 using a pointlike spatial source model, using the nominal PSF of the data set. As a spectral model, we use a log-parabola function (see Equation (B1)). Inspecting the residual significance map for energies between 0.27 TeV and 1 TeV, shown in Figure 10(a), we find that the model based on the nominal PSF is not able to describe the data perfectly, indicating a systematic mismodeling of the PSF. Relative to the total emission, the remaining residuals are small: in a map computed without including a

source model for PKS 2155–304, the peak significance exceeds 250$\sigma$. The steeply falling γ-ray spectrum of PKS 2155–304 (we obtain $\alpha = 3.63 \pm 0.02$) leads to a much smaller excess of γ-ray events at higher energies, meaning that the systematic effect—if still present there—is no longer visible.

To quantify the effect, we fit a Gaussian to the observed and predicted distribution of angular offsets with respect to the position of PKS 2155–304 (see Figure 10(c)), finding that the predicted distribution, based on the nominal PSF, is about 5% too broad. As a verification, we repeat the analysis with a PSF that has been artificially narrowed (by scaling the radial axis) by 5%. Indeed, the residual significance map obtained from this, displayed in Figure 10(b), is compatible with statistical fluctuations only. We conclude that a mismodeling of the PSF at a level of 5% with respect to its width should be taken into account as a systematic uncertainty.

We therefore repeat the analysis of the Tarantula Nebula region, scaling the PSF for each data set such that it becomes narrower/broader by 5%. Affected most by this effect are the measured extensions ($\sigma_{\text{Gauss}}$) of the modeled γ-ray sources, for which we find PSF-related systematic uncertainties of $\pm 0.0025$ (N157B), $\pm 0.0033$ (30 Dor C), and $(\pm_{-0.0035}^{+0.0044})^\circ$ (R136). A comparison with the total systematic errors (see Table 3) shows that the latter are strongly dominated by the PSF uncertainty.

### B.2.2. Residual Hadronic Background

Refitting the normalization of the residual hadronic background model for each of the stacked data sets shows that variations of the order of 0.5% are still possible. We therefore repeat the source modeling with the background model normalization varied by $\pm 0.5\%$ for all data sets simultaneously. Variations larger than this are possible for each individual observation run. These, however, are expected to average out when combining the observations into stacked data sets. The systematic errors resulting from this estimation are always smaller than the statistical errors. The largest effect is observed for the spectrum normalization parameters ($\phi_0$), where we find uncertainties of 0.6% (N157B), 3.5% (30 Dor C), and 4.2% (R136).

### B.2.3. Energy Scale

Lastly, we study the effect of a wrongly calibrated energy scale. To this end, we evaluate the effective area and energy dispersion matrix at energies that are scaled up or down by 10%. Such a miscalibration of the energy scale can arise, for example, from variations of the aerosol level in the atmosphere (Hahn et al. 2014). Similarly to the variation of the residual hadronic background, the energy scale mostly affects the spectrum normalization parameters. We find uncertainties of 10% (N157B), 17% (30 Dor C), and 18% (R136). The energy-scale-related uncertainties are thus of the same magnitude as the statistical errors.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Count profiles for the slice across the target region indicated in Figure 5(b). The black data points, identical in all panels, show the observed number of events in 0.03 wide bins along the slice. The dark gray histogram displays the predicted residual hadronic background, while the light gray histogram denotes the total model prediction including γ-ray source models. (a) After adding N157B to the ROI model. (b) After adding 30 Dor C to the ROI model. (c) After adding R136 to the ROI model. Indicated in the top right corner of each panel is the result of a χ² test between the total model prediction and the observed data.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Evaluation of systematic uncertainties of the PSF using PKS 2155-304 data. (a) Residual significance map obtained with the nominal PSF, in the energy range 0.27-1 TeV. (b) Same but with the PSF made narrower by 5%. In both maps, the green cross marks the position of PKS 2155-304, which is modeled as a pointlike source. (c) Distribution of $\theta^2$, where $\theta$ is the angular offset between the reconstructed direction of each event and the position of PKS 2155-304. The black histogram and data points show the observed data, while the blue histogram displays the expected distribution based on the nominal PSF. The dashed gray and solid blue lines show fits of a Gaussian to the two distributions, respectively. These fits indicate that the predicted distribution is about 5% broader than the observed one (in terms of $\theta$, not $\theta^2$).</figcaption></figure>

## Appendix C
### Fit of Primary Particle Spectra

To fit primary particle spectra, we employ the NAIMA package (Zabalza 2015). Specifically, we use the NaimaS-p spectralModel wrapper class in GAMMAPY, which allows us to fit the models to our binned data sets directly (as opposed to fitting them to flux points derived from these). For both 30 Dor C and R136, we fit a leptonic model in which the γ-ray emission is due to IC emission from CR electrons and a hadronic model in which the γ-rays are produced in interactions of CR nuclei with gas. As the energy spectrum for the primary particles, we use a power-law model (see Equation (B2)), with $E_0 = 3$ TeV as a fixed parameter. For the

leptonic models, we include target photon fields as described in Appendix D. For the hadronic models, we use in the fit an arbitrarily chosen gas density of $1 \text{ cm}^{-3}$. The spectrum normalization parameter ($\phi_0$) is inversely proportional to the gas density, such that the results can easily be scaled to different densities. The hadronic model is furthermore based on the parameterization of γ-ray production cross sections presented in Kafexhiu et al. (2014), which includes a nuclear enhancement factor that is around 1.7 at TeV energies. In view of the LMC exhibiting a lower metallicity than the Milky Way (e.g., Choudhury et al. 2016), this factor might be slightly too large for the environments studied in this work, which would imply slightly larger spectrum normalization parameters. The
<table><caption>**Table 4**
Best-fit Parameters of the Primary Particle Spectrum Models</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td colspan="3">30 Dor C (Leptonic)</td></tr><tr><td>φ<sub>0</sub></td><td>10<sup>34</sup> eV<sup>-1</sup></td><td>10.1 ± 1.8<sub>stat</sub></td></tr><tr><td>Γ</td><td>...</td><td>3.27 ± 0.11<sub>stat</sub></td></tr><tr><td colspan="3">R136 (Leptonic)</td></tr><tr><td>φ<sub>0</sub></td><td>10<sup>34</sup> eV<sup>-1</sup></td><td>7.0 ± 2.4<sub>stat</sub></td></tr><tr><td>Γ</td><td>...</td><td>3.19 ± 0.17<sub>stat</sub></td></tr><tr><td colspan="3">30 Dor C (Hadronic)</td></tr><tr><td>φ<sub>0</sub></td><td>10<sup>37</sup> eV<sup>-1</sup></td><td>5.7 ± 1.1<sub>stat</sub></td></tr><tr><td>Γ</td><td>...</td><td>2.64 ± 0.08<sub>stat</sub></td></tr><tr><td colspan="3">R136 (Hadronic)</td></tr><tr><td>φ<sub>0</sub></td><td>10<sup>37</sup> eV<sup>-1</sup></td><td>4.0 ± 1.5<sub>stat</sub></td></tr><tr><td>Γ</td><td>...</td><td>2.59 ± 0.13<sub>stat</sub></td></tr></tbody></table>

**Note.** Statistical errors are at a 68% confidence level. The normalization φ₀ is specified at an energy of E₀ = 3 TeV. IC target photon fields assumed in the leptonic models are described in Appendix D. For the fit of the hadronic models, an ambient gas density of 1 cm⁻³ has been assumed.

fit is again carried out in three dimensions, with the same spatial models as used in our previous fit (see Appendix B.1). The parameters of the spatial models are essentially identical to those obtained previously and thus not reported again here. The primary particle spectrum model parameters are summarized in Table 4.

## Appendix D
IC Target Photon Fields

For the spectral modeling of the γ-ray emission from 30 Dor C and R136 with IC models (see Section 3 and Appendix C), we need as input an estimate for the target photon fields. Besides the cosmic microwave background, we include infrared and optical photon fields that we match to measurements in the 30 Doradus region taken from Israel et al. (2010) and Meixner et al. (2013). These can approximately be described as two blackbody radiation fields with temperatures of 50 K and 4000 K and energy densities of 2.25 eV cm⁻³ and 0.88 eV cm⁻³, respectively. We note that the estimate for the far-infrared field is roughly consistent with that of H.E.S.S. Collaboration (2012) for the 30 Doradus region.

In addition, we include for each 30 Dor C and R136 a dedicated ultraviolet radiation field that describes the emission from the hot, massive stars in the respective star clusters. For 30 Dor C, we find from Testor et al. (1993) a mean temperature of 34,000 K and follow Härer et al. (2023) to derive an energy density of 10 eV cm⁻³. In the same manner, using the catalog from Brands et al. (2022), we obtain an effective temperature of 45,000 K and an energy density of 20 eV cm⁻³ for R136. Without the ultraviolet fields included, the normalization parameters φ₀ of the leptonic models stated in Table 4 are reduced by about 5%. Their presence thus does not affect the conclusions drawn in this work.

## ORCID iDs

F. Aharonian https://orcid.org/0000-0003-1157-3915

H. Ashkar https://orcid.org/0000-0002-2153-1818

M. Backes https://orcid.org/0000-0002-9326-6400

V. Barbosa Martins https://orcid.org/0000-0002-5085-8828

R. Batzofin https://orcid.org/0000-0002-5797-3386

Y. Becherini https://orcid.org/0000-0002-2115-2930

D. Berge https://orcid.org/0000-0002-2918-1824

K. Bernlöhr https://orcid.org/0000-0001-8065-3252

M. Böttcher https://orcid.org/0000-0002-8434-5692

M. de Bony de Lavergne https://orcid.org/0000-0002-4650-1666

R. Brose https://orcid.org/0000-0002-8312-6930

F. Brun https://orcid.org/0000-0003-0770-9007

S. Casanova https://orcid.org/0000-0002-6144-9122

M. Cerruti https://orcid.org/0000-0001-7891-699X

A. Chen https://orcid.org/0000-0001-6425-5692

G. Cotter https://orcid.org/0000-0002-9975-1829

A. Djannati-Ataï https://orcid.org/0000-0002-4924-1708

G. Fontaine https://orcid.org/0000-0002-6443-5025

S. Funk https://orcid.org/0000-0002-2012-0080

J. F. Glicenstein https://orcid.org/0000-0003-2581-1742

T. L. Holch https://orcid.org/0000-0001-5161-1168

Zhiqiu Huang https://orcid.org/0000-0002-9239-323X

M. Jamrozy https://orcid.org/0000-0002-0870-7778

B. Khélifi https://orcid.org/0000-0001-6876-5577

Nu. Komin https://orcid.org/0000-0003-3280-0582

D. Kostunin https://orcid.org/0000-0002-0487-0076

A. Kundu https://orcid.org/0000-0003-2128-1414

M. Lemoine-Goumard https://orcid.org/0000-0002-4462-3686

J.-P. Lenain https://orcid.org/0000-0001-7284-9220

F. Leuschner https://orcid.org/0000-0001-9037-0272

J. Mackey https://orcid.org/0000-0002-5449-6131

V. Marandon https://orcid.org/0000-0001-9077-4058

G. Martí-Devesa https://orcid.org/0000-0003-0766-6473

R. Marx https://orcid.org/0000-0002-6557-4924

A. Mitchell https://orcid.org/0000-0003-3631-5648

L. Mohrmann https://orcid.org/0000-0002-9667-8654

A. Montanari https://orcid.org/0000-0002-3620-0173

E. Moulin https://orcid.org/0000-0003-4007-0145

J. Niemiec https://orcid.org/0000-0001-6036-8569

S. Ohm https://orcid.org/0000-0002-3474-2243

L. Olivera-Nieto https://orcid.org/0000-0002-9105-0518

M. Ostrowski https://orcid.org/0000-0002-9199-7031

S. Panny https://orcid.org/0000-0001-5770-3805

G. Peron https://orcid.org/0000-0003-3255-0077

G. Pühlhofer https://orcid.org/0000-0003-4632-4644

A. Reimer https://orcid.org/0000-0001-8604-7077

O. Reimer https://orcid.org/0000-0001-6953-1385

B. Reville https://orcid.org/0000-0002-3778-1432

F. Rieger https://orcid.org/0000-0003-1334-2993

G. Rowell https://orcid.org/0000-0002-9516-1581

B. Rudak https://orcid.org/0000-0003-0452-3805

E. Ruiz-Velasco https://orcid.org/0000-0001-6939-7825

V. Sahakian https://orcid.org/0000-0003-1198-0043

A. Santangelo https://orcid.org/0000-0003-4187-9560

M. Sasaki https://orcid.org/0000-0001-5302-1866

F. Schüssler https://orcid.org/0000-0003-1500-6571

H. M. Schutte https://orcid.org/0000-0002-1769-5617

S. Spencer https://orcid.org/0000-0001-5516-1205
S. Steinmassl https://orcid.org/0000-0002-2865-8563

I. Sushch https://orcid.org/0000-0002-2814-1257

A. M. Taylor https://orcid.org/0000-0001-9473-4758

R. Terrier https://orcid.org/0000-0002-8219-4667

N. Tsuji https://orcid.org/0000-0001-7209-9204

C. van Eldik https://orcid.org/0000-0001-9669-645X

S. J. Wagner https://orcid.org/0000-0002-7474-6062

A. Wierzcholska https://orcid.org/0000-0003-4472-7204

M. Zacharias https://orcid.org/0000-0001-5801-3945

A. A. Zdziarski https://orcid.org/0000-0002-0333-2452

## References

Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, NatAs, 3, 561

Ashton, T., Backes, M., Balzer, A., et al. 2020, APh, 118, 102425

Badmaev, D. V., Bykov, A. M., & Kalyashova, M. E. 2024, MNRAS, 527, 3749

Bamba, A., Ueno, M., Nakajima, H., & Koyama, K. 2004, ApJ, 602, 257

Berezinskii, V. S., Bulanov, S. V., Dogiel, V. A., & Ptuskin, V. S. 1990, in Astrophysics of Cosmic Rays, ed. V. L. Ginzburg (Amsterdam: North-Holland)

Bernlöhr, K. 2008, APh, 30, 149

Bestenlehner, J. M., Crowther, P. A., Caballero-Nieves, S. M., et al. 2020, MNRAS, 499, 1918

Brands, S. A., de Koter, A., Bestenlehner, J. M., et al. 2022, A&A, 663, A36

Bykov, A. M., Gladilin, P. E., & Osipov, S. M. 2013, MNRAS, 429, 2755

Cash, W. 1979, ApJ, 228, 939

Cesarsky, C. J., & Montmerle, T. 1983, SSRv, 36, 173

Chen, Y., Wang, Q. D., Gotthelf, E. V., et al. 2006, ApJ, 651, 237

Choudhury, S., Subramaniam, A., & Cole, A. A. 2016, MNRAS, 455, 1855

Chu, Y.-H., & Kennicutt, R. C. 1994, ApJ, 425, 720

Cignoni, M., Sabbi, E., van der Marel, R. P., et al. 2015, ApJ, 811, 76

Crowther, P. A., Caballero-Nieves, S. M., Bostroem, K. A., et al. 2016, MNRAS, 458, 624

Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, MNRAS, 408, 731

CTA Consortium 2023, MNRAS, 523, 5353

de Naurois, M., & Rolland, L. 2009, APh, 32, 231

Deil, C., Donath, A., Terrier, R., et al. 2020, Gammapy, v0.18.2, Zenodo, doi:10.5281/zenodo.4701500

Dembinski, H., Ongmongkolkul, P., Deil, C., et al. 2020, iminuit v1.3.8, Zenodo, doi: 10.5281/zenodo.3949207

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

Eichler, D. 1979, ApJ, 229, 419

Fermi-LAT Collaboration 2016, A&A, 586, A71

Gaensler, B. M., Haverkorn, M., Staveley-Smith, L., et al. 2005, Sci, 307, 1610

Gaustad, J. E., McCullough, P. R., Rosing, W., & Van Buren, D. 2001, PASP, 113, 1326

Gillessen, S. 2004, PhD thesis, Ruprecht-Karls-Universität Heidelberg, https://pure.mpg.de/pubman/faces/View woltemOverviewPage.jsp?paperId=item_919099

Ginzburg, V. L., & Syrovatskii, S. I. 1964, The Origin of Cosmic Rays (London: Macmillan)

Gupta, S., Nath, B. B., Sharma, P., & Eichler, D. 2020, MNRAS, 493, 3159

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

Härer, L. K., Reville, B., Hinton, J., Mohrmann, L., & Vieu, T. 2023, A&A, 671, A4

H.E.S.S. Collaboration 2005, A&A, 430, 865

H.E.S.S. Collaboration 2006, A&A, 457, 899

H.E.S.S. Collaboration 2007, ApJL, 664, L71

H.E.S.S. Collaboration 2010, A&A, 520, A83

H.E.S.S. Collaboration 2012, A&A, 545, L2

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

H.E.S.S. Collaboration 2018, A&A, 610, L17

H.E.S.S. Collaboration 2021, A&A, 655, A7

H.E.S.S. Collaboration, Blackwell, R., Braiding, C., et al. 2022, A&A, 666, A124

Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27

Hunter, J. D. 2007, CSE, 9, 90

Indebetouw, R., Brogan, C., Chen, C. H. R., et al. 2013, ApJ, 774, 73

Israel, F. P., Wall, W. F., Raban, D., et al. 2010, A&A, 519, A67

Johansson, L. E. B., Greve, A., Booth, R. S., et al. 1998, A&A, 331, 857

Kafexhiu, E., Aharonian, F., Taylor, A. M., & Vila, G. S. 2014, PhRvD, 90, 123014

Kalari, V. M., Rubio, M., Elmegreen, B. G., et al. 2018, ApJ, 852, 71

Kavanagh, P. J., Sasaki, M., Bozzetto, L. M., et al. 2015, A&A, 573, A73

Kavanagh, P. J., Vink, J., Sasaki, M., et al. 2019, A&A, 621, A138

Koo, B.-C., & McKee, C. F. 1992, ApJ, 388, 93

LHAASO Collaboration 2023, PhRvL, 131, 151001

Lopez, L. A., Grefenstette, B. W., Auchettl, K., Madsen, K. K., & Castro, D. 2020, ApJ, 893, 144

Lortet, M. C., & Testor, G. 1984, A&A, 139, 330

Lucke, P. B., & Hodge, P. W. 1970, AJ, 75, 171

Massey, P., & Hunter, D. A. 1998, ApJ, 493, 180

Mathewson, D. S., Ford, V. L., Tuohy, I., et al. 1985, ApJS, 58, 197

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

Meixner, M., Panuzzo, P., Roman-Duval, J., et al. 2013, AJ, 146, 62

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

Morlino, G., Blasi, P., Peretti, E., & Cristofari, P. 2021, MNRAS, 504, 6096

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

Parsons, R. D., & Hinton, J. A. 2014, APh, 56, 26

Particle Data Group 2022, PTEP, 2022, 083C01

Pietrzyński, G., Graczyk, D., Gieren, W., et al. 2013, Natur, 495, 76

Price-Whelan, A. M., Lim, P. L., Earl, N., et al. 2022, ApJ, 935, 167

Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123

Robitaille, T. P., Tollerud, E. J., Greenfield, P., et al. 2013, A&A, 558, A33

Sabbi, E., Lennon, D. J., Gieles, M., et al. 2012, ApJL, 754, L37

Sano, H., Yamane, Y., Voisin, F., et al. 2017, ApJ, 843, 61

Sasaki, M., Knies, J., Haberl, F., et al. 2022, A&A, 661, A37

Schneider, F. R. N., Ramírez-Agudelo, O. H., Tramper, F., et al. 2018, A&A, 618, A73

Smith, D. A., & Wang, Q. D. 2004, ApJ, 611, 881

Testor, G., Schild, H., & Lortet, M. C. 1993, A&A, 280, 426

Tibet ASγ Collaboration 2021, PhRvL, 126, 141101

Townsley, L. K., Broos, P. S., Feigelson, E. D., et al. 2006, AJ, 131, 2140

Vieu, T., Gabici, S., Tatischeff, V., & Ravikularaman, S. 2022, MNRAS, 512, 1275

Vieu, T., & Reville, B. 2023, MNRAS, 519, 136

Wang, Q., & Helfand, D. J. 1991, ApJ, 370, 541

Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377

Webb, G. M., Axford, W. I., & Forman, M. A. 1985, ApJ, 298, 684

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

Yamaguchi, H., Bamba, A., & Koyama, K. 2009, PASJ, 61, S175

Yamane, Y., Sano, H., Filipović, M. D., et al. 2021, ApJ, 918, 36

Zabalza, V. 2015, ICRC (The Hague), 236, 922