Abstract
HTML
Reference
Related
PDF
-
As observational sensitivity improves, an increasing number of extended gamma-ray sources have been detected across the sky [1–5]. These large sources exhibit significant interactions with their surrounding environments. At present, most large extended ultrahigh-energy (UHE, photon energy
$ E\geq 100\,\mathrm{TeV} $ ) Galactic sources are associated with pulsar wind nebulae (PWNe) and TeV halos [6–9]. Young PWNe in their early expansion phase exhibit strong broadband emission across multiple wavelengths, although this emission is confined to relatively compact physical regions. As these nebulae evolve, the reverse shock from the supernova remnant (SNR) moves inward toward the pulsar. Owing to the asymmetric nature of the initial explosion and inhomogeneities in the surrounding medium, this process often distorts the morphology of the PWN, resulting in irregular and asymmetric structures, such as those observed in Vela-X [10, 11]. These young PWNe typically display energy-dependent morphological features. For example, in sources such as HESS J1825-137 and HESS J1303-631 [12–14], the emission region appears to shift toward the pulsar with increasing photon energy, reflecting the historical imprint of lower-energy electron-positron pairs left behind by earlier outflows. In addition, the spatial extent of the emission region tends to decrease significantly at higher energies. This behavior is reminiscent of X-ray observations and suggests that advection dominates particle transport in such systems. In contrast, TeV halos such as Geminga, Monogem, and LHAASO J0622+3754 [15, 16] also exhibit large angular extensions. However, the diffuse morphology of these halos, together with the extended particle distributions in the surrounding interstellar medium (ISM), implies that they are powered by middle-aged rather than young pulsars. Despite these observations, the number of well-studied examples of such large, extended gamma-ray sources remains limited. This scarcity hinders a comprehensive understanding of the mechanisms governing particle acceleration, transport dynamics, and evolutionary pathways in PWNe and related systems.Since becoming operational in 2021, the Large High Altitude Air Shower Observatory (LHAASO), one of the most sensitive TeV instruments, has uncovered numerous new faint extended sources. The newly discovered extended gamma-ray source LHAASO J2238+5900, reported in the First LHAASO Catalog [5], has been suggested to be linked to PSR J2238+5903, with a low chance probability of about 0.03%. PSR J2238+5903 is considered the most likely counterpart, given its high spin-down luminosity of approximately
$ 8.9 \times 10^{35}\ \mathrm{erg\ s^{-1}} $ and characteristic age of 26.6 kyr [17]. In the X-ray band, the pulsar was observed with XMM-Newton (Obs.ID=13286), but no extended nebula was detected because of the short exposure time of only 9.84 ks [17]. In the GeV regime, the pulsar exhibits clear pulsed emission detected by Fermi-LAT and is cataloged as 4FGL J2238.5+5903 [18, 19]. In the TeV regime, several observatories have conducted follow-up observations based on the Fermi results. VERITAS accumulated 32.5 hours of exposure while searching for pulsations from PSR J2238+5903 and reported a 1.9σ hint [20]. The MAGIC collaboration conducted a 44-hour observation targeting PSR J2238+5903 to search for a potential PWN, but no significant excess was observed. An upper limit on the differential flux was obtained at the 95% confidence level, approximately$ 10^{-13}\,{\rm{TeV^{-1}\,cm^{-2}\,s^{-1}}}$ at a few TeV [21]. The Tibet$ AS\gamma $ experiment reported faint indications of possible TeV emission, with a 2.5σ hint near the position of the pulsar [22]. The MILAGRO collaboration observed a higher significance of 4.7σ [23], which nevertheless remains insufficient to support a detailed analysis. Additionally, the extended gamma-ray region contains another pulsar, PSR J2240+5832, which has a high spin-down luminosity of$ 2.2 \times 10^{35}\ \mathrm{erg\ s^{-1}} $ and a characteristic age of 144 kyr. However, given its relatively large offset of approximately$ 0.6^\circ $ from the emission center, it is not considered a prime candidate.In this study, we present a detailed analysis of LHAASO J2238+5900 using LHAASO data. The paper is structured as follows: Section 2 describes the instrument and data analysis methods used in this study. Section 3 presents the analysis results for LHAASO J2238+5900. In Section 4, we discuss the origin of the emission and the broadband spectral properties of the source. Finally, Section 5 summarizes the conclusions of this study.
-
As a ground-based cosmic-ray observatory, LHAASO is located on Mount Haizi in Daocheng, Sichuan Province, China, at an altitude of 4,410 m. It is a hybrid observatory consisting of a 1.3
$ \rm{km}^2 $ array (KM2A), a Water Cherenkov Detector Array (WCDA), and wide-field-of-view air Cherenkov telescopes (WFCTA). Each sub-array and telescope can operate independently or in combination to observe air showers, such as PeV events from the Crab Nebula [24]. The sub-arrays are designed to cover distinct energy bands, enabling gamma-ray detection from hundreds of GeV to beyond 1 PeV. More information on the detectors can be found in [25]. After the secondary particles trigger the detectors, we reconstruct the arrival direction, core position, and energy of the showers. Using the parameter$ compactness $ , the rejection power of WCDA can reach up to 99.8% for large showers with 700 hits. With the aid of the muon detectors (MDs), KM2A achieves effective gamma-hadron discrimination, suppressing the hadron survival rate to 0.025% at 100 TeV. The celestial region from$ 0^\circ $ to$ 360^\circ $ in right ascension and$ -20^\circ $ to$ 80^\circ $ in declination is divided into spatial bins of size$ 0.1^\circ \times 0.1^\circ $ , which are filled with the numbers of detected events. The cosmic-ray background events in each bin are estimated using the direct integration method. More details on the reconstruction and performance of WCDA and KM2A can be found in [26, 27].Our analysis uses WCDA data from March 2021 to July 2025 across five bins of the number of triggered units (
$ N_{\text{hit}} $ : 100–200,200–300,300–500,500–800, and$ \geq $ 800), together with KM2A data from December 2019 to July 2025, covering 25 TeV to 1 PeV and logarithmically binned at five bins per decade. A three-dimensional likelihood analysis is performed to model the spatial morphology and the spectral energy distribution (SED) simultaneously. The statistical significance of the source is assessed using the test statistic (TS), defined as:$ \mathrm{TS} = 2 \ln \left(\dfrac{\mathcal{L}_{s+b}}{\mathcal{L}_{b}} \right) $ , where$ \mathcal{L}_{s+b} $ and$ \mathcal{L}_{b} $ denote the maximum likelihood values under the signal-plus-background and background-only (null) hypotheses, respectively. The distribution of TS is asymptotically$ \chi^2 $ , with the number of degrees of freedom given by the difference in the number of free parameters between the two hypotheses. In the case of a point source with a fixed position and spectral shape, for which the normalization is the sole free parameter, the pretrial significance is approximately$ \sqrt{\mathrm{TS}} $ . To eliminate the influence of nearby sources, a large region of interest (ROI) is selected as a$ 20^{\circ}\times10^{\circ} $ rectangle centered on the SNR G106.3+2.7 to simultaneously cover both sources. We iteratively add sources with a power-law spectral shape and a 2D Gaussian spatial model to the region until the log-likelihood value shows no significant improvement, corresponding to a significance threshold of$ \Delta {\rm{TS}}> 25 $ . A subsequent test of the log-parabola model was also performed to examine spectral curvature. Four extended 1LHAASO sources (J2200+5643u, J2228+6100u, J2229+5927u, and J2238+5900), along with a new source spatially associated with SNR G108.2-0.6, are potential contributors to the gamma-ray emission in the ROI. A detailed analysis of the sources other than LHAASO J2238+5900 is beyond the scope of this work and will be presented in a forthcoming paper.After fixing the number of sources, we compare the 2D Gaussian and diffusion models to determine the optimal spatial model for the target source. The diffusion template is taken from [16] and is represented by the equation:
$ f(\theta) \propto \frac{1}{\theta_{d}(\theta + 0.085\theta_{d})}{\rm{exp}}[-1.54(\theta/\theta_{d})^{1.52}], $
(1) where
$ \theta_d = (180^\circ/\pi) \cdot 2\sqrt{D(E_{e})t_{E}}/d $ represents the diffusion size,$ D(E_e) $ is the diffusion coefficient at energy$ E_e $ ,$ t_E $ takes the value of$ \min (t_{\rm cool}, t_{\rm age}) $ , and d is the distance to the pulsar. In addition, a forward-folding method is used, in which simulated excess events are compared with the observed results to derive the differential energy spectrum[28]. -
Between the two tested models, the 2D Gaussian template is favored, yielding a TS value for the full ROI that is approximately 42 higher than that of the diffusion model, with
$ \theta_d=1.49^\circ\pm0.06^\circ $ . Under the Gaussian model, the best-fit extension is$ \sigma=0.54^{\circ}\pm0.01^\circ $ , and the position is (RA=$ 339.48^{\circ}\pm0.04^{\circ} $ , Decl=$ 58.97^{\circ}\pm0.02^{\circ} $ ). Evidently, LHAASO J2238+5900 is in close proximity to PSR J2238+5903, with an angular separation of approximately$ 0.14^{\circ} $ , which is considerably smaller than the radius of our instrument's PSF. We rebinned the original data into three energy intervals: 3–30 TeV, 25–100 TeV, and$ >100 $ TeV. Figure 1 shows the significance maps of LHAASO J2238+5900 in the three energy intervals, with diffuse backgrounds and other sources subtracted. LHAASO J2238+5900 is detected with a significance of$ 7.9\sigma $ above 100 TeV, suggesting that it is a PeVatron candidate.
Figure 1. (color online) Significance maps of the LHAASO J2238+5900 region, divided into three energy bins. The lime-colored plus, lime-colored cross, and red star indicate the positions of PSR J2238+5903, PSR J2240+5832, and LHAASO J2238+5900 (this work), respectively. The white circle in the bottom-right corner indicates the size of the point-spread function (68% containment).
To investigate the position and morphology in greater detail, we performed a 3D likelihood analysis in each of the three energy intervals. The morphological and spectral parameters of the other sources were held fixed, while the extensions of LHAASO J2238+5900 and the normalizations of all sources in the ROI were allowed to vary. The surface-brightness distributions were profiled as a function of angular distance from the best-fit gamma-ray emission centroid, as shown in Figure 2. Within about
$ 1^\circ $ of the centroid, the excess events are primarily attributed to LHAASO J2238+5900. In our study, the gamma-ray emission centroid in the three energy intervals is measured to be offset from PSR J2238+5903 by$ 0.14^{\circ}\pm0.03^{\circ} $ ,$ 0.11^{\circ}\pm0.02^{\circ} $ , and$ 0.17^{\circ}\pm0.06^{\circ} $ , respectively. Considering the statistical uncertainties, the position of the gamma-ray emission does not show a significant variation with energy. In addition, the source size (the Gaussian width$ \sigma_{\text{ext}} $ ) is measured to be$ 0.58^{\circ}\pm0.02^{\circ} $ ,$ 0.50^{\circ}\pm0.02^{\circ} $ , and$ 0.43^{\circ}\pm0.04^{\circ} $ , respectively. The energy-dependent extent, shown in Figure 3, is well described by a power law of the form$ \sigma_{\text{ext}}=(0.71\pm0.06)\times E^{(-0.09\pm0.03)} $ ($ \chi^2/ndf=0.2/1 $ ). However, this represents only a$ 2.7\sigma $ improvement over a constant fit ($ \chi^2/ndf=7.4/2 $ ), suggesting that the data do not strongly favor an energy-dependent model.
Figure 2. (color online) One-dimensional distribution of the gamma-ray emission from LHAASO J2238+5900 after subtraction of the diffuse background and contributions from other sources. The solid line represents the best-fit Gaussian model, and the shaded band indicates the
$ \pm 1\sigma $ statistical uncertainty.
Figure 3. (color online) The measured extension (
$ \sigma_{\text{ext}} $ ) of LHAASO J2238+5900 across three energy intervals. The solid black line indicates the power-law fit, and the red dashed line indicates the result of the constant fit. The energy-dependent power-law model is slightly favored over the constant model at the$ 2.7\sigma $ level.Using the Gaussian template, the energy spectrum in the source region is derived, as shown in Figure 4. The spectrum exhibits a rapid cutoff at high energies and is parameterized as
$ N_0(E/E_0)^{-\Gamma}\exp(-{E/E_{c}}) $ .$ E_0 $ is the reference energy and is set to 30 TeV. A comparison of the fitting results shows a clear preference for the exponential cutoff spectral model ($ \chi^2/ndf=10.7/8 $ ) over the simple power-law model ($ \chi^2/ndf=343.4/9 $ ). The best-fit parameters are$ N_0 = (3.21\pm 0.28)\times 10^{-15}\,{\rm{TeV^{-1}\,cm^{-2}\,s^{-1}}}$ ,$ \Gamma=1.97\pm0.06 $ , and$ E_{\text{c}}=41.0\pm3.5\,\text{TeV} $ . The integrated flux above 1 TeV is$ 2.4\times10^{-12}\,{\rm{cm^{-2}\,s^{-1}}}$ , corresponding to approximately$ 12\% $ of the Crab flux. With the gamma-ray spectrum extending up to$ \sim $ 220 TeV, the inferred parent particle energies are approximately 560 TeV for electrons (via IC scattering [15]) and 2 PeV for protons (via$ \pi^0 $ decay [29]). As noted in [5], the primary systematic uncertainty in the spectrum arises from the atmospheric model used in the Monte Carlo simulations, leading to an effect of about$ 7\% $ on the flux and 0.02 on the spectral index for the KM2A data. For the WCDA data, the flux uncertainty is$ 8\% $ , while the uncertainty in the spectral index depends on the shape of the energy spectrum.
Figure 4. (color online) The spectral energy distribution (SED) of LHAASO J2238+5900 and the corresponding fit results.
-
The most likely counterpart of the gamma-ray emission from LHAASO J2238+5900 is PSR J2238+5903. Distance estimates for this pulsar vary depending on the method used. The Third Fermi-LAT Catalog of Gamma-ray Pulsars (3PC) places an upper limit of 11.3 kpc using the ``DMM'' method, which estimates the distance to a pulsar based on its dispersion measure [30]. In contrast, an alternative approach based on a phenomenological gamma-ray luminosity law suggests a much smaller distance of 2.11 kpc [31]. Recently, FAST detected radio pulsations from PSR J2238+5903, with a distance of approximately 7 kpc estimated using the YMW16 Galactic electron density model [32], and we adopt this distance in this work. Currently, the pulsar lacks deep X-ray observations. In addition, Fermi-LAT observations indicate weak off-pulse emission in this region. We derived the upper limits using the same spatial template as LHAASO J2238+5900 (see Appendix A for details). To investigate the gas distribution in the region, we used CO data from the Milky Way Imaging Scroll Painting (MWISP) survey, a large-scale CO survey led by the Purple Mountain Observatory in China [33]. Several small clouds are observed around the source; however, they are offset from the brightest region of LHAASO J2238+5900 and do not cover the extended emission area (see Appendix B). In addition, the precise distances to these clouds remain unknown, making it unclear whether they are physically associated with the LHAASO source. We also searched the catalog of [34] for cloud information along the line of sight to LHAASO J2238+5900 within an angular extent of
$ 0.5^{\circ} $ . Four clouds were identified at distances of 1.87, 4.97, 6.21, and 6.97 kpc, with relatively low densities peaking at only around 12$ {\rm{cm^{-3}}}$ . The fourth cloud, [MML2017]4246, has a distance consistent with that of PSR J2238+5903 and an average angular radius of$ 0.36^{\circ} $ . However, it is offset by$ \sim 0.5^{\circ} $ from the center of the gamma-ray emission (Figure B2).
Figure B2. (color online) Column density of molecular gas toward LHAASO J2238+5900. The extended region of LHAASO is marked by a solid circle. The lime-colored plus sign and cross indicate the positions of PSR J2238+5903 and PSR J2240+5832, respectively. The light blue dashed circle denotes the molecular cloud [MML2017]4246.
We consider both leptonic and hadronic origins and fit the Fermi-LAT and LHAASO data using the PYTHON package Naima [35–37]. Given the presence of a cutoff in the gamma-ray spectrum, the electron and proton spectra are modeled as power laws with exponential cutoffs. In the leptonic model, the gamma-ray photons are produced through inverse Compton (IC) scattering by electrons. We consider three seed photon fields: the Cosmic Microwave Background (CMB), far-infrared dust emission (FIR), and near-infrared stellar emission (NIR). The energy densities of these fields are
$ U=0.26, 0.32, 0.31\,{\rm{eV\,cm^{-3}}}$ , respectively. The corresponding temperatures are$ T=2.7, 29.0, 5360\,{\rm{K}}$ , as derived from the Interstellar Radiation Field (ISRF) model proposed by [38]. As shown in Figure A1, electrons with a spectral index of 1.55 and an exponential cutoff at 45 TeV can account for the data well. The total energy required for this population is$ 2.6\times10^{47} $ erg. In the hadronic model, we calculate the proton-proton process assuming a maximum density of 12$ {\rm{cm^{-3}}}$ . Although a proton spectrum with an index of 1.2 and a cutoff energy of 115 TeV can reproduce the Fermi-LAT and LHAASO data, it is substantially harder than the typical spectra of protons accelerated in SNRs. This suggests that other proton spectral shapes or multi-zone models may be required. Given the unresolved nature of proton acceleration efficiency in pulsars, our study focuses primarily on the leptonic origin.
Figure A1. (color online) Multiwavelength SED of LHAASO J2238+5900, showing the hadronic (dashed line) and leptonic (solid line) models. The red squares represent the spectral data points from LHAASO, while the green squares indicate the upper-limit points derived from Fermi-LAT, as described in APPENDIX A.
-
Based on a distance of 7 kpc, the corresponding 1–30 TeV gamma-ray luminosity is estimated to be
$ 7.2\times 10^{34}\,{\rm{erg\,s^{-1}}}$ , implying a gamma-ray conversion efficiency of 8%. This value is slightly higher than the typical range for very-high-energy (VHE) PWNe, for which efficiencies are generally between 0% and 7% [39]. It is plausible that the pulsar powers a PWN responsible for the photon production. In general, PWNe exhibit notable energy-dependent evolution, as demonstrated by sources such as HESS J1825-137 and HESS J1303-631 [12–14], in which both the spatial extent of the emission decreases significantly and the emission centroid moves progressively closer to the pulsar with increasing energy. The pulsars PSR J2238+5903 and PSR B1823-13 have comparable ages. However, the centroid of LHAASO J2238+5900 does not vary significantly with energy, and its extension shows only a marginally significant decreasing trend.The extension evolution indicates that particle transport around LHAASO J2238+5900 may be attributable to a diffusive transport mechanism. In a diffusion-dominated scenario, the diffusion coefficient is energy-dependent and is given by
$ D(E)=D_0(E_e/E_0)^{\delta} $ , leading to the expression for the extent size,$ R=\sqrt{2D\tau} $ , where τ is taken to be the minimum of the accelerator age and the cooling timescale of the radiating particles. For synchrotron and IC scattering, the cooling time is assumed to scale as$ \tau\propto1/E_e $ . In the Thomson regime, where$ E_{\gamma}\propto E_{e}^2 $ , we have$ R\propto E_{\gamma}^{(\delta-1)/4} $ (in the Klein–Nishina regime,$ E_{\gamma}\propto E_{e} $ and$ R\propto E_{\gamma}^{(\delta-1)/2} $ ). Typically, the parameter δ ranges from 0 to 1, corresponding to energy-independent diffusion and Bohm diffusion, respectively. From the fit to the energy-dependent extent, we derive$ \delta \approx 0.6 $ for the Thomson regime (and$ \delta\approx0.8 $ for the KN regime). In our scenario, assuming a well-determined pulsar age and a magnetic field of 3 μG, electrons with energies of roughly 30 TeV would have a cooling time comparable to the pulsar age. This implies that at lower energies, τ remains constant, corresponding to$ R\propto E_{\gamma}^{\delta/4} $ . Since our first energy interval lies below 30 TeV, the maximum extension could be larger around this energy, resulting in a steeper fitted line and a slightly smaller value of δ. -
Adopting a distance of 7 kpc for PSR J2238+5903, the physical extent enclosing 90% of the source emission is about 140 pc. Given this large physical scale, LHAASO J2238+5900 may be classified as a halo. Nonetheless, no significant deviations from the 2D Gaussian model are observed when the morphology is fitted with a diffusion template. This does not necessarily argue against the halo interpretation, as the morphology of pulsar halos can be affected by various factors, such as the particle transport mechanism (e.g., isotropic diffusion, ballistic transport, and anisotropic diffusion), non-steady injection, projection effects, and contamination from unresolved sources [40–42]. According to the relation between electrons and inverse Compton scattering photons [15], 20 TeV photons correspond to electrons with an energy of approximately 100 TeV. The derived diffusion coefficient for 100 TeV electrons, estimated at
$ 9.4\times10^{28}(d/7\ \text{kpc})^2\ \text{cm}^2\ \text{s}^{-1} $ , is approximately one order of magnitude higher than those inferred for other well-known TeV halos, including Geminga, PSR B0656+14 [15], and LHAASO J0621+3755 [16]. Since the inferred diffusion coefficient scales quadratically with distance, this discrepancy may be attributed to uncertainty in the distance measurement. Alternatively, it could be explained by the young nature of this halo, in which the diffusion-suppression mechanism has not yet fully developed.An interesting aspect of this system is the relative youth of PSR J2238+5903 compared with other pulsars that power TeV halos. The TeV halos observed thus far are typically associated with middle-aged pulsars, which have likely exited their SNRs and formed distinctive bow-shock PWNe. Electrons then escape from the nebula and begin diffusing through the surrounding ISM. Once the pulsar reaches a distance of roughly 0.68 times the SNR radius, its velocity surpasses the local sound speed. Pulsars such as W44 and the Vela pulsar, which are still embedded within their respective SNRs, are believed to have already developed BSPWN structures, as proposed by [43]. In these cases, the high pulsar velocity or distortion caused by the reverse shock can lead to the dispersal or escape of electrons [44–46]. Search results from the HAWC Collaboration revealed some extended TeV sources surrounding young systems [47]. For instance, PSR J0359+5414, a relatively young and radio-quiet pulsar, is believed to power a large extended TeV halo [48]. Finally, although no emission from the associated SNR has been detected at other wavelengths, the possibility that SNR contamination affects the LHAASO gamma-ray source cannot be ruled out.
-
In this work, we present precise measurements of the newly discovered UHE gamma-ray source LHAASO J2238+5900. A Gaussian fit indicates a spatial extension of
$ 0.54^{\circ}\pm0.01^{\circ} $ . The source is located at an angular distance of$ 0.14^{\circ} $ from PSR J2238+5903. Given its high spin-down luminosity, this pulsar is the most likely counterpart of LHAASO J2238+5900. Multiwavelength observations reveal no significant counterpart emission, and there is no compelling evidence for dense molecular or atomic material surrounding the gamma-ray source. A leptonic origin can naturally account for the observed photons. The high-energy electrons responsible for the inverse Compton scattering could be supplied by either a PWN or a TeV halo. To fully elucidate the transport mechanism, data with improved spatial resolution, such as those provided by Cherenkov telescopes, will be crucial. -
Photon data files with timing analysis results for PSR J2238+5903 are provided in 3PC, covering the period from August 4, 2008, to December 22, 2019 [30]. The data cover the energy range from 20 MeV to 1 TeV and are centered on the pulsar. Events were filtered using a zenith-angle cut (
$ \gt 105^{\circ} $ ) and poor-data-quality flags. In this analysis, we use the Python package$ FermiPy $ to analyze the same data. The ROI was selected as a$ 15^\circ\times15^{\circ} $ region centered on the position of 4FGL J2238.5+5903. The source models were constructed based on the latest Fermi Large Area Telescope Fourth Source Catalog [1]. Additionally, the Galactic and extragalactic diffuse emission models were obtained from gll_iem_v07.fits and iso_P8R3_SOURCE_V3_v1.txt, respectively. Since PSR J2238+5903 is a bright gamma-ray pulsar, we adopt an off-pulse analysis using the same phase intervals (0.375–0.5 and 0.625–1.0) defined in [49]. Therefore, we derived the upper-limit flux using the same spatial template as LHAASO J2238+5900 ($ \sigma=0.5^\circ $ ). The results are shown in Figure A1. -
The MWISP survey covers most of the Galactic plane at low Galactic latitude (
$ |b|<5^{\circ} $ ). It includes data from three CO isotopologue lines:$ {}^{12}\rm{CO(J=1\rightarrow 0}) , {}^{13}\rm CO(J= 1\rightarrow 0) $ , and$ \rm{C{}^{18}O(J=1\rightarrow 0}) $ . The$ {}^{12}\rm{CO}$ line has a noise level of approximately 0.5 K at a velocity resolution of 0.16$ \,{\rm{km\,s^{-1}}}$ , whereas the$ {}^{13}{\rm{CO}}$ and$ {\rm{C^{18}O}}$ lines have a noise level of approximately 0.25 K at 0.17$ \,{\rm{km\,s^{-1}}}$ . In this work, we select a region of$ 4.5^{\circ} \times 4.0^{\circ} $ centered on PSR J2238+5903. Along the line of sight, three distinct velocity intervals,$ (-7\sim-1) $ ,$ (-15\sim-7) $ , and$ (-58\sim -48)\, {\rm{\,km\,s^{-1}}}$ , are above the background level (Figure B1). The cloud's column densities are shown in Figure B2.
Figure B1. (color online)
$ {}^{12}{\rm{CO}}$ spectrum over the selected region.
The Extended Ultrahigh-energy Gamma-Ray Emission in the Vicinity of PSR J2238+5903
- Received Date: 2026-07-23
- Available Online: 2026-12-01
Abstract: We present a comprehensive analysis of the recently discovered TeV gamma-ray source LHAASO J2238+5900. Based on data collected by LHAASO, our fitting results indicate that the source is significantly extended, with an angular extension of $0.54^{\circ}\pm0.01^{\circ}$, and is spatially coincident with the pulsar PSR J2238+5903. Its spectrum is characterized by a power law with an index of $\Gamma=1.97\pm0.06$ and an exponential cutoff at $41.0\pm3.5$ TeV. In addition, the source exhibits a significant signal of 7.9σ above 100 TeV, implying that it is a PeVatron candidate. Although the gamma-ray emission is consistent with a pulsar wind nebula (PWN) scenario, the relatively large extension also permits a halo interpretation, potentially caused by electron-positron pairs escaping from the PWN.
DownLoad: