Possible GeV gamma-ray emission from the pulsar wind nebula in CTA 1

We report a detection of GeV $\gamma$-ray emission potentially originating from the pulsar wind nebula in CTA 1 by analyzing about 15 yr of Fermi Large Area Telescope data. By selecting an energy range from 50 GeV to 1 TeV to remove contamination from the $\gamma$-ray pulsar PSR J0007+7303, we have discovered an extended $\gamma$-ray source with a TS value of $\sim$ 44.94 in the region of CTA 1. The obtained flux is measured to be 6.71 $\pm$ 2.60 $\times$ $10^{-12}$ erg $\mathrm{cm}^{-2}$ $\mathrm{s}^{-1}$ with a spectral index of 1.61 $\pm$ 0.36, which allows for a smooth connection with the flux in the TeV band. CTA 1 is also considered to be associated with 1LHAASO J0007+7303u, which is an Ultra-High-Energy source listed in the recently published catalog of the Large High Altitude Air Shower Observatory. We assume that the radiation originates from the pulsar wind nebula and that its multi-wavelength spectral energy distribution can be explained well with a time-dependent one-zone model.


INTRODUCTION
The detection and characterization of Ultra-High-Energy (UHE; > 0.1 PeV) -ray sources have significant implications for our understanding of the origin and acceleration mechanisms of cosmic rays.The Large High Altitude Air Shower Observatory (LHAASO) is a revolutionary project that aims to explore the high-energy universe with remarkable sensitivity (Cao et al. 2019).The LHAASO collaboration has made a breakthrough by reporting the detection of 12 UHE -ray sources with statistical significance exceeding 7 and a maximum energy of up to 1.4 PeV (Cao et al. 2021), sparking further exploration of UHE gamma-ray sources.Recently, the first catalog of -ray sources detected by LHAASO was published (Cao et al. 2023), containing a total of 90 -ray sources with an extended size smaller than 2 • and a significance of detection at a > 5 significance level.Among these sources, 43 exhibit UHE emission at a > 4 significance level.1LHAASO J0007+7303u is a TeV source in the first LHAASO catalog, detected by both the Water Cherenkov Detector Array (WCDA) and the 1 km 2 array (KM2A) as a UHE source with a significance of TS 100 = 171.6.This source is considered to be associated with the shell-type supernova remnant (SNR) CTA 1 (approximately 0. • 12 away) (Cao et al. 2023).
CTA 1 (G119.5+10.2) was first proposed as a composite SNR through radio observations by Harris & Roberts (1960), and it was first detected in the X-ray band by Seward et al. (1995).In the center of the region, there was a faint compact point source called RX J0007.0+7303, which was considered to be a pulsar candidate powering the nebula (Seward et al. 1995;Slane et al. 1997).Deeper observations of RX J0007.0+7303 were conducted using    −, ASCA satellites (Slane et al. 2004), and the ℎ observatory (Halpern et al. 2004), revealing that it is a point source embedded in a nebula with a jet-like feature.Additionally, the source 3EG ★ E-mail: fangjun@ynu.edu.cnJ0010+7309, which is associated with RX J0007.0+7303, has been reported as a radio-quiet -ray pulsar candidate (Mattox et al. 1996).Although no -ray pulsations had been found, the constant -ray flux, the hard power-law spectrum, and the steepening of the spectrum above approximately 2 GeV, confirmed that the source was a -ray pulsar (Brazier et al. 1998).The Fermi Large Area Telescope (Fermi-LAT) satellite then discovered a -ray pulsar (PSR J0007+7303) in CTA 1, which had a period of 316.86 ms and a period derivative of 3.614 ×10 −13 s s −1 , making it the first pulsar discovered by -rays (Abdo et al. 2008).Subsequently, X-ray pulsations were detected (Caraveo et al. 2010;Lin et al. 2010).However, neither radio nor optical counterparts of RX J0007.0+7303 were detected (Mignani et al. 2013).
The TeV -ray emission from CTA 1 was discovered by Aliu et al. (2013) using the     observatory.The extended source, VER J0006+729, was detected at a 6.5 significance level, with a twodimensional Gaussian of 0.30 • × 0.24 • at a distance of 5 ′ from PSR J0007+7303.Previously, Abdo et al. (2012) reported potentially GeV extended emission after analyzing 2 yr of Fermi-LAT data in the offpulse emission of PSR J0007.0+7303, but they only provided upper limits.On the other hand, Li et al. (2016) using 7 yr of Fermi-LAT data and updated response functions, did not detect GeV emission compatible with CTA 1 or the highly energetic (> 100 GeV) pulsar wind nebula.Instead, they provided an upper limit.In a more recent study, Ackermann et al. (2018) reported GeV emission with a radius of approximately 0.98 • , spatially consistent with the radio emission from CTA 1.However, there was a discrepancy in the flux normalization between the GeV and TeV spectra.In this work, we reanalyze the GeV emission from CTA 1 using about 15 yr of Fermi-LAT data.
First, Section 1 serves as the introduction.Then, in Section 2, we present a detailed analysis of the Fermi-LAT data.Next, Section 3 contains specific discussions related to the analysis.Lastly, Section 4 offers a summary of the paper.

