Open Access
Issue
A&A
Volume 660, April 2022
Article Number A8
Number of page(s) 6
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202142097
Published online 30 March 2022

© M. Breuhaus et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Open Access funding provided by Max Planck Society.

1. Introduction

The High-Altitude Water Cherenkov Gamma-Ray Observatory (HAWC), the Tibet ASγ, and the Large High Altitude Air Shower Observatory (LHAASO) collaborations have reported the detection of multiple Galactic ultra-high energy (UHE; Eγ ≥ 100 TeV) gamma-ray sources (Abeysekara et al. 2020; Amenomori et al. 2019, 2021; Cao et al. 2021a). The latest published data from the partially completed LHAASO observatory (Cao et al. 2021a) reveal 12 sources with significant detections (> 7σ) above several hundred TeV. The completed LHAASO facility will greatly enhance our capability to observe the Galaxy’s most extreme particle accelerators. The unprecedented sensitivity in the UHE regime provided by LHAASO and other planned facilities, such as the Cherenkov Telescope Array (CTA; Cherenkov & Acharya 2019) and the Southern Wide-field Gamma-ray Observatory (SWGO; Albert et al. 2019), is particularly important for efforts to identify Galactic sources of cosmic rays (CRs) with energies higher than 1 PeV. Acceleration to the knee energy in the CR spectrum is a long-standing issue for the supernova-remnant (SNR) CR origin theory (Lagage & Cesarsky 1983; Bell 2013), and though UHE γ rays provide the best method for identifying PeV CR accelerators in our Galaxy, a critical examination of all such sources is warranted.

Hard γ-ray spectra at these energies are generally perceived to favour hadronic emission mechanisms since the Klein-Nishina suppression of inverse Compton (IC) emission at UHEs can, in many circumstances, disfavour hard equilibrium spectra. However, this may not apply if the emission region is moderately photon dominated (e.g. Blumenthal & Gould 1970; Agaronyan & Ambartsumyan 1985; Zdziarski & Krolik 1993). Recently, Breuhaus et al. (2021a; henceforth B21) pointed out that the emission from high-power pulsars may naturally account for hard spectrum UHE sources (see also Di Mauro et al. 2021; Sudoh et al. 2021). Observations reported by the HAWC collaboration already indicate that UHE γ-ray emission in the vicinity of pulsars with high spin-down powers of Ė > 1036 erg s−1 may be common (Albert et al. 2021a). As most of the UHE sources reported to date are spatially coincident or in close proximity to powerful pulsars, disentangling the contribution from competing sources is a necessary step in current and future investigations.

Of the 12 LHAASO-detected sources, Cao et al. (2021a) reported the spectral energy density distributions for a selection of three: LHAASO J2226+6057, LHAASO J1825−1326, and LHAASO J1908+0621. For the last source, a phenomenological fit comparing a hadronic and a leptonic model was presented. For the leptonic fit, the required power-law index of the electron injection spectrum was −1.75, with a cutoff energy of 800 TeV. Such hard injection spectra are not generally expected, at least not within a relativistic shock acceleration scenario. At first sight, this might further favour a hadronic model; however, it was assumed that the magnetic energy density in the emission zone exceeded that of the photon field in their model. The other two sources are already listed in the HAWC 100 TeV source catalogue (Abeysekara et al. 2020), and consistent matches to the UHE emission using standard theoretical injection spectra were presented in B21. It is thus important to explore if these results still apply in light of the improved data in the UHE regime, as well as if such a scenario might also apply for new sources. It is of course possible that the situation is more complex, with more than one source powering the detected emission (Albert et al. 2021b; Crestan et al. 2021). Naively, one might expect that several nearby or overlapping sources accelerating particles to PeV energies is improbable; future high-resolution instruments may be necessary to unambiguously determine if this is indeed the case.

In this paper we revisit the predictions of B21 in light of the latest LHAASO results. We focus on the three sources for which spectral energy distributions (SEDs) are available, two of which overlap with the HAWC sources previously considered in B21. In Sect. 2 we describe the methodology, and Sect. 3 presents the results. We discuss our findings in Sect. 4.

2. Method

