Discovery and timing of pulsar J2016$+$3711 in supernova remnant CTB 87 with FAST

We report on our discovery of the radio pulsar, PSR J2016$+$3711, in supernova remnant (SNR) CTB 87, with a $\sim10.8\sigma$ significance of pulses, which confirms the compact nature of the X-ray point source in CTB 87. It is the first pulsar discovered in SNRs using Five-hundred-meter Aperture Spherical radio Telescope (FAST). Its integrated radio pulse profile can be well described by a single component, with a width at 50% of the peak flux density of about 28.1$^\circ$ and an effective width of about 32.2$^\circ$. The mean flux density at 1.25 GHz is estimated to be about 15.5$\mu$Jy. Combined with the non-detection of the radio pulse at lower frequencies, the radio spectral index of the pulsar is constrained to be $\lesssim 2.3$. We also present the timing solution based on 28 follow-up FAST observations. Our results reveal a period of 50.81 ms, period derivative of $7.2\times 10^{-14}$ s s$^{-1}$, and dispersion measure of 428 pc cm$^{-3}$. The strength of the equatorial surface magnetic dipole magnetic field is inferred to be about $1.9\times10^{12}$ G. Using the ephemeris obtained from the radio observations, we searched Fermi-LAT data for gamma-ray pulsations but detected no pulsed signal. We also searched for radio pulses with FAST toward the X-ray counterpart of the gamma-ray binary HESS J1832$-$093 proximate to SNR G22.7$-$00.2 but found no signal.


INTRODUCTION
Pulsar wind nebulae (PWNe), one of the most important components of supernova remnants (SNRs), are produced by the relativistic winds of rotation-powered neutron stars as they interact with the surroundings.Nonetheless, only about one percent of the known pulsars are found to be directly associated with PWNe and SNRs (e.g.Lyne & Graham-Smith 2012).On the other hand, although the existence of pulsars in PWNe is beyond doubt, about 30% of the known PWNe and PWN candidates have not been confirmed with pulsed emission signals to date (according to Ferrand & Safi-Harb 2012).Finding pulsars in PWNe and SNRs is crucial to studying the pulsar formation and the supernova explosion mechanism, and bridging the gap between the theoretical prediction and the observational results.
SNR CTB 87, appearing located in a superbubble (Liu et al. 2018) and in physical contact with a molecular cloud (e.g.Liu et al. 2018;Kothes et al. 2003), is characterized by the hosted PWN with trailing morphology in X-rays (Matheson et al. 2013).A point-like X-ray source, CXOU J201609.2+371110, has been found in the southeastern head portion of the trailing structure, which was suggested to be the pulsar candidate (Matheson et al. 2013).Several surveys have ★ E-mail: ygchen@nju.edu.cn† E-mail: wangpei@nao.cas.cncovered SNR CTB 87 to search for the radio pulsar.Using the Low Frequency Array (LOFAR), the point-like X-ray source was observed at about 150 MHz with a sensitivity of about 0.4 mJy (Straal & van Leeuwen 2019).Using the Arecibo radio telescope, CTB 87 was fully searched at 430 MHz with a sensitivity of 0.05 mJy (Gorham et al. 1996).Also, it was searched using the 76 m Lovell radio telescope at 606 MHz with sensitivities of 3 mJy (Biggs & Lyne 1996) and 0.16 mJy (Lorimer et al. 1998), respectively.However, no radio pulsation from CTB 87 was ever detected in the previous searches.Proximate to the southwest of SNR G22.7−00.2, an X-ray source (XMMU J183245−0921539) was found consistent with the location of the gamma-ray binary HESS J1832−093 (HESS Collaboration et al. 2015).Recently, near-infrared spectroscopic observations found the companion to be an O6 V star (van Soelen et al. 2024), while the nature of the compact object remains unknown.The X-ray spectrum of XMMU J183245−0921539 is well-fitted by a power-law model with a photon index suggesting a pulsar nature (Mori et al. 2017).
The steep ATCA spectrum of the radio source coincident with the binary is also consistent with the synchrotron radiation from a radio pulsar residing in the binary (Tam et al. 2020).
Although non-detection of radio pulse may be because the beam of the pulsar does not point toward us, it could also result from the instrumental limitations of the surveys.Now, the Five-hundred-meter Aperture Spherical radio Telescope (FAST), with a large effective area, allows us to search for faint pulsars with unprecedented high sensitivity.Actually, a number of faint pulsars have been detected (e.g.Li et al. 2018;Cameron et al. 2020;Zhang et al. 2019;Han et al. 2021) since the first one discovered with FAST (Qian et al. 2019).
We have observed two X-ray point-like sources, CXOU J201609.2+371110 in SNR CTB 87 and XMMU J183245−0921539 proximate to SNR G22.7−00.2,respectively, with FAST.In this paper, we report on our searches for radio pulses toward them and the detection of pulses in CTB 87.We describe the observations and reduction of the data in Section 2. We present our results, the discovery and timing solution of the pulsar in CTB 87, in Section 3. We discuss the properties of the detected pulsar in Section 4 and summarise this paper in Section 5.

