The ultrahighenergy UHE gammaray source 1LHAASO J00077303u is positionally associated with the composite SNR CTA1 that is located at high Galactic Latitude bapprox 10.5circ. This provides a rare opportunity to spatially resolve the component of the pulsar wind nebula PWN and supernova remnant SNR at UHE. This paper conducted a dedicated data analysis of 1LHAASO J00077303u using the data collected from December 2019 to July 2023. This source is well detected with significances of 21sigma and 17sigma at 8100 TeV and gt100 TeV respectively. The corresponding extensions are determined to be 0.23circpm0.03circ and 0.17circpm0.03circ. The emission is proposed to originate from the relativistic electrons and positrons accelerated within the PWN of PSR J00077303. The energy spectrum is well described by a powerlaw with an exponential cutoff function dNdE  42.4pm4.1fracE20rm TeV2.31pm0.11expfracE110pm25rm TeV rm TeV1 cm2 s1in the energy range from 8 TeV to 300 TeV implying a steadystate parent electron spectrum dNedEepropto fracEe100rm TeV3.13pm0.16expfracEe373pm70rm TeV2 at energies above approx 50 rm TeV. The cutoff energy of the electron spectrum is roughly equal to the expected current maximum energy of particles accelerated at the PWN terminal shock. Combining the Xray and gammaray emission the current spaceaveraged magnetic field can be limited to approx 4.5rm mu G. To satisfy the multiwavelength spectrum and the gammaray extensions the transport of relativistic particles within the PWN is likely dominated by the advection process under the freeexpansion phase assumption.
['Astrophysics - High Energy Astrophysical Phenomena']
Deep view of Composite SNR CTA1 with LHAASO in gammarays up to 300 TeV
{'No Header': '. Article .  \nSPECIAL TOPIC:  \nJanuary 2023 Vol. 66 No. 1: 000000  \nhttps: // doi.org / ??', 'Deep view of Composite SNR CTA1 with LHAASO in γ -rays up to 300 TeV': 'LHAASO Collaboration *  \n(The LHAASO Collaboration authors and a ffi liations are listed after the references.)  \nReceived ; accepted  \nThe ultra-high-energy (UHE) gamma-ray source 1LHAASO J0007 + 7303u is positionally associated with the composite SNR CTA1 that is located at high Galactic Latitude b ≈ 10 . 5 · . This provides a rare opportunity to spatially resolve the component of the pulsar wind nebula (PWN) and supernova remnant (SNR) at UHE. This paper conducted a dedicated data analysis of 1LHAASO J0007 + 7303u using the data collected from December 2019 to July 2023. This source is well detected with significances of 21 σ and 17 σ at 8 -100 TeV and > 100 TeV, respectively. The corresponding extensions are determined to be 0.23 · ± 0.03 · and 0.17 · ± 0.03 · . The emission is proposed to originate from the relativistic electrons accelerated within the PWN of PSR J0007 + 7303. The energy spectrum is well described by a power-law with an exponential cuto ff function dN / dE = (42 . 4 ± 4 . 1)( E 20 TeV ) -2 . 31 ± 0 . 11 exp( -E 110 ± 25 TeV ) TeV -1 cm -2 s -1 in the energy range from 8 TeV to 300 TeV, implying a steady-state parent electron spectrum dNe / dEe ∝ ( Ee 100 TeV ) -3 . 13 ± 0 . 16 exp[( -Ee 373 ± 70 TeV ) 2 ] at energies above ≈ 50 TeV. The cuto ff energy of the electron spectrum is roughly equal to the expected current maximum energy of particles accelerated at the PWN terminal shock. Combining the X-ray and gamma-ray emission, the current space-averaged magnetic field can be limited to ≈ 4 . 5 µ G. To satisfy the multi-wavelength spectrum and the γ -ray extensions, the transport of relativistic particles within the PWN is likely dominated by the advection process under the free-expansion phase assumption.  \nPWN, γ -rays , UHE  \nPACS number(s): 97.60.Bw;95.85.Pw;95.85.Ry;  \nCitation:  \nZhen Cao, et al,  \n, Sci. China-Phys. Mech. Astron. 66 , 000000 (2023), https: // doi.org / ??', '1 Introduction': 'CTA1 is a famous composite supernova remnant (SNR), characterized by a shell SNR and a central pulsar wind nebula (PWN). It was discovered in the radio band and first proposed as a SNR by Harris & Roberts (1960). The detailed radio observation further revealed a northwest-incomplete shell with a diameter ≈ 1 . 8 · and a bridge-like structure towards the center, with a kinematic distance of ≈ 1 . 4 kpc according to the associated HI shell (Pineault et al., 1993, 1997; Sun et al., 2011). The first X-ray detection in the direction of CTA1 was made by ROSATSeward et al. (1995), which revealed both thermal emission from the outer shock-heated region and a non-thermal component from the central region. Subsequent studies of the central X-ray source (RX J0007.0 + 7302) with ASCA, XMM-Newton, Chandra and Suzaku observatory have resolved the X-ray emission', 'Physics, Mechanics & Astronomy': "into a point-like source likely corresponding to the central pulsar, and a di ff use nebula, implying a X-ray PWNe driven by the central active pulsar (Halpern et al., 2004; Lin et al., 2012; Slane et al., 1997, 2004).  \nThe point-like source has been firmly identified as a pulsar (PSR J0007 + 7303), owing to the successful search for pulsation in GeV band (Abdo et al., 2008) and X-ray band (Caraveo et al., 2010; Lin et al., 2010). PSR J0007 + 7303, with a period of ≈ 316 ms, has a su ffi cient spin-down power ˙ E = 4 . 5 × 10 35 erg s -1 to drive the X-ray PWN. In addition, the corresponding characteristic age τ c is ≈ 13 . 9 kyr,consistent with the SNR age derived by Slane et al. (2004). In GeV regime, searches for the extended emission associated with radio shell and / or X-ray nebulae are ongoing (Abdo et al., 2012; Acero et al., 2013; Li et al., 2016). The latest results of the E > 1 GeV extended emission searching shows a ≈ 0 . 98 · disk-shaped source, potentially associated with the SNR due to its overlap with the radio emission contours (Ackermann et al., 2018). In the E > 50 GeV energy range, Zhou et al. (2024) report a possible ∼ 0 . 4 · extended γ -ray emission from the PWN.  \nComposite SNR CTA1 is associated with an extended TeV γ -ray source, VER J0006 + 729, discovered by VERITAS experiment in the energy range 0 . 6 TeV -17 . 8 TeV (Aliu et al., 2013). The TeV source VER J0006 + 729 shows an elliptical Gaussian morphology with 1 σ angular extension of 0 . 30 · along the major axis and 0 . 24 · along the minor axis, with the orientation of the major axis of ≈ 17 . 8 · west of north, centered at the position of the pulsar PSR J0007 + 7303. Due to the positional and morphological coincidence of the TeV emission with the X-ray PWN, the TeV source VER J0006 + 729 is discussed as TeV PWN in studies such as Aliu et al. (2013) and Torres et al. (2014a). Generally, PWNe are clouds of magnetised plasma that are created inside SNRs by the highly relativistic outflow ('wind') of a pulsar, observed via synchrotron emission produced when energetic electrons interact with the magnetic field, or through inverse Compton (IC) radiation generated by the scattering of electrons o ff the background photon fields (see Gaensler & Slane (2006) and Mitchell & Gelfand (2022), for a comprehensive review). The latter mechanism is believed to be relevant to the TeV γ -ray emission of the PWN. Alternatively, hadronic mechanisms could also be the origin of TeV emission, although unidentified by observation, in which case relativistic hadrons collide with the ambient medium, producing TeV emission through π 0 decay. The study of UHE γ -ray emission is crucial in determining whether the γ -ray emission is predominantly from hadronic processes or the IC process, as IC emission above 100 TeV energies undergoes suppression due to the Klein-Nishina e ff ect.  \nIndeed, the vast majority of Galactic TeV emittiers have been identified in the population of PWNe (H. E. S. S. Collaboration et al., 2018). Among them, CTA1 holds a special status as it is located at a relatively high Galactic Latitude b ≈ 10 . 5 · , experiencing minimal influence from Galactic di ff use emission and nearby TeV sources. This positioning makes it an excellent candidate for accurate measuring characteristics of TeV emission, serving as a valuable object for probing the physical model of TeV PWNe and distinguishing between leptonic and hadronic processes. Thanks to the wide field view, high sensitivity, and broad energy range of the Large High Altitude Air Shower Observatory (LHAASO) (Cao, 2010), our understanding of CTA1 has made impressive progress. Using approximately 2 years of data, the LHAASO collaboration has reported an extended TeV source, 1LHAASO J0007 + 7303u, which is tentatively shaped by a Gaussian with 1 σ angular extension of ∼ 0 . 2 · at energies E > 25 TeV and is positionally coincident with the X-ray / TeV PWN in Composite SNR CTA1. Particularly noteworthy is the significant detection of the ultra-high-energy (UHE, E > 100 TeV) γ -ray emission from the source 1LHAASO J0007 + 7303u, at a 13 σ confidence level (Cao et al., 2023).  \nThis paper presents a deep observation of CTA1 with LHAASO using approximately 3 years of data, and proposes some possible mechanisms for the UHE emission through dedicated data analysis and discussion. Section 2 provides a brief introduction to the LHAASO detector array and the analysis method. In Section 3, we report the analysis results based on LHAASO data. A discussion of the multi-wavelength observation and a model for CTA1 is presented in Section 4. Finally, we summarize the conclusions in Section 5.", '2 LHAASO Data Analysis': 'LHAASO,constructed on Mountain Haizi, in Sichuan province, China, is a complex extensive air shower (EAS) array with a high sensitivity ( ∼ 1% CU 1 year) and a wide field of view (FOV, ∼ 2.24 Sr for the maximum zenith angle of 50 · ) for CRs and γ -rays. It consists of three subarrays, i.e., Water Cherenkov Detector Array (WCDA), Kilomiter Square Array (KM2A) and Wide Field-of-view Cherenkov Telescope Array (WFCTA). By Combining WCDA and KM2A, LHAASO can cover the energy range from ∼ 1 TeV to > 1 PeV for γ -ray observation, with the low-energy threshold depending on the zenith angle distribution to some extent. The performance of WCDA and KM2A has been studied in detail employing the Monte Carlo simulations (Cao et al., 2024a), and calibrated using the measurements of Crab Nebula as a standard candle (Aharonian et al., 2021a,b).', '2.1 Data selection and binning': "The KM2A events utilized in this study were taken during the period from 17th December 2019 to 31st July 2023. The data quality control system and the long-term performance monitoring of KM2A data can be found in (Cao et al., 2024b). After rigorous data selection, the total live time amounts to 1216 days. These events were further reduced and reconstructed, according to the selection criteria and reconstructing methods outlined in Aharonian et al. (2021a). We applied a γ -ray / background discrimination cut to select out all γ -like events, with the survival fraction of cosmic ray background events at approximately 4 × 10 -4 at 50 TeV energy compared to the CR observation. For the analysis presented in this paper, only KM2A γ -like events with zenith angles less than 50 degrees and reconstructed energies ( Erec ) above 40 TeV were included. The low-energy threshold is dictated by the higher threshold for the higher zenith angle of the source. At the LHAASO site, the CTA1 (Decl. = 73 · ) culminates at a zenith angle θ = 43.7 · and lies at zenith angles θ < 50 · for 6.2 hours per sidereal day. We binned our selected data into 5 logarithmically spaced Erec bins per decade. Within each Erec bin, we generated an 'on map' by filling events into a grid of the spatial pixels with dimensions of ∆ R . A . × ∆ Decl . = 0 . 1 · × 0 . 1 · based on the reconstructed direction. The number of isotropic background events, still predominantly cosmic-ray events, in each spatial pixel was estimated using the 'direct integral method' (Fleysher et al., 2004), from which we derived the 'background map'.  \nThe WCDA data used in this study were collected from 5th March 2021 to 31st July 2023, covering a total of 796 days of live time. To obtain the available γ -like event set for source analysis, we implemented quality cuts, direction reconstruction, and γ -ray / background discrimination, as detailed in Aharonian et al. (2021b). The number of hit ( Nhit ) for each event was selected as a shower energy estimator. We divided our data into 5 intervals of Nhit , namely [100,200], [200,300], [300,500], [500,800], and [800, 2000]. Within each Nhit bin, we generated the 'on map' and 'backgound map' following the aforementioned procedure for KM2A data.  \nA region of interest (ROI) specific to the CTA1 was defined as a 6 · × 6 · square region centered at the position of pulsar PSR J0007 + 7303 (R.A. = 1.757 · ; Decl. = 73.052 · ). We selected the 'on map' and 'background map' within ROI and for all energy bins to following analysis. In total, there are 14 energy bins represented by Nhit bins for WCDA data and Erec bins for KM2A data. The angular and energy resolution are dependent on the zenith angle and the energy. Considering the position of pulsar PSR J0007 + 7303, the distribution of zenith angle is constant, and thus the median energy in each energy bin slightly varies with the spectral shape of the source.", '2.2 Maximum Likelihood Analysis and Statistic Test': "We utilized a maximum likelihood fit using parametric spatial and spectral models to determine the statistical significance of γ -ray emissions detected by LHAASO. This process involved generating expected source counts map ('src map') of γ -ray signals in our ROI through forward-folded models with the detector response. We then compared the sum of the model counts listed in 'src map' and the background events listed in 'background map' with the actual observed counts ( Nobs ) listed in 'on map'. Given a model θ with spatial and spectral parameters, we maximized the likelihood of the model comparing with Nobs as follows:  \nln L ( Nobs | θ ) = Nbins X j = 0 lnP ( N j obs | θ ) (1)  \nwhere Nbins is the number of the data bins in the ROI and the analysis energy, P is the Poisson probability of detecting N j obs events in each bin given the model parameters θ .  \nThe ratio of the maximum likelihood defines a test statistic (TS):  \nTS = 2 ln( L ( θ 1) / L ( θ 0)) (2)  \nwhere θ 0 and θ 1 represent the parameters of the model of the null and alternative hypotheses, respectively. Based on Wilks' Theorem, the TS follows a χ 2 n distribution, where n is the degrees of freedom (dof.) derived from the di ff erence in the number of free parameters between the models. Typically, detection significance is determined by comparing the likelihood of the background-only model (null hypothesis) with that of the signal plus background model (alternative hypothesis). Unless otherwise mentioned, the TS utilized in this study is for detection significance. We employed the di ff erence in detection significance  \n( ∆ TS) to compare the spatial and / or spectral models of the source. It is important to note that the di ff erence in TS cannot be used to quantitatively determine the preferred model when the models are not nested.Alternatively, the Akaike information criterion test (AIC, Akaike, 1974) can be considered. The AIC is defined as AIC = 2 k -2 ln L , where k is the number of parameters in the model. In this context, the best hypothesis is considered to be the one that minimizes the AIC. A qualitative strength of evidence rules to assess the significance of a model is based on the di ff erence in AIC ( ∆ AIC) values between the two models. If ∆ AIC > 5, then it is considered strong evidence against the model with higher AIC and ∆ AIC > 10 is considered as decisive evidence against the model with higher AIC (Kulkarni & Desai, 2017; Liddle, 2007).", '3.1 Morphology': 'We evaluated the TS value for a point-like signal with a power-law spectral shape dN / dE ∼ E -2 . 7 at each pixel in the map. In this case, the significance is represented by √ TS . The √ TS map is displayed within a 2 . 5 · × 2 . 5 · square area centered at the position of pulsar PSR J0007 + 7303 in the energy range of 8 TeV < E < 100 TeV and E > 100 TeV, respectively, as illustrated in Figure 1. In both the 8 TeV < E < 100 TeV and E > 100 TeV energy ranges, an obvious γ -ray emission excess structure is observed in the inner region of the radio shell.  \nFigure 1 The significance maps of CTA1 in the energy ranges of 8 TeV < E < 100 TeV (a) and E > 100 TeV (b). Overlaid are the GB6 image countour (4850MHz) in gray. Cyan stars represent the position of the pulsar PSR J0007 + 7303.  \n<!-- image -->  \nWe studied the morphology of the γ -ray emission in the energy range of 8 TeV < E < 100 TeV and E > 100 TeV, respectively. Due to the impact of angular resolution, the actual morphology of the γ -ray emission in CTA1 should be addressed by using a forward-folded approach. We first explored the routine and mathematical-empirical geometrical models, such as point, Gaussian and disk. Motivated by the observations of VERITAS in the energy range of 0 . 6 TeV < E < 17 . 7 TeV, we also examined the elliptical Gaussian model. However, due to the limited angular resolution of LHAASO, which is 3-10 times larger than that of Image Atmosphere Cherenkove Telescope (IACT), we encountered convergence issues in the case of all parameters being free. To address this, we maintained the rotation angle fixed. For each energy range, the spectrum of the γ -ray emission in CTA1 region was initially modeled by a simple power-law for the 3D maximum likelihood analysis. Fit results are detailed in Table 1.  \nComparing the Gaussian model with the point model, we find that the Gaussian model improved the fit with a significance of 6.7 σ ( ∆ TS = 45 . 4) at energies 8-100 TeV and 3.5 σ ( ∆ TS = 13 . 2) at energies > 100 TeV, respectively. These results confirm that the TeV source detected by LHAASO is indeed extended. Using an elliptical gaussian intensity distribution instead of a Gaussian profile only marginally improves the fit likelihood by 1 . 6 σ ( ∆ TS = 2 . 8 ) at energies 8-100 TeV and by a negligible amount at energies > 100 TeV. This suggests that LHAASO cannot claim an obviously asymmetric structure as VERITAS has previously reported. To quantitatively compare the Gaussian and disk models, we evaluated their AIC values since these models  \nTable 1 Morphological Fit Results for CTA1.  \nP 95 is the statistical positional uncertainty at a 95% confidence level. The Major Axis is the radius of 100% flux region for Disk model and the radius of 39% flux region for Gaussian model. The AIC values are subtracted from that of the point model for a clear comparison. For the elliptical Gaussian model, we fixed the rotation angle (from west to north 17.8 · ) determined from VERITAS observations to ensure fitting convergence.  \nTable 2 Spectral Fitting Results for CTA1 above 8 TeV  \nNote: PL stands for the power-law model defined by dN / dE = N 0( E / E 0) -α , PLC represents the power-law with exponential cuto ff defined by dN / dE = N 0 × ( E / E 0) -α exp [ -( E / Ec )], and LP represents the log-parabola model defined by F ( E ) = N 0( E / E 0) -( α + β log 10( E / E 0)) . For the LP model, the latter value listed in α column is the β parameter. The AIC value is subtracted from that of the PL spectrum for a clear comparison.  \nare not nested. Due to a small di ff erence in AIC ( ∆ AIC < 4), we are unable to clearly distinguish between these two extended models. For the subsequent analysis, we tentatively consider the Gaussian profile as our benchmark spatial model. With this choice, the significance is estimated to be 21 σ and 17 σ at 8 -100 TeV and > 100 TeV, respectively, according to the TS values shown in Table 1.', '3.2 Spectrum': 'For the spectral analysis of CTA1, we performed a maximum likelihood fitting in the energy range above 8 TeV, using a Gaussian model described above. We compared three di ff erent spectral shapes for CTA1: a simple power-law (PL), a powerlaw with exponential cuto ff (PLC), and a log-parabola (LP). As shown in Table 2, the addition of a curvature in the spectrum is significantly improve the fit ( > 6 σ ). Comparing the LP and PLC spectrum, the best spectral model is PLC with an improvement in modeling with ∆ AIC = 8.4. Therefore, we can conclude that the best-fit spectral model is PLC model with an index of 2.31 ± 0.13 and a cuto ff energy of 110 ± 25 TeV. The integral energy flux above 8 TeV is F > 8 TeV ≈ 4 . 9 × 10 -12 erg cm -2 s -1 . We determined the spectral points by performing a maximum likelihood analysis in each Nhit or Erec bin considering the fixed PLC spectral shape. Upper limits on flux at a 95% confidence level were derived when CTA1 had TS < 4 (2 σ ) in a given bin. The di ff erential photon spectrum is shown in Figure 2, with spectral data points listed in Table 3.  \nUsing the PLC spectrum presented in Table 2, the energy of the gamma-ray events from the CTA1 region were re-estimated event by event. We detailed the events with reconstruction energy above 250 TeV considering more stringent exclusion of cosmic ray background. Four gamma-like events with reconstructed energy larger than 250 TeV are listed in Table 4. The  \nFigure 2 LHAASOdi ff erential γ -ray spectrum of the PWN CTA1. The red line and the orange butterfly represent the best-fit PLC model and its uncertainties, respectively.  \n<!-- image -->  \nmaximum energy of the events is ≈ 300 TeV.  \nTable 3 Di ff erential Flux Measurements of CTA1Note: Mid-Energy represents the median energy corresponding to di ff erent bins. The error is statistic. The systematic uncertainties of the flux are estimated to 8%.', '3.3 Systematic Uncertainties': 'To assess the robustness of our results we performed a number of systematic checks similar to those employed in Cao et al. (2024c). The pointing error were estimated to 0 . 04 · for WCDA and KM2A data. A conservative estimation of the systematic error in Gaussian extension measurement was on the order of 0 . 08 · . The systematic uncertainties in flux were estimated to be 8%, which arise from the atmospheric model used in the Monte Carlo simulations. Furthermore, we assessed the uncertainties arising from the GDE by incorporating a dust template into our fitting. The impact on the flux was found to be less than 2%.  \nTable 4 Photons with Erec > 250 TeV for CTA1  \nNote: Erec represents reconstructed energy based the PLC spectral shape. Ne is the number of electromagnetic particles. N µ is the number of muons detected in the region with a distance farther than 15 m from the core of the shower.', '4.1 The Nature of the γ -ray emission detected by LHAASO': "To provide a clear identification of the γ -ray emission in the energy range from 8 TeV to > 100 TeV detected by LHAASO, we compare its position and extension to observations in other bands, as shown in Figure 3. The γ -ray emission detected by LHAASO can be firmly associated with the X-ray PWN powered by the young pulsar PSR J0007 + 7303, as the position of the γ -ray emission coincides with that of X-ray PWN, while being clearly distinguished from the SNR radio shell. The nonthermal X-ray emission has an extension of approximately 0 . 3 · in ASCA observation (Slane et al., 1997) and 0 . 16 · in Suzaku observation (Lin et al., 2012). Due to the constraints of the small field of view of the X-ray telescope and the influence of thermal X-rays, it is possible that the size of the X-ray PWN is more extensive. Conservatively, we consider the lower limit (LL) size of the X-ray PWN at 0 . 3 · , and thus the 39% flux radius of the X-ray PWN is larger than 0 . 15 · assuming a Gaussian X-ray flux distribution. There is no conflict in extension between the X-ray and γ -ray observations. We favor a one-zone leptonic model in which the same population of electrons contributes to both the X-ray and γ -ray emission. Additionally, evidence for leptonic process includes possible energy-dependent morphology detection. The cooling of energetic electrons as they are transported away from the pulsar results in the size of the emission shrinks and its center moves closer to the pulsar at higher energies. This behavior is supported by observations from Fermi-LAT, VERITAS, and LHAASO, with a confidence level of 2.5 σ .  \nIt is widely believed that the γ -ray radiation in the PWNe is produced in the leptonic scenario. However, as the IC emission above 100 TeV energies undergoes suppression due to the Klein-Nishina e ff ect, the hadronic emission is possibly identified in UHE band. Indeed, Crab spectrum detected by LHAASO implies the possible contribution of the hadronic process due to the hardening at PeV energies(Liu & Wang, 2021). In the case of CTA1, we specifically investigate whether the UHE photons are dominated by the π 0 -decay gamma rays generated by the CR protons interacting with the interstellar medium (ISM). To investigate the environment around the pulsar, Mart'ın et al. (2016) conducted a dedicated study based on a survey data of the intensity of CO line provided in Dame et al. (2001). Most of the data in the velocity range -26 km s -1 < v < -6 km s -1 , which corresponds to the distance of CTA1, is tagged as noise. Due to the absence of a molecular cloud, we disfavor the UHE γ -ray emission originating from protons escaping a shock around CTA1 SNR and illuminating the molecular cloud. Additionally, specific regions were assigned a density of 0.017 and 0.037 cm -3 based on the X-ray absorption observations along the line of sight of CTA1, indicating a lower density in the vicinity of CTA1 (Slane et al., 1997). Assuming an average density of 0 . 1 cm -3 , which is the upper limit from the ISM density of the Galaxy in the direction of CTA1, we can conservatively estimate a total energy of proton with energies > 1 PeV to be Wp > 1 × 10 47 ergs based on the UHE γ -ray flux observed. The energy budget is larger than the total energies from the pulsar spin-down power, roughly estimated by ˙ E τ c ≈ 1 × 10 46 ergs. we can also exclude that the CTA1 PWN accelerates a su ffi cient number of PeV protons to dominate the UHE γ -ray emission observed.", '4.2 The broad-band spectrum: implications for the magnetic field and particle spectrum': "In the leptonic scenario, the broad-band nonthermal emission is dominated by two mechanisms: (1) synchrotron radiation for the emission from radio to X-ray, and (2) IC radiation for emission in the γ -ray band. In the one-zone model, the same population of relativistic electrons can produce the X-ray and γ -ray emission by interacting with the ambient magnetic field and the soft photon field, respectively.  \nThe γ -ray emission is contributed by IC scattering of relativistic electrons interacting with several target photon fields,  \nFigure 3 Multi-wavelength observations in CTA1 region. The circles represent 39% flux region. The dashed circle represents the lower limit of the extension. The size of the X-ray PWN has not been determined to date. The non-thermal X-ray emission is detected by ASCA in the 0 . 3 · region around the pulsar. Conservatively, the lower limit radius of the 39% X-ray flux region is estimated to 0 . 15 · , assuming a Gaussian flux distribution.  \n<!-- image -->  \nincluding the cosmic microwave background (CMB), infrared radiation field and star-light photon field. In the composite SNR CTA1 region, almost all of the VHE and UHE γ -ray fluxes come from the IC scattering with CMB, since the contributions from interactions with other photon fields are suppressed due to their low number densities and the Klein-Nishina e ff ect (Mart'ın et al., 2016; Torres et al., 2014a). Because the property of the CMB is well quantified, the γ -ray emission in LHAASO band could provide precise information about the spectrum of parent electrons. At energies above a few tens of TeV, the energy of the upscattered CMB photon E γ and that of the parent electron Ee could be linked through the following simple relation (LHAASO Collaboration et al., 2021),  \nEe ≃ 0 . 85( E γ/ 300 TeV) 0 . 77 PeV . (3)  \nThe energy range of γ -rays detected by LHAASO is from ≈ 8 TeV to ≈ 300 TeV, which corresponds to the parent electron energy ranging from ≈ 50 TeV to ≈ 850 TeV. We assume that the electron spectrum follows a steady-state power-law distribution terminated by a super-exponential cuto ff , i.e. E -α e exp[ -( Ee / Ee , c ) 2 ]. The spectral index and cuto ff energy of parent electron in the energy band 50 TeV ≲ Ee ≲ 880 TeV can be constrained to α = 3 . 13 ± 0 . 16 and Ee , c = 373 ± 70 TeV, respectively, by fitting to the LHAASO data as shown in Figure 4.  \nThe mean energies of the synchrotron ( Esyn ) and IC photons produced by the same population of electrons in the ambient magnetic field and CMB are related by  \nEsyn = 4 . 6( E γ/ 50 TeV) 1 . 5 ( B / 4 . 5 µ G ) keV . (4)  \nThe energy range of the X-rays detected by ASCA is Esyn = 0 . 5 -10 KeV, which corresponds to the γ -ray energy range of E γ ≈ 11 -83 TeV considering a magnetic filed of 4.5 µ G. This implies that the observed X-ray and γ -ray emission from  \nCTA1 are roughly produced by the population of electrons in the same energy range, in which the electron spectrum can be well constrained by LHAASO. Utilizing both the ASCA and LHAASO observations of the PWN, we can roughly constrain the current space-averaged magnetic field strength B to ∼ 4.5 µ G for the PWN CTA1, as shown in Figure 4. The deduced magnetic field is at the same level with the interstellar magnetic field. Such low magnetic field has also been reported in other PWNe by analysing their synchrotron spectra, like HESS J1825-137 and HESS J1809-193. The magnetic field in young PWNe decreases with time during the free-expansion phase, due to the expansion of the volume of PWNe. Therefore, the magnetic field inside PWNe is expected to reach its minimum at the end of the free-expansion phase. We also note that although the magnetic field in PWN CTA1 is as low as the interstellar magnetic field, this does not necessarily imply that the turbulence is also similar. The possible stronger turbulence in the PWN could prevent particles from escaping to the ISM. It's also worth noting that the non-thermal X-ray flux is not extracted from the entire TeV PWN region detected by LHAASO, but rather from a specific 0 . 3 · region within it. A more extended non-thermal X-ray survey around the entire TeV PWN is needed to provide a precise estimation of the current space-averaged magnetic field strength.  \nIn the PWN CTA1, the energy loss of the electrons is primarily dominated by synchrotron radiation, with a cooling timescale given by tc , syn ≈ 13 . 0( Ee / 40 TeV) -1 ( B / 5 µ G ) -2 kyr. Considering a time-averaged magnetic field strength of B ∼ 5 µ G, the cooling time roughly equals the pulsar's characteristic age of 13.9 kyr. Electrons with an energy greater than ≈ 40 TeV are expected to be in the fast cooling regime. Assuming that particle escape and the time evolution of energy injection rate can be neglected, the expected spectral index of the injected electrons can be approximated as p ∼ α -1 = 2 . 13 ± 0 . 16 in the energy range from ≈ 40 TeV to ≈ 373 TeV, considering the spectral steeping as being due to only the synchrotron loss. Conversely, electrons with Ee ≲ 40 TeV are in the slow cooling regime, where the spectral index of the integrated electrons remains almost unchanged with respect to that of the injected spectrum. Therefore, a spectral break is expected in the steady-state electron spectrum at an energy around 40 TeV, which is indicated by the possible break observed in the γ -ray flux around the transition energy ( ∼ 8 TeV) between the VERITAS and LHAASO energy range (see Figure 4). In our time-dependent model, we take into account the decrease in energy injection rate and magnetic field strength over time, as well as the potential escape of particles. Through a detailed numerical simulation of the evolution of the particle distribution (see appendix), we find that an intrinsic injection spectrum with an index of p = 2 . 2 could indeed result in a present electron spectrum with index of α ≈ 3 . 11 in the energy range from ≈ 40 TeV to ≈ 350 TeV, and of α ≈ 2 . 18 below ≈ 40 TeV. This is consistent with our previous simple explanation of the parent electron spectral properties of the PWN CTA1.  \nIn the case of standard shock acceleration in the Bohm di ff usion regime, a simple analytical presentation of the electron spectrum at shock over the entire energy range follows the form of a super-exponential cuto ff , i.e., Q ( Ee ) ∝ E -p e exp( -E 2 e / E 2 e , max ) (Zirakashvili & Aharonian, 2007). Considering a compression ratio of 3-4, the index p is estimated to be 2-2.25, which is consistent with the spectral index of the injected electrons inferred by the LHAASO observations. The maximum energy ( Ee , max ) of accelerated electrons is limited by the requirement of confinement of the particles inside the termination shock. This requires that the Larmor radius RL is smaller than the termination shock radius Rs , i.e., RL = ε Rs , where ε is the so-called containment factor with the value of 0 < ε < 1. The Larmor radius of the electron with energy Ee is  \nRL = Ee / ( eBs ) , (5)  \nwhere e is the electron charge, Bs is the post-shock magnetic field strength, defined as (Kennel & Coroniti, 1984)  \nBs ∼ (3( η BL ( t ) / c ) 0 . 5 ) / Rs , (6)  \nwhere L ( t ) is the spin-down luminosity of the associated pulsar at time t , η B is the magnetic energy fraction. Using Equation (5) and (6) and the condition RL = ε Rs , the maximum energy of the accelerated electrons is given by  \nE e , max = 3 e ε p η B L ( t ) / c . (7)  \nWe can see that the maximum energy of electrons evolves as the spin-down power with time. Taking into account the spindown power ˙ E = 4 . 5 × 10 35 erg / s at current time, the maximum energy can be estimated as Ee , max ≈ 430( ε 0 . 16 )( η B 0 . 6 ) 0 . 5 TeV. The cuto ff energy of the steady-state electron spectrum, approximately 373 TeV, should be the time-integrated value over the last few thousand years. However, the electrons injected in earlier stages experienced more energy losses, the time-integrated average maximum energy is expected to be close to the maximum energy of the recently injected electrons, which indicates that the maximum energy of electrons injected at current time could be estimated by the maximum energy inferred from our observation. Moreover, in our time-dependent modeling, the cuto ff energy of approximately 350 TeV in the accumulated electron spectrum  \n(see Fig A1 in the Appendex) indeed roughly coincides with the maximum energy of 430 TeV in the currently injected electron spectrum. This implies that the cuto ff energy derived from the steady-state electron spectrum could be a good estimation of the current particle acceleration ability of the termination shock in PWNe similar to CTA1.  \nFigure 4 The di ff erential energy spectrum of the PWN CTA1. The red line represents the expectation from a one-zone leptonic model, assuming a magnetic field strength of 4.5 µ G and considering only CMB target photons. The electron spectrum follows a broken power law distribution, where dNe / dEe ∝ E -3 . 13 e exp[ -( Ee / 373 TeV) 2 ] for energies E > 40 TeV and dNe / dEe ∝ E -2 . 13 e for energies E < 40 TeV. The time-dependent model (pure advection scenario) is referenced in the Appendix.  \n<!-- image -->", '4.3 The Morphology: implications for propagation mechanism': "Relativistic particles accelerated at the termination shock propagate in PWN mainly by two transport mechanisms, i.e. advection and di ff usion. To investigate the e ff ects of particle advection and di ff usion in the PWN, we performed detailed simulations of the spectrum and morphology of CTA1 for each transport scenario. The description of our model can be referred to the Appendix.  \nWe first tune the parameters for both pure advection and pure di ff usion scenarios to explain the multi-band non-thermal spectrum from the PWN CTA1. The parameters used to fit the flux data are almost similar in these two cases, except for those related to particle transport (see Table A1 in the Appendix). As a result, the calculated fluxes in the two cases are also similar, so we only plot the result in advection case for convenience (see Figure 4).  \nThe morphology of the emission is an interplay between particle transport and energy losses in the PWN. Figure 5 shows the energy dependence of the γ -ray extension for di ff erent transport mechanisms. In our model, the γ -ray extension at a given energy refers to the 39% flux size, in which the γ -ray flux is 39% of the total flux in the PWN. It is clear that overall the extension of the γ -ray emission decreases with increasing photon energies in both pure advection and di ff usion scenarios (for Kolmogorov turbulence), as the electrons responsible for producing higher energy γ -ray photons have a shorter lifetime, hence a shorter travelling distance. We can also see that the energy dependence of the size is not strictly a power-law in the whole energy range, because the electrons responsible for the lower energy γ -rays (below few TeV) are not completely cooled at present age. The calculated size in di ff usion scenario is less dependent on energy than that in advection scenario (the curve is flatter), particularly in the higher energy range. This can be explained by the fact that higher energy electrons di ff use faster, which compensates the e ff ect of stronger energy losses on the propagation distance.  \nFigure 5 shows that the size measured by Fermi-LAT, VERITAS and LHAASO could be well explained by the advection scenario, however, the extension calculated in the di ff usion scenario is much smaller. The reason for this di ff erence is that in the advection scenario, particles propagate faster in the inner nebula than in the outer region, causing them to accumulate at a larger radius. While in the case of di ff usion, the particle density decreases roughly exponentially with increasing distance from  \nthe central source, resulting in a more compact distribution of particles. The di ff erent impacts of particle transport on particle distribution are illustrated in Figure A2 in the Appendix. For the di ff usion scenario, it is unlikely to fit the data by simply increasing the di ff usion coe ffi cient, since larger di ff usion coe ffi cient means greater escape losses, which will make it hard to explain the observed flux. Moreover, when particles start to escape the PWN, further increasing di ff usion coe ffi cient would have very limited impact on the extension, because escaped particles are assumed not to contribute any emission outside the PWN. In fact, recent works suggest that escaped particles also produce IC emissions by interacting with the CMB photons, an e ff ect which we did not consider here (Martin et al., 2024). However, the similar size of the X-ray and γ -ray emission region may indicate that particle escape is not important for PWN CTA1, if the magnetic field is confined in the PWN.  \nIn our simulation, the radius of PWN is a crucial factor which can a ff ect the size of the emission. The PWN radius is mainly dependent on the mass of SN ejecta, the density of ISM, the initial spin-down luminosity and the age of the system. For simplicity, we assume that the evolution of PWN CTA1 is still in the free-expansion phase, and choose reasonable parameters to make the PWN radius as large as possible. For the chosen parameters (see Table A1 in the Appendix), the PWN radius at present time is about 11.3 pc (A larger radius requires extreme parameters). The largest 39% flux size of the γ -ray emission corresponding to this PWN radius in the pure di ff usion scenario is around 0.1 · , which is significantly smaller than the observed size (see Figure 5). On the other hand, we can estimate the upper limit of the PWN radius by the fact that it should be smaller than the corresponding SNR radius, which is ∼ 20.4 pc infered from the radio observation (Mart'ın et al., 2016). Even for such radius, the expected 39% flux size in the di ff usion scenario is ∼ 0.2 · , which still does not fit the data as well as the advection scenario. Therefore, the pure di ff usion scenario is unlikely to explain the observed PWN size, and we suggest that the particle transport is dominated by advection in the nebula of CTA1.  \nThe radius of PWN can also be influenced by its evolutionary stage. While we have assumed free-expansion phase in our calculation, it remains unclear whether the PWN CTA1 is in the free-expansion phase or the compression phase. Its relatively large age ( ∼ 10 kyr) compared to the typical young PWNe indicates it may have already passed the free-expansion phase. Nonetheless, the large PWN radius ( ∼ 10 pc indicated by observations) seems not to favor compression by the reverse shock. The attempt to fit the CTA1 spectrum using the reverberation model by Mart'ın et al. (2016) also failed due to the incapability of fitting the VERITAS data. This may imply that the CTA1 is in an intermediate phase, where the PWN is currently under transition from the free-expansion phase to the compression phase. There is another possibility that some parts of the PWN are in free-expansion phase, while other parts in compression phase, due to di ff erent density of the ISM in the northwest and southeast direction to the PWN, as suggested by Mart'ın et al. (2016). If the PWN CTA1 is indeed in this intermediate phase, our modelling based on the assumption of free-expansion phase may still applies to some extent, since the compression e ff ect of the reverse shock may be relatively small at that time. Finally, a more sophisticated PWN model is needed to study particle transport during the reverberation phase. Additionally, future detailed observations in radio and X-rays may help distinguish between the di ff erent evolutionary stages of the PWN CTA1.", '5 Conclusions': 'Using about three years of KM2A data and two years of WCDA data, we have a deep view of γ -rays emission from composite SNR CTA1 in energy range from ≈ 8 TeV to ≈ 300 TeV. The LHAASO source is detected with a significance of 21 σ in the energy band 8-100 TeV and 17 σ at energies above 100 TeV. The source is a significant extended source with a 39% flux radius of ≈ 0 . 23 · at energy band 8-100 TeV and of ≈ 0 . 17 · at energies above 100 TeV, based on a Gaussian profile. The best spectral model is PLC with dN / dE = (42 . 4 ± 4 . 1)( E 20 TeV ) -2 . 31 ± 0 . 11 exp( -E 110 ± 25 TeV ) TeV -1 cm -2 s -1 in the energy range from ≈ 8 TeV to ≈ 300 TeV. The integral energy flux above 8 TeV, F > 8 TeV ≈ 4 . 9 × 10 -12 erg cm -2 s -1 . Given that the emission is extended, with a centroid near the pulsar PSR J0007 + 7303 and its X-ray PWN, we confirm that the γ -ray emission detected by LHAASO is physically associated with the PWN in the composite SNR CTA1. Additionally, the plausible morphological evolution in the γ -ray emission can be confirmed by combining Fermi-LAT, VERITAS and LHAASO observations, further implying the γ -ray leptonic origin. We also investigate whether or not the UHE γ -ray emission is from a hadronic process. Due to the absence of a molecular cloud and the lower density of the ISM, we can almost rule out the UHE γ -ray emission being dominated by a hadronic process. In the leptonic scenario, we can accurately determine the electron spectrum above 50 TeV energies. The steady-state electron spectrum index of approximately 3.13 implies that the intrinsic injected electrons index should be around 2.13, which is consistent with the expected acceleration mechanism at the terminal shock of PWNe. The cuto ff energy might be treated as a good indicator of the maximum energy of the current injected spectrum. Combined with X-ray observations,  \nFigure 5 The γ -ray extensions of the PWN CTA1. The simulated extension corresponds to the 39% flux size, which means that the flux inside this region is 39% of the total flux in the PWN. The black and grey dashed lines represent the extension calculated by the pure advection and pure di ff usion scenarios.  \n<!-- image -->  \nthe current space-averaged magnetic field of the CTA1 PWN can be roughly limited to ∼ 4 . 5 µ G. Under the assumption of free-expansion phase, we favor advection dominating the particle transportation in the CTA1 PWN, in order to satisfy the multi-wavelength spectrum and the γ -ray extensions.', 'Acknowledgements': 'We would like to thank all sta ff members who work at the LHAASO site above 4400 meters above sea level year-round to maintain the detector and keep the water recycling system, electricity power supply and other components of the experiment operating smoothly. We are grateful to Chengdu Management Committee of Tianfu New Area for the constant financial support for research with LHAASO data. We deeply appreciate the computing and data service support provided by the National High Energy Physics Data Center for the data analysis in this paper. This research work is also supported by the following grants: in China by the National Natural Science Foundation of China NSFC No. 12393851, No.12393854, No.12393852, No.12393853, No.12022502, No.12205314, No.12105301, No.12261160362, No.12105294, No.U1931201, No.2024NSFJQ0060, and in Thailand by the National Science and Technology Development Agency (NSTDA) and the National Research Council of Thailand (NRCT) under the High-Potential Research Team Grant Program (N42A650868)..', 'AUTHOR CONTRIBUTIONS': 'S.Q. Xi, S.Z. Chen and B. Li led the drafting of the text. Y.Z. Li and S.Q. Xi conducted the data analysis, and Y. Huang performed the cross-check for the data analysis. B. Li and S.Q. Xi led the modeling e ff orts. S.Z. Chen and S.C. Hu provided the sky map and contributed to the background estimation. Zhen Cao, the spokesperson of the LHAASO Collaboration and the principal investigator of the LHAASO project. All other authors participated in data analysis, including detector calibration, data processing, event reconstruction, data quality checks, and various simulations, and provided valuable comments on the manuscript.  \nτ  \n0', 'Appendix': 'To obtain the particle spatial and energy distribution in PWN, we need to solve the advection-di ff usion equation (assuming spherical symmetry):  \n∂ N e( r , E e , t ) ∂ t = Q e( E e , t ) + 1 r 2 ∂ ∂ r " r 2 D ( E e , t ) ∂ N e( r , E e , t ) ∂ r # -1 r 2 ∂ ∂ r h r 2 V ( r ) N e( r , E e , t ) i + ∂ ∂ E e h ˙ E e N e( r , E e , t ) i , (A1)  \nwhere N e( r , E e , t ) is the di ff erential number density of electrons and positrons. The first term of the equation is the source term. The injection spectrum of electrons is assumed to follow a broken power-law with supper exponential cuto ff distribution. It\'s a good approximation to express the spectrum as (LHAASO Collaboration et al., 2021)  \nQ e( E e , t ) = Q 0( t ) E -α 2 e h 1 + ( E e / E b) -( α 2 -α 1) i -1 exp h -GLYPH<0> E e / E e , max( t ) GLYPH<1> 2 i , (A2)  \nwhere Q 0( t ) is the normalization factor, which can be determined by solving L ( t ) = R E e , max E e , min E e Q ( E e , t ) d E e. E e , min and E e , max are the minimum and maximum energy of injected electrons. The maximum energy of electrons accelerated at the termination shock could be estimated by Equation (7). E b is the break energy. α 1 and α 2 are the spectral indexes of the lower energy and higher energy part of the electron spectrum.  \nAssuming a braking index n = 3, the spin-down luminosity L ( t ) of the pulsar evolves with time as (Gaensler & Slane, 2006)  \nL ( t ) = ˙ E GLYPH<16> 1 + t age /τ 0 GLYPH<17> 2 (1 + t /τ 0) 2 , (A3)  \nwhere ˙ E is the pulsar\'s spin-down power at its current age t age  \nt age = τ c GLYPH<20> 1 -( P 0 P ) 2 GLYPH<21> , (A4)  \nwith τ c = P / 2 ˙ P the characteristic age of pulsar. P 0 and P are the initial and current spin periods, respectively. The initial spin-down timescale of the pulsar is defined as  \n=  \nτ  \nc(  \nP  \n0  \nP  \n)  \n.  \n(A5)  \nThe second term on the right hand side of Equation (A1) is the di ff usion term. We assume that the di ff usion coe ffi cient is homogeneous in the PWN, while it has time and energy dependence of D ( E e , t ) = D 0( t ) ( E e / 1 TeV) δ , with D 0( t ) being the normalization factor at 1 TeV and δ = 1 / 3 for the Kolmogorov turbulence. D 0( t ) is assumed to be inversely proportional to the magnetic field B ( t ). The evolution of magnetic field B ( t ) in PWN could be calculated by (Pacini & Salvati, 1973)  \ndW B( t ) dt = η B L ( t ) -W B( t ) R pwn( t ) dR pwn( t ) dt , (A6)  \nwhere W B = B 2 R 3 pwn / 6 is the total magnetic energy in PWN and η B the fraction of the spin-down luminosity converted into magnetic energy.  \nThe third term describes particle advection by the flow downstream the termination shock. The bulk velocity of the downstream plasma flow is supposed to decreases with radius as V ( r ) = V 0( r / R ts) -β , where R ts( t ) is the radius of termination shock and β (0 ≤ β ≤ 2) describes the spatial dependence of advection velocity.  \nThe last term represents the energy losses of relativistic electrons due to synchrotron radiation, IC scattering and adiabatic cooling, which are given by  \n˙ E e = 4 3 σ T γ 2       U B + X Ui GLYPH<0> 1 + 4 γϵ 0 , i GLYPH<1> 3 / 2       + E e 3 1 r 2 ∂ ∂ r h r 2 V ( r ) i , (A7)  \nwhere U B = B 2 / 8 π is the energy density of magnetic field and Ui is the energy density of the i -th component of the interstellar radiation field (ISRF). Following Torres et al. (2014b), we consider three components for the radiation field in the direction  \n2  \nof CTA1: the cosmic microwave background (CMB, T = 2.73 K, U = 0.25 eV cm -3 ), far-infrared radiation field ( T = 70 K, U = 0.1 eV cm -3 ), and near-infrared radiation field ( T = 5000 K, U = 0.1 eV cm -3 ). ϵ 0 , i is the average photon energy of the radiation field, which equals to 2 . 8 k B Ti for blackbody / greybody radiation with Ti being the temperature.  \nEquation (A1) is solved numerically in this work. The computational domain is the region between the termination shock and the outer boundary of the PWN. The evolution of the radii of termination shock and PWN ( R pwn( t ) and R ts( t )) could be obtained using the method given by Gelfand et al. (2009). At the termination shock, the inner boundary condition requires that the number of particles injected from the upper stream should equal to that transported downstream by advection and di ff usion (Peng et al., 2022). At the outer boundary, a free escape condition is imposed to simulate the particle escape from PWN. The model described above is then applied to the nebula of CTA1. The parameters of the model used to simulate the the particle distribution and emission from the nebula of CTA1 are shown in Table A1.  \nTable A1 Parameters used to fit the SED and γ -ray extensions of the PWN CTA1 in the pure advection and pure di ff usion scenarios.  \nThe time-integrated electron spectrum calculated by our model is shown in Figure A1. In order to compare the numerical result with our measurements, we use a broken power-law with an exponential cuto ff function to fit the spectrum from 0.5 TeV  \nto 1 PeV, in which the break energy is fixed at 40 TeV. The spatial distribution of electrons in di ff erent transport scenarios are plotted in Figure A2.  \nFigure A1 The spectral energy distribution of the electrons in PWN CTA1.  \n<!-- image -->  \n<!-- image -->  \nFigure A2 Particle number densities as a function of the distance from the termination shock in di ff erent transport scenarios. The parameters used to plot are given in Table A1. The particle distribution for the advection case (a) is nearly constant for the flow velocity profile V ∝ r -1 . 9 . The bump near the PWN outer boundary is due to the higher pulsar spin-down power in the early times. In the di ff usion case (b), the particle distribution decreases nearly exponentially with the radius.  \n<!-- image -->', 'References': "Abdo, A. A., Ackermann, M., Atwood, W. B., et al. 2008, Science, 322, 1218, doi: 10.1126 / science.1165572 Abdo, A. A., Wood, K. S., DeCesar, M. E., et al. 2012, ApJ, 744, 146, doi: 10.1088 / 0004-637X / 744 / 2 / 146 Acero, F., Ackermann, M., Ajello, M., et al. 2013, ApJ, 773, 77, doi: 10.1088 / 0004-637X / 773 / 1 / 77 Ackermann, M., Ajello, M., Baldini, L., et al. 2018, ApJS, 237, 32, doi: 10.3847 / 1538-4365 / aacdf7 Aharonian, F., An, Q., Axikegu, et al. 2021a, Chinese Physics C, 45, 025002, doi: 10.1088 / 1674-1137 / abd01b -. 2021b, Chinese Physics C, 45, 085002, doi: 10.1088 / 1674-1137 / ac041b Akaike, H. 1974, IEEE Transactions on Automatic Control, 19, 716  \nAliu, E., Archambault, S., Arlen, T., et al. 2013, ApJ, 764, 38, doi: 10.1088 / 0004-637X / 764 / 1 / 38  \nCao, Z. 2010, Chinese Physics C, 34, 249, doi: 10.1088 / 1674-1137 / 34 / 2 / 018  \nCao, Z., Aharonian, F., An, Q., et al. 2023, arXiv e-prints, arXiv:2305.17030, doi: 10.48550 / arXiv.2305.17030  \n- -. 2024a, Radiation Detection Technology and Methods, doi: 10.1007 / s41605-024-00467-8  \nCao, Z., Aharonian, F., Axikegu, et al. 2024b, arXiv e-prints, arXiv:2405.11826, doi: 10.48550  \n/  \narXiv.2405.11826  \n- Cao, Z., Aharonian, F., An, Q., et al. 2024c, ApJS, 271, 25, doi: 10.3847 / 1538-4365 / acfd29\n- Caraveo, P. A., De Luca, A., Marelli, M., et al. 2010, ApJL, 725, L6, doi: 10.1088 / 2041-8205 / 725 / 1 / L6\n- Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792, doi: 10.1086 / 318388\n- Fleysher, R., Fleysher, L., Nemethy, P., Mincer, A. I., & Haines, T. J. 2004, ApJ, 603, 355, doi: 10.1086 / 381384\n- Gaensler, B. M., & Slane, P. O. 2006, ARA&A, 44, 17, doi: 10.1146 / annurev.astro.44.051905.092528\n- Gelfand, J. D., Slane, P. O., & Zhang, W. 2009, ApJ, 703, 2051, doi: 10.1088 / 0004-637X / 703 / 2 / 2051\n- H. E. S. S. Collaboration, Abdalla, H., Abramowski, A., et al. 2018, A&A, 612, A2, doi: 10.1051 / 0004-6361 / 201629377\n- Halpern, J. P., Gotthelf, E. V., Camilo, F., Helfand, D. J., & Ransom, S. M. 2004, ApJ, 612, 398, doi: 10.1086 / 422409\n- Harris, D. E., & Roberts, J. A. 1960, PASP, 72, 237, doi: 10.1086 / 127521\n- Kennel, C. F., & Coroniti, F. V. 1984, ApJ, 283, 710, doi: 10.1086 / 162357  \nKulkarni, S., & Desai, S. 2017, Ap&SS, 362, 70, doi: 10.1007 / s10509-017-3047-6  \nLHAASO Collaboration, Cao, Z., Aharonian, F., et al. 2021, Science, 373, 425, doi: 10.1126 / science.abg5137  \nLi, J., Torres, D. F., de O˜na Wilhelmi, E., Rea, N., & Martin, J. 2016, ApJ, 831, 19, doi: 10.3847  \n/  \n0004-637X  \n/  \n831  \n/  \n1  \n/  \n19  \nLiddle, A. R. 2007, MNRAS, 377, L74, doi: 10.1111  \n/  \nj.1745-3933.2007.00306.x  \n- Lin, L. C. C., Huang, R. H. H., Takata, J., et al. 2010, ApJL, 725, L1, doi: 10.1088 / 2041-8205 / 725 / 1 / L1\n- Lin, L. C. C., Takata, J., Kong, A. K. H., et al. 2012, MNRAS, 426, 2283, doi: 10.1111 / j.1365-2966.2012.21832.x\n- Liu, R.-Y., & Wang, X.-Y. 2021, ApJ, 922, 221, doi: 10.3847 / 1538-4357 / ac2ba0\n- Mart'ın, J., Torres, D. F., & Pedaletti, G. 2016, MNRAS, 459, 3868, doi: 10.1093 / mnras / stw684\n- Martin, P., de Guillebon, L., Collard, E., et al. 2024, arXiv e-prints, arXiv:2407.07583, doi: 10.48550 / arXiv.2407.07583\n- Mitchell, A. M. W., & Gelfand, J. 2022, in Handbook of X-ray and Gamma-ray Astrophysics, 61, doi: 10.1007 / 978-981-16-4544-0 157-1\n- Pacini, F., & Salvati, M. 1973, ApJ, 186, 249, doi: 10.1086 / 152495\n- Peng, Q.-Y., Bao, B.-W., Lu, F.-W., & Zhang, L. 2022, ApJ, 926, 7, doi: 10.3847 / 1538-4357 / ac4161\n- Pineault, S., Landecker, T. L., Madore, B., & Gaumont-Guay, S. 1993, AJ, 105, 1060, doi: 10.1086 / 116493\n- Pineault, S., Landecker, T. L., Swerdlyk, C. M., & Reich, W. 1997, A&A, 324, 1152\n- Seward, F. D., Schmidt, B., & Slane, P. 1995, ApJ, 453, 284, doi: 10.1086 / 176388\n- Slane, P., Seward, F. D., Bandiera, R., Torii, K., & Tsunemi, H. 1997, ApJ, 485, 221, doi: 10.1086 / 304416\n- Slane, P., Zimmerman, E. R., Hughes, J. P., et al. 2004, ApJ, 601, 1045, doi: 10.1086 / 380498  \nSun, X. H., Reich, W., Wang, C., Han, J. L., & Reich, P. 2011, A&A, 535, A64, doi: 10.1051 / 0004-6361 / 201117679  \nTorres, D. F., Cillis, A., Mart'ın, J., & de O˜na Wilhelmi, E. 2014a, Journal of High Energy Astrophysics, 1, 31, doi: 10.1016 / j. jheap.2014.02.001  \n- -. 2014b, Journal of High Energy Astrophysics, 1, 31, doi: 10.1016 / j.jheap.2014.02.001  \nZhou, L., Wu, K., Gong, Y., & Fang, J. 2024, MNRAS, doi: 10.1093 / mnras / stae720  \nZirakashvili, V. N., & Aharonian, F. 2007, A&A, 465, 695, doi: 10.1051 / 0004-6361:20066494  \nZhen Cao 1 , 2 , 3 , F. Aharonian 4 , 5 , Axikegu 6 , Y.X. Bai 1 , 3 , Y.W. Bao 7 , D. Bastieri 8 , X.J. Bi 1 , 2 , 3 , Y.J. Bi 1 , 3 , W. Bian 9 , A.V. Bukevich 10 , Q. Cao 11 , W.Y. Cao 12 , Zhe Cao 13 , 12 , J. Chang 14 , J.F. Chang 1 , 3 , 13 , A.M. Chen 9 , E.S. Chen 1 , 2 , 3 , H.X. Chen 15 , Liang Chen 16 , Lin Chen 6 , Long Chen 6 , M.J. Chen 1 , 3 , M.L. Chen 1 , 3 , 13 , Q.H. Chen 6 , S. Chen 17 , S.H. Chen 1 , 2 , 3 , S.Z. Chen 1 , 3 , T.L. Chen 18 , Y. Chen 7 , N. Cheng 1 , 3 , Y.D. Cheng 1 , 2 , 3 , M.C. Chu 19 , M.Y. Cui 14 , S.W. Cui 11 , X.H. Cui 20 , Y.D. Cui 21 , B.Z. Dai 17 , H.L. Dai 1 , 3 , 13 , Z.G. Dai 12 , Danzengluobu 18 , X.Q. Dong 1 , 2 , 3 , K.K. Duan 14 , J.H. Fan 8 , Y.Z. Fan 14 , J. Fang 17 , J.H. Fang 15 , K. Fang 1 , 3 , C.F. Feng 22 , H. Feng 1 , L. Feng 14 , S.H. Feng 1 , 3 , X.T. Feng 22 , Y. Feng 15 , Y.L. Feng 18 , S. Gabici 23 , B. Gao 1 , 3 , C.D. Gao 22 , Q. Gao 18 , W. Gao 1 , 3 , W.K. Gao 1 , 2 , 3 , M.M. Ge 17 , T.T. Ge 21 , L.S. Geng 1 , 3 , G. Giacinti 9 , G.H. Gong 24 , Q.B. Gou 1 , 3 , M.H. Gu 1 , 3 , 13 , F.L. Guo 16 , J. Guo 24 , X.L. Guo 6 , Y.Q. Guo 1 , 3 , Y.Y. Guo 14 , Y.A. Han 25 , O.A. Hannuksela 19 , M. Hasan 1 , 2 , 3 , H.H. He 1 , 2 , 3 , H.N. He 14 , J.Y. He 14 , Y. He 6 , Y.K. Hor 21 , B.W. Hou 1 , 2 , 3 , C. Hou 1 , 3 , X. Hou 26 , H.B. Hu 1 , 2 , 3 , Q. Hu 12 , 14 , S.C. Hu 1 , 3 , 27 , C. Huang 7 , D.H. Huang 6 , T.Q. Huang 1 , 3 , W.J. Huang 21 , X.T. Huang 22 , X.Y. Huang 14 , Y. Huang 1 , 2 , 3 , Y.Y. Huang 7 , X.L. Ji 1 , 3 , 13 , H.Y. Jia 6 , K. Jia 22 , H.B. Jiang 1 , 3 , K. Jiang 13 , 12 , X.W. Jiang 1 , 3 , Z.J. Jiang 17 , M. Jin 6 , M.M. Kang 28 , I. Karpikov 10 , D. Khangulyan 1 , 3 , D. Kuleshov 10 , K. Kurinov 10 , B.B. Li 11 , C.M. Li 7 , Cheng Li 13 , 12 , Cong Li 1 , 3 , D. Li 1 , 2 , 3 , F. Li 1 , 3 , 13 , H.B. Li 1 , 3 , H.C. Li 1 , 3 , Jian Li 12 , Jie Li 1 , 3 , 13 , K. Li 1 , 3 , S.D. Li 16 , 2 , W.L. Li 22 , W.L. Li 9 , X.R. Li 1 , 3 , Xin Li 13 , 12 , Y.Z. Li 1 , 2 , 3 , Zhe Li 1 , 3 , Zhuo Li 29 , E.W. Liang 30 , Y.F. Liang 30 , S.J. Lin 21 , B. Liu 12 , C. Liu 1 , 3 , D. Liu 22 , D.B. Liu 9 , H. Liu 6 , H.D. Liu 25 , J. Liu 1 , 3 , J.L. Liu 1 , 3 , M.Y. Liu 18 , R.Y. Liu 7 , S.M. Liu 6 , W. Liu 1 , 3 , Y. Liu 8 , Y.N. Liu 24 , Q. Luo 21 , Y. Luo 9 , H.K. Lv 1 , 3 , B.Q. Ma 29 , L.L. Ma 1 , 3 , X.H. Ma 1 , 3 , J.R. Mao 26 , Z. Min 1 , 3 , W. Mitthumsiri 31 , H.J. Mu 25 , Y.C. Nan 1 , 3 , A. Neronov 23 , K.C.Y. Ng 19 , L.J. Ou 8 , P. Pattarakijwanich 31 , Z.Y. Pei 8 , J.C. Qi 1 , 2 , 3 , M.Y. Qi 1 , 3 , B.Q. Qiao 1 , 3 , J.J. Qin 12 , A. Raza 1 , 2 , 3 , D. Ru ff olo 31 , A. S'aiz 31 , M. Saeed 1 , 2 , 3 , D. Semikoz 23 , L. Shao 11 , O. Shchegolev 10 , 32 , X.D. Sheng 1 , 3 , F.W. Shu 33 , H.C. Song 29 , Yu.V. Stenkin 10 , 32 , V. Stepanov 10 , Y. Su 14 , D.X. Sun 12 , 14 , Q.N. Sun 6 , X.N. Sun 30 , Z.B. Sun 34 , J. Takata 35 , P.H.T. Tam 21 , Q.W. Tang 33 , R. Tang 9 , Z.B. Tang 13 , 12 , W.W. Tian 2 , 20 , L.H. Wan 21 , C. Wang 34 , C.B. Wang 6 , G.W. Wang 12 , H.G. Wang 8 , H.H. Wang 21 , J.C. Wang 26 , Kai Wang 7 , Kai Wang 35 , L.P. Wang 1 , 2 , 3 , L.Y. Wang 1 , 3 , P.H. Wang 6 , R. Wang 22 , W. Wang 21 , X.G. Wang 30 , X.Y. Wang 7 , Y. Wang 6 , Y.D. Wang 1 , 3 , Y.J. Wang 1 , 3 , Z.H. Wang 28 , Z.X. Wang 17 , Zhen Wang 9 , Zheng Wang 1 , 3 , 13 , D.M. Wei 14 , J.J. Wei 14 , Y.J. Wei 1 , 2 , 3 , T. Wen 17 , C.Y. Wu 1 , 3 , H.R. Wu 1 , 3 , Q.W. Wu 35 , S. Wu 1 , 3 , X.F. Wu 14 , Y.S. Wu 12 , S.Q. Xi 1 , 3 , J. Xia 12 , 14 , G.M. Xiang 16 , 2 , D.X. Xiao 11 , G. Xiao 1 , 3 , Y.L. Xin 6 , Y. Xing 16 , D.R. Xiong 26 , Z. Xiong 1 , 2 , 3 , D.L. Xu 9 , R.F. Xu 1 , 2 , 3 , R.X. Xu 29 , W.L. Xu 28 , L. Xue 22 , D.H. Yan 17 , J.Z. Yan 14 , T. Yan 1 , 3 , C.W. Yang 28 , C.Y. Yang 26 , F. Yang 11 , F.F. Yang 1 , 3 , 13 , L.L. Yang 21 , M.J. Yang 1 , 3 , R.Z. Yang 12 , W.X. Yang 8 , Y.H. Yao 1 , 3 , Z.G. Yao 1 , 3 , L.Q. Yin 1 , 3 , N. Yin 22 , X.H. You 1 , 3 , Z.Y. You 1 , 3 , Y.H. Yu 12 , Q. Yuan 14 , H. Yue 1 , 2 , 3 , H.D. Zeng 14 , T.X. Zeng 1 , 3 , 13 , W. Zeng 17 , M. Zha 1 , 3 , B.B. Zhang 7 , F. Zhang 6 , H. Zhang 9 , H.M. Zhang 7 , H.Y. Zhang 17 , J.L. Zhang 20 , Li Zhang 17 , P.F. Zhang 17 , P.P. Zhang 12 , 14 , R. Zhang 14 , S.B. Zhang 2 , 20 , S.R. Zhang 11 , S.S. Zhang 1 , 3 , X. Zhang 7 , X.P. Zhang 1 , 3 , Y.F. Zhang 6 , Yi Zhang 1 , 14 , Yong Zhang 1 , 3 , B. Zhao 6 , J. Zhao 1 , 3 , L. Zhao 13 , 12 , L.Z. Zhao 11 , S.P. Zhao 14 , X.H. Zhao 26 , F. Zheng 34 , W.J. Zhong 7 , B. Zhou 1 , 3 , H. Zhou 9 , J.N. Zhou 16 , M. Zhou 33 , P. Zhou 7 , R. Zhou 28 , X.X. Zhou 1 , 2 , 3 , X.X. Zhou 6 , B.Y. Zhu 12 , 14 , C.G. Zhu 22 , F.R. Zhu 6 , H. Zhu 20 , K.J. Zhu 1 , 2 , 3 , 13 , Y.C. Zou 35 , X. Zuo 1 , 3 , B. Li 7 (The LHAASO Collaboration)  \n- 1 Key Laboratory of Particle Astrophysics & Experimental Physics Division & Computing Center, Institute of High Energy Physics, Chinese Academy of Sciences, 100049 Beijing, China\n- 2 University of Chinese Academy of Sciences, 100049 Beijing, China\n- 3 TIANFU Cosmic Ray Research Center, Chengdu, Sichuan, China\n- 4 Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, 2 Dublin, Ireland\n- 5 Max-Planck-Institut for Nuclear Physics, P.O. Box 103980, 69029 Heidelberg, Germany\n- 6 School of Physical Science and Technology & School of Information Science and Technology, Southwest Jiaotong University, 610031 Chengdu, Sichuan, China\n- 7 School of Astronomy and Space Science, Nanjing University, 210023 Nanjing, Jiangsu, China\n- 8 Center for Astrophysics, Guangzhou University, 510006 Guangzhou, Guangdong, China\n- 9 Tsung-Dao Lee Institute & School of Physics and Astronomy, Shanghai Jiao Tong University, 200240 Shanghai, China\n- 10 Institute for Nuclear Research of Russian Academy of Sciences, 117312 Moscow, Russia\n- 11 Hebei Normal University, 050024 Shijiazhuang, Hebei, China\n- 12 University of Science and Technology of China, 230026 Hefei, Anhui, China\n- 13 State Key Laboratory of Particle Detection and Electronics, China\n- 14\n- Key Laboratory of Dark Matter and Space Astronomy & Key Laboratory of Radio Astronomy, Purple Mountain Observatory,\n- Chinese Academy of Sciences, 210023 Nanjing, Jiangsu, China\n- 15 Research Center for Astronomical Computing, Zhejiang Laboratory, 311121 Hangzhou, Zhejiang, China\n- 16 Key Laboratory for Research in Galaxies and Cosmology, Shanghai Astronomical Observatory, Chinese Academy of Sciences, 200030 Shanghai, China  \n- 17 School of Physics and Astronomy, Yunnan University, 650091 Kunming, Yunnan, China\n- 18 Key Laboratory of Cosmic Rays (Tibet University), Ministry of Education, 850000 Lhasa, Tibet, China\n- 19 Department of Physics, The Chinese University of Hong Kong, Shatin, New Territories, Hong Kong, China\n- 20 Key Laboratory of Radio Astronomy and Technology, National Astronomical Observatories, Chinese Academy of Sciences, 100101 Beijing, China  \n21  \nSchool of Physics and Astronomy (Zhuhai) & School of Physics (Guangzhou) & Sino-French Institute of Nuclear Engineer-  \n- ing and Technology (Zhuhai), Sun Yat-sen University, 519000 Zhuhai & 510275 Guangzhou, Guangdong, China\n- 22 Institute of Frontier and Interdisciplinary Science, Shandong University, 266237 Qingdao, Shandong, China\n- 23 APC, Universit'e Paris Cit'e, CNRS / IN2P3, CEA / IRFU, Observatoire de Paris, 119 75205 Paris, France\n- 24 Department of Engineering Physics & Department of Astronomy, Tsinghua University, 100084 Beijing, China\n- 25 School of Physics and Microelectronics, Zhengzhou University, 450001 Zhengzhou, Henan, China\n- 26 Yunnan Observatories, Chinese Academy of Sciences, 650216 Kunming, Yunnan, China\n- 27 China Center of Advanced Science and Technology, Beijing 100190, China\n- 28 College of Physics, Sichuan University, 610065 Chengdu, Sichuan, China\n- 29 School of Physics, Peking University, 100871 Beijing, China\n- 30 Guangxi Key Laboratory for Relativistic Astrophysics, School of Physical Science and Technology, Guangxi University,\n- 530004 Nanning, Guangxi, China\n- 31 Department of Physics, Faculty of Science, Mahidol University, Bangkok 10400, Thailand\n- 32 Moscow Institute of Physics and Technology, 141700 Moscow, Russia\n- 33 Center for Relativistic Astrophysics and High Energy Physics, School of Physics and Materials Science & Institute of Space Science and Technology, Nanchang University, 330031 Nanchang, Jiangxi, China\n- 34 National Space Science Center, Chinese Academy of Sciences, 100190 Beijing, China\n- 35 School of Physics, Huazhong University of Science and Technology, Wuhan 430074, Hubei, China"}