To calculate the emission, we used a simple single-zone treatment. In this model, accelerated electrons are injected over a time frame that corresponds to the characteristic age of the associated pulsar into a region with a constant homogeneous magnetic field and isotropic radiation fields. The injection rate is assumed to follow a time-independent (isotropic) power-law distribution with an exponential cutoff:

(1)

In other words, we ignore the evolution of the source and restrict our attention to the highest energy electrons, for which the cooling time is less than or close to the dynamical age. We note that the impact of considering a finite source age rather than an equilibrium is modest in most cases, but it improves the fit and modifies the best-fit injection index in the case of younger sources and wide spectral coverage.

The total electron spectrum is calculated by taking energy losses due to both IC and synchrotron emission into account. On timescales shorter than the dynamical time, adiabatic losses are assumed to be negligible. The solution can be obtained using a standard single-zone approach (see e.g. Atoyan & Aharonian 1999), which has been implemented into the open-source code GAMERA (Hahn 2015). Particle escape is neglected, a valid assumption if the particles diffuse sufficiently slowly and cool fast. The energy is therefore deposited close to the source, and the spectral data points obtained cover the entire emission region. We can check a posteriori that these hold.

For the radiation fields, we used the large-scale Galactic emission model of Popescu et al. (2017) together with the cosmic microwave background (CMB) at the respective source locations, assuming the sources are located at the position of the pulsar counterparts at their nominal distances. We retained the option to include a local enhancement factor of the diffuse photon field due to nearby photon sources or the location in spiral arms (see B21). Between source and detection at Earth, γ rays are absorbed due to pair production with mostly far-infrared (FIR) Galactic radiation fields and the CMB, the latter being important for energies above ∼100 TeV (e.g. Gould & Rephaeli 1978; Vernetto & Lipari 2016). In the worst case, for the source LHAASO J1908+0621 associated with the pulsar PSR 1907+0632, the transmissivity at 1 PeV is 67%.

Cao et al. (2021a) provide a list of potential particle acceleration sites, including SNRs and pulsars, associated with the detected UHE sources. We focus here on the possible contribution from pulsars alone. For LHAASO J2226+6057 there is only one known high-power pulsar in the region of interest, while for the other sources two powerful pulsars are located within 1° of the measured positions. The key physical parameters of these pulsars are provided in Table 1, and the characteristic spin-down age is taken as a proxy for the true age of the pulsar. In each case we then fitted the normalisation, the injection index, α, and the cutoff energy, Ecut (within physically motivated bounds) to the data for a fixed magnetic field for each pulsar association. The errors in determining Ecut are found to be highly asymmetric. We therefore estimated the 95% confidence interval separately using the probability density function of the χ2 distribution and its dependence on Ecut.

Table 1.

Adopted properties of the pulsars with possible association with the three LHAASO sources (Cao et al. 2021a, and references therein).

3. Results

We applied the method described above to the sources LHAASO J2226+6057, LHAASO J1825-1326, and LHAASO J1908+0621, for which the spectra are provided in Cao et al. (2021a).

3.1. LHAASO J2226+6057

This source can be associated with the SNR G106.3+2.7 as well as the pulsar J2229+6114 and its nebula (the ‘Boomerang Nebula’; e.g. Kothes et al. 2006). The centroid of the 100 TeV γ-ray source was also shown by the Tibet ASγ Collaboration (2021) to be spatially coincident with a molecular cloud, which might favour a hadronic origin. However, given the mature age of the SNR (≲104 years), rather unique conditions are necessary to allow ongoing PeV acceleration at the SNR shock. We consider here an alternative scenario where the electrons and positrons are accelerated as part of the energy dissipation process of the pulsar wind, for example via shock acceleration (e.g. Giacinti & Kirk 2018).

