# Very-high-energy $\gamma$-ray emission from young massive star clusters in
the Large Magellanic Cloud

F. Aharonian  Dublin Institute for Advanced Studies, 31 Fitzwilliam Place,
Dublin 2, Ireland  Max-Planck-Institut für Kernphysik, Saupfercheckweg 1,
69117 Heidelberg, Germany  Yerevan State University, 1 Alek Manukyan St,
Yerevan 0025, Armenia  F. Ait Benkhali  Landessternwarte, Universität
Heidelberg, Königstuhl, D 69117 Heidelberg, Germany  J. Aschersleben  Kapteyn
Astronomical Institute, University of Groningen, Landleven 12, 9747 AD
Groningen, The Netherlands  H. Ashkar Laboratoire Leprince-Ringuet, École
Polytechnique, CNRS, Institut Polytechnique de Paris, F-91128 Palaiseau,
France  M. Backes University of Namibia, Department of Physics, Private Bag
13301, Windhoek 10005, Namibia  Centre for Space Research, North-West
University, Potchefstroom 2520, South Africa  V. Barbosa Martins Deutsches
Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany  R.
Batzofin Institut für Physik und Astronomie, Universität Potsdam, Karl-
Liebknecht-Strasse 24/25, D 14476 Potsdam, Germany  Y. Becherini Université de
Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France  Department
of Physics and Electrical Engineering, Linnaeus University, 351 95 Växjö,
Sweden  D. Berge Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738
Zeuthen, Germany  Institut für Physik, Humboldt-Universität zu Berlin,
Newtonstr. 15, D 12489 Berlin, Germany  K. Bernlöhr Max-Planck-Institut für
Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany  M. Böttcher Centre
for Space Research, North-West University, Potchefstroom 2520, South Africa
J. Bolmont  Sorbonne Université, CNRS/IN2P3, Laboratoire de Physique Nucléaire
et de Hautes Energies, LPNHE, 4 place Jussieu, 75005 Paris, France  M. de Bony
de Lavergne IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France
J. Borowska  Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr.
15, D 12489 Berlin, Germany  R. Brose Dublin Institute for Advanced Studies,
31 Fitzwilliam Place, Dublin 2, Ireland  A. Brown  University of Oxford,
Department of Physics, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH,
UK  F. Brun IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France
B. Bruno  Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre
for Astroparticle Physics, Nikolaus-Fiebiger-Str. 2, 91058 Erlangen, Germany
C. Burger-Scheidlin  Dublin Institute for Advanced Studies, 31 Fitzwilliam
Place, Dublin 2, Ireland  S. Casanova Instytut Fizyki Ja̧drowej PAN, ul.
Radzikowskiego 152, 31-342 Kraków, Poland  J. Celic  Friedrich-Alexander-
Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics,
Nikolaus-Fiebiger-Str. 2, 91058 Erlangen, Germany  M. Cerruti Université de
Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France  T. Chand
Centre for Space Research, North-West University, Potchefstroom 2520, South
Africa  S. Chandra  Centre for Space Research, North-West University,
Potchefstroom 2520, South Africa  A. Chen School of Physics, University of the
Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, Johannesburg, 2050 South
Africa  J. Chibueze  Centre for Space Research, North-West University,
Potchefstroom 2520, South Africa  O. Chibueze  Centre for Space Research,
North-West University, Potchefstroom 2520, South Africa  G. Cotter University
of Oxford, Department of Physics, Denys Wilkinson Building, Keble Road, Oxford
OX1 3RH, UK  P. Cristofari  Laboratoire Univers et Théories, Observatoire de
Paris, Université PSL, CNRS, Université Paris Cité, 5 Pl. Jules Janssen, 92190
Meudon, France  J. Devin  Laboratoire Univers et Particules de Montpellier,
Université Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, F-34095
Montpellier Cedex 5, France  A. Djannati-Ataï Université de Paris, CNRS,
Astroparticule et Cosmologie, F-75013 Paris, France  J. Djuvsland  Max-Planck-
Institut für Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany  A.
Dmytriiev  Centre for Space Research, North-West University, Potchefstroom
2520, South Africa  K. Egberts  Institut für Physik und Astronomie,
Universität Potsdam, Karl-Liebknecht-Strasse 24/25, D 14476 Potsdam, Germany
S. Einecke  School of Physical Sciences, University of Adelaide, Adelaide
5005, Australia  K. Feijen  Université de Paris, CNRS, Astroparticule et
Cosmologie, F-75013 Paris, France  M. Filipovic  School of Science, Western
Sydney University, Locked Bag 1797, Penrith South DC, NSW 2751, Australia  G.
Fontaine Laboratoire Leprince-Ringuet, École Polytechnique, CNRS, Institut
Polytechnique de Paris, F-91128 Palaiseau, France  S. Funk Friedrich-
Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle
Physics, Nikolaus-Fiebiger-Str. 2, 91058 Erlangen, Germany  S. Gabici
Université de Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France
Y.A. Gallant  Laboratoire Univers et Particules de Montpellier, Université
Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, F-34095 Montpellier
Cedex 5, France  J.F. Glicenstein IRFU, CEA, Université Paris-Saclay, F-91191
Gif-sur-Yvette, France  J. Glombitza  Friedrich-Alexander-Universität
Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Nikolaus-
Fiebiger-Str. 2, 91058 Erlangen, Germany  G. Grolleron  Sorbonne Université,
CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE, 4
place Jussieu, 75005 Paris, France  L. Haerer  Max-Planck-Institut für
Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany  B. Heß  Institut für
Astronomie und Astrophysik, Universität Tübingen, Sand 1, D 72076 Tübingen,
Germany  J.A. Hinton  Max-Planck-Institut für Kernphysik, Saupfercheckweg 1,
69117 Heidelberg, Germany  W. Hofmann  Max-Planck-Institut für Kernphysik,
Saupfercheckweg 1, 69117 Heidelberg, Germany  T. L. Holch Deutsches
Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany  D. Horns
Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, D
22761 Hamburg, Germany  Zhiqiu Huang Max-Planck-Institut für Kernphysik,
Saupfercheckweg 1, 69117 Heidelberg, Germany  M. Jamrozy Obserwatorium
Astronomiczne, Uniwersytet Jagielloński, ul. Orla 171, 30-244 Kraków, Poland
F. Jankowsky  Landessternwarte, Universität Heidelberg, Königstuhl, D 69117
Heidelberg, Germany  I. Jung-Richardt  Friedrich-Alexander-Universität
Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Nikolaus-
Fiebiger-Str. 2, 91058 Erlangen, Germany  E. Kasai  University of Namibia,
Department of Physics, Private Bag 13301, Windhoek 10005, Namibia  K.
Katarzyński  Institute of Astronomy, Faculty of Physics, Astronomy and
Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun,
Poland  R. Khatoon  Centre for Space Research, North-West University,
Potchefstroom 2520, South Africa  B. Khélifi Université de Paris, CNRS,
Astroparticule et Cosmologie, F-75013 Paris, France  W. Kluźniak  Nicolaus
Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18,
00-716 Warsaw, Poland  Nu. Komin Corresponding author: contact.hess@hess-
experiment.eu School of Physics, University of the Witwatersrand, 1 Jan Smuts
Avenue, Braamfontein, Johannesburg, 2050 South Africa  K. Kosack  IRFU, CEA,
Université Paris-Saclay, F-91191 Gif-sur-Yvette, France  D. Kostunin Deutsches
Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany  A. Kundu
Centre for Space Research, North-West University, Potchefstroom 2520, South
Africa  R.G. Lang  Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen
Centre for Astroparticle Physics, Nikolaus-Fiebiger-Str. 2, 91058 Erlangen,
Germany  S. Le Stum  Aix Marseille Université, CNRS/IN2P3, CPPM, Marseille,
France  A. Lemière  Université de Paris, CNRS, Astroparticule et Cosmologie,
F-75013 Paris, France  M. Lemoine-Goumard Université Bordeaux, CNRS, LP2I
Bordeaux, UMR 5797, F-33170 Gradignan, France  J.-P. Lenain Sorbonne
Université, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes
Energies, LPNHE, 4 place Jussieu, 75005 Paris, France  F. Leuschner Institut
für Astronomie und Astrophysik, Universität Tübingen, Sand 1, D 72076
Tübingen, Germany  J. Mackey Dublin Institute for Advanced Studies, 31
Fitzwilliam Place, Dublin 2, Ireland  V. Marandon IRFU, CEA, Université Paris-
Saclay, F-91191 Gif-sur-Yvette, France  G. Martí-Devesa Universität Innsbruck,
Institut für Astro- und Teilchenphysik, Technikerstraße 25, 6020 Innsbruck,
Austria  R. Marx Landessternwarte, Universität Heidelberg, Königstuhl, D 69117
Heidelberg, Germany  A. Mehta  Deutsches Elektronen-Synchrotron DESY,
Platanenallee 6, 15738 Zeuthen, Germany  A. Mitchell Friedrich-Alexander-
Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics,
Nikolaus-Fiebiger-Str. 2, 91058 Erlangen, Germany  R. Moderski  Nicolaus
Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18,
00-716 Warsaw, Poland  M.O. Moghadam  Institut für Physik und Astronomie,
Universität Potsdam, Karl-Liebknecht-Strasse 24/25, D 14476 Potsdam, Germany
L. Mohrmann Corresponding author<EMAIL_ADDRESS>Max-Planck-
Institut für Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany  A.
Montanari Landessternwarte, Universität Heidelberg, Königstuhl, D 69117
Heidelberg, Germany  E. Moulin IRFU, CEA, Université Paris-Saclay, F-91191
Gif-sur-Yvette, France  M. de Naurois  Laboratoire Leprince-Ringuet, École
Polytechnique, CNRS, Institut Polytechnique de Paris, F-91128 Palaiseau,
France  J. Niemiec Instytut Fizyki Ja̧drowej PAN, ul. Radzikowskiego 152,
31-342 Kraków, Poland  S. Ohm Deutsches Elektronen-Synchrotron DESY,
Platanenallee 6, 15738 Zeuthen, Germany  L. Olivera-Nieto Max-Planck-Institut
für Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany  E. de Ona
Wilhelmi  Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738
Zeuthen, Germany  M. Ostrowski Obserwatorium Astronomiczne, Uniwersytet
Jagielloński, ul. Orla 171, 30-244 Kraków, Poland  S. Panny Universität
Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstraße 25, 6020
Innsbruck, Austria  U. Pensec  Sorbonne Université, CNRS/IN2P3, Laboratoire de
Physique Nucléaire et de Hautes Energies, LPNHE, 4 place Jussieu, 75005 Paris,
France  G. Peron  Université de Paris, CNRS, Astroparticule et Cosmologie,
F-75013 Paris, France  G. Pühlhofer Institut für Astronomie und Astrophysik,
Universität Tübingen, Sand 1, D 72076 Tübingen, Germany  A. Quirrenbach
Landessternwarte, Universität Heidelberg, Königstuhl, D 69117 Heidelberg,
Germany  S. Ravikularaman  Université de Paris, CNRS, Astroparticule et
Cosmologie, F-75013 Paris, France  Max-Planck-Institut für Kernphysik,
Saupfercheckweg 1, 69117 Heidelberg, Germany  M. Regeard  Université de Paris,
CNRS, Astroparticule et Cosmologie, F-75013 Paris, France  A. Reimer
Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstraße
25, 6020 Innsbruck, Austria  O. Reimer Universität Innsbruck, Institut für
Astro- und Teilchenphysik, Technikerstraße 25, 6020 Innsbruck, Austria  H. Ren
Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, 69117 Heidelberg,
Germany  M. Renaud  Laboratoire Univers et Particules de Montpellier,
Université Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, F-34095
Montpellier Cedex 5, France  B. Reville Max-Planck-Institut für Kernphysik,
Saupfercheckweg 1, 69117 Heidelberg, Germany  F. Rieger  Max-Planck-Institut
für Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany  G. Rowell School
of Physical Sciences, University of Adelaide, Adelaide 5005, Australia  B.
Rudak Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul.
Bartycka 18, 00-716 Warsaw, Poland  E. Ruiz-Velasco Max-Planck-Institut für
Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany  K. Sabri
Laboratoire Univers et Particules de Montpellier, Université Montpellier,
CNRS/IN2P3, CC 72, Place Eugène Bataillon, F-34095 Montpellier Cedex 5, France
V. Sahakian Yerevan Physics Institute, 2 Alikhanian Brothers St., 0036
Yerevan, Armenia  H. Salzmann  Institut für Astronomie und Astrophysik,
Universität Tübingen, Sand 1, D 72076 Tübingen, Germany  A. Santangelo
Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, D 72076
Tübingen, Germany  M. Sasaki Friedrich-Alexander-Universität Erlangen-
Nürnberg, Erlangen Centre for Astroparticle Physics, Nikolaus-Fiebiger-Str. 2,
91058 Erlangen, Germany  J. Schäfer  Friedrich-Alexander-Universität Erlangen-
Nürnberg, Erlangen Centre for Astroparticle Physics, Nikolaus-Fiebiger-Str. 2,
91058 Erlangen, Germany  F. Schüssler IRFU, CEA, Université Paris-Saclay,
F-91191 Gif-sur-Yvette, France  H.M. Schutte Centre for Space Research, North-
West University, Potchefstroom 2520, South Africa  H. Sol  Laboratoire Univers
et Théories, Observatoire de Paris, Université PSL, CNRS, Université Paris
Cité, 5 Pl. Jules Janssen, 92190 Meudon, France  S. Spencer Friedrich-
Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle
Physics, Nikolaus-Fiebiger-Str. 2, 91058 Erlangen, Germany  Ł. Stawarz
Obserwatorium Astronomiczne, Uniwersytet Jagielloński, ul. Orla 171, 30-244
Kraków, Poland  S. Steinmassl Max-Planck-Institut für Kernphysik,
Saupfercheckweg 1, 69117 Heidelberg, Germany  C. Steppa  Institut für Physik
und Astronomie, Universität Potsdam, Karl-Liebknecht-Strasse 24/25, D 14476
Potsdam, Germany  K. Streil  Friedrich-Alexander-Universität Erlangen-
Nürnberg, Erlangen Centre for Astroparticle Physics, Nikolaus-Fiebiger-Str. 2,
91058 Erlangen, Germany  I. Sushch Centre for Space Research, North-West
University, Potchefstroom 2520, South Africa  A.M. Taylor Deutsches
Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany  R.
Terrier Université de Paris, CNRS, Astroparticule et Cosmologie, F-75013
Paris, France  M. Tsirou  Deutsches Elektronen-Synchrotron DESY, Platanenallee
6, 15738 Zeuthen, Germany  N. Tsuji RIKEN, 2-1 Hirosawa, Wako, Saitama
351-0198, Japan  C. van Eldik Friedrich-Alexander-Universität Erlangen-
Nürnberg, Erlangen Centre for Astroparticle Physics, Nikolaus-Fiebiger-Str. 2,
91058 Erlangen, Germany  M. Vecchi  Kapteyn Astronomical Institute, University
of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands  C. Venter
Centre for Space Research, North-West University, Potchefstroom 2520, South
Africa  J. Vink  GRAPPA, Anton Pannekoek Institute for Astronomy, University
of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands  S.J.
Wagner Landessternwarte, Universität Heidelberg, Königstuhl, D 69117
Heidelberg, Germany  R. White  Max-Planck-Institut für Kernphysik,
Saupfercheckweg 1, 69117 Heidelberg, Germany  A. Wierzcholska Instytut Fizyki
Ja̧drowej PAN, ul. Radzikowskiego 152, 31-342 Kraków, Poland  M. Zacharias
Landessternwarte, Universität Heidelberg, Königstuhl, D 69117 Heidelberg,
Germany  Centre for Space Research, North-West University, Potchefstroom 2520,
South Africa  A.A. Zdziarski Nicolaus Copernicus Astronomical Center, Polish
Academy of Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland  A. Zech
Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS,
Université Paris Cité, 5 Pl. Jules Janssen, 92190 Meudon, France  N. Żywucka
Centre for Space Research, North-West University, Potchefstroom 2520, South
Africa