OBSERVATION AND DATA REDUCTION
Our radio observations (PI: Y. Chen) were carried out with FAST during 2019 July 15-16 in the "risk-sharing" open session, for a total time of about 1 hr for each source.The two target positions are centred at (20 h 16 m 09 s .2,+37 • 11 ′ 10 ′′ .5,J2000) and (18 h 32 m 45 s .0,−9 • 21 ′ 53 ′′ .9,J2000) (denoted as positions P CTB87 and P G22.7 hereafter) for the point-like X-ray sources in SNR CTB 87 and close to SNR G22.7−00.2,respectively (see Fig. 1; Matheson et al. 2013;HESS Collaboration et al. 2015).In the observations, the central beam of the 19-beam 1.05-1.45GHz receiver (Jiang et al. 2019) was used, with the tracking mode applied.The half-power beamwidth of the telescope is about 3.2 ′ at 1.25 GHz and the pointing accuracy of the telescope is better than 16 ′′ .The digital back-end provided a bandwidth of 400 MHz with 4096 frequency channels and a time resolution of 49.152 s.The data were recorded in the search mode of PSRFITS format.
We used the PulsaR Exploration and Search TOolkit (PRESTO) software (Ransom 2001) to reduce the data.The data were first cleaned from the radio frequency interference (RFI) and then were used to search for periodic signals.The dispersion search ranges of the targets (see Table 1) were limited based on the dispersion measure (DM) predicted by the electron density model (Yao et al. 2017).The data were de-dispersed in steps of 0.1 pc cm −3 below DM = 230 pc cm −3 , steps of 0.3 pc cm −3 for 230 pc cm −3 < DM < 446 pc cm −3 , steps of 0.5 pc cm −3 for 446 pc cm −3 < DM < 800 pc cm −3 and steps of 1 pc cm −3 above 800 pc cm −3 , which were chosen using the ddplan routine in PRESTO.We selected candidates of pulse signals with PRESTO-reported significance above 4.The selected candidates were then folded and inspected by eyes to clarify them as either RFI or possible radio pulse.We also searched for single pulses using the single_pulse_search.py routine in PRESTO.Yet no positive candidate of single pulse has been found toward any of the targets.
After the radio pulses from the CTB 87 PWN were discovered in the "risk-sharing" observation in 2019 (see §3.1 below), 28 followup timing observations were then carried out until 2023, with the same receiver and back-end setup as for the discovery observation.A summary of the follow-up observations used in the timing analysis is given in Table 2.The pulse time of arrivals (ToAs) was measured for each observation by cross-correlating the pulse profile against a high-S/N template, which was obtained by fitting a set of Gaussian curves to the best detection.We carried out the subsequent timing analysis of the ToAs with the TEMPO pulsar timing package (Nice et al. 2015) 1 to obtain more accurate measurements of the timing parameters.DM was measured using ToAs from multifrequency subbands.As shown in Table 2, the gap between some observations can be 1 month or even much longer, we thus used the DRACULA script (Freire & Ridolfi 2018) 2 to help derive phase-connected timing solution for the pulsar.

Detection of radio pulses from PSR J2016+3711 and timing
We did not detect radio pulses toward P G22.7 3 but detected radio pulses toward P CTB87 (Fig. 2), with a significance of ∼ 10.8 and a reduced  2 of 4.3.We will refer to this pulsar as PSR J2016+3711 hereafter.The period () of the detected pulses in the 2019 observation is 50.806263(65)ms, and the DM of J2016+3711 is 429.5(5.0)pc cm −3 .We note that the period is somewhat similar to that estimated from the X-ray observation of the host PWN (e.g. ≈ 80 ms in Matheson et al. 2013 or 65 ms in Guest et al. 2020).
Although unlikely, the radio pulses could potentially arise from sideband effect of the single-dish telescope.To clarify if the detected radio pulses are a result of the effect, we searched for published known pulsars in a region centred at position P CTB87 with an angular radius of 5 • using the ATNF database (Manchester et al. 2005) 4 , the Pulsar ALFA Survey Project5 , the FAST Galactic Plane Pulsar Snapshot survey (Han et al. 2021; Zhou et al. 2023) 6 , and the Commensal Radio Astronomy FAST Survey7 .With the DM limited in range 429.5±100 pc cm −3 , the results are shown in Table 3.There is one pulsar, J2022+3842, which has similar DM and period as those of the radio pulses we detected.We folded our data with the period of J2022+3842, yet no radio pulse was detected.Therefore, the radio pulsation we detected is from the newly discovered pulsar at position P CTB87 in SNR CTB 87, i.e., J2016+3711.
In Table 4, we present the timing solution for J2016+3711.We only obtained a timing solution with several JUMPs left (see Table 2).One possible explanation is the discontinuity in the pulsar's frequency or its differential, e.g.glitches, which are commonly observed in young pulsars (Lyne & Graham-Smith 2012).The fitted position of J2016+3711 is spatially coincident with the X-ray point source detected by Chandra (Matheson et al. 2013).The postfit timing residuals with all of the ToAs obtained by the timing solution are shown in Fig. 3.

