3HWC J0631+107/LHAASO J0631+1040: A TeV Halo Powered by the Pulsar J0631+1036?

PSR J0631+1036, discovered by Zepka et al. (1996), is a middle-aged one having spin period P ≃ 0.288 s, characteristic age τ_c ≃ 43.6 kyr, and spin-down luminosity $\dot{E}\simeq 1.7\times {10}^{35}$ erg s⁻¹. Based on the new electron-density model for the Galaxy (Yao et al. 2017), its distance is D ≃ 2.1 kpc given in the Australia Telescope National Facility Pulsar Catalogue (Manchester et al. 2005). In X-rays, observational studies have not detected the pulsar down to ≃4.9 × 10³⁰(D/2.1)² erg s⁻¹ (in 0.5–2.0 keV; Becker & Truemper 1997; Kennea et al. 2002).

This seemingly "normal" pulsar, along with several tens of others, has been selected as one likely associated with the Galactic very-high-energy (VHE; >100 GeV) TeV sources, namely J0631+107 reported by the High-Altitude Water Cherenkov (HAWC) Observatory in the third HAWC catalog (3HWC; Albert et al. 2020) and J0631+1040 by the Large High Altitude Air Shower Observatory (LHAASO; Cao et al. 2019) in the The First LHAASO Catalog of Gamma-Ray Sources (Cao et al. 2023). The reasons for finding associations between VHE sources and pulsars are the following. First, pulsars with τ_c ≤ 100 kyr can have pulsar wind nebulae (PWNs), which are considered as one primary type of the Galactic TeV sources (e.g., H.E.S.S. Collaboration et al. 2018a, 2018b). Second, as inspired by the detection of the extended TeV emissions around two nearby pulsars, Geminga and Monogem (Abeysekara et al. 2017a), a new type of TeV sources, the so-called TeV halos that are powered by middle-aged pulsars, has been proposed (Linden et al. 2017; Linden & Buckman 2018). Third, among more than 100 sources detected in recent years with the VHE facilities, mainly the High Energy Spectroscopy System (HESS; H.E.S.S. Collaboration et al. 2018a), the HAWC, and the LHAASO, a significant number of the sources do not have typical known counterparts, i.e., the PWNs, the supernova remnants (SNRs), or other types of VHE emitters (H.E.S.S. Collaboration et al. 2018a; Albert et al. 2020; Cao et al. 2023).

Given these and intrigued by the third one described above, we have been carrying out multiwavelength studies of the TeV sources that do not have obvious counterparts at other energy bands (Xing et al. 2022; Zheng et al. 2023). In our studies, we noted that 3HWC J0631+107 (hereafter J0631+107), also LHAASO J0631+1040 given the positional match between these two sources, has a clean field in high energies. There are no known PWNs or SNRs found such as in the TeV online catalog (TeVCat; Wakely & Horan 2008), and PSR J0631+1036 is the only notable source based on the position matches (e.g., Cao et al. 2023). Interestingly, this pulsar has bright GeV γ-ray emission, detected with the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope (Fermi) from the early observations (Weltevrede et al. 2010). We thus conducted detailed analysis of the Fermi LAT data for the pulsar. The analysis and results are described below in Section 2, based on which we argue that J0631+107 is likely a TeV halo powered by PSR J0631+1036; the related discussion is presented in Section 3.

2.1. LAT Data and Source Model

2.2. Timing Analysis of PSR J0631+1036

PSR J0631+1036 is bright in the LAT energy band and located in a clean field, as shown in the left panel of Figure 1, a test statistic (TS) map calculated for the source region from the whole LAT data (Section 2.3.1). The pulsar is the only 4FGL-DR3 source in the 2° × 2° TS map. Also seen is its positional match to J0631+107.

Zoom In Zoom Out Reset image size

In order to check if there are other sources hiding in the bright emission of the pulsar, we worked to obtain its pulsed emission through timing analysis. On the first try, we folded the photons within a 6° radius (∼size of the point-spread function of LAT at 100 MeV ⁴ ) aperture centered at the pulsar according to the ephemeris given in the LAT Gamma-ray Pulsar Timing Models Database ⁵ (Ray et al. 2011), but no clear pulse profile over the ≃14.5 yr long data could be obtained.