Using the pulsar parameters and photon field associated with the pulsar’s position (see Table 1), fits based on the LHAASO data alone are not well constrained. While theoretical predictions for relativistic shock acceleration in general disfavour injection power-law indices α < 2 (e.g. Sironi et al. 2015), equally good fits are possible over a broad range of α and cutoff energies. Figure 1 compares the LHAASO data with an exemplary model fit for B = 3 μG and a fixed injection index of α = 2.2. The resulting normalisation parameter is 0 = 9 ± 2 × 1033 erg−1 s−1, and the cutoff energy is at a 95% confidence level. The cutoff energy is comfortably below the maximum potential drop of the pulsar (see for example Amato & Olmi 2021 or B21 Appendix A). The resulting flux corresponds to a fraction η = 0.13% of the pulsar’s spin-down power injected into electrons above 1 TeV, or ≈1% for electrons above 1 GeV assuming no break in the injected spectrum. PSR J2229+6114 is therefore easily able to provide the necessary power input. Multi-wavelength models of the source have been explored by other authors, Liu et al. (2020) and Yu et al. (2022), where their best fits correspond to a 4 μG magnetic field, with α = 2.4 − 2.5. It should be noted that in both studies an infrared (IR) photon field of 0.3 eV cm−3 was adopted, which differs from the model of Popescu et al. (2017) used here (see Table 1).

thumbnail Fig. 1.

Data and models for LHAASO J2226+6057. The orange curve is a model fit for a magnetic field of B = 3 μG and a fixed power-law index of the injection spectrum of α = 2.2, leaving 0 and Ecut as free parameters. The resulting values were 0 = 9 × 1033 ergxs s−1 and Ecut = 420 TeV. The dotted line shows the emission without absorption.

An injection spectrum of α = 2.2 is also consistent with the radio spectrum, Sν ∝ ν−0.59, measured from the Boomerang Nebula. Kothes et al. (2006) interpreted the break at ν ≈ 5 GHz as an indication that this is in fact the cooled spectrum, though this would require a magnetic field strength of 2.6 mG. Recent work has questioned whether such intense fields can exist on large scales (Liu et al. 2020). However, Kothes et al. (2006) themselves interpret the radio emission as originating in the ‘crushed’, part of the nebula due to a reverse shock resulting from the SNR interaction with dense H I material in the north-east direction. In this scenario, the termination shock in the south-west region would provide more favourable conditions not only for acceleration but also for escape into the low-density cavity evacuated by the SNR. This would also be consistent with the observed offset between the emission centres of the Boomerang Nebula and the LHAASO emission, which are separated by ≈0.4° to the south-west. While a detailed numerical model is needed to explore the complex transport expected in such conditions, it is possible to make rough estimates of the expected diffusion coefficient to account for the radial extent of the LHAASO emission. We define a radial diffusion coefficient Dr = R2/4τ, where R is the radius of the observed emission and τ is the minimum of the system age and the total (energy-dependent) cooling time. Adopting a half-angle of ≈1° (see Cao et al. 2021a, Fig. 1) at a distance of 0.8 kpc, we find, for example, at a representative energy of 100 TeV, Dr = 1.5 × 1027 cm2 s−1. While this is smaller than the Galactic average of , it is consistent with the Bohm limit for a 3 μG field, DBohm(100 TeV) = 1.1 × 1027 cm2 s−1.

LHAASO J2226+6057 is a source of considerable interest (e.g. Fujita et al. 2021; Ge et al. 2021), with new multi-wavelength data forthcoming. Such broadband coverage will be essential for discriminating different interpretations regarding the origins of the UHE emission.

3.2. LHAASO J1908+0621

This source overlaps with the well-known very-high energy (VHE; e.g. Aharonian et al. 2009) and UHE source eHWC J1907+063 (Abeysekara et al. 2020). However, since two powerful pulsars are spatially coincident with the LHAASO source, we must consider them separately. We performed a joint fit using the new data from the LHAASO observatory together with the HAWC data. The result for B = 3 μG associated with the pulsar PSR 1907+0602 is shown in Fig. 2 along with the data from LHAASO (Cao et al. 2021a) and the data from HAWC (Abeysekara et al. 2020). We also show the systematic and statistical error band from HAWC.

thumbnail Fig. 2.

Model fit for LHAASO J1908+0621 and eHWC J1907+063 together with the corresponding data and the model from B21 for eHWC J1907+063. Red data points show the data from the LHAASO collaboration (Cao et al. 2021a) and the blue triangles data from the HAWC collaboration for eHWC J1907+063 (Abeysekara et al. 2020). The blue band shows the systematic and statistical error for eHWC J1907+063 (Abeysekara et al. 2020) and the dashed-dotted blue line the corresponding model from B21 extended to higher energies and corrected for absorption. The orange curve represents the joint fit to the data from HAWC and LHAASO for B = 3 μG, resulting in 0 = 1 × 1036 erg−1 s−1 and α = 2.66, assuming an association with PSR 1907+0602. Ecut was only constrained by the fit to be larger than 1.4 PeV and was set to equal the Hillas limit of 4.2 PeV. The solid green line shows the result from the fit that includes the radio upper limits with the B-field as a free parameter but a fixed value for α = 2.2. The dotted lines in the corresponding colours show the emission without absorption.

