Two Candidate Pulsar TeV Halos Identified from Property-similarity Studies

Teraelectronvolt halos have been suggested to be a common phenomenon associated with middle-aged pulsars. Based on our recent work on the middle-aged pulsar J0631+1036, which is the only known source positionally coincident with a hard teraelectronvolt γ-ray source and likely powers the latter as a teraelectronvolt halo, we select three candidate teraelectronvolt halos from the first Large High Altitude Air Shower Observatory (LHAASO) catalog of γ-ray sources. The corresponding pulsars, given by the positional coincidences and property similarities, are PSR J1958+2846, PSR J2028+3332, and PSR J1849-0001. We analyze the gigaelectronvolt γ-ray data obtained with the Large Area Telescope on board the Fermi Gamma-ray Space Telescope for the first two pulsars, as the last is γ-ray quiet. We remove the pulsed emissions of the pulsars from the source regions from timing analysis, and determine that there are no residual gigaelectronvolt emissions in the regions as any possible counterparts to the teraelectronvolt sources. Considering the previous observational results for the source regions and comparing the two pulsars to Geminga (and Monogem), the LHAASO-detected teraelectronvolt sources are likely the pulsars’ respective teraelectronvolt halos. We find that the candidate and identified teraelectronvolt halos, including that of PSR J1849-0001, have luminosities at 50 TeV (estimated from the differential fluxes) approximately proportional to the spin-down energy Ė of the pulsars, and the ratios of the former to the latter are ∼6 × 10−4.

In our recent studies of the region of a middle-aged pulsar J0631+1036 (Zheng et al. 2023), which has a high positional coincidence with a TeV source likely detected with the High-Altitude Water Cherenkov (HAWC) Observatory as 3HWC J0631+107 (Albert et al. 2020) and with the Large High Altitude Air Shower Observatory (LHAASO; Cao et al. 2019) as 1LHAASO J0631+1040 (Cao et al. 2023), no GeV γ-ray emission was found at the region in the data obtained with the Large Area Telescope (LAT) onboard the Fermi Gamma-ray Space Telescope (Fermi).The non-detection, which sets a constraint on the existence of a PWN associated with the pulsar, strongly suggests the TeV source as a TeV halo powered by the pulsar, as PWNe detected at TeV energies appear to have soft emissions (H.E. S. S. Collaboration et al. 2018b) and most of them can be detectable at GeV energies with Fermi LAT (see, e.g., Smith et al. 2023;Zheng et al., in preparation).This possibility is supported by the great similarities in the X-ray properties of the pulsar with Geminga and the properties of the TeV source with the TeV halo of Geminga.
Following the studies, we have checked for potential TeV halos among the sources in the current VHE source catalogs, which include that reported by the High Energy Spectroscopy System (HESS) Galactic plane survey (H.E. S. S. Collaboration et al. 2018a).As argued in Zheng et al. (2023), sources were selected if they have hard TeV emission (i.e., possibly having an energy spectrum peaking around ∼25 TeV) and do not have obvious supernova remnant (SNR) or PWN counterparts, but have a high positional coincidence with a pulsar.We found three such sources, HESS J1849−000 (or 1LHAASO J1848−0001u), 1LHAASO J1959+2846u, and 1LHAASO J2028+3352 that are in possible association with PSR J1849−0001, PSR J1958+2846, and PSR J2028+3332 (hereafter J1849, J1958, and J2028, respectively).However, the first one J1849 is an X-ray pulsar that is both radio and γ-ray quiet (Gotthelf et al. 2011;Bogdanov et al. 2019).We thus conducted analysis of the Fermi LAT data similar to that in Zheng et al. (2023) for J1958 and J2028.These two pulsars are γray bright, first discovered from the Fermi LAT observations (Abdo et al. 2009a;Pletsch et al. 2012), while J2028 is radio quiet (Grießmeier et al. 2021 and references therein).
In this paper, we report on the results from our analysis for the two pulsars.The analysis and results are presented below in Section 2. Based on the results and properties of the TeV sources (1LHAASO J1959+2846u and 1LHAASO J2028+3352) and pulsars, we suggest that the TeV sources are likely TeV halos of the pulsars.This discussion, including a summary for the possible properties of the TeV halos, is presented in Section 3.