We then changed to use the method fully described in Xing et al. (2022). In essence, we divided the data into sets of 200 days, and assigned pulse phases to the photons according to the ephemeris in the Database (Ray et al. 2011) by using the Fermi TEMPO2 plug-in (Edwards et al. 2006; Hobbs et al. 2006). We were able to obtain empirical Fourier template profiles before and after MJD 56770, generate the times of arrival (TOAs) for each set of ∼200 days data, and obtain timing solutions by fitting the TOAs with high-order frequency derivatives. We could not extended the timing solutions to times longer than MJD 57930, probably due to the glitches of the pulsar at MJD 58341 and 58352 (Lower et al. 2021; Basu et al. 2022).

With the two timing solutions, the photons during the two time periods were folded respectively. The two pulse profiles had a phase mismatch of ≃0.3075 (which was directly read off from the profiles because of the clear pulse shape). After applying the phase shift to the photons of the second time period, the pulse profile over nearly 9 yr was obtained (Figure 2). Based on the pulse profile, we defined phase 0.0625–0.5625 as the on-pulse phase range and phases 0.0–0.0625 and 0.5625–1.0 as the off-pulse phase ranges.

Zoom In Zoom Out Reset image size

2.3. Likelihood and Spectral Analysis

2.3.1. Whole and On-pulse Data

We performed standard binned likelihood analysis to the 0.1–500 GeV LAT data during the whole ∼14.5 yr time period and the on-pulse phase range during the ∼9 yr time period. The spectral parameters of the sources within 5° from the pulsar in the source model were set free, while the spectral parameters of the other sources were fixed at the values given in 4FGL-DR3. In addition, the normalizations of the two background components were set free.

For the emission at the pulsar's position, in the whole data or those of the on-pulse phase, we used a subexponentially cutoff power-law model (Abdollahi et al. 2020), $\displaystyle \frac{{dN}}{{dE}}={N}_{0}{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-{\rm{\Gamma }}-\tfrac{d}{2}\mathrm{ln}\left(\tfrac{E}{{E}_{0}}\right)-\displaystyle \frac{{db}}{6}{\mathrm{ln}}^{2}\left(\displaystyle \frac{E}{{E}_{0}}\right)-\displaystyle \frac{{{db}}^{2}}{24}{\mathrm{ln}}^{3}\left(\tfrac{E}{{E}_{0}}\right)}$, where Γ and d are the photon index and the local curvature at E₀ respectively, and b is a measure of the shape of the exponential cutoff. We fixed b = 2/3, which is such set for most of the γ-ray pulsars in the LAT catalogs (e.g., Abdollahi et al. 2020, 2022).

The likelihood analysis results, including the TS values, are provided in Table 1. The parameter values of the pulsar are consistent with those given in 4FGL-DR3. A 0.1–500 GeV TS map of a 2° × 2° region centered at the pulsar was calculated and is shown in the left panel of Figure 1.

Table 1. Binned Likelihood Analysis Results

Phase Range (src)	F0.1−500/10−8	Γ	d	TS	$\mathrm{log}(L)$
	(photons s−1 cm−2)
Catalog values	⋯	1.85 ± 0.05	0.54 ± 0.07	2000	⋯
Whole data	3.01 ± 0.23	2.02 ± 0.03	0.47 ± 0.05	2570	⋯
On-pulse	2.23 ± 0.19	1.96 ± 0.04	0.62 ± 0.07	2279	⋯
Off-pulse (PS1)	1.35 ± 0.29	2.66 ± 0.13	⋯	60	763784.9
Off-pulse (2PS-PS1)	0.92 ± 0.41	2.63 ± 0.17	⋯	29	763786.9
(2PS-PS2)	0.45 ± 0.36	2.53 ± 0.25	⋯	9	⋯

Download table as:  ASCIITypeset image

We extracted the γ-ray spectrum of PSR J0631+1036 in the on-pulse phase data. The spectral bands were set as 10 evenly divided in logarithm from 0.1 to 500 GeV. In the analysis of obtaining fluxes in the bands, the spectral normalizations of the sources within 5° of the pulsar were set as free parameters, while all the other parameters of the sources were fixed at the values obtained from the above binned likelihood analysis. For the spectral data points, we kept those with TS ≥ 4 and derived 95% flux upper limits otherwise. The obtained spectrum is shown in Figure 3.

Zoom In Zoom Out Reset image size

2.3.2. Off-pulse Data

With the off-pulse phase ranges obtained above (Figure 2), we examined the source region by first calculating a TS map in 0.1–500 GeV from the off-pulse data. No source emission could be detected at the pulsar's position, with a 95% flux upper limit of ∼10⁻¹³ erg cm⁻² s⁻¹ (assuming Γ = 2 in 0.1–500 GeV; see Figure 3). However residual emission (TS ∼ 50) southeast of the pulsar is seen (middle panel of Figure 1).