Data Selection
We analyzed about 15 yr (from 2008 August 4 to 2023 May 24) of data toward CTA 1 in a 10 • radius centered on PSR J0007+7303 (logged as 4FGL J0007.0+7303 in 4FGL-DR3; Abdollahi et al. 2022) with the available software Fermipy (Wood et al. 2017).The newest response functions P8R3_SOURCE_V3 and the Pass 8 events with "source" event class (evtype = 3 and evclass = 128) was selected, while the expression of "(DATA_QUAL > 0)&&(LAT_CONFIG==1)" was used to filter a good time interval.The maximum zenith angle was fixed to 90 • to limit contamination of -rays from the Earth's limb.To search for potential -rays from the PWN, we chose an energy band of 50 GeV−1 TeV to suppress the contamination of PSR J0007+7303.
The sources in the Data Release 3 of the fourth Fermi-LAT source catalog (4FGL-DR3; Abdollahi et al. 2022) were considered in the analysis.The isotropic -ray background emission (iso_P8R3_SOURCE_V3_v1) and the diffuse Galactic interstellar emission (gll_iem_v07) were used to construct the background model.In the analysis that follows, the normalization and spectral parameters are free for the sources within 3 • in the background model.The normalization parameters of the isotropic and Galactic components are also left free.

Spatial Analysis
To investigate the possibility of detecting the radiation originating from PWN in CTA 1, we first generated a test statistic (TS) map centred on 4FGL J0007.0+7303(R.A. = 1.• 768, Decl.= 73.• 051; Abdollahi et al. 2022).In order to minimize contamination from PSR J0007+7303, we set an energy threshold of 50 GeV.We performed a fitting procedure for the background sources and subtracted 4FGL J0007.0+7303 to create the excess TS map in the energy range of 50 GeV to 1 TeV.Here, 4FGL J0007.0+7303, which corresponds to PSR J0007+7303, is a point source listed in the 4FGL catalog with a PLSuperExpCutoff4 spectrum (Abdollahi et al. 2022).
In Figure 1, a significant excess of GeV -ray radiation is evident at the observed positions of the other detectors.Based on LHAASO's observations (Cao et al. 2023), 1LHAASO J0007+7303u was detected by KM2A as an extended source located at R.A. = 1 • .91 and Decl.= 73 • .07 with  39 ∼ 0 • .17and by WCDA as a point source located at R.A. = 1 • .48 and Decl.= 73 • .15 with a 95% error radius of ∼ 0 • .22.The locations can be seen in Figure 1, denoted by cyan and green circles.Meanwhile, magenta ellipse of 0. • 30 × 0. • 24 indicates the VERITAS's observation (Aliu et al. 2013).To determine the best morphological model describing the excess radiation, we utilized    to calculate the best-fit location based on the initial position with the maximum TS value.The best-fit location for the excess radiation was R.A. = 1 • .76 and Decl.= 73 • .02,with a 1 error of 0 • .06,marked in Figure 1 as a black cross and labelled as SrcA for convenience.
We first tested the power-law spectrum point-source model, with a spectral index of 2.0, at the best-fit position.Spatial extension models such as the uniform disk and the 2D Gaussian models were also applied to judge the best-fit model.We applied both models to the extension analysis of SrcA using the .programme in   (Wood et al. 2017).As a result, SrcA obtained a   ext value of 27.61 in the 2D Gaussian template and 28.33 in the uniform disk template.Here, the   ext value is defined as   ext = 2 (ln L ext − ln L ps ), where the ln L ext and ln L ps represent the maximum likelihood values of the extended and point-like model, re-spectively (Lande et al. 2012).Therefore, we will chose the extended model to describe SrcA as   ext ⩾ 16 (a significance level of 4 ).
We calculated the delta log-likelihood value, defined as ln L ext − ln L ps , for various radii in order to determine the optimal extended model.The distribution of delta log-likelihood values for the various templates is shown in Figure 2. Comparing the two templates, we found that the disk template had a higher   ext value and less uncertainty.As a result, we selected the disk template to analyze the excess -ray radiation in this region.Using the best-fit positions and extension of the different extended models, we performed a likelihood analysis again.During the fitting process, the spectrum type of power-law / =  0 (/ 0 ) − was used.The results of the fitting, which including the point-source model's, are conveniently presented in Table 1.Acero et al. (2013) collected a large number of PWNe detected by Fermi-LAT and studied the correlations between different energy bands, establishing new constraints.Using 45 months of Fermi-LAT data, they derived an upper limit of 6.9 × 10 −12 erg cm −2 s −1 on the flux of CTA 1 between the 10−316 GeV band.Li et al. (2016) performed a detailed analysis of the radiation properties of PSR J0007+7303 during the off-peak and the on-peak phase intervals.However, no point-like or extended -ray emission between 10 and 300 GeV was detected during the off-peak phase of PSR J0007+7303, resulting an upper limit of 6.5×10 −12 erg cm −2 s −1 instead.According to The Third Fermi Large Area Telescope Catalog of Gamma-ray Pulsars (Smith et al. 2023), the spectral energy distribution (SED) of PSR J0007+7303 peaks at 2.44 ± 0.03 GeV.We assumed that the contribution of pulsar to GeV flux can be described by a PLEC4 model (For details, see Abdollahi et al. 2022).As shown in Figure 3, the flux of the pulsar exponentially decays in the GeV band.When selecting the energy range from 50 GeV to 1 TeV, an extended source in the region overlapping observations from other telescopes exhibits GeV -ray emission with a hard spectral index of ∼ 1.61.This indicates that the high-energy radiation may originate from the PWN in CTA 1.This suggests that above 50 GeV, the pulsar no longer contributes significantly to the -ray emission.Therefore, to obtain higher photon statistics, the utilization of pulsar gating is unnecessary.