LAT data and source model
We used the γ-ray data collected with Fermi LAT (Atwood et al. 2009).The two regions of interest (RoIs) were set to have a size of 15 • × 15 • with each centered at the positions of J1958 and J2028 respectively.The events in each RoI in the energy range of 0.1-500 GeV over the time period of from 2008 August 04 15:43:36 (UTC) to 2023 February 16 00:00:00 (UTC; approximately 14.5 yr) were selected from the latest Fermi Pass 8 database.Following the recommendations of the LAT team 1 , we excluded the events with quality flags of 'bad' and zenith angles ≥ 90 • .
1 http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/We used the Fermi LAT 12-year source catalog (4FGL-DR3; Abdollahi et al. 2022) to construct the source models.The sources within 15 • of the pulsar targets in the catalog were included in the source models, and their spectral forms provided by the catalog were used.[Since the Fermi LAT Fourth source catalog was just updated (4FGL-DR4; Ballet et al. 2023) while we were finishing this paper, we checked for possible differences between the two releases.For the RoIs, there are no significant differences, and for the regions of approximately 5 • around the pulsars, there were no differences.]The spectral model files gll iem v07.fits and iso P8R3 SOURCE V3 v1.txt were used as the background Galactic and extragalactic diffuse emissions respectively.

Timing analysis
The two pulsars are bright in the LAT γ-ray band.We selected photons within an 1 • radius circular region centered at each of them for the purpose of constructing their pulse profiles and thus determining the offpulse phase ranges.

PSR J1958+2846
We first tested to fold the photons with the ephemeris given in LAT Gamma-ray Pulsar Timing Models (GPTM) database2 (Ray et al. 2011).No clear pulse profile over the 14.5-yr time period could be obtained.
We then used the method described in Xing et al. (2022), in which the photons during MJD 54682-56540 (the similar time period to that in LAT GPTM database) were folded according to the known ephemeris by using the Fermi TEMPO2 plug-in (Edwards et al. 2006;Hobbs et al. 2006).An empirical template profile was constructed, based on which the times of arrival (TOAs) for as many as possible sets of ∼ 200 day LAT data were generated.The template and TOAs were obtained using the maximum likelihood method described in Ray et al. (2011).We fitted the TOAs with TEMPO2 by adding high-order frequency derivatives and updated the ephemeris.However this ephemeris, whose main parameters f , f 1 , and f 2 are given in Table 1, could only cover the data before MJD 56540.
For the data after MJD 56540, we were able to find two template profiles during MJD 56540-56740 and MJD 56740-59540, while the same method as the above was used.The updated ephemerides are given in Table 1.
We folded the photons in the three time periods according to the ephemerides respectively to construct the  56740-59540  3.44349004(5)  −2.521(1)  1.9(3)  J2028  54682-59991 5.65906764938(2) −0.1555548(9) 0.000353(9) The frequency epoch is MJD 55555 for both pulsars.pulse profiles.There were phase shifts between them, ≃ 0.4375 between the first and second and ≃ 0.125 between the first and the third.After making the corrections for the phase shifts, a pulse profile over nearly 13.3-yr was obtained (left panel of Figure 1).Based on the profile, we defined phase 0.0-0.625 as the onpulse phase range and phase 0.625-1.0as the offpulse phase range.

PSR J2028+3332
For this pulsar, the selected photons could be easily folded according to the ephemeris given in LAT GPTM database.A pulse profile over the whole LAT data time period was obtained, while the ephemeris used is given in Table 1.Based on the pulse profile, which is shown in the right panel of Figure 1, we defined phase 0.0-0.5625and 0.9375-1.0 as the onpulse phase ranges and phase 0.5625-0.9375as the offpulse phase range.

