High-energy studies of the 3HWC J1954+286 region: likely Gamma-ray detection of the supernova remnant G65.1+0.6

We carry out high-energy studies of the region of the Galactic TeV source 3HWC J1954+286, whose location coincides with those of PSR~J1954+2836 and supernova remnant (SNR) G65.1+0.6. Analyzing the GeV $\gamma$-ray data obtained with the Large Area Telescope (LAT) onboard {\it the Fermi Gamma-ray Space Telescope}, we are able to separate the pulsar's emission from that of the region. Excess power-law--like emission of a $\sim 6\sigma$ significance level at the region is found, for which we explain as arising from the SNR~G65.1+0.6. Given the low-significance detection, either a hadronic or a leptonic model can provide a fit to the power-law spectrum. Considering the properties of the pulsar and the SNR, we discuss the possible origin of the TeV source, and suggest that it is likely the TeV halo associated with the pulsar.


INTRODUCTION
Recent very-high energy (VHE; >100 GeV) surveys have shown that our Galaxy may contain many TeV sources, as approximately 100 such sources have thus far been detected. Notably there are 78 sources reported by the High Energy Spectroscopy System (HESS) Galactic plane survey (H. E. S. S. Collaboration et al. 2018a) to be located along the latitudes of ≤3 • within longitudes from 250 • to 65 • , and there are 65 sources (a few are common detections with the HESS) reported by the High-Altitude Water Cherenkov (HAWC) Observatory in the third HAWC catalog (3HWC; Albert et al. 2020). While sources like the supernova remnants (SNRs) and pulsar wind nebulae (PWNe) are considered as primary VHE emitters, mostly due to leptonic processes of high-energy particles produced in them, a significant fraction of the reported TeV sources can not be associated with any potential VHE sources (H. E. S. S. Collaboration et al. 2018a; Albert et al. 2020). Thus a question is raised about the nature of the TeV sources without known counterparts: whether there are un-revealed SNRs/PWNe or they are some other types of VHE sources. One possibility, recently invoked by the detections of extended TeV emissions around the Geminga and Monogem pulsars (Abeysekara et al. 2017a), is that pulsars could commonly produce (detectable) TeV halos (Linden et al. 2017) through the inverse-Compton scattering process of ∼10 TeV electrons/positrons emitted from pulsars. Indeed, a few of the TeV sources do appear to be positionally coincident with relatively energetic pulsars (H. E. S. S. Collaboration et al. 2018b; Albert et al. 2020), although it should be noted that the positions of the sources determined from VHE observations generally have large uncertainties, ∼0. • 1.
For the purpose of investigating the nature of the TeV sources not in obvious associations, we took advantage of the available different high-energy data, in particular the all-sky survey ones from the Large Area Telescope (LAT) onboard the Fermi Gamma-ray Space Telescope (Fermi), and have carried out systematic studies of the TeV sources by using a multi-wavelength study approach. In this paper, we report our study of the region of 3HWC J1954+286.
For this VHE source, initially excess TeV emission with a 4.3σ significance was found with Milagro at the position of the Fermi LAT source 0FL J1954.4+2838 .
A TeV source named 2HWC J1955+285 was listed in the second HAWC catalog in this region with a Test Statistic (TS) value of 25.4 (Abeysekara et al. 2017b), and this source's TS value was increased to be 48.3 in 3HWC (with an 1σ positional uncertainty of 0. • 14). Furthermore, a source was detected to have the maximum emission energy of 0.42 PeV with the Kilometer Square Array (KM2A) of the Large High Altitude Air Shower Observatory (LHAASO; named LHAASO J1956+2845; Cao et al. 2021), although its position is 0. • 33 away from the 3HWC position. Therefore the nature of this TeV source, which even has a possible PeV counterpart, is extremely interesting.
As pointed out by Abeysekara et al. (2017b), a young pulsar, PSR J1954+2836, is known to be located within the 1σ error circle of 3HWC J1954+286. Its spindown luminosityĖ ≃ 1.0 × 10 36 erg s −1 , and it was detected in γ rays with Fermi LAT (Saz Parkinson et al. 2010), having a γ-ray luminosity of ∼ 4.1 × 10 34 erg s −1 (where a dispersion-measurement distance of 1.96 kpc is used). It can be noted that there is a SNR, G65.1+0.6, also in this source region ). This SNR was discovered in the radio band with low surface brightness (Landecker et al. 1990;Kothes et al. 2006). It has a large SNR-like shell with a size of ∼90 ′ ×50 ′ (Landecker et al. 1990;Tian & Leahy 2006). A distance of 9.2 kpc was derived for this SNR based on the HI ob- Figure 2. Phase-connected pulse profile (top) and twodimensional phaseogram (bottom) in 32 phase bins obtained for PSR J1954+2836 during the whole LAT data time period. The dotted and dashed lines mark the bridge and offpulse phase ranges respectively. Two spin cycles are shown for clarity. servations, and a Sedov age of 4-14×10 4 yr was estimated (Tian & Leahy 2006).
Below we first describe our analysis of the Fermi LAT data in Section 2 for the region, in which we were able to remove the emission of the pulsar. We also analyzed archival X-ray data covering the region. The results from the analyses are presented and discussed in Section 3.  Fermi Pass 8 database was used. Following the recommendations of the LAT team 1 , the events with quality flags of 'bad' and zenith angles larger than 90 degrees were excluded. A source model was constructed on the basis of the recently released Fermi LAT 10-year source catalog (4FGL-DR2; Abdollahi et al. 2020;Ballet et al. 2020). The sources in 4FGL-DR2 within a 20 degree radius circular region from PSR J1954+2836 were included in the source model. Their spectral forms provided in the catalog were used. The background Galactic and extragalactic diffuse spectral models (gll iem v07.fits and iso P8R3 SOURCE V3 v1.txt respectively) were also included in the source model.