Spectral Analysis
Based on the uniform disk model with an extension of  68 = 0.41 • ±0.03 • , we obtained a global fitted flux of (6.71 ± 2.60) × 10 −12 erg cm −2 s −1 with a spectral index of 1.61 ± 0.36 and a TS value of ∼ 44.94 calculated by the likelihood analysis in the energy range of 50 GeV−1 TeV.We divided the data from 50 GeV to 1 TeV band into 3 logarithmically equal energy bins and performed the same likelihood analysis for each bin.The results are displayed as blue dots in Figure 3, with specific values listed in Table 2.In the Fermi High-Latitude Extended Sources Catalog (Ackermann et al. 2018), a significantly extended (∼ 0. • 98) -ray source (FHES J0006.7+7314) was reported, which is thought to be associated with the CTA 1.However, there was a mismatch in flux normalization observed between the two spectra.In our work, the derived flux does not exceed the upper limit defined in previous work (Acero et al. 2013;Li et al. 2016).The obtained hard spectral index value of ∼ 1.61 agrees with the typical hard spectral index observed in PWNe by Fermi-LAT (Acero et al. 2013).Furthermore, as shown in Figure 4, the GeV SED could smoothly connect to the TeV one.(Cao et al. 2023).Magenta ellipse of 0. • 30 × 0. • 24 indicates the VERITAS's observation (Aliu et al. 2013).The yellow circle indicates the best-fit extension radius of the GeV excess under the uniform disk model.The white contours correspond to the radio observation of the SNR shell (Pineault et al. 1997).The map was smoothed using a Gaussian function with a kernel radius of 0.