The resulting fit parameters for a PSR 1907+0602 association are: 0 = (10 ± 2) × 10635 erg−1 s−1, α = 2.66 ± 0.03. The best fit cutoff is Ecut = 10 PeV, but any value above 1.4 PeV is consistent with the data at a 95% confidence level. The Hillas limit for a pulsar with the spin-down power of PSR 1907+0602 is 4.2 PeV. In Fig. 2 this value was chosen for Ecut. This scenario requires ∼40% of the spin-down power of the pulsar in electrons above 1 TeV. To avoid exceeding the available energy budget, electrons should be injected at a minimum energy of 240 GeV. However, if a break in the injection spectrum below 1 TeV exists, this condition can be relaxed. Data in the GeV to TeV regime can help to detect possible spectral breaks, and evidence for such a break exists (Crestan et al. 2021). For an association with PSR 1907+0632, the values 0 = (4.2 ± 0.5) × 1036 erg−1 s−1 and α = 2.79 ± 0.03 are obtained. The cutoff energy is constrained to be > 2 PeV at a 95% confidence level. In addition to being in tension with the Hillas limit of 1.8 PeV, approximately seven times the available pulsar power must be transferred to electrons above 1 TeV to account for the emission. Therefore, it is highly unlikely that PSR 1907+0632 alone can account for this source, even allowing for evolution of the injected power over the cooling timescale.

We also compared the fit here to the equilibrium model from B21 shown in Fig. 2. The curve, developed initially to fit just the HAWC data, provides a reasonable match to the UHE emission detected by LHAASO. Both models are nearly identical at low energies but diverge slightly at higher energies due to the higher cutoff energy in the new model. A more gradual turnover at the highest energies is needed to match the slightly reduced flux values of the LHAASO data relative to the HAWC data in the overlapping energy range. As shown previously (B21), data from the Infrared Astronomical Satellite (IRAS) survey exclude enhancement in the FIR radiation fields in the range 60 μm to 100 μm, and the radiation fields cannot be very different from the average Galactic field. This leaves the magnetic field as the only uncertain environmental parameter. The fit parameters for larger magnetic field values result in smaller α and smaller Ecut. Insisting that α ≳ 2.0 implies an upper limit on the magnetic field strength of ∼8 μG, with a resulting cutoff energy of at a 95% confidence level. However, additional physical effects, such as energy-dependent confinement, could in principle allow injection spectra harder than α = 2.0 to occur and should not be ruled out a priori.

For B = 3 μG the cutoff energy is constrained to be very close to 1 PeV. In scenarios with larger B-fields, the value for Ecut is reduced, but also less well constrained. Even for B = 8 μG, a cutoff energy of 1 PeV is still consistent within errors.

A detection of radio emission from the source region can help to constrain the magnetic field, but no emission on the scales of the LHAASO source extent has been detected so far. Duvidovich et al. (2020) calculated upper limits at 1.5 GHz and 6 GHz based on the non-detection of 1 mJy beam−1 (10 μJy beam−1) for a beamsize of 51.1″ × 39.5″ (12.2″ × 8.6″). With an estimated radial size of 1° (Cao et al. 2021a, Fig. 1), these limits translate into upper limits of 3.0 × 10−13 erg cm−2 s−1 (2.3 × 10−13 erg cm−2 s−1). The model is not necessarily consistent with these limits, but if electrons are indeed injected only above 240 GeV to fulfil the energy requirement, or a break in the injection spectrum below TeV energies exists, the emission is below the upper limit. Models with lower minimum injection energy and without a break fulfilling the upper limit exist as well: In a different approach, we required consistency with the radio limits and performed a fit with the B-field as a free parameter but with a fixed value for α = 2.2. This leads to B = (5.7 ± 0.3) μG, 0 = (2.8 ± 0.2) × 1035 erg−1 s−1, and . In this scenario, 30% of the spin-down power is injected into electrons above 1 TeV and a minimum injection energy of 10 GeV is required to be consistent with energy constraints.