To further check the residual emission, a TS map in 1–500 GeV was also obtained (right panel of Figure 1). As can be seen, there seemingly are two sources. We ran gtfindsrc in Fermitools to determine their positions. The obtained best-fit positions are (R.A., decl.) = (9803, 1059) and (R.A., decl.) = (9834, 1052) for point source 1 (PS1) and 2 (PS2), respectively, with the 1σ nominal uncertainties of 011 and 013 (indicated in Figure 1). By including PS1 or PS1+PS2 in the source model, we performed the likelihood analysis to the off-pulse data. We found that PS2 was not significantly detected, indicated by the likelihood value being similar to that when only PS1 was in the source model (see Table 1). We extracted the spectrum of PS1 (Figure 3), which could be fitted with a power law with Γ ≃ 2.66).

We conducted analysis of the LAT data for PSR J0631+1036, because of its possible association with the TeV source J0631+107 and the absence of PWN/SNR-like counterparts in the source region. By timing the pulsar, we were able to removed its pulsed emission in a ∼9 yr time period of the data. No off-pulse emission was detected at the pulsar's position. Residual emission, PS1, was seen approximately 016 east of the pulsar. The emission was soft, mostly detectable at ≲1 GeV (Figure 3). We checked the SIMBAD database for possible counterparts to it, but no obvious ones (particularly in radio or X-rays) were found within its error circle. The nature of PS1 is thus not clear. Given the positional offset and its soft emission, it is not likely in association with the pulsar or the TeV source.

Then it is straightforward to note the similarities of PSR J0631+1036 to Geminga. They have similar P values and are both γ-ray bright, while the former's τ_c is younger by a factor of ∼8 and $\dot{E}$ higher by a factor of ∼5. Given these and our analysis results for the field, we thus argue that J0631+107 is likely the TeV halo of PSR J0631+1036. In Figure 4, we compare this pulsar to Geminga. The latter's X-ray, γ-ray, and TeV halo fluxes, shown in the figure, are scaled to the distance (2.1 kpc) of the former, where the nominal distance 250 pc is used for Geminga (Manchester et al. 2005). As can be seen, the X-ray upper limit on PSR J0631+1036 (Kennea et al. 2002) is approximately consistent with the X-ray fluxes of Geminga or its PWN (Posselt et al. 2017), where the fluxes of all the components of the latter's PWN are added together as the flux in 0.3–8 keV. Because the TeV halo of Geminga is extended with fine structures (Abeysekara et al. 2017a), we adopt the flux measurement from the second HAWC catalog (Abeysekara et al. 2017b), in which a 2° extension was used. The scaled flux at 7 TeV is ∼7 × 10⁻¹⁶ TeV⁻¹ cm⁻¹ s⁻¹, approximately ∼1/6 of that of J0631+107. Interestingly, the ratio is similar to that in $\dot{E}$ of Geminga to PSR J0631+1036. The size of the TeV halo of PSR J0631+1036 would be roughly 024 by taking that of Geminga as a standard (Linden et al. 2017), smaller than the upper limit of 030 set by the LHAASO (Cao et al. 2023).

Zoom In Zoom Out Reset image size

We further note that the emission of J0631+107 is hard, as LHAASO detected it in 25–100 TeV with Γ ≃ 3.3 but did not detect it in 1–25 TeV (Figure 3). This spectral feature is similar to that of Geminga's TeV halo, indicated by the LHAASO detection of Γ ≃ 1.7 in 1–25 TeV and Γ ≃ 3.7 in ≥25 TeV (i.e., the spectrum likely peaks around ∼25 TeV). This type of spectra is harder than those of PWNs, since the latter have Γ ≳ 2 in 1–10 TeV and thus some of them can be detected with Fermi LAT (H.E.S.S. Collaboration et al. 2018b and references therein); indeed, part of the purpose of this work was to search for a PWN in the off-pulse data. Hopefully with more data collected with LHAASO in the near future, the similarity of the spectrum of J0631+107 to that of the Geminga's TeV halo can be established, and thus firmly confirm the TeV halo nature of J0631+107.

This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This research is supported by the Basic Research Program of Yunnan Province (No. 202201AS070005), the National Natural Science Foundation of China (12273033), and the Original Innovation Program of the Chinese Academy of Sciences (E085021002).