Likelihood and spectral analysis 2.3.1. Onpulse data
We performed standard binned likelihood analysis to the 0.1-500 GeV LAT data during the onpulse phase ranges of two pulsars determined above.The source models described in Section 2.1 were used.The sources within 5 • from each of the pulsars were set to have free spectral parameters, while for the other sources in the 10 1 10 3 10 5 10 7 10 9 Energy (M(v) source models, their spectral parameters were fixed at the values given in the catalog.The background normalizations were set as the free parameters.
For fitting the emissions of the pulsars, we used a PLSuperExpCutoff4 (PLSEC) model (Abdollahi et al. 2022)), where Γ is the photon index, d the local curvature at E 0 , and b a measure of the shape of the exponential cutoff and was fixed at 2/3 (a characteristic value used for the γ-ray pulsars in 4FGL-DR3).From the likelihood analysis, we obtained Γ = 2.12 ± 0.01 and d = 0.63 ± 0.03 (Γ = 2.140 ± 0.004 and d = 1.109 ± 0.007) for J1958 (J2028).These values are close to those given in 4FGL-DR3.These results, as well as the corresponding TS values, are provided in Table 2.
We also obtained the flux measurements of the γ-ray emissions of J1958 and J2028 from their onpulse data.The number of the energy bins was set to be 10 by evenly dividing the energy range of from 0.1 to 500 GeV in logarithm.Fluxes were obtained from the maximum likelihood analysis of the data in each energy bin.In this analysis, the spectral normalizations of the sources within 5 • from each pulsar were set as free parameters, and all other parameters of the sources in the source models were fixed at the values obtained from the above binned likelihood analysis.When the TS value of a spectral data point was less than 4, we used the 95% flux upper limit instead.The obtained spectra are shown in Figure 2.

Offpulse data
We performed standard binned likelihood analysis to the 0.1-500 GeV LAT data during the offpulse phase ranges of the two pulsars.The parameter setups were the same as those in the analysis of the onpulse data (Section 2.3.1).We assumed a power law for any emission at the position of each pulsar, dN/dE = N 0 (E/1 GeV) −Γ .No emissions were detected.When we fixed Γ = 2, the resulting TS values were ∼0 (Table 2).
To verify the non-detection results in the offpuse data, we calculated 0.1-500 GeV TS maps for the two pulsar regions.PWNe or TeV halos may show extended weak emission at GeV energies (e.g., Di Mauro et al. 2021).However, as shown by the TS maps (Figure 3), no residual emissions are seen in each of the pulsar regions.We also tested to assume a uniform-disk model (with a radius value varied upto 1 • ) centered at each of the pulsars and perform binned likelihood analysis.No extended emissions were detected, which confirms the visual inspection results of the TS maps.

DISCUSSION
By selecting TeV sources sharing the similarities with the TeV halo of Geminga, i.e., having hard emission and being positionally coincident with only a middleaged pulsar, we have found 1LHAASO J1959+2846u and 1LHAASO J2028+3352 as possible TeV halos.The corresponding pulsars J1958 and J2028 both have bright γ-ray emissions in the Fermi LAT GeV band.We conducted timing analysis of the LAT data to remove the pulsed emissions of the pulsars from the source regions, and found that no GeV γ-ray emissions were detected in each of the regions.The non-detections largely constrain the existence of PWNe associated with the pulsars, as we note that the spectral upper limits of ∼10 −13 erg s −1 cm −2 (Figure 2) set a luminosity limit on a PWN at a distance of 2 kpc down to ∼10 32 erg s −1 , lower than those of the known PWNe or candidate PWNe (e.g., Ackermann et al. 2011).Combining the non-detections of any GeV emissions with the hard emission property of the two TeV sources, it is likely that they are TeV halos of the two pulsars.
Below we discuss the different detection results for the source regions of the two pulsars in Section 3.1 & 3.2, which help strengthen their likely associations with the two TeV sources.In Section 3.3, we show the detailed similarities of the two pulsars and their presumed TeV halos with the other pulsar TeV halos, in which we include the pulsar J1849.