Systematic Uncertainties Analysis
The systematic errors in the Fermi-LAT data are considered from two main origins: imperfect modelling of the Galactic diffuse emission and uncertainties regarding the calibration of the effective area.
For the first aspect, we estimated the systematic error by repeated analyses using the Galactic diffuse model with artificially altered normalization of ±6% (Abdo et al. 2010(Abdo et al. , 2013)), and for the second aspect, we used the bracketing Aeff method 1 (Ackermann et al. 2012).The specific systematic error values we listed in

DISCUSSION
The hadronic origin of VER J0006+729 has been extensively investigated.According to the findings of Martín et al. (2016), the mass of the molecular clouds is insufficient to explain the observed TeV emission.Additionally, Aliu et al. (2013) found that the much smaller extent of VER J0006+729 compared to the supernova remnant (SNR) contradicts the SNR shell hypothesis.Therefore, we favour the leptonic PWN origin scenario for VER J0006+729.
The SED of PWNe exhibit a bimodal structure, with the lowenergy component generated through synchrotron from radio to Xray band, while the high-energy component in the -ray band is produced by inverse Compton scattering off the soft photon field.10 12 10 9 10 6 10 3 10 0 10 3 10 6 10 9 Photon Energy (MeV)   (Aliu et al. 2013;Giacani et al. 2013) and the X-ray data (Lin et al. 2012) are shown in this figure as red and orange dots, respectively.The fluxes in GeV -rays with the statistic error (this work) and in TeV -rays with     (Aliu et al. 2013) and LHAASO (Cao et al. 2023) are also shown in this figure.Aliu et al. (2013) used the 1.4 GHz image from Pineault et al. (1997) to estimate the radio flux within a 20 arcmin radius around the pulsar, and here we adopt this flux as an upper limit for the radio emission.Giacani et al. (2013) also presented new high angular resolution and high sensitivity radio observations toward PSR J0007+7303 at 1.5 GHz with the Jansky Very Large Array, and set an upper limit considering a size for the nebula of 20 arcmin in radius.In the X-ray band, Lin et al. (2012) detected a ∼10 arcmin extended feature with  that may correspond to the bow shock of the nebula.In the TeV -ray band, the TeV source VER J0006+729 associated with CTA 1 was detected by     (Aliu et al. 2013), which has a photon spectrum described by a power-law spectrum / =  0 (/3  ) − with a differential spectral index of  = 2.2±0.2  ±0.3  and the nor- CTA 1 is also thought to be associated with a UHE source 1LHAASO J0007+7303u catalogued in the first LHAASO catalog (Cao et al. 2023), with a flux of (5.01±1.11)×10−13   −1  −2  −1 detected by WCDA at 3 TeV and a flux of (3.41 ± 0.27) × 10 −16   −1  −2  −1 detected by KM2A at 50 TeV.The aforementioned observational data, as well as the Fermi-LAT data analyzed in this work, are plotted in Fig. 4.
We considerd a time-dependent one-zone model for the radiative evolution of a PWN (For details, see Torres et al. 2014;H. E. S. S. Collaboration et al. 2018).The spin-down power of an energetic pul-sar is continuously transferred to high-energy particles and magnetic field in the PWN, and the spin-down power L(t) evolves in time as (Gaensler & Slane 2006) where  0 () is the initial luminosity, the breaking index n is assumed to be 3.0 here.And the initial spin-down time-scale  0 of the pulsar is (Gaensler & Slane 2006) where   = /2  is the characteristic age can be derived from the period and the period derivative.+ (, ). (3) We assumed that the injection rate of particles (, ) follow a broken power-law distribution, i.e., (, ) =  0 () in which   is the break Lorentz factor, and  1 and  2 are the spectral indices.
The expansion behaviour of the PWN in the SNR environment is determined by the age of the system t, the initial spin-down timescale  0 , and the reverse-shock interaction time    , which as follows (Mayer et al. 2012;H. E. S. S. Collaboration et al. 2018): The magnetic field strength evolves with time can be calculated by (Tanaka & Takahara 2010;Martín et al. 2012) The young pulsar PSR J0007+7303 in CTA 1 has a period of ∼316 ms, a spin-down power of ∼4.5 × 10 35 erg s −1 and a characteristic age of   = 13 kyr (Abdo et al. 2012), which is energetic to power the nebula (Slane et al. 2004).We adopt a distance of 1.4 kpc in this paper (Pineault et al. 1993).The interstellar radiation field including CMB, the near-infrared (NIR) background with  NIR = 25.0K and  NIR = 0.3 eV cm −3 and far-infrared (FIR) background with  FIR = 3000 K and  FIR = 0.6 eV cm −3 were calculated using GALPROP code (Porter et al. 2006) by Zhu et al. (2018).We assumed that the nebula has an age of 10 kyr, and that the corresponding  0 and  0 are 3900 yr and 2.78 × 10 36 erg s −1 .For the particle spectrum, we assumed a broken power-law with  1 = 1.8,  2 = 2.0 and   = 8 × 10 4 , and the resulting SED at t age is shown in Fig. 4. With the magnetic energy ratio  = 0.4, the magnetic field strength in the nebula is 3.25 G.
As the PWNe evolve continuously inside the SNR, the expansion of the PWNe results in a gradual decrease in magnetic field strength (Torres et al. 2014).In the model proposed by Gelfand et al. (2009), as the PWNe evolve to about 10 kyr, the magnetic field strength decreases to a few G.In this work, we obtained a relatively low magnetic field strength of 3.25 G.Similar low values have also been adopted for several PWNe, such as HESS J1826−130 (Burgess et al. 2022) As shown in Fig. 4, with these reasonable parameters, the observed SED can be well reproduced.In addition, we present in Fig. 5 the evolution of the SED at different ages from 1 to 14 kyr, which is able to fit the observation data points well when the nebula age is 10 kyr.Finally, the fitting parameters we used for the model are summarized in Table 3.
It is well known that PWNe are a prominent population of TeV sources in the Milky Way2 as observed by various instruments in recent years (Park & VERITAS Collaboration 2015;H. E. S. S. Collaboration et al. 2018;Albert et al. 2020).In particular, the detection of photons with energies up to ∼ 1 PeV from the Crab Nebula implies the presence of PeV electron accelerators (Pevatrons) (Lhaaso Collaboration et al. 2021).Observations from LHAASO have revealed numerous hundred TeV sources in the Milky Way (Cao et al. 2021(Cao et al. , 2023)), a significant fraction of which can be attributed to PWNe (Wu et al. 2023;Tian et al. 2023;Xia et al. 2023) and are plausible candidates for PeVatrons (de Oña Wilhelmi et al. 2022).The detection of UHE -rays from CTA 1, suggests its potential as a PeVatron.Future observations in the UHE band can provide further insights into the energy distribution of particles within CTA 1.