The diffusion coefficient for B = 3 μG at 100 TeV for a source size of 1° is 1.3 × 1028 cm2 s−1, which is an order of magnitude above the value for Bohm diffusion (1.1 × 1027 cm2 s−1). For the B = 5.7 μG used in the second model, the diffusion coefficient is 3.9 × 1028 cm2 s−1, again far above the corresponding Bohm diffusion coefficient of 5.8 × 1026 cm2 s−1.

3.3. LHAASO J1825−1326

The environment of this source is equally complex. As with LHAASO J1908+0621, it also possesses a HAWC counterpart, and a similar joint fit was performed. The magnetic field in the electron cooling zone was again initially fixed to be 3 μG. Two candidate pulsars are potentially associated with LHAASO J1825−1326: PSR J1826-1334 and PSR J1826-1256. For an association with PSR J1826-1334, the fit gives 0 = (5 ± 2) × 1035 erg−1 s−1 and α = 2.23 ± 0.09, and Ecut is constrained to at a 95% confidence level. This model corresponds to η = 50% and is shown in Fig. 3. A minimum injection energy of 90 GeV is required to not exceed the available power input from the pulsar, but a spectral break below TeV energies relaxes this condition. The cutoff energy of the model is consistent with the Hillas limit of 4.2 PeV. The γ-ray emission from the fit for an association with PSR J1826−1256 is nearly identical for the energy range shown, and the resulting fit parameters in this case are 0 = (4 ± 2) × 1035 erg−1 s−1, α = 2.39 ± 0.09, and at a 95% confidence level, with corresponding efficiency η = 20%. A minimum injection energy of 25 GeV is required if no spectral break below 1 TeV exists. The cutoff energy is below the Hillas limit of 4.7 PeV.

thumbnail Fig. 3.

Model fit for LHAASO J1825−1326 and eHWC J1825−134 together with the corresponding data and the model from B21 for eHWC J1825−134. Red data points show the data from the LHAASO collaboration (Cao et al. 2021a) and the black triangles data from the HAWC collaboration for eHWC J1825−134. The blue band shows the systematic and statistical error for eHWC J1825−134 (Abeysekara et al. 2020) and the dashed-dotted blue line the corresponding model from B21 extended to higher energies and corrected for absorption. The orange curve represents the joint fit to the data from HAWC and LHAASO for B = 3 μG, resulting in 0 = 5 × 1035 erg−1 s−1, α = 2.23, and Ecut = 390 TeV. The fit result for the association with PSR J1826−1256 is nearly identical to the orange curve. The dotted lines in the corresponding colours show the emission without absorption.

The model from B21 for the HAWC source eHWC J1825−134 is depicted in Fig. 3 and closely follows the model from the joint fit. While the B-field was the same in both cases, the radiation fields were not: In B21 enhancement factors of 3 and 5 were used for the association with PSR J1826−1334 and PSR J1826−1256, respectively. Due to the presence of a cooling break in the new model, no enhancement is needed. Fits for a higher B-field value decrease the resulting value of α. For example, for B = 5 μG and an association with PSR J1826−1334, the resulting parameters are 0 = (1.8 ± 0.8) × 1035 erg−1 s−1, α = 1.9 ± 0.1, and . An enhancement in the radiation field will increase the resulting value for α. As shown in B21, enhancements of the IR and ultraviolet (UV) fields up to a factor of 16 are compatible with upper limits from IRAS data, and therefore a variety of models for larger B-field values are viable as long as the source-specific radiation fields are not constrained further.

As for LHAASO J1908+0621, no radio emission on the scales of the LHAASO source extent has been detected so far. Duvidovich et al. (2019) calculated an upper limit at 1.4ĠHz based on the non-detection of 0.24 mJy beam−1 for a beamsize of 9.24″ × 6.43″. With the estimated radial size of 1° given by Cao et al. (2021a), this limit translates into an upper limit of 2.3 × 10−12 erg cm−2 s−1. The models developed above fulfil this constraint. The diffusion coefficients at 100 TeV are 2.3 × 1028 cm2 s−1 (PSR J1826−1334) and 6.0 × 1027 cm2 s−1 (PSR J1826−1256), both of which are above the Bohm diffusion case of 1.1 × 1027 cm2 s−1.