PSR J1958+2846
As shown in the left panel of Figure 3, in addition to 1LHAASO J1959+2846u that is detected with an extension size of 0. • 29 in the energy range of ≥25 TeV and has not been detected in 1-25 TeV, the Milagro gamma-ray observatory (MGRO) had a 4σ detection at the position of J1958 in 1-100 TeV (Abdo et al. 2009b), and HAWC and LHAASO respectively reported two nearby sources 3HWC J1957+291 (Albert et al. 2020) and LHAASO J1956+2845 (Cao et al. 2021).The MGRO detection is consistent with the LHAASO results (Figure 2), but the latter two have separation distances of 0. • 54 and 0. • 74 respectively.We suspect that the latter two are individual sources, not likely in association, which can be verified by near-future updating observations.In addition, the MAGIC telescopes searched for the PWN of J1958 (Fernandez-Barral et al. 2017) and VERITAS telescopes searched for pulsed γray emission of J1958 (Archer et al. 2019), both at the energy of 100 GeV, but no emission was detected (see Figure 2 for their upper limits).Finally, we note that there are two candidate SNRs identified from an optical Hα survey (Sabin et al. 2013), which positionally overlap with the pulsar or 1LHAASO J1959+2846u.However very limited information is available for these two candidate SNRs (e.g., Green 2019).We suggest they are not likely the counterpart to the TeV source given the non-detection of any GeV γ-ray emissions in the offpulse data.

PSR J2028+3332
The field of J2028 is rather clean, as 1LHAASO J2028+3352, having a positional error circle of 0. • 86 and an extension size of 1. • 7, is the only known VHE source reported (Figure 3).The source has similar hard emission, detected in 25-100 TeV but not in 1-25 TeV.There was a reported upper limit on the pulsar from MGRO observations (Abdo et al. 2009b), and the upper limit is lower than that of the LHAASO measurements (Figure 2).This inconsistency is not understood, probably due to the simple assumption of a point source with a fixed power-law spectral form in the MGRO data analysis.We note that there are two Fermi LAT sources also detected in the error circle, 4FGL J2025.3+3341 and 4FGL J2027.0+3343.Given in the LAT catalog, the first is identified as a blazar, thus not likely the counterpart to 1LHAASO J2028+3352, and the second is an unknown-type source with its emission described with a power law and is faint (TS∼40).The nature of this second GeV source remains to be investigated, helping clarify whether its presence is purely because of a coincidence due to the large error circle of the TeV source.
In addition, since J2028 is radio quiet, there is no distance estimation for the source.We obtained its phase-averaged γ-ray flux, ≃5.1×10 −11 erg cm −2 s −1 in 0.1-500 GeV.Given its spin-down energy Ė = 3.5 × 10 34 erg s −1 , the distance should be ≤2.4 kpc.If we further consider an efficiency of 10% for converting Ė to the γ-ray emission (Smith et al. 2023), the distance would be ≃0.76 kpc.Therefore the pulsar is likely close, which can explain the large extension size of 1LHAASO J2028+3352 if we consider the latter as the putative TeV halo.