###### 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 high-energy cosmic rays, potentially
beyond PeV energies. Here, we report the detection of very-high-energy
$\gamma$-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 $\gamma$-ray emission from 30 Dor C, the only
superbubble detected at TeV energies presently. The $\gamma$-ray luminosity
above $0.5\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$ of both sources is
$(2-3)\times 10^{35}\,\mathrm{erg}\,\mathrm{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 $\gamma$-ray emission from each source is extended with a
significance of $>3\sigma$ and a Gaussian width of about
$30\text{\,}\mathrm{p}\mathrm{c}$. For 30 Dor C, a connection between the
$\gamma$-ray emission and the nonthermal X-ray emission appears likely.
Different interpretations of the $\gamma$-ray signal from R136 are discussed.

Young star clusters (1833) — Massive stars (732) — Large Magellanic Cloud
(903) — Gamma-ray astronomy (628)

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

tablenum SIXtablenum

##  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 $\gamma$-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$\gamma$ 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 $\gamma$-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 $\gamma$-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\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$) $\gamma$-ray domain.

In fact, H.E.S.S. Collaboration et al. (2022) have recently been able to
associate the VHE $\gamma$-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
$\gamma$-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 $\gamma$-ray emission with the star cluster is not always firm. The
detection and observation of further YMCs with VHE $\gamma$-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 $\gamma$-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 $\gamma$-rays so far111The $\gamma$-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. (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
$12\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ (e.g. Massey & Hunter, 1998;
Brands et al., 2022), R136 is also relatively young, implying that only 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 paper, we report the discovery of VHE $\gamma$-ray emission from the
direction of R136 with the High Energy Stereoscopic System (H.E.S.S.),
achieved through a multicomponent, likelihood-based modeling of the spatial
and spectral distribution of $\gamma$-ray-like events. Moreover, from the same
analysis, we provide updated results on the SB 30 Dor C. Throughout the paper,
we assume a distance to the LMC of $50\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$
(Pietrzyński et al., 2013).

##  2 Data Analysis

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

H.E.S.S. has collected an extensive data set on the Tarantula Nebula region.
Observations are taken in individual “runs” of usually
$28\text{\,}\mathrm{min}$ duration. Here, we have analyzed a total of 794 runs
taken with all four of the initial telescopes, corresponding to an observation
time of $\approx$$360\text{\,}\mathrm{h}$. Of these, 301 runs
($\approx$$138\text{\,}\mathrm{h}$) have been taken between 2004 December 30
and 2013 February 7, while the remaining 493 runs
($\approx$$222\text{\,}\mathrm{h}$) date to between 2017 October 14 and 2022
February 2.222The 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\text{\,}{\mathrm{m}}^{2}$
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).
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 $\gamma$-rays proceeds via measuring the Cherenkov radiation
that is emitted by secondary particles in the extensive air shower that is
launched when the primary $\gamma$-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.5^{\circ}$. Owing to the relatively strict
selection cuts, the achieved energy threshold for the data set is
$0.5\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$. 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 Gammapy v0.18.2 (Donath et al., 2023; Deil et al., 2020), employing a
three-dimensional binned likelihood fit (Mattox et al., 1996) in a
$5^{\circ}\times 5^{\circ}$ region of interest (ROI), with spatial pixels of
size $0.02^{\circ}\times 0.02^{\circ}$ (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
$\gamma$-ray emission, following an iterative procedure in which we continue
to add sources to the ROI model until no significant $\gamma$-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_{\mathrm{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).

Figure 1:  $\gamma$-ray flux maps of the target region of the analysis. The
maps show the $\gamma$-ray flux $F$ in units of
$10^{-8}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}\,\mathrm{sr}^{-1}$, integrated
above an energy of $0.5\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$, assuming a
power-law spectrum with index $-2.5$. Smoothing with a top-hat kernel of
$0.07^{\circ}$ radius has been applied. (a) Entire emission. (b) Residual
emission after subtraction of the emission from N 157B predicted by the best-
fit ROI model. Dashed blue lines show flux contours at $(2.5/7.5/12.5)\times
10^{-8}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}\,\mathrm{sr}^{-1}$ in panel _a_ and
$(1.5/3)\times 10^{-8}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}\,\mathrm{sr}^{-1}$
in panel _b_. Pixels with a negative excess after background subtraction are
clipped at zero. The light green contour lines denote H$\alpha$ emission as
inferred by SHASSA (Gaustad et al., 2001).