The properties of the integrated radio pulse profile
The integrated radio pulse profiles are unique signatures that differ from pulsar to pulsar (e.g.Lyne & Graham-Smith 2012).The integrated radio pulse of the pulsar in CTB 87, J2016+3711, can be well described by a single component with a Gaussian profile (Fig. 4).The width of the profile of a pulsar is usually measured at 10% and 50%  , where  S/N = 10.8 is the signal-to-noise ratio of the radio pulse.The 50% duty cycle of the radio pulse is about (7.8 ± 0.8)%.
Also, we can estimate the effective pulse width  eff , which is defined as the width of a boxcar-like pulse with the same energy and amplitude of the real pulse (e.g.Kramer et al. 1998).The value of  eff is estimated to be 32.2 • ±1.7 • , which is equivalent to 4.5±0.2ms in units of time.

The flux density of PSR J2016+3711
The flux density of the radio pulses at 1.25 GHz can be estimated as (e.g.Dewey et al. 1985) where  ≈ 1.5 is a factor taking into account the losses and system imperfections,   = 2 the polarisation channels, Δ = 3.4 × 10 8 Hz the effective bandwidth of the data,  int = 3470 s the on-source integrating time,  the radio pulsation period of the pulsar,  sys ≈ 25 K (Jiang et al. 2019) the system temperature of the telescope,  sky ∼ 50 K (Kothes et al. 2020) the averaged background temperature of the sky toward the pulsar at ∼ 1.25 GHz, and  ≈ 16 K/Jy (Jiang et al. 2019) the gain of the telescope.The flux density at 1.25 GHz is estimated to be 15.5 ± 0.7Jy, which indicates the faint nature of the pulsar's radio emission and explains why it has remained undetected in early surveys.Notably, there are only 91 pulsars (< 3%) in the ATNF catalog (v1.71) (Manchester et al. 2005) with a flux density less than 16 Jy at 1.4 GHz, among which over 75% were discovered by FAST (Han et al. 2021).
The flux densities of pulsars are known to have steep spectra (e.g.Sieber 1973), with   ∝  −  , where   is the flux density of the pulsar at frequency , and  the spectral index.Here we crudely constrain the spectral index of PSR J2016+3711, based on the estimated flux density at 1.25 GHz and the non-detection at lower frequencies in early surveys.
The pulsar candidate in SNR CTB 87 had been searched using the Low Frequency Array (LOFAR) at about 150 MHz with a sensitivity of about 0.4 mJy (Straal & van Leeuwen 2019), the Arecibo radio telescope at 430 MHz with a sensitivity of 0.05 mJy (Gorham et al. 1996), and the 76 m Lovell radio telescope at 606 MHz with a sensitivity of 3 mJy (Biggs & Lyne 1996) and 0.16 mJy (Lorimer et al. 1998), yet no radio pulse had been detected.The 5 values of the sensitivities are here used as the upper limits of the flux densities.4 as a function of the observation date for PSR J2016+3711. a range of ∼ 1.4-1.8(e.g.Lorimer et al. 1995;Toscano et al. 1998;Maron et al. 2000;Bates et al. 2013).
It is noteworthy that the confinement is valid only if the spectral property of the pulsar can be well described by a single power law.Actually, many young pulsars have broken power-law spectra, and some even show a turn-over at low frequencies (e.g.Bilous et al. 2020;Jankowski et al. 2018;Murphy et al. 2017).If there is a similar case for pulsar J2016+3711, the spectral property of it would be more complex than derived above.