SUMMARY
In this paper, we analyzed the GeV -ray emission toward CTA 1 using about 15 yr of Fermi-LAT data.Under the hypothesis of a uniform disk model with an extension of  68 = 0.41 • ± 0.03 • , we obtained a global flux of (6.71 ± 2.60) × 10 −12 erg cm −2 s −1 with a spectral index of 1.61 ± 0.36 and a TS value of ∼ 44.94 in the 50 GeV−1 TeV energy band.The GeV emission could potentially be from the PWN in the CTA 1, wherein the source is powered by PSR J0007+7303 and emits -rays through inverse Compton scattering.The SED could be reasonably reproduced by a time-dependent onezone model with a broken power-law spectrum.CTA 1 has also been detected by     and LHAASO instruments.The UHE -ray emission detected by LHAASO implies the potential for particle acceleration to the PeV range.Further in-depth observations and analyses are needed to clarify the maximum energy of the particles in the PWN.

Figure 1 .
Figure 1.The excess TS map in the 50 GeV−1 TeV energy band centred on 4FGL J0007.0+7303 with each pixel of 0. • 01.All the background sources are removed in the map.The red, black and yellow crosses indicate the position of 4FGL J0007.0+7303, the best-fit position of the point-source model from this work and other 4FGL sources, respectively.The cyan circle (a radius of ∼ 0 • .17)shows the 95% statistic upper limits of the extension detected by WCDA, and the green circle (a radius of ∼ 0 • .22)shows the 39% containment radius of the 2D-Gaussian model detected by KM2A(Cao et al. 2023).Magenta ellipse of 0. • 30 × 0. • 24 indicates the VERITAS's observation(Aliu et al. 2013).The yellow circle indicates the best-fit extension radius of the GeV excess under the uniform disk model.The white contours correspond to the radio observation of the SNR shell(Pineault et al. 1997).The map was smoothed using a Gaussian function with a kernel radius of 0. • 06.