##  3 Results

We show a $\gamma$-ray flux map of the Tarantula Nebula region in Fig. 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 $\gamma$-ray
source in this region is the pulsar wind nebula (PWN) N 157B, 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 $\gamma$-ray
emission, adding source components – beginning with a model for N 157B – 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 paper. This is illustrated in Fig. 1(b), where
the emission from N 157B 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

_N 157B_ — We do not investigate N 157B in detail in this paper, but note that
our analysis yields a nonzero extension of
$\sigma_{\mathrm{Gauss}}^{\mathrm{N\,157B}}=(0.82\pm 0.20_{\mathrm{stat}}\pm
0.18_{\mathrm{sys}})^{\prime}$ for this source. We caution, however, that the
extended source model is preferred over a pointlike model by only
$\approx$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 N 157B is less extended than suggested
by our best fit. We therefore do not claim an extension and provide an upper
limit of $\sigma_{\mathrm{Gauss}}^{\mathrm{N\,157B}}<1.14^{\prime}$ (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 N 157B. The obtained spectrum for
N 157B, shown in Fig. 7 in Appendix B, is compatible with our previously
published result.

_30 Dor C_ — For the first time, we find HESS J0535$-$691, associated with 30
Dor C, to be an extended $\gamma$-ray source. The best-fit Gaussian width is
$\sigma_{\mathrm{Gauss}}^{\mathrm{30\,Dor\,C}}=(1.91\pm
0.40_{\mathrm{stat}}\pm 0.20_{\mathrm{sys}})^{\prime}$, which corresponds to
$(27.8\pm 5.8_{\mathrm{stat}}\pm 2.9_{\mathrm{sys}})\,\mathrm{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 Fig.
2. The best-fit position deviates by $1.1^{\prime}$ 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_{\mathrm{30\,Dor\,C}}=-2.57\pm 0.09_{\mathrm{stat}}$.

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
$\alpha$–$\zeta$ found by Lortet & Testor (1984). The pink contour line
outlines the SB 30 Dor C as visible in X-rays
($12\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$) 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 $\gamma$-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 12CO(1–0) emission measured with the Atacama Large Millimeter/submillimeter
Array (Yamane et al., 2021).

_R136_ — Lastly, HESS J0538$-$691 is detected as a new $\gamma$-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^{\prime\prime}$
(see Fig. 3). Because there is no other plausible counterpart, we associate
HESS J0538$-$691 with R136.333The 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 Sect. 4. Similarly to 30 Dor C, we
find a preference (3.1$\sigma$) for an extended source, with a Gaussian width
of $\sigma_{\mathrm{Gauss}}^{\mathrm{R136}}=(2.30\pm
0.54_{\mathrm{stat}}\,{}_{-0.21}^{+0.26}|_{\mathrm{sys}})^{\prime}$, or
$(33.5\pm
7.9_{\mathrm{stat}}\,{}_{-3.1}^{+3.8}|_{\mathrm{sys}})\,\mathrm{pc}$. For the
energy spectrum, we find a power-law spectral index
$\Gamma_{\mathrm{R136}}=-2.54\pm 0.15_{\mathrm{stat}}$.

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
$\gamma$-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 12CO(1–0) emission
(Johansson et al., 1998), while the light red contours indicate the positions
of dense “knots” of molecular gas identified through 12CO(2–1) emission
(Kalari et al., 2018).

###  3.2 Spectral results and energy requirements

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

Figure 4:  $\gamma$-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\,(\mathrm{d}N/\mathrm{d}E)\,\cdot\,4\pi d^{2}$, with $d$ the
distance to the source. For comparison, the $\gamma$-ray luminosity of the YMC
Westerlund 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).

Fig. 4 also shows the $\gamma$-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_{\mathrm{p}}^{\mathrm{30\,Dor\,C}}$$\,=\,$$-2.64\pm
0.08_{\mathrm{stat}}$ and
$\Gamma_{\mathrm{p}}^{\mathrm{R136}}$$\,=\,$$-2.59\pm 0.13_{\mathrm{stat}}$.
Assuming an extrapolation of the particle spectrum to
$1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ with the same spectral index, this
would imply a total energy requirement for protons of
$W_{\mathrm{p}}^{\mathrm{30\,Dor\,C}}\approx$2.1\text{\times}{10}^{53}$\,(n/$1\text{\,}{\mathrm{cm}}^{-3}$)^{-1}\,$\mathrm{e}\mathrm{r}\mathrm{g}$$
and
$W_{\mathrm{p}}^{\mathrm{R136}}\approx$1.1\text{\times}{10}^{53}$\,(n/$1\text{\,}{\mathrm{cm}}^{-3}$)^{-1}\,$\mathrm{e}\mathrm{r}\mathrm{g}$$,
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
$\gamma$-ray flux in the $110\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ 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
$\sim$$1\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$. In this case, one obtains a
lower requirement of
$W_{\mathrm{p}}^{\mathrm{30\,Dor\,C}}\approx$1.4\text{\times}{10}^{52}$\,(n/$1\text{\,}{\mathrm{cm}}^{-3}$)^{-1}\,$\mathrm{e}\mathrm{r}\mathrm{g}$$
(or, for comparison, although no limit from _Fermi_ -LAT is available,
$W_{\mathrm{p}}^{\mathrm{R136}}\approx$9.7\text{\times}{10}^{51}$\,(n/$1\text{\,}{\mathrm{cm}}^{-3}$)^{-1}\,$\mathrm{e}\mathrm{r}\mathrm{g}$$).
These values can be treated as lower limits for the required energy in
protons.

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_{\mathrm{e}}^{\mathrm{30\,Dor\,C}}=-3.27\pm
0.11_{\mathrm{stat}}$ and $\Gamma_{\mathrm{e}}^{\mathrm{R136}}=-3.19\pm
0.17_{\mathrm{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\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$, and
taking into account the energy-dependent cooling due to IC scattering, we find
a minimum required injection power for electrons of
$L_{\mathrm{e}}^{\mathrm{30\,Dor\,C}}\approx$6.4\text{\times}{10}^{36}\text{\,}\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$
for 30 Dor C and
$L_{\mathrm{e}}^{\mathrm{R136}}\approx$3.9\text{\times}{10}^{36}\text{\,}\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$
for R136. In the presence of a magnetic field of
$B=$5\text{\,}\mathrm{\SIUnitSymbolMicro G}$$, a value roughly representative
for the LMC as a whole (Gaensler et al., 2005), additional synchrotron losses
lead to larger requirements of
$L_{\mathrm{e}}^{\mathrm{30\,Dor\,C}}\approx$8.5\text{\times}{10}^{36}\text{\,}\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$
and
$L_{\mathrm{e}}^{\mathrm{R136}}\approx$5.3\text{\times}{10}^{36}\text{\,}\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$,
respectively. For $B=$15\text{\,}\mathrm{\SIUnitSymbolMicro G}$$, as derived
by H.E.S.S. Collaboration (2015) within a leptonic model for 30 Dor C, the
corresponding values are
$L_{\mathrm{e}}^{\mathrm{30\,Dor\,C}}\approx$1.5\text{\times}{10}^{37}\text{\,}\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$
and
$L_{\mathrm{e}}^{\mathrm{R136}}\approx$9.5\text{\times}{10}^{36}\text{\,}\mathrm{erg}\text{\,}{\mathrm{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
$\gamma$-ray emission and the YMCs LH 90444While 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). and R136. 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
$\gamma$-ray emission from the newly detected source HESS J0538$-$691 will be
discussed at the end of this section.

We provide in Table 4 an overview of the relevant properties of the YMCs and
the surrounding medium, as well as results derived from the H.E.S.S.
$\gamma$-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.

Table 1: YMC properties.

Property | LH 90 | R136 | Ref.
---|---|---|---
Cluster age [$\mathrm{Myr}$] | 4 | 1.5 | 1,2
Wind poweraaAveraged over the cluster lifetime. [${10}^{38}\text{\,}\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$] | 1.5 | 10 | 3,4
Wind velocity [$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$] | 3000 | 3000 | 5,6
Average ISM density [${\mathrm{cm}}^{-3}$] | 100 | 100 | 7–11
Magnetic field [$\mathrm{\SIUnitSymbolMicro G}$] | 15 | 15bbNo literature estimate available. | 12
SB radius [pc] | 74 | 56 |
Termination shock radius [pc] | 7.9 | 8.7 |
2D Gaussian width [pc] | 27.8 | 33.5 |
68% containment radius [pc] | 42.0 | 50.5 |
Spectral index | $-2.57$ | $-2.54$ |
FluxccThe integrated flux and luminosity are given above an energy of $0.5\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$, which is the energy threshold of the H.E.S.S. analysis. [${10}^{-13}\text{\,}{\mathrm{cm}}^{-2}\text{\,}{\mathrm{s}}^{-1}$] | 4.8 | 3.6 |
LuminosityccThe integrated flux and luminosity are given above an energy of $0.5\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$, which is the energy threshold of the H.E.S.S. analysis. [${10}^{35}\text{\,}\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$] | 2.9 | 2.2 |
Req. power ($pp$) [${10}^{36}\text{\,}\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$] | 1.1 | 2.0 |
Req. power (IC) [${10}^{37}\text{\,}\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$] | 1.5 | 0.95 |

Note. — 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 $\gamma$-ray
emission measured with H.E.S.S. (cf. Sect. 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\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$ (cf. Sect. 3),
considered no escape or radiation losses, and assumed a continuous injection
over the lifetime of the respective cluster.

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.

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
($\gtrsim$3000\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$) shock, which
cannot be the forward shock of the expanding SB, as H$\alpha$ observations
indicate a shell expanding at only
$\lesssim$100\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$. 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
($\sim${10}^{-3}\text{\,}{\mathrm{cm}}^{-3}$$) 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\text{\,}\mathrm{\SIUnitSymbolMicro G}$; H.E.S.S. Collaboration 2015).
Based on its total mass ($\sim$$$2.2\text{\times}{10}^{4}$\,M_{\odot}$;
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$\alpha$ 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 (cf. Johansson et al. 1998; Kalari et al. 2018 and Fig. 3). For the
purpose of our discussion, we nevertheless compute and list in Table 4 the
expected size of a putative SB around R136 and its termination shock.

It is interesting that in terms of their $\gamma$-ray emission, 30 Dor C and
R136 look similar: both appear extended with a Gaussian width of around
$30\text{\,}\mathrm{p}\mathrm{c}$ 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 $\gamma$-ray emission must originate from different
processes or be created at different sites.

In both cases, the $\gamma$-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
$\gamma$-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 $\gamma$-ray emission. In that case, on the other hand, one would expect
the centroid of the $\gamma$-ray emission to coincide with the positions of
the densest gas clouds (see Figs. 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
$\gamma$-ray emission (see Table 4) 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 $\sim$10%. While this would
be a surprisingly large efficiency for a leptonic accelerator, we stress that
the wind power listed in Table 4
($1.5\text{\times}{10}^{38}\text{\,}\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$)
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
$\sim$5\text{\times}{10}^{38}\text{\,}\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$,
which alleviates the requirements slightly. Nevertheless, it appears
challenging in this scenario to entirely explain the $\gamma$-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 $\gamma$-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
$\sim$1% in the leptonic case.

Considering, on the other hand, a hadronic origin of the $\gamma$-ray
emission, we deduce a minimum ratio between required power in protons and
power provided by the cluster wind of $\sim$0.7% for 30 Dor C and of
$\sim$0.2% for R136. We assume in this case that the $\gamma$-ray emission
originates from interactions of CRs in the dense gas clouds that surround the
respective SBs (see Figs. 2 and 3); this is in line with the extension of the
$\gamma$-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
$\mathcal{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_{\mathrm{e}}\approx-3.2$) are consistent with a synchrotron-cooled
injection spectrum of $\Gamma_{\mathrm{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_{\mathrm{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 $\gamma$-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 $\gamma$-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 $\sim${10}^{6}\text{\,}\mathrm{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 $\gamma$-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 $\gamma$-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 $\sim${10}^{37}\text{\,}\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$$
(cf. Sect. 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 $\gamma$-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 $\gamma$-ray sources. The brightest one, HESS
J0537$-$691, is associated with the PWN N 157B. The other two sources can be
associated with YMCs and/or their surrounding SBs. Both sources are extremely
luminous in $\gamma$-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 $\gamma$-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 $\gamma$-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 $\gamma$-ray source associated
with the YMC R136, labeled HESS J0538$-$691\. This source is similar in terms
of spatial extension and $\gamma$-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
$\gamma$-ray emission appear viable.

The detection of $\gamma$-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 $\gamma$-ray
emission and its interpretation. Our analysis thus provides crucial
information to understand the ability of YMCs to accelerate CRs better.

## Acknowledgements

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.

##  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
$\mathrm{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).

Table 2: Summary of analysis settings.

Setting | Value
---|---
ROI center (J2000) | R.A. $5^{\mathrm{h}}35^{\mathrm{m}}28.25^{\mathrm{s}}$
Decl. $-69^{\circ}16^{\prime}13.08^{\prime\prime}$
ROI size | $5^{\circ}\times 5^{\circ}$
Spatial pixel size | $0.02^{\circ}\times 0.02^{\circ}$
Maximum offset angle | $1.5^{\circ}$
Energy binning | 16 bins decade-1
Energy range |  $0.5\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$–$100\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$

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

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 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 N
132D (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^{\prime}=R\cdot(E/$1\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$)^{-\delta}$.
We show in Fig. 5 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 Fig. 5 shows
that this is very nearly the case for this analysis, indicating that an
excellent description of the hadronic background has been achieved. In Fig. 5,
we provide a zoom-in of the target region of this work, with the location of
three detected $\gamma$-ray sources indicated (cf. Sect. 3 and Appendix B).

Figure 5:  (a) Significance map for the entire ROI, above an energy of
$0.5\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$, after adjustment of the hadronic
background model. The map has been smoothed with a top-hat kernel of
$0.07^{\circ}$ radius, as indicated in the bottom right corner. The dark gray
lines enclose regions that have been excluded in the fit of the background,
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 $\gamma$-ray excess profiles have been computed (cf. Fig.
9).

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.

##  Appendix B Source Fitting

In this appendix, we detail the modeling of the $\gamma$-ray sources in the
Tarantula Nebula region as displayed in Fig. 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.

Table 3: Best-fit parameters of the $\gamma$-ray source models.

Parameter | Unit | Value
---|---|---
N 157B / HESS J0537$-$691
R.A. | deg |  $84.4394\pm 0.0048_{\mathrm{stat}}$ ($5^{\mathrm{h}}37^{\mathrm{m}}45.5^{\mathrm{s}}\pm 1.1^{\mathrm{s}}_{\mathrm{stat}}$)
Decl. | deg |  $-69.1713\pm 0.0016_{\mathrm{stat}}$ ($-69^{\circ}10^{\prime}17^{\prime\prime}\pm 6_{\mathrm{stat}}^{\prime\prime}$)
$\sigma_{\mathrm{Gauss}}$ | deg | $0.0137\pm 0.0033_{\mathrm{stat}}\pm 0.0030_{\mathrm{sys}}$
$\phi_{0}$ | $10^{-13}\,\mathrm{TeV}^{-1}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}$ | $8.69\pm 0.56_{\mathrm{stat}}\pm 0.85_{\mathrm{sys}}$
$\alpha$ | — | $2.03\pm 0.07_{\mathrm{stat}}\pm 0.08_{\mathrm{sys}}$
$\beta$ | — | $0.311\pm 0.037_{\mathrm{stat}}$
30 Dor C / HESS J0535$-$691
R.A. | deg |  $84.021\pm 0.018_{\mathrm{stat}}$ ($5^{\mathrm{h}}36^{\mathrm{m}}5.0^{\mathrm{s}}\pm 4.3^{\mathrm{s}}_{\mathrm{stat}}$)
Decl. | deg |  $-69.197\pm 0.006_{\mathrm{stat}}$ ($-69^{\circ}11^{\prime}49^{\prime\prime}\pm 22_{\mathrm{stat}}^{\prime\prime}$)
$\sigma_{\mathrm{Gauss}}$ | deg | $0.0319\pm 0.0066_{\mathrm{stat}}\pm 0.0034_{\mathrm{sys}}$
$\phi_{0}$ | $10^{-13}\,\mathrm{TeV}^{-1}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}$ | $2.54\pm 0.37_{\mathrm{stat}}\,{}_{-0.40}^{+0.44}|_{\mathrm{sys}}$
$\Gamma$ | — | $2.57\pm 0.09_{\mathrm{stat}}$
R136 / HESS J0538$-$691
R.A. | deg |  $84.692\pm 0.038_{\mathrm{stat}}$ ($5^{\mathrm{h}}38^{\mathrm{m}}46^{\mathrm{s}}\pm 9^{\mathrm{s}}_{\mathrm{stat}}$)
Decl. | deg |  $-69.103\pm 0.013_{\mathrm{stat}}$ ($-69^{\circ}06^{\prime}11^{\prime\prime}\pm 47_{\mathrm{stat}}^{\prime\prime}$)
$\sigma_{\mathrm{Gauss}}$ | deg | $0.0384\pm 0.0090_{\mathrm{stat}}\,{}_{-0.0037}^{+0.0045}|_{\mathrm{sys}}$
$\phi_{0}$ | $10^{-13}\,\mathrm{TeV}^{-1}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}$ | $1.90\pm 0.58_{\mathrm{stat}}\,{}_{-0.38}^{+0.45}|_{\mathrm{sys}}$
$\Gamma$ | — | $2.54\pm 0.15_{\mathrm{stat}}$

Note. — Coordinates are given for the epoch J2000. The flux normalization
$\phi_{0}$ is specified at a reference energy of
$E_{0}=$1\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$. None of the sources
exhibit significant emission above $30\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$
(cf. Figs. 4 and 7). Statistical errors are at a 68% confidence level.
Systematic errors have been derived as explained in Sect. 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^{\prime\prime}$–$20^{\prime\prime}$ (Gillessen, 2004).

###  B.1 Modeling Procedure and Results

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

Figure 6:  Maps of residual significance of the analysis target region. (a)
After the background model adjustment (cf. Appendix A). (b) After adding N
157B 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.

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

$\frac{\mathrm{d}N}{\mathrm{d}E}=\phi_{0}\left(\frac{E}{E_{0}}\right)^{-\alpha-\beta\ln(E/E_{0})}\,,$
(B1)

where the reference energy $E_{0}=$1\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$
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 Fig. 8, the
best-fit position of the model is in excellent agreement with the center of
the PWN as seen in X-rays.

Figure 7:  Spectral energy distribution for N 157B (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).  Figure 8:
Chandra X-ray image of N 157B (Chen et al., 2006), with properties of the
best-fit $\gamma$-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.

After the inclusion of N 157B, the largest excess emission is visible around
30 Dor C (cf. Fig. 6(b)), which we add as a second source to our model. This
leads to an improvement in the fit statistic of $\mathrm{TS}\approx 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
model the energy spectrum with a power-law function,

$\frac{\mathrm{d}N}{\mathrm{d}E}=\phi_{0}\left(\frac{E}{E_{0}}\right)^{-\Gamma}$
(B2)

(with $E_{0}=$1\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$ again fixed), and
find it to be in good agreement with our earlier measurement (H.E.S.S.
Collaboration, 2015), as shown in Fig. 7. A log-parabola model for 30 Dor C
improves the fit only marginally ($\mathrm{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 $\mathrm{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 Fig. 7.

With N 157B, 30 Dor C, and R136 included in the ROI model, the residual
significance map (see Fig. 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 Fig. 9 one-dimensional profiles along the slice displayed in Fig. 5. 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.

Figure 9:  Count profiles for the slice across the target region indicated in
Fig. 5. The black data points, identical in all panels, show the observed
number of events in $0.03^{\circ}$ 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 $\gamma$-ray
source models. (a) After adding N 157B 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 $\chi^{2}$ test between
the total model prediction and the observed data.

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 N 157B with a
pointlike one. Surprisingly, given the relatively small statistical error of
$0.0033^{\circ}$ (at a 68% confidence level) on the measured extension of
$\sigma_{\mathrm{Gauss}}^{\mathrm{N\,157B}}=0.0137^{\circ}$, we find the
pointlike model to be disfavored by only $\approx$1.3$\sigma$. We attribute
this to the presence of the two other $\gamma$-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 N 157B to flatten toward
smaller values. Despite the marginal preference for the Gaussian spatial model
for N 157B, 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 N 157B is modeled as
a pointlike source, we obtain
$\sigma_{\mathrm{Gauss}}^{\mathrm{30\,Dor\,C}}=(0.0334\pm
0.0066_{\mathrm{stat}})^{\circ}$ and
$\sigma_{\mathrm{Gauss}}^{\mathrm{R136}}=(0.0431\pm
0.0089_{\mathrm{stat}})^{\circ}$, which is compatible with our default result
(cf. Table 3) within the uncertainties. Interpreted as an upper limit, our
analysis yields $\sigma_{\mathrm{Gauss}}^{\mathrm{N\,157B}}<0.0190^{\circ}$ at
a 95% confidence level (statistical uncertainties only).

We have also tested a model in which the Gaussian spatial model for N 157B is
allowed to be elongated. While this improves the fit quality considerably when
only N 157B is included as a $\gamma$-ray source ($\mathrm{TS}\approx 9.1$),
this is no longer the case when 30 Dor C and R136 are included as well
($\mathrm{TS}\approx 0.5$). Conversely, the model with all three sources
included is strongly preferred over that with only N 157B, even when modeled
as elongated ($\mathrm{TS}\approx 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 $\mathrm{TS}\approx 50$), implying that also just the
emission around R136 cannot be explained by allowing the model for N 157B to
be elongated.

###  B.2 Estimation of Systematic Uncertainties

In this section, we explain our procedure for estimating the systematic
uncertainties on the $\gamma$-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 order of
$10^{\prime\prime}$–$20^{\prime\prime}$ (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 $\gamma$-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 $\approx$$274.2\text{\,}\mathrm{h}$ of observation
time) on PKS 2155$-$304 taken between 2004 July 14 and 2012 November 10 above
an energy threshold of $0.27\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$.
Observations taken during the very strong outburst of PKS 2155$-$304 in July
2006 (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 (cf. Eq. B1). Inspecting the residual significance map for energies
between $0.27\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$ and
$1\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$, shown in Fig. 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 $\gamma$-ray spectrum of PKS 2155$-$304 (we
obtain $\alpha=3.63\pm 0.02$) leads to a much smaller excess of $\gamma$-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 Fig. 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
Fig. 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.

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.271\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$. (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}$).

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_{\mathrm{Gauss}}$) of
the modeled $\gamma$-ray sources, for which we find PSF-related systematic
uncertainties of $\pm 0.0025^{\circ}$ (N 157B), $\pm 0.0033^{\circ}$ (30 Dor
C), and $(_{-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% (N 157B), 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% (N 157B), 17% (30 Dor C), and 18%
(R136). The energy-scale-related uncertainties are thus of the same magnitude
as the statistical errors.

##  Appendix C Fit of Primary Particle Spectra

To fit primary particle spectra, we employ the naima package (Zabalza, 2015).
Specifically, we use the NaimaSpectralModel 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 $\gamma$-ray emission is due to IC
emission from CR electrons and a hadronic model in which the $\gamma$-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 Eq. B2), with
$E_{0}=$3\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$ 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{\,}{\mathrm{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 $\gamma$-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 fit is again carried out in three dimensions, with the same
spatial models as used in our previous fit (cf. Sect. 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.

Table 4: Best-fit parameters of the primary particle spectrum models.

Parameter | Unit | Value
---|---|---
30 Dor C (Leptonic)
$\phi_{0}$ | $10^{34}\,\mathrm{eV}^{-1}$ | $10.1\pm 1.8_{\mathrm{stat}}$
$\Gamma$ | — | $3.27\pm 0.11_{\mathrm{stat}}$
R136 (Leptonic)
$\phi_{0}$ | $10^{34}\,\mathrm{eV}^{-1}$ | $7.0\pm 2.4_{\mathrm{stat}}$
$\Gamma$ | — | $3.19\pm 0.17_{\mathrm{stat}}$
30 Dor C (Hadronic)
$\phi_{0}$ | $10^{37}\,\mathrm{eV}^{-1}$ | $5.7\pm 1.1_{\mathrm{stat}}$
$\Gamma$ | — | $2.64\pm 0.08_{\mathrm{stat}}$
R136 (Hadronic)
$\phi_{0}$ | $10^{37}\,\mathrm{eV}^{-1}$ | $4.0\pm 1.5_{\mathrm{stat}}$
$\Gamma$ | — | $2.59\pm 0.13_{\mathrm{stat}}$

Note. — Statistical errors are at a 68% confidence level. The normalization
$\phi_{0}$ is specified at an energy of
$E_{0}=$3\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$. 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\text{\,}{\mathrm{cm}}^{-3}$ has
been assumed.

##  Appendix D IC Target Photon Fields

For the spectral modeling of the $\gamma$-ray emission from 30 Dor C and R136
with IC models (see Sect. 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\text{\,}\mathrm{K}$ and $4000\text{\,}\mathrm{K}$ and
energy densities of $2.25\text{\,}\mathrm{eV}\text{\,}{\mathrm{cm}}^{-3}$ and
$0.88\text{\,}\mathrm{eV}\text{\,}{\mathrm{cm}}^{-3}$, 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\text{\,}\mathrm{K}$ and follow Härer et al.
(2023) to derive an energy density of
$10\text{\,}\mathrm{eV}\text{\,}{\mathrm{cm}}^{-3}$. In the same manner, using
the catalog from Brands et al. (2022), we obtain an effective temperature of
$45\,000\text{\,}\mathrm{K}$ and an energy density of
$20\text{\,}\mathrm{eV}\text{\,}{\mathrm{cm}}^{-3}$ for R136. Without the
ultraviolet fields included, the normalization parameters $\phi_{0}$ 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.

## References

  * Aharonian et al. (2019) Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, Nat. Astron., 3, 561, doi: 10.1038/s41550-019-0724-0
  * Ashton et al. (2020) Ashton, T., Backes, M., Balzer, A., et al. 2020, Astropart. Phys., 118, 102425, doi: 10.1016/j.astropartphys.2019.102425
  * Badmaev et al. (2024) Badmaev, D. V., Bykov, A. M., & Kalyashova, M. E. 2024, MNRAS, 527, 3749, doi: 10.1093/mnras/stad3389
  * Bamba et al. (2004) Bamba, A., Ueno, M., Nakajima, H., & Koyama, K. 2004, ApJ, 602, 257, doi: 10.1086/380957
  * Berezinskii et al. (1990) Berezinskii, V. S., Bulanov, S. V., Dogiel, V. A., & Ptuskin, V. S. 1990, Astrophysics of Cosmic Rays, ed. V. L. Ginzburg (North-Holland)
  * Bernlöhr (2008) Bernlöhr, K. 2008, Astropart. Phys., 30, 149, doi: 10.1016/j.astropartphys.2008.07.009
  * Bestenlehner et al. (2020) Bestenlehner, J. M., Crowther, P. A., Caballero-Nieves, S. M., et al. 2020, MNRAS, 499, 1918, doi: 10.1093/mnras/staa2801
  * Brands et al. (2022) Brands, S. A., de Koter, A., Bestenlehner, J. M., et al. 2022, A&A, 663, A36, doi: 10.1051/0004-6361/202142742
  * Bykov et al. (2013) Bykov, A. M., Gladilin, P. E., & Osipov, S. M. 2013, MNRAS, 429, 2755, doi: 10.1093/mnras/sts553
  * Cash (1979) Cash, W. 1979, ApJ, 228, 939, doi: 10.1086/156922
  * Cesarsky & Montmerle (1983) Cesarsky, C. J., & Montmerle, T. 1983, Space Sci. Rev., 36, 173, doi: 10.1007/BF00167503
  * Chen et al. (2006) Chen, Y., Wang, Q. D., Gotthelf, E. V., et al. 2006, ApJ, 651, 237, doi: 10.1086/507017
  * Choudhury et al. (2016) Choudhury, S., Subramaniam, A., & Cole, A. A. 2016, MNRAS, 455, 1855, doi: 10.1093/mnras/stv2414
  * Chu & Kennicutt (1994) Chu, Y.-H., & Kennicutt, R. C. 1994, ApJ, 425, 720, doi: 10.1086/174017
  * Cignoni et al. (2015) Cignoni, M., Sabbi, E., van der Marel, R. P., et al. 2015, ApJ, 811, 76, doi: 10.1088/0004-637X/811/2/76
  * H.E.S.S. Collaboration (2005) H.E.S.S. Collaboration. 2005, A&A, 430, 865, doi: 10.1051/0004-6361:20041853
  * H.E.S.S. Collaboration (2006) —. 2006, A&A, 457, 899, doi: 10.1051/0004-6361:20065351
  * H.E.S.S. Collaboration (2007) —. 2007, ApJ, 664, L71, doi: 10.1086/520635
  * H.E.S.S. Collaboration (2010) —. 2010, A&A, 520, A83, doi: 10.1051/0004-6361/201014484
  * H.E.S.S. Collaboration (2012) —. 2012, A&A, 545, L2, doi: 10.1051/0004-6361/201219906
  * H.E.S.S. Collaboration (2015) —. 2015, Science, 347, 406, doi: 10.1126/science.1261313
  * H.E.S.S. Collaboration (2018) —. 2018, A&A, 610, L17, doi: 10.1051/0004-6361/201732426
  * H.E.S.S. Collaboration (2021) —. 2021, A&A, 655, A7, doi: 10.1051/0004-6361/202141486
  * H.E.S.S. Collaboration et al. (2022) H.E.S.S. Collaboration, Blackwell, R., Braiding, C., et al. 2022, A&A, 666, A124, doi: 10.1051/0004-6361/202244323
  * Crowther et al. (2010) Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, MNRAS, 408, 731, doi: 10.1111/j.1365-2966.2010.17167.x
  * Crowther et al. (2016) Crowther, P. A., Caballero-Nieves, S. M., Bostroem, K. A., et al. 2016, MNRAS, 458, 624, doi: 10.1093/mnras/stw273
  * CTA Consortium (2023) CTA Consortium. 2023, MNRAS, 523, 5353, doi: 10.1093/mnras/stad1576
  * de Naurois & Rolland (2009) de Naurois, M., & Rolland, L. 2009, Astropart. Phys., 32, 231, doi: 10.1016/j.astropartphys.2009.09.001
  * Deil et al. (2020) Deil, C., Donath, A., Terrier, R., et al. 2020, Gammapy, v0.18.2, Zenodo, doi: 10.5281/zenodo.4701500
  * Dembinski et al. (2020) Dembinski, H., Ongmongkolkul, P., Deil, C., et al. 2020, Zenodo, doi: 10.5281/zenodo.3949207
  * Donath et al. (2023) Donath, A., Terrier, R., Remy, Q., et al. 2023, A&A, 678, A157, doi: 10.1051/0004-6361/202346488
  * Eichler (1979) Eichler, D. 1979, ApJ, 229, 419, doi: 10.1086/156969
  * _Fermi_ -LAT Collaboration (2016) _Fermi_ -LAT Collaboration. 2016, A&A, 586, A71, doi: 10.1051/0004-6361/201526920
  * Gaensler et al. (2005) Gaensler, B. M., Haverkorn, M., Staveley-Smith, L., et al. 2005, Science, 307, 1610, doi: 10.1126/science.11088
  * Gaustad et al. (2001) Gaustad, J. E., McCullough, P. R., Rosing, W., & Van Buren, D. 2001, PASP, 113, 1326, doi: 10.1086/323969
  * Gillessen (2004) Gillessen, S. 2004, PhD thesis.  https://pure.mpg.de/pubman/faces/ViewItemOverviewPage.jsp?itemId=item_919099
  * Ginzburg & Syrovatskii (1964) Ginzburg, V. L., & Syrovatskii, S. I. 1964, The Origin of Cosmic Rays (Macmillan)
  * Gupta et al. (2020) Gupta, S., Nath, B. B., Sharma, P., & Eichler, D. 2020, MNRAS, 493, 3159, doi: 10.1093/mnras/staa286
  * Hahn et al. (2014) Hahn, J., de los Reyes, R., Bernlöhr, K., et al. 2014, Astropart. Phys., 54, 25, doi: 10.1016/j.astropartphys.2013.10.003
  * Härer et al. (2023) Härer, L. K., Reville, B., Hinton, J., Mohrmann, L., & Vieu, T. 2023, A&A, 671, A4, doi: 10.1051/0004-6361/202245444
  * Høg et al. (2000) Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27
  * Hunter (2007) Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90, doi: 10.1109/MCSE.2007.55
  * Indebetouw et al. (2013) Indebetouw, R., Brogan, C., Chen, C. H. R., et al. 2013, ApJ, 774, 73, doi: 10.1088/0004-637X/774/1/73
  * Israel et al. (2010) Israel, F. P., Wall, W. F., Raban, D., et al. 2010, A&A, 519, A67, doi: 10.1051/0004-6361/201014073
  * Johansson et al. (1998) Johansson, L. E. B., Greve, A., Booth, R. S., et al. 1998, A&A, 331, 857
  * Kafexhiu et al. (2014) Kafexhiu, E., Aharonian, F., Taylor, A. M., & Vila, G. S. 2014, Phys. Rev. D, 90, 123014, doi: 10.1103/PhysRevD.90.123014
  * Kalari et al. (2018) Kalari, V. M., Rubio, M., Elmegreen, B. G., et al. 2018, ApJ, 852, 71, doi: 10.3847/1538-4357/aa99dc
  * Kavanagh et al. (2015) Kavanagh, P. J., Sasaki, M., Bozzetto, L. M., et al. 2015, A&A, 573, A73, doi: 10.1051/0004-6361/201424354
  * Kavanagh et al. (2019) Kavanagh, P. J., Vink, J., Sasaki, M., et al. 2019, A&A, 621, A138, doi: 10.1051/0004-6361/201833659
  * Koo & McKee (1992) Koo, B.-C., & McKee, C. F. 1992, ApJ, 388, 93, doi: 10.1086/171132
  * LHAASO Collaboration (2023) LHAASO Collaboration. 2023, Phys. Rev. Lett., 131, 151001, doi: 10.1103/PhysRevLett.131.151001
  * Lopez et al. (2020) Lopez, L. A., Grefenstette, B. W., Auchettl, K., Madsen, K. K., & Castro, D. 2020, ApJ, 893, 144, doi: 10.3847/1538-4357/ab8232
  * Lortet & Testor (1984) Lortet, M. C., & Testor, G. 1984, A&A, 139, 330
  * Lucke & Hodge (1970) Lucke, P. B., & Hodge, P. W. 1970, AJ, 75, 171, doi: 10.1086/110959
  * Massey & Hunter (1998) Massey, P., & Hunter, D. A. 1998, ApJ, 493, 180, doi: 10.1086/305126
  * Mathewson et al. (1985) Mathewson, D. S., Ford, V. L., Tuohy, I., et al. 1985, ApJS, 58, 197, doi: 10.1086/191037
  * Mattox et al. (1996) Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396, doi: 10.1086/177068
  * Meixner et al. (2013) Meixner, M., Panuzzo, P., Roman-Duval, J., et al. 2013, AJ, 146, 62, doi: 10.1088/0004-6256/146/3/62
  * Mohrmann et al. (2019) Mohrmann, L., Specovius, A., Tiziani, D., et al. 2019, A&A, 632, A72, doi: 10.1051/0004-6361/201936452
  * Morlino et al. (2021) Morlino, G., Blasi, P., Peretti, E., & Cristofari, P. 2021, MNRAS, 504, 6096, doi: 10.1093/mnras/stab690
  * Ohm et al. (2009) Ohm, S., van Eldik, C., & Egberts, K. 2009, Astropart. Phys., 31, 383, doi: 10.1016/j.astropartphys.2009.04.001
  * Parsons & Hinton (2014) Parsons, R. D., & Hinton, J. A. 2014, Astropart. Phys., 56, 26, doi: 10.1016/j.astropartphys.2014.03.002
  * Particle Data Group (2022) Particle Data Group. 2022, Prog. Theor. Exp. Phys., 2022, 083C01, doi: 10.1093/ptep/ptac097
  * Pietrzyński et al. (2013) Pietrzyński, G., Graczyk, D., Gieren, W., et al. 2013, Nature, 495, 76, doi: 10.1038/nature11878
  * Price-Whelan et al. (2022) Price-Whelan, A. M., Lim, P. L., Earl, N., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Robitaille et al. (2013) Robitaille, T. P., Tollerud, E. J., Greenfield, P., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Sabbi et al. (2012) Sabbi, E., Lennon, D. J., Gieles, M., et al. 2012, ApJ, 754, L37, doi: 10.1088/2041-8205/754/2/L37
  * Sano et al. (2017) Sano, H., Yamane, Y., Voisin, F., et al. 2017, ApJ, 843, 61, doi: 10.3847/1538-4357/aa73e0
  * Sasaki et al. (2022) Sasaki, M., Knies, J., Haberl, F., et al. 2022, A&A, 661, A37, doi: 10.1051/0004-6361/202141054
  * Schneider et al. (2018) Schneider, F. R. N., Ramírez-Agudelo, O. H., Tramper, F., et al. 2018, A&A, 618, A73, doi: 10.1051/0004-6361/201833433
  * Smith & Wang (2004) Smith, D. A., & Wang, Q. D. 2004, ApJ, 611, 881, doi: 10.1086/422181
  * Testor et al. (1993) Testor, G., Schild, H., & Lortet, M. C. 1993, A&A, 280, 426
  * Tibet AS$\gamma$ Collaboration (2021) Tibet AS$\gamma$ Collaboration. 2021, Phys. Rev. Lett., 126, 141101, doi: 10.1103/PhysRevLett.126.141101
  * Townsley et al. (2006) Townsley, L. K., Broos, P. S., Feigelson, E. D., et al. 2006, AJ, 131, 2140, doi: 10.1086/500532
  * Vieu et al. (2022) Vieu, T., Gabici, S., Tatischeff, V., & Ravikularaman, S. 2022, MNRAS, 512, 1275, doi: 10.1093/mnras/stac543
  * Vieu & Reville (2023) Vieu, T., & Reville, B. 2023, MNRAS, 519, 136, doi: 10.1093/mnras/stac3469
  * Wang & Helfand (1991) Wang, Q., & Helfand, D. J. 1991, ApJ, 370, 541, doi: 10.1086/169840
  * Weaver et al. (1977) Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377, doi: 10.1086/155692
  * Webb et al. (1985) Webb, G. M., Axford, W. I., & Forman, M. A. 1985, ApJ, 298, 684, doi: 10.1086/163652
  * Wilks (1938) Wilks, S. S. 1938, Ann. Math. Stat., 9, 60, doi: 10.1214/aoms/1177732360
  * Yamaguchi et al. (2009) Yamaguchi, H., Bamba, A., & Koyama, K. 2009, PASJ, 61, S175, doi: 10.1093/pasj/61.sp1.S175
  * Yamane et al. (2021) Yamane, Y., Sano, H., Filipović, M. D., et al. 2021, ApJ, 918, 36, doi: 10.3847/1538-4357/ac0adb
  * Zabalza (2015) Zabalza, V. 2015, in Proc. 34th Int. Cosmic Ray Conf. (ICRC2015), 922, doi: 10.22323/1.236.0922