Timing analysis of PSR J1954+2836
PSR J1954+2836 is bright in the LAT γ-ray band, and so we first worked to separate its emission from other sources in the region through timing analysis. Following Saz Parkinson et al. (2010), we selected photons within a 1 • radius circular region centered at the pulsar in order to construct its pulse profile. We tested to fold them with the ephemeris given in the second Fermi LAT catalog of γ-ray pulsars (Abdo et al. 2013; see also Table 1, with the first set of the parameters from the online material of Abdo et al. 2013). No clear pulse profile over the 13-yr-long time period could be obtained.
We changed to assign pulse phases to the photons during MJD 54682-55800 (the same time period as that in 1 http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/ Abdo et al. 2013) according to the known ephemeris, using the Fermi TEMPO2 plugin Hobbs et al. 2006). An empirical Fourier template profile was extracted. Using this template, we were able to generate the time-of-arrivals (TOAs) for each sets of ∼200-day data over the whole LAT data time period. The template and TOAs were obtained using the maximum likelihood method described in Ray et al. (2011). We fitted the TOAs with TEMPO2. The known ephemeris has three timing parameters: f , f 1, and f 2 (Table 1). We updated the ephemeris by adding highorder frequency derivatives, while we tried including as many 200-day datasets as possible. However an updated ephemeris (its f , f 1, and f 2 are given in Table 1) only describing the TOAs before MJD ∼57500 was obtained.
For the data after MJD ∼57500, we were able to find a template profile during MJD 58300-59470 (where the latter date is the end time of the whole data), by searching around the pulsar catalog ephemeris values. Same as the above by adding high-order frequency derivatives, we obtained an ephemeris that can describe the TOAs after MJD 57500 (the main parameters are given in Table 1).
We folded the LAT photons before and after MJD 57500 according to the two ephemerides respectively. The folded pulse profiles are shown in Figure 1. We cross-correlated the two profiles in the Fourier domain and obtained a phase shift of ≃0.35 (Figure 1). Applying the phase shift to the second part of the photons, the pulse profile over the whole data period was obtained. In Figure 2, the pulse profile and two-   Figure 3) removed and all catalog sources removed, and the right panel shows the region in 0.5-500 GeV band with all catalog sources removed. The cyan contour is an approximate radio shape of the SNR G65.1+0.6 (obtained from the SNRcat: http://snrcat.physics.umanitoba.ca/SNRrecord.php?id=G065.1p00.6), and the black and white dashed lines mark the 2σ error circles centered at the best-fit positions for PS1 and PS2 (cf., Section 2.1.4). Other marks are the same as those in Figure 3.
dimensional phaseogram in 32 phase bins are plotted. Based on the profile, we defined phase 0.0-0.19 and 0.47-0.66 as the onpulse phase ranges, phase 0.19-0.47 as a bridge phase range, and phase 0.66-1.0 as the offpulse phase range.

Likelihood analysis for the onpulse and bridge data
We performed standard binned likelihood analysis to the 0.1-500 GeV LAT data during the onpulse and bridge phase ranges of the pulsar determined above. The RoI and the source model, set in Section 2.1.1, were used. The sources in our source model within 5 degrees from the pulsar were set to have free spectral parameters; for the other sources, their spectral parameters were fixed at the values given in the catalog. The background normalizations were set as the free parameters.
For the emission at the pulsar's position in the two phase ranges, which is presumably dominated by that from the pulsar (Figure 3), we used a sub-exponentially cutoff power-law model, , where Γ and E c are the photon index and cutoff energy respectively, and b is a measure of the shape of the exponential cutoff. We fixed b to be 2/3, which is a characteristic value used for the γ-ray pulsars in 4FGL-DR2. We obtained Γ = 1.42 ± 0.01 and E c = 1.60 ± 0.03 (Γ = 1.29 ± 0.06 and E c = 0.96 ± 0.07) for the onpulse (bridge) phase range. The values are close to those given in 4FGL-DR2 (Γ = 1.43 ± 0.06, E c = 1.4 ± 0.1 GeV for the whole data). These results, as well as the corresponding Test Statistic (TS) values, are provided in Table 2.
In order to show the appearances of the pulsar and other detected sources during the different phase ranges of the pulsar, we calculated the 0.1-500 GeV TS maps of a 3 • × 3 • region centered at the pulsar ( Figure 3). All catalog sources in the source model were kept in the TS maps. The pulsar was dominantly bright during the onpulse phase range, and in the bridge and offpulse phase ranges, a nearby source, 4FGL J1958.7+2846, appeared brighter.

Likelihood analysis for the offpulse data
For the offpulse phase range, if we still assumed the sub-exponentially cutoff model for the 0.1-500 GeV γray emission at the position of PSR J1954+2836, the cutoff could not be significantly determined. If we instead used a simple power law dN/dE = N 0 E −Γ , we would obtain Γ = 2.6 ± 0.1, which is significantly larger than those in the onpulse and bridge phase ranges.
In order to examine the region, we calculated 0.1-500 GeV TS maps of a 3 • × 3 • region centered at the pulsar. The left and middle panels of Figure 4 respectively show the region with all catalog sources kept (except 4FGL J1958.7+2846) and removed; the nearby source 4FGL J1958.7+2846 (cf., Figure 3) was removed in the left panel such that the other sources in the region can be relatively well seen.
The TS map shown in the middle panel of Figure 4 indicates that the residual γ-ray emission is obviously not centered at the pulsar; moreover, there appear to likely have two point sources. We ran gtfindsrc in Fermitools to determine their positions. The obtained best-fit positions are (R.A., Decl.) = (298. • 44, 28. • 38) and (R.A., Decl.) = (298. • 4, 28. • 9) for point source 1 (PS1) and 2 (PS2), respectively, with the 2σ nominal uncertainties of 0. • 07 and 0. • 2 (see Figure 4). The position for PS1 seems Figure 5. γ-ray spectra and corresponding 0.1-500 GeV spectral models of PSR J1954+2836 during the onpulse and bridge phase ranges (red and green respectively), and of PS1 during the offpulse phase range (black and light gray area). The dark gray area is the HAWC spectrum for 3HWC J1954+286 (Albert et al. 2020), and the gray bar is the flux measurement for LHAASO J1956+2845 at the energy of 100 TeV (Cao et al. 2021). The thick solid lines and long-dashed lines mark the upper limits obtained with the VERITAS from the point-source search and extended-source search, respectively (the spectral indices of 2 and 2.4 were considered; see Abeysekara et al. 2018). to be off from the highest-TS location. However as we checked TS maps at different higher energy bands such as 0.3-500, 0.5-500, and 1-500 GeV, we found that the PS1 position is right at the locations with the highest TS values in them. The 0.5-500 GeV TS map is shown as an example in the right panel of Figure 4. The localization result obtained with 0.1-500 GeV data is thus likely correct, as the large Point Spread Function (PSF) of LAT (possibly resulting in contamination from nearby sources) as well as the strong Galactic background in the low energy range of <300 MeV can cause the off-peak appearance. We also noted that in the 0.5-500 GeV TS map, PS2 is nearly invisible, indicating that it had very soft emission and likely is a separate source not in association with PS1.
As the location of PS1 is away from the pulsar but within the SNR G65.1+0.6 (Figure 4), we considered it as a possible counterpart to the SNR and re-performed the likelihood analysis to the LAT data during the offpulse phase range, in which we removed PSR J1954+2836 from the source model and included PS1. A simple power law was assumed for the source, and we obtained Γ = 2.4 ± 0.1 (the details are given in Table 2). We noted that the radio flux density contour of the SNR is much extended and contains PS2, and so we tested to use the radio contour as a template. The likelihood analysis was conducted by including this template in the source model. The likelihood value L ext was compared to that L 2PS when PS1 and PS2 were considered as two point sources or that L PS1 when only PS1 was considered. The resulting results were −2 × log(L ext /L 2PS ) ∼23 or −2 × log(L ext /L PS1 ) ∼16. As either value is significantly large, suggesting that no extended γ-ray emission was preferred over the individual point sources.

Spectral and variability analysis
We extracted the γ-ray spectra of PSR J1954+2836 during the onpulse and bridge phase ranges and that of PS1 during the offpulse phase range. The maximum likelihood analysis of the LAT data in 10 evenly divided energy bands in logarithm from 0.05 to 500 GeV was conducted. In this analysis, the spectral normalizations of the sources within 5 degrees from each target were set as free parameters, while all the other parameters of the sources were fixed at the values obtained from the maximum likelihood analysis in the above section. For the results, we kept the spectral data points when the fluxes are >2 times larger than the uncertainties and derived 95% flux upper limits otherwise. The obtained spectra are shown in Figure 5, and the values are given in Table 3.
In order to fully study the γ-ray emission properties of PS1, we also searched for its long-term variability. We calculated the variability index TS var with 80 time bins (each bin was constructed from 60-day data) in the energy range of 0.1-500 GeV, following the procedure introduced in Nolan et al. (2012). If its fluxes were consistent with a constant, TS var would be distributed as a χ 2 function with 79 degrees of freedom. The computed TS var value was 96.9, lower than the threshold value of 111.1 (at a 99% confidence level) for identifying a variable source. Therefore PS1 had no significant long-term variability.

Chandra X-ray data analysis
Chandra X-ray telescope observed PSR J1954+2836 on 2011 December 26 (Obs. ID 12148) with an effective exposure time of 9.84 ks. The Advanced CCD Imaging Spectrometer array was used in the observation to image the source region. We processed the data with the software CIAO v4.13 and the Chandra Calibration Database (CALDB v4.9.5). Task wavdetect was used for detecting sources in the field, but no sources were found at the pulsar's position or in the nearby regions. The upper limit on the count rate at the position is 4.56×10 −4 cts s −1 . We estimated an unabsorbed model flux based on the count-rate upper limit. Assuming a power law with a photon index of 1.5 and a hydrogen column density of N H ∼ 1.1×10 22 cm −2 towards the source position (HI4PI Collaboration et al. 2016), the flux upper limit of ∼ 8.9 × 10 −15 erg s −1 cm −2 in 0.5-7.0 keV is given by Chandra tool PIMMS.

RESULTS AND DISCUSSION
To investigate the nature of the VHE source 3HWC J1954+286, we analyzed 13 years of Fermi LAT Pass 8 data for the source's region. By determining the onpulse and bridge phase ranges of PSR J1954+2836 and using the different phase-range LAT data accordingly, we were able to separate the pulsar's emission from that of the region. The pulsar's emission in the onpulse and bridge phase ranges can be well fitted with a typical spectral shape of pulsars, a sub-exponentially cutoff power law (Abdo et al. 2013;Abdollahi et al. 2020). After removing the pulsar's emission component (i.e., by using the data only during the offpulse phase range), excess emission has been found in the region. There are likely two point sources, and their determined positions are significantly away from that of the pulsar. Also from the likelihood analysis of the data in 0.1-500 GeV band, we found that the emission of each of the two sources is consistent with being described with a power law (see Table 2). Given the positional differences and spectral types, the two sources very likely are not emission from the pulsar in the offpulse phase range.
For the HAWC TeV source 3HWC J1954+286, both PS1 and PS2 appear away from it, having 2.3σ and 1.9σ positional differences respectively. As PS2 was faint (Table 2) and limited property information was extracted for it, below we mainly discuss the more significantly detected PS1 (∼6σ detection significance). We have obtained its spectrum. The extension of the spectrum according to its power-law index would predict approximately an order of magnitude lower TeV fluxes than those detected by HAWC ( Figure 5). Based on these facts, the excess emission is not likely the GeV counterpart to the TeV source either. We note that the Very Energetic Radiation Imaging Telescope Array System (VERITAS) also observed the region (Abeysekara et al. 2018), and the flux upper limits derived from the observation for a point source have reached the upper bound of the extension of the PS1's spectrum ( Figure 5), which potentially provides a constraint on any TeV emission of PS1. An additional note is that the HAWC source would be extended, such that the VERITAS upper limits become consistent with the HAWC flux measurements ( Figure 5; see discussion in Abeysekara et al. 2018).
One possibility is if the excess emission could be the PWN associated with the pulsar. However, if the excess emission is the PWN, we would have detected it in X-rays. According to theoretical calculations as well as observational evidence given in Zhu et al. (2018), a pulsar like PSR J1954+2836 would have an X-ray-to-γ-ray luminosity ratio of 5 × 10 −3 . The source PS1 has a γ-ray luminosity of L γ ≃ 4.5 × 10 33 (d/1.96 kpc) 2 erg s −1 , where the pulsar's dispersion-measurement distance is assumed. This γ-ray luminosity would predict an X-ray luminosity of 2 × 10 31 erg s −1 , or 4 × 10 −14 erg s −1 cm −2 (at 1.96 kpc). Comparing this value to the flux upper limit set by the Chandra observation, we would see some evidence for an X-ray source at the pulsar's location in the X-ray image. Therefore the non-detection of any X-ray sources in the sensitive Chandra observation does not support it as the PWN. In addition, the power-law-like spectrum of PS1 also does not support the PWN origin since the PWNe generally have a γ-ray spectrum similar to those of pulsars (Ackermann et al. 2011), although the γ-ray efficiency L γ /Ė ≃ 4.5 × 10 −3 would be consistent with the range of γ-ray efficiency values of the known PWNe (Ackermann et al. 2011).
The possibility remaining to be considered is that the excess emission PS1 would be the γ-ray counterpart to the SNR G65.1+0.6, as their positions are coincident with each other. Below in Section 3.1, we discuss possible scenarios to explain the excess γ-ray emission as arising from the SNR. In the end the section 3.2, we discuss the speculation if the TeV source could be related to the SNR or the pulsar in the region given their properties.

Emission scenarios for the SNR G65.1+0.6
At the SNR's distance of 9.2 kpc, the 0.1-500 GeV luminosity was ≃ 1×10 35 erg s −1 , which is in the range of middle-aged SNRs (∼10 35 − 10 36 erg s −1 ; Acero et al. 2016). By interacting with dense clouds, middle-aged SNRs have γ-ray emission dominated by that arising from the hadronic process. They may have flux peaks at ∼1 GeV (e.g. W44, Giuliani et al. 2011;HB 21, Pivato et al. 2013;IC 443, Ackermann et al. 2013). For several low-significance ones (e.g. MSH 17−39, Castro et al. 2013;Kes 27, Xing et al. 2015), γray emission could only be significantly detected above ∼1 GeV band. The γ-ray spectrum we have obtained is similar in this respect, although there are only three flux data points and the uncertainties are large.
Below we first tested if the γ-ray spectrum could be fitted with a hadronic model. As a leptonic model is often considered for γ-ray emission of SNRs, we also tested the leptonic scenario. We simply assumed that the particles accelerated in the SNR have a power-law form with a high-energy cutoff: where, i = e, p, α i and E i,c are the power-law index and high-energy cutoff, respectively. The normalization A i is determined by the total energy (W i ) in particles with energies above 1 GeV. We used a python package naima (Zabalza 2015), a code for computation of non-thermal radiation from relativistic particle populations based on emcee (Foreman-Mackey et al. 2013), to explore the parameter space. For the hadronic model, the γ-ray emission is produced by the proton-proton inelastic collision. Due to lack of constraints on the cutoff energy, we fixed E p,c = 1 PeV in our calculation. Assuming the target density n 0 = 1 cm −3 , we obtained α p = 2.83 +0.43 −0.40 and W p = 1.14 +0.65 −0.40 × 10 51 erg. The resulting best-fit spectral energy distribution (SED) with 1σ confidence interval is shown in the left panel of Figure 6. The fitted proton index is consistent with that of the interacting SNRs (Acero et al. 2016). Although the energy in protons exceeds the typical explosion energy, it is acceptable if the target density is greater than 1 cm −3 (or we may also suspect that the distance of the SNR would be significantly smaller than 9.2 kpc).
The leptonic scenario considers the inverse Compton (IC) process of high-energy electrons upper-scattering the seed photons in the region. The cosmic microwave background as well as the infrared component of the interstellar radiation field (with a temperature of 35 K and an energy density of 0.2 eV cm −3 ; Porter et al. 2006;Shibata et al. 2011) were included in the calculation.
The same population of electrons also produce emission from radio to X-rays via the synchrotron process with the average magnetic field B. The resulting bestfit SED with 1σ confidence interval is displayed in the right panel of Figure 6. The parameters (with 1σ uncertainties) are α e = 1.82 +0.11 −0.13 , E e,c = 152 +202 −68 GeV, B = 2.4 +2.4 −1.0 µG, and W e = 6.0 +7.9 −3.7 × 10 49 erg. For the SNR with an age of ∼ 10 5 yr (Tian & Leahy 2006), the shock velocity is expected to decrease below ∼100 km s −1 . Thus the cooling-limited maximum energy of electrons is ∼ 7(B/2.4 µG) −0.5 TeV (e.g., Zhang & Liu 2019). Considering the higher magnetic field in the early evolution stage, the fitted cutoff energy of electrons is in line with expectations.
From the fitting, we found that the SED can be described with reasonable parameters of the SNR in either a hadronic or a leptonic scenario. The fitting results thus support the association of the γ-ray source PS1 with the SNR, in addition to their positional coincidence.

Origin of the TeV emission?
As both HAWC and LHAASO KM2A have detected TeV and upto PeV emissions respectively in the region of SNR G65.1+0.6 and PSR J1954+2836 (while it can be noted that the two VHE positions have a ∼ 2σ difference) and there are no other potential VHE sources nearby, it is conceivable to think that the VHE emissions would be associated with either the SNR or the pulsar.
Considering that the GeV emission we have detected as arising from the SNR G65.1+0.6, the extension of the simple power-law GeV spectrum would predict lower than the detected TeV-PeV fluxes. Also the GeV emission appears away from the HAWC position (with a 2.3σ difference). If the VHE emissions are associated with the SNR is thus highly questionable. While our leptonic modeling does not require the association as the model fit drops steeply in the VHE range ( cf., Figure 6; we note that unless the one-zone assumption in our model is changed), the hadronic model fit is constrained by the VERITAS upper limits. As shown in the left panel of Figure 6, the GeV emission could not be associated with the HAWC source; otherwise the VHE emission from SNR G65.1+0.6 would have been detected by the VERITAS observation. In any case the picture for this region hopefully could be clarified with more and deeper VHE observations, in particular with the LHAASO Water Cherenkov Detector Array (WCDA; (Bai et al. 2019)) or the Cherenkov Telescope Array (CTA) north (Cherenkov Telescope Array Consortium et al. 2019). The former has a sensitivity 2 times that of HAWC and the latter would observe the region deeper than VER-ITAS and potentially detect the VHE emission from SNR G65.1+0.6 if its emission is due to the hadronic process.
As PSR J1954+2836 is located right within the 1σ error circle of 3HWC J1954+286, their association could be a more likely possibility. This pulsar was identified to have pulsed emission detected in >10 GeV band (Ajello et al. 2017; see also Figure 6 and Table 3 in this work). We check the pulse-phase distribution of high energy ≥ 12 GeV photons in the 0. • 1 circular region around the pulsar (note that the LAT's 68% containment angle at 12 GeV is ∼ 0. • 1). As shown in Figure 7, their pulse phases match the onpulse phase ranges of the pulsar well. As a comparison in the bridge and offpulse phase ranges, few such photons were seen. In addition, the maximum energy of the photons is ∼ 300 GeV (Figure 7). The results do suggest that the pulsar at least have emission upto ∼100 GeV. Search for pulsed VHE emission from this pulsar could provide interesting results, as whether young, energetic pulsars would have VHE emission similar to the Crab (VERITAS Collaboration et al. 2011;Aleksić et al. 2012)  On the other hand as first pointed out by Linden et al. (2017), pulsars could generally have TeV halos such as the cases of the Geminga and Monogem pulsars (Abeysekara et al. 2017a). Using the TeV halo of Geminga as a typical case, estimation about the TeV halo of PSR J1954+2836 can be made (Linden et al. 2017): its 7 TeV flux and spatial extension would be ∼ 1.4 × 10 −14 (d/1.96 kpc) −2 TeV −1 cm −2 s −1 and ∼ 0. • 2, respectively, where the distance of 0.19 kpc and the 7 TeV flux of 4.9 × 10 −14 TeV −1 cm −2 s −1 for the Geminga are used. The estimated 7 TeV flux is approximately 2 times the HAWC measurement and 3HWC J1954+286 is given as a point source (Albert et al. 2020), but considering the simple scaling estimation and the uncertainty on the pulsar's distance, the differences may not be significant. Given the considerations discussed above, we thus suggest that the VHE emissions would more likely be the TeV halo associated with the pulsar.  (7) MJD 54682-57500 10.7863425017(4) -2.46225(5) 0.11 (7) MJD 57500-59470 10.786340(3) -2.41(6) -7(9) Note: the frequency epoch is MJD 55225.