As shown above, reasonable leptonic single-zone models exist for all sources. Due to the uncertainty in the environmental conditions, considerable flexibility in the fitting models is permitted. In general, however, larger magnetic field values require harder injection spectra, while enhanced radiation fields have the opposite effect. To constrain parameters and to confirm or rule out different scenarios, it is therefore necessary to either have more information about the specific environmental conditions or include multi-wavelength data in the modelling1. This requires a detailed source-specific investigation. For each source, multiple possible counterparts exist, and taking different sources into account might be necessary, especially at energies below several TeV. In all models, the cutoff energies were constrained to be between a few hundred TeV to more than 1 PeV in some cases. Nearby pulsars can supply sufficient power to account for the observed flux levels, and therefore the nebulae of high spin-down power pulsars may frequently house particles with energies in excess of 100 TeV. The estimated diffusion coefficients are orders of magnitudes below the interstellar CR diffusion coefficients but similar to estimates for the region around Geminga (Abeysekara et al. 2017). Such low diffusion coefficients could be caused by self-confinement (Evoli et al. 2018).

4. Conclusion

In this paper we have shown that the emission from LHAASO J2226+6057, LHAASO J1825−1326, and LHAASO J1908+0621 is consistent with simple one-zone models of IC emission provided that, in the case of LHAASO J2226+6057 and LHAASO J1825−1326, the magnetic fields in these systems are somewhat lower than typical Galactic average values or local radiation field energy densities are somewhat higher than the large-scale average. For each source, there exists at least one plausibly associated pulsar providing the necessary power to account for the observations above 1 TeV. This fact, coupled with the validity of the IC emission scenario as demonstrated here, can be seen as strong evidence that high spin-down power pulsars routinely accelerate electrons to hundreds of TeV and even PeV energies. The detection of PeV γ rays from the Crab nebula, recently reported by LHAASO, supports this claim (Cao et al. 2021b). Since all such neutron stars are thought to launch an ultra-relativistic wind, such extreme accelerators may be surprisingly common, provided they satisfy the associated Hillas condition (Hillas 1984). However, as discussed in B21, special conditions may be necessary for these sources to produce detectable fluxes.

Regions with strong radiation fields and/or low magnetic fields exist in many parts of the Galaxy (e.g. Breuhaus et al. 2021b), providing suitable conditions for hard IC emission at energies above 100 TeV. Radiation-dominated sources may be particularly prevalent close to the Galactic centre region (e.g. Hinton & Aharonian 2007), where the energy density in the diffuse photon field increases rapidly with decreasing galactocentric radius (Popescu et al. 2017). Future detection of new UHE γ-ray sources is anticipated, with many new results expected from the LHAASO observatory. Future southern hemisphere observatories, such as CTA South and SWGO, which are better positioned to observe the Galactic plane and central region, are highly desirable in this effort.

The search for a firm identification of the sources of PeV protons and nuclei is ongoing. To confirm or rule out different scenarios, it is crucial to have detailed information about the environmental conditions and to take multi-wavelength data into account. In some cases, it may be necessary to consider multi-source models. VHE γ-ray spectral information alone is, in most cases, not sufficient to discriminate between leptonic and hadronic scenarios, but we note that spatially resolved VHE-UHE emission might well be, in particular if the energy-dependent morphology implied by the cooling-limited IC picture presented here can be demonstrated. The CTA project is particularly important in this regard (Cherenkov & Acharya 2019).


1

Noting that great care is needed when combining data on different angular scales and/or variable surface brightness sensitivity to avoid reaching inappropriate conclusions.

Acknowledgments

The authors thank F. Aharonian for valuable comments. For the numerical calculations, we made use of the open source GAMERA code (Hahn 2015). We also used the python packages SciPy (Virtanen et al. 2020), NumPy (Harris et al. 2020), Astropy (Astropy Collaboration 2013, 2018) and Matplotlib (Hunter 2007).