Figure 4 .
Figure 4. Broad-band SED of the PWN model at t age = 10 kyr with observed data.The purple dotted line represents the synchrotron radiation, and the red, blue, orange and green dotted lines are for the inverse Compton scatterings off the synchrotron photons, CMB and FIR and NIR background, respectively.The total SED is shown by the black solid line.The radio data(Aliu et al. 2013;Giacani et al. 2013) and the X-ray data(Lin et al. 2012) are shown in this figure as red and orange dots, respectively.The fluxes in GeV -rays with the statistic error (this work) and in TeV -rays with    (Aliu et al. 2013) and LHAASO(Cao et al. 2023) are also shown in this figure.

Table 1 .
(Aliu et al. 2013)of CTA 1.The black dashed line indicates the energy spectrum of the PSR J0007+7303 described by a PLEC4 model.The blue dots with 1 statistic error and the systematic error in brown represents the results of the Fermi-LAT data analysis in this work (between the 50 GeV−1 TeV energy band), while the best-fit spectrum is shown by the blue lines and the TS value in each energy bin is represented by the gray-shaded region.The red line indicates the flux upper limit in the 10−300 GeV energy range(Li et al. 2016).The cyan and green dots with statistical uncertainties represent the results of LHAASO's observations detected by WCDA and KM2A, respectively(Cao et al. 2023).The magenta points from    's observation indicate the spectrum of VER J0006+729(Aliu et al. 2013).Spatial Distribution Analysis of SrcA Using Different Models in the 50-1000 GeV Energy Range Figure 2. Extension likelihood curves for uniform disk model (left), and 2D Gaussian (right).The best-fit for the extension is represented by the blue vertical line, and the 1 uncertainty is shown as a purple band.

Table 2 .
The Energy Flux of SrcA with Fermi-LAT Data in the 50 GeV−1 TeV Energy Band.The energy fluxes of SrcA with 1 statistic error in each energy bin are given, and the second one indicates systematic error.

Table 3 .
Fitting Parameters of the PWN model