The properties of PSR J2016+3711
The observed period  ≈ 50.81 ms and first period derivative  ≈ 7.25×10 −14 s s −1 are crucial to derive some basic physical properties of PSR J2016+3711.The characteristic age of the pulsar can be obtained through  c = /[( − 1) ] ≈ 11.1 kyr for a dipole radiator with the braking index  = 3.Interestingly, the previous inferences using the empirical relation between the X-ray luminosity of PWNe  The signal can be seen in every observation by folding the data with the timing solution.
* The JUMP positions.2020) with a distance 6.1 kpc adopted, are very close to this estimate.The strength of the equatorial surface dipole magnetic field, can be inferred as (Gaensler & Slane 2006) G, where  45 is the moment of inertia of the pulsar in units of 10 45 g cm −2 , and  10 is the radius of the pulsar in units of 10 km.Also, the spin-down luminosity of the pulsar is estimated to be  = 4 2  45 / 3 ≈ 2.2 × 10 37 erg s −1 , which can be used to explore the PWN scenario for the observed TeV gamma-ray emission (Abeysekara et al. 2018).
Using DM = 428.0pc cm −3 and the electron density model (Yao et al. 2017), the distance to the pulsar is estimated as ∼ 13.3 kpc, which is about twice the estimate obtained from extinction-distance relation (6.1±0.9 kpc) (Kothes et al. 2003).This may be because the line of sight will penetrate the Orion spur, along which many objects (e.g.H II regions) should contribute significant DMs.Furthermore, we note that the distances to many pulsars estimated from the electron density model could be overestimated for similar reasons (e.g.Yao et al. 2017;Han et al. 2021).

Searching for a gamma-ray counterpart of the radio pulsation
Pulsars have long been suggested to emit gamma-rays (e.g.Romani & Yadigaroglu 1995).Actually, about 294 pulsars have been detected with pulsed gamma-ray emission to date (Smith et al. 2023).
To search for the possible gamma-ray pulsation, we collected the Fermi-LAT Pass 8 data toward a region that is centred at the pulsar with a radius of 1 • , from 2021-08-06 00:00:00 (UTC) to 2023-09-09 00:00:00 (UTC).The data collected span over two years and cover the same MJD range as the timing analysis.We analysed the data with the standard software Fermitools8 version 2.0.8 released on 2021 January 20.We only selected events with zenith angle < 90 • to filter out the background gamma-rays from the Earth's limb.To optimise the signal-to-noise ratio, the analysis was restricted to the energy range 500 MeV-10 GeV.We folded the Fermi-LAT data with the ephemeris (parameter file) obtained from the radio timing results using the TEMPO2 software (Edwards et al. 2006) and checked the results using H-test.Yet, no clear signature of pulsed profile was found.We also performed the gtpsearch routine of Fermi Science Tools software packages (v10r0p5) to search for pulsation frequencies near the radio pulsation frequency.We adopted the spin frequency, the first and second spin frequency derivatives at 58871 MJD listed in Table 4 for the search.Still, there was no gamma-ray pulsation found.
Nonetheless, the existence of pulsed gamma-ray emission could not thus be excluded because of the limited statistics and the lack of a precisely determined pulsar ephemeris.More follow-up radio observations spanning over years would be helpful to obtain a more accurate timing solution, which could then be used to fold the gammaray data and search for the pulsation.

Orientation of the emitting regions
The Chandra X-ray observation of the CTB87 PWN displays a jets+torus structure (Matheson et al. 2013).The projected elliptical shape of the torus suggests an inclination angle ∼ 55 • for the equatorial wind.Similar to that in the G54.1+0.3PWN, the torus is one-sidedly brightened and can be explained with a Doppler boosting of sub-relativistic downstream flow (Lu et al. 2002).The jets, aligned with the pulsar spin axis, are at an angle ∼ 55 • with the line of sight.
The radio pulse profile of PSR J2016+3711 is narrow ( 50 ≈ 28 • ,  eff ≈ 32 • ), without broad wings, which is different from the pulse profile of PSR J1930+1852 in the G54.1+0.3PWN with broad wings (Camilo et al. 2002).This may indicate that either the pulsar radio beam starting near the magnetic polar cap is intrinsically narrow, or the line of sight sweeps just across a small segment of a broad beam.But, as a reference, the estimates of the beam radius from the relation obtained by Kuz'min et al. (1984)  = 5 •  −1/2 ∼ 22 • or the relation for the case of magnetic inclination angle close to 90 • (Lyne & Manchester 1988)  = 6.5 •  −1/3 ∼ 17 • appear not in favor of the broad beam possibility.The detection of the radio pulsation of PSR J2016+3711 implies no small magnetic inclination angle (between the pulsar's spin, defined by the jets, and the rotating magnetic dipole axis, surrounded by the radio beam that sweeps past the Earth).If the gamma-ray pulsation is indeed undetectable, there could be two possibilities.One could be that the pairs are not accelerated to the gamma-ray emitting energies.The other could be that the direction angles of gamma-ray emission from the outer magnetosphere are offset from the magnetic dipole axis, avoiding the Earth.