References

  1. Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017, Science, 358, 911 [NASA ADS] [CrossRef] [Google Scholar]
  2. Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2020, Phys. Rev. Lett., 124, 021102 [NASA ADS] [CrossRef] [Google Scholar]
  3. Agaronyan, F. A., & Ambartsumyan, A. S. 1985, Astrophysics, 23, 650 [Google Scholar]
  4. Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009, A&A, 499, 723 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Albert, A., Alfaro, R., Ashkar, H., et al. 2019, ArXiv e-prints [arXiv:1902.08429] [Google Scholar]
  6. Albert, A., Alfaro, R., Alvarez, C., et al. 2021a, ApJ, 911, L27 [NASA ADS] [CrossRef] [Google Scholar]
  7. Albert, A., Alfaro, R., Alvarez, C., et al. 2021b, ApJ, 907, L30 [NASA ADS] [CrossRef] [Google Scholar]
  8. Amato, E., & Olmi, B. 2021, Universe, 7, 448 [NASA ADS] [CrossRef] [Google Scholar]
  9. Amenomori, M., Bao, Y. W., Bi, X. J., et al. 2019, Phys. Rev. Lett., 123, 051101 [NASA ADS] [CrossRef] [Google Scholar]
  10. Amenomori, M., Bao, Y. W., Bi, X. J., et al. 2021, Phys. Rev. Lett., 127, 031102 [NASA ADS] [CrossRef] [Google Scholar]
  11. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [Google Scholar]
  12. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  13. Atoyan, A. M., & Aharonian, F. A. 1999, MNRAS, 302, 253 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bell, A. R. 2013, Astropart. Phys., 43, 56 [NASA ADS] [CrossRef] [Google Scholar]
  15. Blumenthal, G. R., & Gould, R. J. 1970, Rev. Mod. Phys., 42, 237 [Google Scholar]
  16. Breuhaus, M., Hahn, J., Romoli, C., et al. 2021a, ApJ, 908, L49 [NASA ADS] [CrossRef] [Google Scholar]
  17. Breuhaus, M., Hahn, J., Romoli, C., et al. 2021b, Int. Cosmic Ray Conf., 37 [Google Scholar]
  18. Cao, Z., Aharonian, F., An, Q., et al. 2021a, Nature, 594, 33 [NASA ADS] [CrossRef] [Google Scholar]
  19. Cao, Z., Aharonian, F., An, Q., et al. 2021b, Science, 373, 425 [Google Scholar]
  20. Cherenkov, T. A. C., Acharya, B. S., et al. 2019, Science with the Cherenkov Telescope Array [Google Scholar]
  21. Crestan, S., Giuliani, A., Mereghetti, S., et al. 2021, MNRAS, 505, 2309 [NASA ADS] [CrossRef] [Google Scholar]
  22. Di Mauro, M., Manconi, S., Negro, M., & Donato, F. 2021, Phys. Rev. D, 104, 103002 [NASA ADS] [CrossRef] [Google Scholar]
  23. Duvidovich, L., Giacani, E., Castelletti, G., Petriella, A., & Supán, L. 2019, A&A, 623, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Duvidovich, L., Petriella, A., & Giacani, E. 2020, MNRAS, 491, 5732 [NASA ADS] [CrossRef] [Google Scholar]
  25. Evoli, C., Linden, T., & Morlino, G. 2018, Phys. Rev. D, 98, 063017 [NASA ADS] [CrossRef] [Google Scholar]
  26. Fujita, Y., Bamba, A., Nobukawa, K. K., & Matsumoto, H. 2021, ApJ, 912, 133 [NASA ADS] [CrossRef] [Google Scholar]
  27. Ge, C., Liu, R.-Y., Niu, S., Chen, Y., & Wang, X.-Y. 2021, The Innovation, 2, 100118 [NASA ADS] [CrossRef] [Google Scholar]
  28. Giacinti, G., & Kirk, J. G. 2018, ApJ, 863, 18 [NASA ADS] [CrossRef] [Google Scholar]
  29. Gould, R. J., & Rephaeli, Y. 1978, ApJ, 225, 318 [NASA ADS] [CrossRef] [Google Scholar]
  30. Hahn, J. 2015, Int. Cosmic R. Conf., 34, 917 [NASA ADS] [Google Scholar]
  31. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
  32. Hillas, A. M. 1984, ARA&A, 22, 425 [Google Scholar]
  33. Hinton, J. A., & Aharonian, F. A. 2007, ApJ, 657, 302 [NASA ADS] [CrossRef] [Google Scholar]
  34. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  35. Kothes, R., Reich, W., & Uyanıker, B. 2006, ApJ, 638, 225 [NASA ADS] [CrossRef] [Google Scholar]
  36. Lagage, P. O., & Cesarsky, C. J. 1983, A&A, 125, 249 [NASA ADS] [Google Scholar]
  37. Liu, S., Zeng, H., Xin, Y., & Zhu, H. 2020, ApJ, 897, L34 [NASA ADS] [CrossRef] [Google Scholar]
  38. Popescu, C. C., Yang, R., Tuffs, R. J., et al. 2017, MNRAS, 470, 2539 [NASA ADS] [CrossRef] [Google Scholar]
  39. Sironi, L., Keshet, U., & Lemoine, M. 2015, Space Sci. Rev., 191, 519 [CrossRef] [Google Scholar]
  40. Sudoh, T., Linden, T., & Hooper, D. 2021, JCAP, 8, 010 [CrossRef] [Google Scholar]
  41. Tibet ASγ Collaboration (Amenomori, M., et al.) 2021, Nat. Astron., 5, 460 [NASA ADS] [CrossRef] [Google Scholar]
  42. Vernetto, S., & Lipari, P. 2016, Phys. Rev. D, 94, 063009 [NASA ADS] [CrossRef] [Google Scholar]
  43. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  44. Yu, H., Wu, K., Wen, L., & Fang, J. 2022, New Astron., 90, 101669 [NASA ADS] [CrossRef] [Google Scholar]
  45. Zdziarski, A. A., & Krolik, J. H. 1993, ApJ, 409, L33 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Adopted properties of the pulsars with possible association with the three LHAASO sources (Cao et al. 2021a, and references therein).