Property comparison
We collected the properties of these two pulsars as well as Geminga, Monogem, and PSR J0631+1036, which are given in Table 3.In addition, although J1849 is γ-ray quiet and we did not conduct analysis for it, we also list it in the table because the corresponding TeV source 1LHAASO J1848-0001u shows similar hard emission.These pulsars have spin periods of 0.2-0.4s except J1849 whose value is only ≃0.04 s The characteristic ages spread from 20 to 600 kyr.They are mostly X-ray faint and have faint (or non-detectable) X-ray PWNe also except J1849 (considering its assumed distance of ∼7 kpc; Gotthelf et al. 2011).Since they are all detected by LHAASO in 25-100 TeV, we also list their fluxes at 50 TeV (Cao et al. 2023) in Table 3.
As noted in Zheng et al. (2023), the putative TeV halo of PSR J0631+1036 has a flux proportional to Ė of the pulsar when comparing to the corresponding values of Geminga (where the HAWC measurements were used).We thus plot the luminosity values at 50 TeV of the (candidate) TeV halos given by the LHAASO detections versus Ė of the pulsars in the left panel of Figure 4, where for the source distances (PSR J0631+1036 and J1958) given by the radio timing we assign 30% uncertainties, for J2028 the upper limit is 2.4 kpc and we also assign a 30% uncertainty on the possible distance of 0.76 kpc, and for PSR J1849 its suggested distance of 7 kpc is assumed to have a 30% uncertainty as well.As can be seen, there is a possible correlation between the 50 TeV luminosities (L 50TeV ) and Ė, although the distances suffer large uncertainties.We test to use the four pulsars with relatively certain properties, Geminga, Monogem, PSR J0631+1036, and J1958 and find a correlation of L 50TeV ≃ (2.5 ± 1.7) Ė0.90±0.01 , where the Markov Chain Monte Carlo (MCMC) code emcee (Foreman-Mackey et al. 2013) was used.The index value 0.90 is close to 1, which does suggest the fluxes of the TeV halos are proportional to Ė.It is also interesting to note that J1849 is along this correlation line.
Another way to show this possible correlation is to plot the ratios of L 50TeV to Ė.In the right panel of Figure 4, we plot the ratios versus ages (τ ) of the pulsars.Despite the large ranges mainly caused by the distance uncertainties, we test to fit the data points of the four pulsars (given above) with a constant c and a function of aτ b kyr (where τ kyr is τ in units of kyr), and obtain c = (5.9± 1.1) × 10 −4 and a = (1.2+2.3 −0.8 ) × 10 −3 , b = −0.18± 0.26, where for the latter the MCMC method was used.As for the best-fits, the corresponding values of χ 2 over degree of freedom (DoF) are 0.4/3 and 1.5/2 respectively, the latter does not provide a better fit and thus the brightnesses of TeV halos are likely not in any close relation with ages of pulsars.We conclude that the ratios are approximately in a range of 4.8-7.0× 10 −4 .
From the comparison analysis, we may summarize the general properties of pulsar TeV halos and thus the method to identify them.The TeV sources should have hard emission with spectrum peaks 25 TeV, which can be used to be differentiated among PWNe (Zheng et al. 2023).The source regions are rather clean, without any potential SNR or PWN counterparts known at energies of X-rays or γ-rays, while middle-aged pulsars can be found positionally coincident with the source regions.The pulsars are generally γ-ray bright and X-ray faint, and may have faint X-ray PWNe.Then if considering the relationship built based on the limited cases (i.e., Figure 4), the TeV halos would have fluxes proportional

Figure 1 .
Figure 1.Phase-connected pulse profile (top) and two-dimensional phaseogram (bottom) constructed for J1958 (left) and J2028 (right).Two spin cycles are shown for clarity.The onpulse and offpulse phase ranges are marked by dashed lines.

Figure 3 .
Figure 3. TS maps of the J1958 (left) and J2028 (right) regions in 0.1-500 GeV calculated for the offpulse phase ranges of the two pulsars.In the left panel, the positional error circle and extension of 1LHASSO J1959+2846u (Cao et al. 2023) are marked by white solid and dashed circles respectively, which are coincident with J1958's position (red diamond).Also shown are the positional error circles of 3HWC J1957+291 (pink dashed; Albert et al. 2020) and LHAASO J1956+2845 (red dashed; Cao et al. 2021).In addition, two candidate SNRs, G65.8−0.5 and G66.0−0.0 (marked by yellow dashed and cyan dashed ellipses respectively), reported in Sabin et al. (2013) are located in the region as well.In the right panel, the positional error circle and extension of 1HASSO J2028+3352 are marked as white solid and dashed circles respectively.The pulsar J2028, along with two Fermi LAT sources, are located in the error circle.

Table 1 .
Timing parameters derived for PSR J1958+2846 and PSR J2028+3332 (Fernandez-Barral et al. 2017)ots) of J1958 (left) and J2028 (right) in 0.1-500 GeV obtained from their onpulse data and the corresponding best-fit PLSEC models (dashed curves).From the offpulse data, the upper limits derived by assuming a power law with Γ = 2 are shown as the red lines.The fluxes or flux upper limits from different VHE observations of the source regions are also shown, which include in the left panel, 1LHAASO J1959+2846u measured with the Kilometer Squared Array (KM2A) and Water Cherenkov Detector Array (WCDA;Cao et al. 2023), MGRO J1958+2848(Abdo et al. 2009b), VERITAS upper limits on J1958(Archer et al. 2019), and MAGIC upper limits on the PWN of J1958(Fernandez-Barral et al. 2017), and in the right panel 1LHAASO J2028+3352 and MGRO upper limit on J2028.

Table 2 .
Binned Likelihood Analysis Results Fluxes in energy 0.1-500 GeV are in units of photon s −1 cm −2 † Fluxes are 95% upper limits derived with Γ = 2 assumed *