SUMMARY
We have performed a FAST observation to search for radio-pulsed signals toward two X-ray point-like sources in SNR CTB 87 and proximate to SNR G22.7−0.02,respectively.We discovered radio pulses from the former with a significance of ∼ 10.8, which is thus revealed to be a pulsar and named J2016+3711.This undoubtedly confirms the compact nature of the X-ray point source in SNR CTB 87.This is the first pulsar discovered in SNRs using FAST thanks to its unprecedented sensitivities.The integrated profile of the pulses can be well described by a single component, with a  50 value of about 28.1 • .Combining the estimated mean flux density at 1.25 GHz (≈ 15.5Jy) with the non-detection of the radio pulse at lower frequencies, the spectral index of the pulsar is constrained to be ≲ 2.3.In the 28 follow-up timing observations with FAST, we identified PSR J2016+3711 and measured a period of 50.8087377(1) ms, period derivative of 7.2 × 10 −14 s s −1 , and DM of 428.0(1) pc cm −3 .Additionally, the strength of the equatorial surface magnetic dipole magnetic field, characteristic age and spin down luminosity of the pulsar are derived to be about 1.9 ×10 12 G, 11.1 kyr, and 2.2 ×10 37 erg s −1 , respectively.Based on the radio timing results, we also searched for gamma-ray pulsation with the Fermi-LAT data but found no pulsation.from NSFC under grant 12273010.YLY acknowledges support from CAS "Light of West China" Program.This work made use of data from FAST, a Chinese national mega-science facility built and operated by the National Astronomical Observatories, Chinese Academy of Sciences.

Figure 1 .
Figure 1.(Left) Dual-wavelength image of SNR CTB 87: NVSS 1.4 GHz radio continuum image in red, and Chandra X-ray (0.5-7 keV) image in blue.The white cross indicates the location of the X-ray point-like source (Matheson et al. 2013).The white circle shows the beam size of FAST at 1.25 GHz toward position P CTB87 .(Right) Tri-wavelength image of the southwestern part of SNR G22.7−0.02:VGPS 1.4 GHz radio continuum image in red, the XMM-Newton X-ray (0.5-6 keV) image in blue, and the H.E.S.S. gamma-ray (> 1 TeV) integral flux image in green.The white cross and circle are similar to those in the left panel, but for position P G22.7 .The radio continuum images in both panels and the gamma-ray image in the right panel are re-gridded with a 4 ′′ × 4 ′′ resolution.

Fig. 5
Fig.5shows the fitting of the flux density of PSR J2016+3711 using   = 15.5(/1.25 GHz) −  Jy.It is shown that the spectral index of the pulsar should not be larger than ∼ 2.3.It appears not to contradict the typical spectral indices of radio pulsars, which are in

Figure 2 .Figure 3 .
Figure 2. Radio pulses detected toward position P CTB 87 in SNR CTB 87.The pulsed signal is present throughout the entire 1 hour of the observation within the bandwidth of the data (1.05-1.45GHz).The DM is 429.5(5.0)pc cm −3 and the period of the radio pulses is 50.806263(65)ms.

Figure 4 .
Figure 4. Integrated pulse profile at 1.25 GHz obtained from the FAST observation.The blue dashed line represents the single component fitted with a Gaussian curve.

Figure 5 .
Figure 5. Flux density of PSR J2016+3711 at 150 MHz (Straal & van Leeuwen 2019), 430 MHz (Gorham et al. 1996), 606 MHz (Lorimer et al. 1998) and 1250 MHz.The 5 values of the flux densities are shown as upper limits at 150 MHz, 430 MHz, and 606 MHz.The solid line represents the flux density distribution with spectral index  equal to 2.3.

Table 1 .
Parameters of the targets Kothes et al. 2003   Su et al. 2014of the peak flux density, as denoted by  10 and  50 , respectively.Because the significance of the pulses detected here is low, we only estimate the value of  50 , which is about 28.1 • , with an uncertainty of  50 =  50 /( √ 2 ln(2) S/N ) = 2.7 •

Table 2 .
Summary of the follow-up observations made with FAST for PSR J2016+3711

Table 3 .
Known pulsars near the P CTB87 position in 5 • region