All Figures

thumbnail Fig. 1.

Data and models for LHAASO J2226+6057. The orange curve is a model fit for a magnetic field of B = 3 μG and a fixed power-law index of the injection spectrum of α = 2.2, leaving 0 and Ecut as free parameters. The resulting values were 0 = 9 × 1033 ergxs s−1 and Ecut = 420 TeV. The dotted line shows the emission without absorption.

In the text
thumbnail Fig. 2.

Model fit for LHAASO J1908+0621 and eHWC J1907+063 together with the corresponding data and the model from B21 for eHWC J1907+063. Red data points show the data from the LHAASO collaboration (Cao et al. 2021a) and the blue triangles data from the HAWC collaboration for eHWC J1907+063 (Abeysekara et al. 2020). The blue band shows the systematic and statistical error for eHWC J1907+063 (Abeysekara et al. 2020) and the dashed-dotted blue line the corresponding model from B21 extended to higher energies and corrected for absorption. The orange curve represents the joint fit to the data from HAWC and LHAASO for B = 3 μG, resulting in 0 = 1 × 1036 erg−1 s−1 and α = 2.66, assuming an association with PSR 1907+0602. Ecut was only constrained by the fit to be larger than 1.4 PeV and was set to equal the Hillas limit of 4.2 PeV. The solid green line shows the result from the fit that includes the radio upper limits with the B-field as a free parameter but a fixed value for α = 2.2. The dotted lines in the corresponding colours show the emission without absorption.

In the text
thumbnail Fig. 3.

Model fit for LHAASO J1825−1326 and eHWC J1825−134 together with the corresponding data and the model from B21 for eHWC J1825−134. Red data points show the data from the LHAASO collaboration (Cao et al. 2021a) and the black triangles data from the HAWC collaboration for eHWC J1825−134. The blue band shows the systematic and statistical error for eHWC J1825−134 (Abeysekara et al. 2020) and the dashed-dotted blue line the corresponding model from B21 extended to higher energies and corrected for absorption. The orange curve represents the joint fit to the data from HAWC and LHAASO for B = 3 μG, resulting in 0 = 5 × 1035 erg−1 s−1, α = 2.23, and Ecut = 390 TeV. The fit result for the association with PSR J1826−1256 is nearly identical to the orange curve. The dotted lines in the corresponding colours show the emission without absorption.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.