On the Potential of Bright, Young Pulsars to Power Ultrahigh Gamma-Ray Sources

On the Potential of Bright, Young Pulsars to Power Ultrahigh Gamma-Ray Sources Emma de Oña Wilhelmi 1 , Rubén López-Coto 2,3 , Elena Amato 4,5 , and Felix Aharonian 6,7 1 Deutsches Elektronen Synchrotron DESY, D-15738 Zeuthen, Germany; [email protected] 2 Istituto Nazionale di Fisica Nucleare, Sezione di Padova, I-35131 Padova, Italy; [email protected] 3 Instituto de Astrofísica de Andalucía, CSIC, E-18080 Granada, Spain 4 INAF, Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy; [email protected] 5 Dipartimento di Fisica e Astronomia, Università di Firenze, Via Sansone 1, I-50019 Sesto Fiorentino (FI), Italy 6 Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin, Ireland; [email protected] 7 Max-Planck-Institut für Kernphysik, P.O. Box 103980, D-69029 Heidelberg, Germany Received 2022 February 7; revised 2022 April 5; accepted 2022 April 12; published 2022 April 28 Abstract The recent discovery of a new population of ultrahigh-energy gamma-ray sources with spectra extending beyond 100 TeV revealed the presence of Galactic PeVatrons—cosmic-ray factories accelerating particles to PeV energies. These sources, except for the one associated with the Crab Nebula, are not yet identified. With an extension of 1° or more, most of them contain several potential counterparts, including supernova remnants, young stellar clusters, and pulsar wind nebulae (PWNe), which can perform as PeVatrons and thus power the surrounding diffuse ultrahigh-energy gamma-ray structures. In the case of PWNe, gamma-rays are produced by electrons, accelerated at the pulsar wind termination shock, through the inverse Compton scattering of 2.7 K cosmic microwave background (CMB)radiation. The high conversion efficiency of pulsar rotational power to relativistic electrons, combined with the short cooling timescales, allow gamma-ray luminosities up to the level of  ~ g L E0.1 . The pulsar spin-down luminosity,  E , also determines the absolute maximum energy of individual photons:  » g E E0.9 Pe V ,max 36 0.65 . This fundamental constraint dominates over the condition set by synchrotron energy losses of electrons for young PWNe with typical magnetic field of ≈100 μG with  - E 10 erg s 37 1. We discuss the implications of E γ,max by comparing it with the highest-energy photons reported by LHAASO from a dozen of ultrahigh-energy sources. Whenever a PWN origin of the emission is possible, we use the LHAASO measurements to set upper limits on the nebular magnetic field. Unified Astronomy Thesaurus concepts: Gamma-ray sources (633);Pulsars (1306);Rotation powered pulsars (1408);Pulsar wind nebulae (2215);High energy astrophysics (739) 1. Introduction The recent discovery of the LHAASO Collaboration (Cao et al. 2021), reporting the detection of a dozen sources with particle spectra reaching PeV (1 PeV =10 15 eV)energies, represents a major step toward the identification of the nature of the sources known as PeV accelerators, or PeVatrons. These observations are complemented in the gamma-ray regime by those from instruments sensitive in the 100 GeV–100 TeV energy range, in particular by water Cherenkov instruments like HAWC and Tibet ASγ(Albert et al. 2020; Amenomori et al. 2021), and Imaging Atmospheric Cherenkov Telescope (IACT) arrays such as H.E.S.S., MAGIC, and VERITAS (Aleksićet al. 2015; Park 2015; H.E.S.S. Collaboration et al. 2018).In general, IACTs provide a superior angular resolution that can be used to localize the emission regions more accurately and identify the accelerator type. The majority of the sources reported by LHAASO in Cao et al. (2021)are described by an extended gamma-ray emission as large as ∼1 o . These sources have, in almost all cases, a sub-100 TeV counterpart. Interestingly, in this sub-100 TeV regime, more than 30% of the sources detected in the Galactic plane have been associated with pulsar wind nebulae (PWNe; Wakely & Horan 2008; H.E. S.S. Collaboration et al. 2018). These associations are based on spatial correlation with energetic pulsars, and spectral– morphological features connecting the usually extended (0°.2)gamma-ray emission with the pulsars. Pulsars (or PWNe)are the only identified source class in which PeV particles have been detected: the Crab Nebula, associated with the young, very energetic pulsar PSR B0531+21, shows a synchrotron (steady and flaring)spectrum in the GeV regime that corresponds to PeV electrons (Abdo et al. 2010,2011), and it has also been recently detected up to an energy of 1.1 PeV by the LHAASO experiment (LHAASO Collaboration et al. 2021). Nevertheless, pulsars appear to be close to the absolute theoretical limit in terms of acceleration rate (Aharonian 1995; LHAASO Collaboration et al. 2021). In the following, we investigate the capability of energetic pulsars to power PeVatrons. 2. Pulsars as Effective PeV Accelerators PWNe have been recognized as one of the most efficient electron factories in our galaxy (H.E.S.S. Collaboration et al. 2018). They are powered by energetic pulsars, which inject ultrarelativistic electrons and positrons in their magnetosphere. These particles form a cold ultrarelativistic wind, expanding with bulk Lorentz factor Γin the range 10 4 –10 7 , until reaching the termination shock (TS; Rees & Gunn 1974; Kennel & Coroniti 1984). At the shock, particles are believed to be accelerated to multi-TeV energies, inflating a nonthermal nebula that constitutes the plerion (Amato 2020). A large fraction of the pulsar spin-down power,  E , is radiated in the The Astrophysical Journal Letters, 930:L2 (6pp), 2022 May 1 https://doi.org/10.3847/2041-8213/ac66cf © 2022. The Author(s). Published by the American Astronomical Society. Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s)and the title of the work, journal citation and DOI. 1 very high energy regime via inverse Compton (IC)scattering, resulting in a power-law spectrum that can extend up to at least a few tens of TeV (Aharonian et al. 1997; Gaensler & Slane 2006; de Jager & Djannati-Ataï 2009; Amato & Olmi 2021). At these energies, the scattering occurs mostly in a deep Klein–Nishina regime, where electrons lose most of their energy in a single scattering event and the maximum energy observed in photons roughly coincides with the maximum energy to which the electrons are accelerated. Radiation losses beyond hundreds of TeV, into the PeV regime, are very rapid and demand an extremely efficient acceleration rate. Indeed the connection between pulsars and PeVatrons provides important constraints, which stand on first principles, independently of more sophisticated modeling. First, the absolute maximum energy the particles can reach, E max , depends ultimately on the maximum potential drop between the pulsar and infinity, ()  F=Ec PSR 12 , with cthe speed of light. Since the particle acceleration, regardless of the acceleration mechanism, is always carried out by the electric field |E|, this maximum energy is related to the electric potential associated with this field. Therefore, E max can be defined in terms of the maximum potential drop and has as an absolute maximum at the value ()  = E qE c max 12 , where qis the charge of the particle. This is equivalent to saying that the maximum energy of the particle depends on the size of the accelerator, which in the case of PWNe is the size of the TS, as deduced from X-ray observations of young, well-studied systems; thus, ∣∣=E E q max R TS . If the wind is described as an ideal magnetohydrodynamic flow (Kennel & Coroniti 1984; Porth et al. 2014; Olmi & Bucciantini 2019), the electric field strength cannot exceed the magnetic one, B TS .Defining the ratio between the two as η e (which is η e 1 in such ideal conditions), the former equation can be written as E max =qη e B TS R TS , which is the well-known Hillas criterion (Hillas 1984). The magnetic energy density (defined as B TS 2/8π)can be expressed as a fraction η B of the pulsar wind energy flux: () ()  php = BE Rc 84;1 TS 2 B TS 2 thus, () ( ) ()   h hm = = - - BREc RE 2 25 G, 2 TS B12 TS 112 B 12 0.1 1 36 12 where R 0.1 is the termination shock radius in units of 0.1pc and  E 36 is in units of 10 36 erg s −1 . Using this expression in the expression of E max described above results in the following limit to the maximum energy of the accelerated particles, independently of whether they are electrons, positrons, or protons: ()  hh»EE2PeV. 3 e max B 12 36 12 The fraction of pulsar wind energy flux transferred to magnetic field is constrained, by energy conservation, to be η B 1. The expression above, which is independent of the particle species, can be used to derive an absolute maximum to the energy to which particles can be accelerated. For electrons in the multi-TeV regime, the most relevant target for IC scattering is the 2.7 K CMB photons. The photon and electron energy can then be related using the approximation presented by LHAASO Collaboration et al. (2021): g E E2.15 e,15 0.77 PeV (with E γ,15 in units of 10 15 eV, or PeV), which provides an accuracy better than 10% above 30 TeV. Thus, Equation (3)can be written as ()  hh» g EE0.9 PeV. 4 eB ,max 1.3 0.65 36 0.65 The above expression provides a direct link between the observed maximum energy in the gamma-ray spectrum and the spin-down power of the pulsar. It can be immediately derived that only very energetic pulsars with at least ⪆  E 103 6 erg s −1 could power the observed PeV gamma-rays. Additionally, a second condition is required to shine in gamma-rays up to PeVs: the acceleration rate τ acc should also overcome the radiative losses of the parent electrons. The first can be expressed as a function of the magnetic field at the TS as τ acc =E e /(η e eB TS c). The overall cooling time, which can be written as ()ttt=+ - 11 loss sync IC 1, is dominated, even for a few μG magnetic field, by synchrotron losses, above a few hundreds of TeV. In the Klein–Nishina limit, the IC cooling time of electrons in 2.7 K CMB depends only on the electron energy as tE10 e IC 12 ,15 0.7 s(Khangulyan et al. 2014). The former can be compared with the synchrotron time t´- - - EB410 e sync 9,15 15 2s, where B −5 is the magnetic field in units of 10 μG. The condition τ acc =τ syn results in the following expression for the maximum energy of the electron population: () // h»- - EB20 PeV. 5 e e,max 12 5 12 This corresponds to a maximum energy of the photons produced by IC scattering on the CMB of ()h» g- - EB5PeV, 6 e ,max 0.65 5 0.65 or, in terms of pulsar  E and TS radius, using Equation (2): ()  hh» g-- ERE2.7 PeV. 7 e ,max 0.65 B0.33 0.1 0.65 36 0.33 The comparison between Equations (3)and (5)shows that radiation losses pose the most serious challenge to reach PeV energies only for young, energetic pulsars (  E 1037 erg s −1 ), with a magnetic field at the wind TS in the 100 μG range, such as the Crab Nebula. 8 This does not mean that the Crab is a poor accelerator. On the contrary, as discussed in LHAASO Collaboration et al. (2021), the Crab Nebula, with an estimated magnetic field of ≈112 μG, requires an acceleration rate corresponding to η e ≈0.16 to reach the energies observed (E γ,max =1.1 PeV). These large values of η e are impressive, being several orders of magnitude larger than those inferred for other powerful accelerators, such as, e.g., supernova remnants (SNRs), for which typically η e ≈10 −3 (Malkov & Drury 2001). Even though larger values of η e 1, are possible during, e.g., gamma-ray flares (Bühler & Blandford 2014), the large value derived in Crab only reflects the relativistic nature of the acceleration, which in such shocks is far from being understood (see, e.g., Amato 2020; Amato & Olmi 2021). The constraints so far discussed are related to the maximum photon energy observed. Additional information can be obtained by comparing the IC luminosity of the PeV sources and the total power injected by the pulsar in its surrounding. 8 Note that in scenarios where acceleration and radiation do not occur in the same region, the absolute maximum energy is still determined by the potential drop. 2 The Astrophysical Journal Letters, 930:L2 (6pp), 2022 May 1 de Oña Wilhelmi et al. The total energy in electrons responsible for IC radiation can be derived using the gamma-ray observations as W e,γ =L γ τ IC . W e,γ cannot exceed W e,PSR , the total energy made available by the pulsar in the form of gamma-ray-emitting electrons. Since the lifetimes of the latter are determined by losses, one can write  gt= W E e,PSR eff loss, with γ eff the fraction of  E converted into gamma-ray-emitting electrons: ⎜⎟ ⎛ ⎝ ⎞ ⎠ ()()   gt t =+ =+ g g - - L E L EEB 1 10 1 260 , 8 e eff IC syn 4,32 36 ,15 1.7 5 2 which also depends strongly on the magnetic field. The maximum allowed conversion efficiency from rotational power to a gamma-ray one (γ eff =1)results in another boundary to pulsars as ultrahigh-energy sources. 3. Comparison with Pulsars in the Region of Interest Twelve ultrahigh-energy gamma-ray sources were reported by Cao et al. (2021), with a spectral energy distribution extending up to more than 100 TeV, one of them being associated with the Crab Nebula. The maximum energy in photons E γ,LHAASO and electrons E e LHAASO , derived from the LHAASO observations, are listed in Table 1. The latest was derived from the photon energy using the expression in Khangulyan et al. (2014). The majority of these sources show a diffuse gamma-ray structure, with angular extensions up to 1 o . This extended structure makes the association with the PeVatron accelerator complex. To explore the possibility of an association of the ultrahigh-energy sources with pulsars, we searched for relatively young (τ<10 6 yr), energetic ( E /d >10 kpc 234 erg s −1 kpc 2 ,or > E 1036 erg s −1 when the distance is unknown)pulsars in the ATNF catalog, 9 located within 1 o around the position of the LHAASO sources. For each of the LHAASO sources, we found at least one pulsar (two in some cases)that could potentially be linked to it, except for the source J2108+5157, for which no bright pulsar is found in the vicinity. The selected pulsars and their properties are listed in Table 1. The LHAASO source associated with the Crab Nebula has been described in detail in LHAASO Collaboration et al. (2021). To evaluate the potential of an associated pulsar to power the LHAASO sources, we estimate the maximum energy to which particles can be accelerated by such a pulsar, assuming η B =η e =1 in Equation (3). The results are summarized in Figure 1, and the maximum energies of the particles and corresponding gamma-ray E γ,max , in the case of IC scattering of the CMB, are listed in Table 2. In Figure 1we place the potential associations in the  E – E e,max plane and compare them with theoretical predictions based on Equation (3)for h h eB 12 ranging from 0.01 to 1. Out of all pulsars possibly associated with the 12 LHAASO sources, only in the Crab pulsar is the maximum energy limited by the radiation losses, while for all other pulsars the most relevant constraint will come from saturation of the full available potential drop (Equation (3)). The upper limit to the maximum electron (and photon)energy, using η e =1 and 100μG, is marked with a dotted horizontal line in Figure 1. Above the red line, the particle flow would require values of η B =1 and η e >1, and would demand nonideal mechanisms (see Amato & Olmi 2021 for a review). Only one pair of an ultrahigh-energy source and energetic pulsar (LHAASO J2032+4102 /J2032+4127)lies above the absolute maximum, resulting in an impossible connection between the two (if the spin-down power of the pulsar is correct within a factor of ∼4). We also marked in Figure 1two remarkable pulsars with vertical blue lines: Geminga and the Crab twin, N157B, located in the Magellanic Cloud. With a moderated spin-down luminosity of 3.26 ×10 34 erg s −1 , but located at a small distance of 250 pc, the Geminga Nebula is a prime target for LHAASO, given its large size of ∼2° (Abeysekara et al. 2017). N157B is, on the contrary, the furthest gamma-ray PWN detected (H.E.S.S. Collaboration Table 1 LHAASO Ultrahigh-energy Sources, Together with the Bright, Young Pulsars Located within 1 o of the LHAASO Source and Their Characteristics LHAASO Source Pulsar Separation  E Age Distance Flux 100TeV E γLHAASO E e LHAASO (deg)×10 36 (erg s −1 )(kyr)(kpc)(c.u.)(PeV)(PeV) J1825-1326 J1826-1256 0.51 3.6 14.4 1.55 3.57 0.42 1.06 B1823-13 0.16 2.8 21.4 3.61 3.57 0.42 1.06 J1839-0545 J1837-0604 0.61 2.0 33.8 4.77 0.70 0.21 0.63 J1838-0537 0.25 6.0 4.9 L0.70 0.21 0.63 J1843-0338 J1841-0345 0.37 0.3 55.9 3.78 0.73 0.26 0.74 J1844-0346 0.41 4.2 11.6 L0.73 0.26 0.74 J1849-0003 J1849-0001 0.10 9.8 43.1 L0.74 0.35 0.92 J1908+0621 J1907+0602 0.32 2.8 19.5 2.37 1.36 0.44 1.10 J1907+0631 0.33 0.5 11.3 3.40 1.36 0.44 1.10 J1929+1745 J1925+1720 0.94 0.9 115.0 5.06 0.38 0.71 1.60 J1928+1746 0.07 1.6 82.6 4.34 0.38 0.71 1.60 J1956+2845 J1954+2836 0.44 1.0 69.4 1.96 0.41 0.42 1.06 J1958+2846 0.54 0.3 21.7 1.95 0.41 0.42 1.06 J2018+3651 J2021+3651 0.42 3.4 17.2 1.80 0.50 0.27 0.75 J2032+4102 J2032+4127 0.41 0.1 201.0 1.33 0.54 1.42 2.79 J2108+5157 J2226+6057 J2229+6114 0.38 22.0 10.5 3.00 1.05 0.57 1.35 Note. The two rightmost columns display the maximum energy quoted by Cao et al. (2021)and its corresponding energy in electrons, using the formulation in Khangulyan et al. (2014). 9 http://www.atnf.csiro.au/research/pulsar/psrcat/ 3 The Astrophysical Journal Letters, 930:L2 (6pp), 2022 May 1 de Oña Wilhelmi et al. et al. 2012), but its large spin-down luminosity (=´ E 4.9 1038 erg s −1 )and similarity to the Crab Nebula makes it also an interesting source to understand the contribution of pulsars to the PeV sky. We can also use the spectral parameters of the PeV sources, in particular E γ,max and gamma-ray luminosity to impose an upper limit on the magnetic field. Constraints are provided by the fact that synchrotron losses should not forbid acceleration up to E e,max (Equation (5)) and the energy input from the pulsar should be sufficient to power the gamma-ray source (Equation (8)). We found that, in general, the latter constraint is much stronger and requires that the magnetic field cannot exceed a few tens of μGauss, which agrees with the typical values derived from very high energy observations in the TeV regime (see, e.g., H.E.S.S. Collaboration et al. 2018). Despite these low constraints in the magnetic field, the Larmour radius of the electrons with the highest energies is still in agreement with the typical size of the TS, defined by the balance between the wind pressure and the one from the surrounding medium (see, e.g., Kargaltsev & Pavlov 2008). Figure 1. Maximum electron energy derived from the LHAASO spectra vs. spin-down power of the co-located pulsars. The right y-axis shows the corresponding gamma-ray energy. The colored area shows the values for η e η B 1/2 ranging from 0.01 to 1, with the red line indicating the limiting value corresponding to maximally efficient acceleration η e =1 and η B =1. The dotted black line marks the upper limit to the maximum energy for young pulsars with a large magnetic field of 100 μG. The blue dashed horizontal lines show the predicted values for PWNe associated with Geminga and N157B. Table 2 LHAASO Ultrahigh-energy Sources and Putative Associated Pulsars, with the Corresponding Constraints on the Maximum Energy and Magnetic Field LHAASO Source Pulsar E γ,max E max B max (PeV)(PeV)(μG) J1825-1326 J1826-1256 2.06 3.79 38 B1823-13 1.77 3.35 14 J1839-0545 J1837-0604 1.44 2.83 33 J1838-0537 2.78 4.90 ?100 J1843-0338 J1841-0345 0.41 1.04 12 J1844-0346 2.25 4.10 ?100 J1849-0003 J1849-0001 3.71 6.26 ?100 J1908+0621 J1907+0602 1.77 3.35 30 J1907+0631 0.63 1.46 9 J1929+1745 J1925+1720 0.91 1.95 9 J1928+1746 1.26 2.53 14 J1956+2845 J1954+2836 0.94 2.00 37 J1958+2846 0.47 1.17 22 J2018+3651 J2021+3651 1.99 3.69 102 J2032+4102 J2032+4127 0.28 0.77 7 J2108+5157 J2226+6057 J2229+6114 5.89 9.38 64 4 The Astrophysical Journal Letters, 930:L2 (6pp), 2022 May 1 de Oña Wilhelmi et al. 4. Concluding Remarks We derived the absolute maximum energy that can be accelerated by pulsars, obtained from the maximum potential drop available, without further assumptions beyond ideal MHD flow. This maximum energy can now be confronted with observational results as those recently published by LHAASO. The extreme energies reached in the sources detected by LHAASO provide direct information about current particle acceleration, given the fast cooling time involved beyond hundreds of TeV. At these energies, the upscattering of the 2.7 K CMB radiation dominates the observed gamma-ray radiation in PWNe, providing a powerful diagnostic tool. Additionally, these multi-TeV electrons propagating in the magnetized nebulae should also power an X-ray nebula, visible at a few keV. The detection of such an extended nebula is challenging for pointing X-ray instruments like XMMNewton or Chandra (Kargaltsev & Pavlov 2008; Liu et al. 2019). However, the new X-ray satellite eROSITA, sensitive to X-rays in the energy range of 0.3–11 keV and with a wide field of view of 0.81 o , is optimal to constrain the X-ray counterpart. The expected sensitivities achieved by eROSITA for extended sources in the energy range of 0.5–2 keV are 1.1 ×10 −13 erg cm −2 s −1 for the first all-sky survey (eRASS:1), and 3.4 ×10 −14 erg cm −2 s −1 for the 4 yr all-sky survey (eRASS:8). With fluxes at 100 TeV ranging from ∼5×10 −13 erg cm -2 s −1 to ∼5×10 −12 erg cm -2 s −1 ,the X-ray counterpart, assuming a magnetic field as low as 3 μG, should be expected with fluxes above ∼2×10 −12 erg cm -2 s −1 in the 0.5–2 keV band, an order of magnitude larger than the eROSITA sensitivity. These numbers should be taken with caution, since the surface brightness might not be homogeneous across the large TeV source region. From the 11 sources considered, two sources stand out, and different accelerators and/or gamma-ray production mechanisms should be investigated: LHAASO J2108+5157 and LHAASO J2032+4102. The first is extensively discussed in The LHAASO Collaboration (2021)and it is found to be pointlike, within the angular resolution of LHAASO for this analysis (0°.26).The closest pulsar is ∼3°away, which at 2–3 kpc corresponds to more than 100 pc away. No counterpart has been found in the TeV regime either. The second source is co-located with a pulsar (PSR J2032+4127)in an interacting binary system (Abeysekara et al. 2018; Williamson 2019). The above considerations only apply to the isolated pulsar; however, the mixing of the two winds could in principle lead to different conclusions. The pulsar powers a compact ∼0°.2 gamma-ray nebula (Aharonian et al. 2002; Konopelko et al. 2007;Albertetal.2008;Aliuetal. 2014), which can only partially be connected to the very extended ultrahigh-energy source. The system is located at the heart of the Cygnus cocoon, a bright GeV and TeV extended diffuse emission, which has also been connected with several individual sources, including the massive stellar cluster Cygnus OB (Ackermann et al. 2011; Aharonian et al. 2019; Abeysekara et al. 2021). The potential connection between the Cygnus cocoon and the LHAASO source opens interesting prospects for stellar clusters as contributors of ultrahigh-energy particles (Bykov et al. 2020). Further information regarding the morphology of these sources should provide crucial insight into the origin of the emission. Indeed, electrons and positrons at these energies undergo fast losses due to synchrotron radiation, and might appear as compact, subdegree regions. However, for low enough (a few ∼μG)magnetic fields, and fast enough transport, these electrons could still fill up a volume larger than a few tens of parsecs, which would match the large extension observed, if located close enough to us (a few kiloparsecs). Alternatively, these electrons might have escaped into the interstellar medium, filling up a halo where particles are essentially free from their parent PWN (López-Coto et al. 2022). This is particularly relevant for nebulae like HESS J1825-137, for which a clear energy-dependent morphology has been established in the TeV regime (H.E.S.S. Collaboration et al., 2019). Observations above 100~TeV should provide a clear picture of the radiative cooling and propagation of electrons. Note that regardless of the propagation regime, the maximum energy to which a particle can be accelerated is always limited by the equations derived here. Gamma-ray images at different energies should also serve as a test-bench for other effects involving ballistic and diffusive propagation, which could play an important role in the observed morphology (Prosekin et al. 2015). EdOW acknowledges the support of the Alexander von Humboldt foundation and DESY (Zeuthen), a member of the Helmholtz Association HGF. R.L.-C. acknowledges the financial support of the European Union Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement 754496-FELLINI. R.L.-C. also acknowledges financial support from the State Agency for Research of the Spanish MCIU through the Centre of Excellence Severo Ochoa award to the Instituto de Astrofísica de Andalucía (SEV-20170709). E.A. acknowledges support from ASI-INAF under grant 2017-14-H.0, and from INAF under grants INAF Mainstream 2018 and PRIN-INAF 2019. ORCID iDs Emma de Oña Wilhelmi https://orcid.org/0000-00025401-0744 Rubén López-Coto https://orcid.org/0000-0002-3882-9477 Elena Amato https://orcid.org/0000-0002-9881-8112 Felix Aharonian https://orcid.org/0000-0003-1157-3915 References Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ,708, 1254 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, Sci,331, 739 Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017, Sci,358, 911 Abeysekara, A. U., Benbow, W., Bird, R., et al. 2018, ApJL,867, L19 Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2021, NatAs,5, 465 Ackermann, M., Ajello, M., Allafort, A., et al. 2011, Sci,334, 1103 Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, NatAs,3, 561 Aharonian, F., Akhperjanian, A., Beilicke, M., et al. 2002, A&A,393, L37 Aharonian, F. A. 1995, NuPhS,39, 193 Aharonian, F. A., Atoyan, A. M., & Kifune, T. 1997, MNRAS,291, 162 Albert, A., Alfaro, R., Alvarez, C., et al. 2020, ApJ,905, 76 Albert, J., Aliu, E., Anderhub, H., et al. 2008, ApJL,675, L25 Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2015, JHEAp,5, 30 Aliu, E., Aune, T., Behera, B., et al. 2014, ApJ,783, 16 Amato, E. 2020, in Proc. of Science, High Energy Phenomena in Relativistic Outflows VII (HEPRO VII)- Relativistic outflows from galactic sources 345 (Trieste: SISSA),033 Amato, E., & Olmi, B. 2021, Univ,7, 448 Amenomori, M., Bao, Y. W., Bi, X. J., et al. 2021, PhRvL,126, 141101 Bühler, R., & Blandford, R. 2014, RPPh,77, 066901 Bykov, A. M., Marcowith, A., Amato, E., et al. 2020, SSRv,216, 42 Cao, Z., Aharonian, F. A., An, Q., et al. 2021, Natur,594, 33 de Jager, O. C., & Djannati-Ataï, O. C. 2009, in Astrophysics and Space Science Library, ed. W. Becker, 357 (Berlin: Springer), 451 Gaensler, B. M., & Slane, P. O. 2006, ARA&A,44, 17 H.E.S.S. Collaboration, Abramowski, Acero, A., et al. 2012, A&A, 545, L2 H.E.S.S. Collaboration, Abdalla, Abramowski, H., et al. 2018, A&A, 612, A1 5 The Astrophysical Journal Letters, 930:L2 (6pp), 2022 May 1 de Oña Wilhelmi et al. H. E. S. S. Collaboration, Abdalla, H., Aharonian, F., et al. 2019, A&A, 621, 116 Hillas, A. M. 1984, ARA&A,22, 425 Kargaltsev, O., & Pavlov, G. G. 2008, in AIP Conf. Ser. 983, 40 Yr of Pulsars: Millisecond Pulsars, Magnetars and More, C. Bassa, ed. C. Bassa et al. (San Francisco, CA: AIP),40 Kennel, C. F., & Coroniti, F. V. 1984, ApJ,283, 710 Khangulyan, D., Aharonian, F. A., & Kelner, S. R. 2014, ApJ,783, 100 Konopelko, A., Atkins, R. W., Blaylock, G., et al. 2007, ApJ,658, 1062 LHAASO Collaboration, Cao, Aharonian, Z., et al. 2021, Sci,373, 425 Liu, R.-Y., Ge, C., Sun, X.-N., & Wang, X.-Y. 2019, ApJ,875, 149 López-Coto, R., de Oña Wilhelmi, E., Aharonian, F., Amato, E., & Hinton, J. 2022, NatAs,6, 199 Malkov, M. A., & Drury, L. O. 2001, RPPh,64, 429 Olmi, B., & Bucciantini, N. 2019, MNRAS,488, 5690 Park, N. & VERITAS Collaboration 2015, in Proc. of ICRC, Performance of the VERITAS experiment 34 (Trieste: SISSA) Porth, O., Komissarov, S. S., & Keppens, R. 2014, MNRAS,438, 278 Prosekins, A. Y., Kelner, S. R., & Aharonian, F. A. 2015, PhRvD,92, 083003 Rees, M. J., & Gunn, J. E. 1974, MNRAS,167, 1 Wakely, S. P., & Horan, D. 2008, Proc. of ICRC, TeVCat: An online catalog for Very High Energy Gamma-Ray Astronomy, 30 ed. R. Caballero et al. (Mexico City: Universidad Nacional Autónoma de México),1341 Williamson, T. 2019, arXiv:1908.04165 Zhen, C., Aharonian, F., An, Q., et al. 2021, ApJ,919, 22 6 The Astrophysical Journal Letters, 930:L2 (6pp), 2022 May 1 de Oña Wilhelmi et al.