Rotational Evolution of the Slowest Radio Pulsar, PSR J0250+5854

More than 2600 pulsars have been observed since the first discovery of a pulsar⁷ (Hewish et al. 1968; Manchester et al. 2005). They are divided into different groups according to their observational properties (Harding 2013). The evolutionary connections between different groups of neutron stars are widely discussed with different ways. The rotating radio transients are scattered throughout the range of normal pulsars on the P–$\dot{P}$ diagram, but with relatively larger characteristic ages. They are taken as old pulsars due to the extreme nulling behavior in terms of the radiation (Zhang et al. 2007). The high-energy bursts of young high magnetic field pulsars look very similar to the emission of magnetars (Gavriil et al. 2008; Göğüş et al. 2016). Possible evolutionary routes from high-$\dot{P}$ rotation-powered pulsars to magnetars are studied by modifying braking mechanisms (Espinoza et al. 2011; Liu et al. 2014). X-ray dim isolated neutron stars (XDINSs for short) are located just below the magnetars on the P–$\dot{P}$ diagram and are identified as old magnetars according to the studies of magnetothermal evolution (Viganò et al. 2013).

Our traditional understandings about the relationship between different groups of pulsars may be updated by the discovery of long-period pulsars. For PSR J0250+5854, the slowest rotating radio pulsar detected to date (Tan et al. 2018), its period and period derivative are 23.5 s and $2.7\times {10}^{-14}\,{\rm{s}}/{\rm{s}}$, which derives a dipole-magnetic field of 5.1 × 10¹³ G at the pole and a characteristic age of 13.7 Myr. On the one hand, PSR J0250+5854 is located under the conventional definition of a death line where radio emission is expected to stop, which challenges traditional knowledge about pulsar radio emission. On the other hand, though the spin period of PSR J0250+5854 is similar to that of magnetars or XDINSs, it is hard to make a connection between them due to the lack of detection of high-energy emission.

To investigate the link between PSR J0250+5854 and other neutron stars, we respectively apply spin-down models of magnetospheric evolution and magnetic field decay to simulate its possible rotational evolutions and to explain its unusual observations. Besides, possible observational properties under these two cases are predicted.

2.1. Spin-down Torque

It is generally accepted that a neutron star should be surrounded by a plasma-filled magnetosphere. There are particles accelerated in the magnetosphere. The radiation of particles will generate the observed pulsar emission. The out-flowing particle winds will take away the rotational energy to spin down the pulsar. Meanwhile, the presence of plasma in the magnetosphere will also generate a torque to align the magnetic axis and the rotational axis (Michel & Goldwire 1970; Philippov et al. 2014). For an oblique rotator, the evolution equation is (Michel & Goldwire 1970; Philippov et al. 2014)

$\begin{eqnarray}&&I\displaystyle \frac{d{\rm{\Omega }}}{{dt}}={\boldsymbol{K}},\end{eqnarray} \tag{ 1 }$

where I is the moment of inertia and a typical value of 10⁴⁵ g cm² is taken in the following calculations, Ω is the angular velocity, and ${\boldsymbol{K}}$ is the torque working on the pulsar. For a spherical system, Equation (1) can be expressed as

$\begin{eqnarray}&&I\displaystyle \frac{d{\rm{\Omega }}}{{dt}}={K}_{\mathrm{spin} \mbox{-} \mathrm{down}},\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&I{\rm{\Omega }}\displaystyle \frac{d\alpha }{{dt}}=-{K}_{\mathrm{alignment}},\end{eqnarray} \tag{ 3 }$

where α is the angle between the magnetic axis and the rotational axis (i.e., the inclination angle), ${K}_{\mathrm{spin} \mbox{-} \mathrm{down}}$ is the torque to spin down the pulsar, and K_alignment is the torque to align the pulsar.

According to the magnetohydrodynamic (MHD) simulations of pulsar magnetospheres (Spitkovsky 2006), the spin-down torque can be expressed as

$\begin{eqnarray}&&{K}_{\mathrm{spin} \mbox{-} \mathrm{down}}=-{k}_{0}\displaystyle \frac{{\mu }^{2}{{\rm{\Omega }}}^{3}}{{c}^{3}}({\sin }^{2}\alpha +{k}_{1}),\end{eqnarray} \tag{ 4 }$

where $\mu =(1/2){{BR}}^{3}$ is the magnetic dipole moment (B is the polar magnetic field and R is the radius of the pulsar), and c is the speed of light. Generally, k₀ is a numerical factor and k₁ represents the effect of out-flowing particles and is model dependent. Comparisons between different MHD simulations were discussed in Tong & Kou (2017). Here, we take k₀ ≈ k₁ ≈ 1 for simplicity. As pulsar is spinning down, the ability of pair production will decrease (Ruderman & Sutherland 1975). For pulsars near the death line, the effect of death must be considered (Zhang et al. 2000; Contopoulos & Spitkovsky 2006). Contopoulos & Spitkovsky (2006) applied the effect of pulsar death to simulate the rotational evolution of pulsars. We employ their treatment of pulsar death, and Equation (2) could be expressed as

$\begin{eqnarray}I\displaystyle \frac{d{\rm{\Omega }}}{{dt}}=-\displaystyle \frac{{\mu }^{2}}{{c}^{3}}{{\rm{\Omega }}}^{3}\left\{\begin{array}{ll}{\sin }^{2}\alpha +\left(1-\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{death}}}{{\rm{\Omega }}}\right) & \mathrm{if}\ {\rm{\Omega }}\gt {{\rm{\Omega }}}_{\mathrm{death}}\\ {\sin }^{2}\alpha & \mathrm{if}\ {\rm{\Omega }}\leqslant {{\rm{\Omega }}}_{\mathrm{death}}\end{array}\right.,\end{eqnarray} \tag{ 5 }$

where ${{\rm{\Omega }}}_{\mathrm{death}}=2\pi /{P}_{\mathrm{death}}$ is the angular velocity when a pulsar dies, and the death period is defined as (Contopoulos & Spitkovsky 2006; Tong & Xu 2012)

$\begin{eqnarray}&&{P}_{\mathrm{death}}=0.885{\left(\displaystyle \frac{B}{{10}^{12}{\rm{G}}}\right)}^{1/2}{\left(\displaystyle \frac{{V}_{\mathrm{gap}}}{{10}^{13}{\rm{V}}}\right)}^{-1/2}\,{\rm{s}},\end{eqnarray} \tag{ 6 }$

where V_gap is the potential drop of the acceleration gap. Equation (5) means that the pulsar is braked by the combination of magneto-dipole radiation and particle wind before its death, but only brakes from the magneto-dipole radiation after its death. Compared with the previous work of Contopoulos & Spitkovsky (2006), a factor of ${\cos }^{2}\alpha$ is omitted here. As discussed in Tong & Kou (2017), it is a weighting factor between the magnetic dipole radiation and the magnetospheric particles. Besides, considering the result of the magnetospheric simulations (Spitkovsky 2006; Li et al. 2012), the ${\cos }^{2}\alpha$ factor may not appear in the particle wind component.

As the two components of the spin-down torque are independent of the inclination angle when α = 0° and α = 90°, the alignment torque and the spin-down torque can be related as ${K}_{\mathrm{alignment}}\,=[{K}_{\mathrm{spin} \mbox{-} \mathrm{down}}(0^\circ )-{K}_{\mathrm{spin} \mbox{-} \mathrm{down}}(90^\circ )]\sin \alpha \cos \alpha$ (Philippov et al. 2014). Hence, Equation (3) could be expressed as

$\begin{eqnarray}&&I{\rm{\Omega }}\displaystyle \frac{d\alpha }{{dt}}=-\displaystyle \frac{{\mu }^{2}}{{c}^{3}}{{\rm{\Omega }}}^{3}\sin \alpha \cos \alpha .\end{eqnarray} \tag{ 7 }$

2.2. Spin Down of PSR J0250+5854

As PSR J0250+5854 is located just below the death line on the P–$\dot{P}$ diagram, its radio emission may tend to stop, and its period may be much close to the death period, ${P}_{\mathrm{obs}}\leqslant {P}_{\mathrm{death}}$. A magnetic field of $B\geqslant 1.276\times {10}^{12}{P}_{\mathrm{obs}}^{2}({V}_{\mathrm{gap}}/{10}^{13}\,{\rm{V}})\,{\rm{G}}$ can be calculated from Equation (6). Hence, an upper limit on the inclination angle of $\alpha =3\buildrel{\circ}\over{.} 4{({V}_{\mathrm{gap}}/{10}^{13}{\rm{V}})}^{-1}$ can be calculated from Equation (5). Generally, ${V}_{\mathrm{gap}}={10}^{13}\,{\rm{V}}$ is taken for normal pulsars (Contopoulos & Spitkovsky 2006), and the corresponding minimum magnetic field and maximum inclination angle are 7.06 × 10¹⁴ G and 34, respectively. Because PSR J0250+5854 is still a radio-loud pulsar, an inclination angle of 2° is assumed. The magnetic field could be calculated from Equations (5) and (6), which is about 7.1 × 10¹⁴ G.⁸

The coupled evolutions of inclination angle and rotation can be calculated from Equations (5) and (7). Choosing the current period and magnetic field, the evolutionary paths for different assumed ages are the same, but with different beginning points. To get an initial period smaller than 20 ms,⁹ an age of 2.7 × 10⁵ yr is assumed, and the corresponding initial period and inclination angle are 15 ms and 885, respectively. The evolutions of the inclination angle and rotation of PSR J0250+5854 are shown as solid lines in Figures 1 and 2. The solid triangle in Figure 1 represents the inclination angle of 2° at its age of 2.7 × 10⁵ yr. Meanwhile, the solid triangles in Figure 2 represent the periods and period derivatives at ages of 10 yr, 100 yr, 10³ yr, 10⁴ yr, 10⁵ yr, and 10⁶ yr, respectively. The death period and age are about 23.6 s and 2.1 × 10⁶ yr, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Different acceleration models yield the potential drops of the order of 10¹² V. A constant potential model of 3 × 10¹² V was taken for an approximate calculation (Ruderman & Sutherland 1975; Xu & Qiao 2001). Here, we discuss a relatively extreme case with V_gap = 10¹² V. Similarly, a minimum magnetic field of 7.06 × 10¹³ G and a maximum inclination angle of 36° can be calculated from Equations (5) and (6). Lower magnetic fields could be calculated by giving larger inclination angles. However, the evolutionary trend stays almost the same. Considering the alignment of the inclination angle, as well as the pulsar's position on the P–$\dot{P}$ diagram, a small inclination angle is assumed for PSR J0250+5854. A magnetic field of 10¹⁴ G will be calculated if an inclination angle of 7° is assumed. Different ages used just result in different starting points of the evolutionary line. Similarly, to get an initial period smaller than 20 ms, an age of 5.4 × 10⁶ yr is assumed here. Coupled evolutions of the rotation and the inclination angle can be calculated from Equations (5) and (7). The initial inclination angle and rotational period are about 878 and 18 ms, respectively. Evolutions of the inclination angle and rotation of PSR J0250+5854 are shown as dashed lines in Figures 1 and 2. The hollow triangle in Figure 1 represents its present inclination angle and age. The hollow triangles in Figure 2 are the periods and period derivatives at ages of 10 yr, 100 yr, 10³ yr, 10⁴ yr, 10⁵ yr, 10⁶ yr, 10⁷ yr, and 10⁸ yr, respectively. The death period and death age are about 28 s and 1.4 × 10⁸ yr, respectively.

As we can see from Figure 1, the inclination angle decreases as the pulsar ages, which means that the magnetic axis and the rotational axis tend to align. In the P–$\dot{P}$ diagram (Figure 2), the pulsar first evolves to the right under the combined effect of particle winds and magneto-dipole radiation, and then turns down under the effect of pulsar "death." It is possible to predict that PSR J0250+5854 should be an old high magnetic field pulsar or magnetar, which is on the death edge. As the magnetospheric particles are exhausted, its radio radiation will tend to stop. That may be why nulling pulses are observed (Zhang et al. 2007; Kou & Tong 2015). Besides, it may have a relatively small inclination angle due to the alignment of the rotational and magnetic axes (Philippov et al. 2014; Tong & Kou 2017). After the death point, the pulsar will only be spun-down by magneto-dipole radiation, and the effective magnetic field should be $B\sin \alpha$. From the calculation results, the magnetic field of PSR J0250+5854 is still high, of the order of 10¹⁴ G. Its true age is much smaller than its characteristic age. It is still possible to observe magnetar-like emissions (Perna & Pons 2011).

Besides, as we can see from evolution tracks in Figure 2, the rotational period increases slowly after the death point because of the decrease in spin-down torque. The period derivative at the death point is ${\dot{P}}_{\mathrm{death}}=5\times {10}^{-16}{({P}_{\mathrm{death}}/1{\rm{s}})}^{3}{({V}_{\mathrm{gap}}/{10}^{13}{\rm{V}})}^{2}{\sin }^{2}\alpha$ (Contopoulos & Spitkovsky 2006; Tong & Xu 2012). For a pulsar with B = 10¹⁴ G and ${V}_{\mathrm{gap}}={10}^{13}\,{\rm{V}}$, its maximum inclination angle is 34, and its expected period and period derivative are 8.8 s and 1.2 × 10⁻¹⁵ s/s, respectively. We could hence predict that, for similar long-period pulsars, most of them should have a relatively small period derivative, i.e., <10⁻¹⁵ s/s.

In this case, rotational evolution of the pulsar under the effect of magnetic field decay will be discussed. We also take the same assumption with k₀ ≈ k₁ ≈ 1 in the spin-down torque, and pulsar rotational evolution should be (Spitkovsky 2006)

$\begin{eqnarray}&&I\dot{{\rm{\Omega }}}=-\displaystyle \frac{{\mu }^{2}}{{c}^{3}}{{\rm{\Omega }}}^{3}({\sin }^{2}\alpha +1)=-\displaystyle \frac{{B}^{2}{R}^{6}}{4{c}^{3}}{{\rm{\Omega }}}^{3}({\sin }^{2}\alpha +1).\end{eqnarray} \tag{ 8 }$

Because magnetars are generally taken to be powered by a high magnetic field, magnetic field decay is naturally and widely discussed (Goldreich & Reisenegger 1992; Heyl & Kulkarni 1998; Viganò et al. 2013; Igoshev & Popov 2018). For an isolated neutron star, numerical calculations of magnetic field decay can be generally expressed in a simple power law with decay index β, $\tfrac{\partial B}{\partial t}\propto -{B}^{1+\beta }$ (Colpi et al. 2000; Dall'Osso et al. 2012). The resulting changes are only quantitative for different forms of magnetic field decay. The magnetic energy may be dissipated by Hall/ohmic diffusion in the stellar crust, and an analytical decay could be expressed as (Aguilera et al. 2008; Popov & Turolla 2012)

$\begin{eqnarray}&&B(t)=\displaystyle \frac{{B}_{0}\exp (-t/{\tau }_{{\rm{o}}})}{1+({\tau }_{{\rm{o}}}/{\tau }_{{\rm{H}}})[1-\exp (-t/{\tau }_{{\rm{o}}})]}+{B}_{\mathrm{fin}},\end{eqnarray} \tag{ 9 }$

where B₀ is the initial field, B_fin is the relic field, and τ_o and τ_H are the characteristic timescale of ohmic and Hall decay, respectively.

3.1. Spin Down of PSR J0250+5854

Assuming an initial magnetic field of B₀ = 4 × 10¹⁵ G, together with ${\tau }_{{\rm{o}}}={10}^{6}\,\mathrm{yr}$, ${\tau }_{{\rm{H}}}={10}^{3}\,\mathrm{yr}$, and ${B}_{\mathrm{fin}}=8\times {10}^{12}\,{\rm{G}}$, the magnetic field evolution can be calculated from Equation (9),¹⁰ which is shown in Figure 3.

Zoom In Zoom Out Reset image size

Given a typical inclination angle of 45°, together with an initial period of 0.01 s, the rotational evolution of PSR J0250+5854 could be calculated from Equations (8) and (9), which is shown in Figure 4. The solid triangles represent its periods and period derivatives at ages of 10 yr, 100 yr, 10³ yr, 10⁴ yr, 10⁵ yr, 10⁶ yr, 10⁷ yr and 10⁸ yr, respectively. The age and present magnetic field are about 1.4 × 10⁵ yr and 3.4 × 10¹³ G, respectively. As we can see from Figure 3, the magnetic field decays sharply at ages between 10³ and 10⁶ yr, which will result in a decrease in the spin-down torque. Correspondingly, the rotational evolution line quickly turns down after age about 10³ yr. As shown in Figure 4, the evolution line passes through the region of SGR 1806-20, a so-called high magnetic field magnetar with dipole-magnetic field about 4 × 10¹⁵ G at the pole (which is chosen as the initial magnetic field). It is possible to predict that the precursor of PSR J0250+5854 is a magnetar. Due to the decay of magnetic field, the magnetic-field-powered magnetospheric activity will become weaker and tend to stop. From the calculation results, it may be hard to observe magnetar-like burst (Perna & Pons 2011).

Zoom In Zoom Out Reset image size

As the magnetic field tends to a constant value after 10⁶ yr, the pulsar will spin down under a constant spin-down torque, and $\dot{P}=4.0\times {10}^{-14}/P$ in the case of α = 45°, according to Equation (8). In this case, for similar long-period pulsars, most of them should have relatively higher period derivatives, i.e., ∼10⁻¹⁵ s/s.

Several cases could be computed by varying the initial magnetic field, B₀; the Hall timescale, τ_H; as well as the inclination angle, α. The evolution tracks are roughly the same, and the results only change quantitatively. The aim of this section is to probe the possible evolution path of the long-period pulsar under the existing parameter space. Hence, a relatively high initial magnetic field is chosen here.

The pulsar will evolve along isodynamic lines under the pure magneto-dipole radiation model. It is hard to explain observations of the pulsar with such a hypothesis. There should not be period clustering among pulsars in this case. However, it would be different when considering physical braking torque. Colpi et al. (2000) discussed the period clustering of the anomalous X-ray pulsars by considering coupled evolution of the thermal and magnetic field. The pulsar clustering could also be identified under the two spin-down models discussed in this paper and in terms of the period derivative.

The pulsar death line is defined as a threshold when pair production ceases, and it is model dependent (Ruderman & Sutherland 1975; Chen & Ruderman 1993; Zhang et al. 2000; Zhou et al. 2017). For a specific pulsar, its death depends on its magnetic field and the potential drop of the acceleration gap (Equation (6)). Therefore, it is not strange to find pulsars under the death lines or outside the death valley. As discussed in Kou & Tong (2015), the particle density may decrease as pulsars spin down. Hence, for pulsars on the death edge, their radio emissions tend to stop, and they may be observed in abnormal ways, such as nulling, intermittent pulsar, or ever rotational radio transients (Zhang et al. 2007; Karako-Argaman et al. 2015).

The evolution of the inclination angle is also widely studied; statistical studies show that long-period pulsars have relatively small inclination angles (Lyne & Manchester 1988; Tauris & Manchester 1998). The same conclusion are performed by a statistical study of the pulse width of radio pulsars with interpulse emission (Maciesiak & Gil 2011). In a 3D magnetosphere, the magnetosphere may also generate a torque to align the rotational axis and magnetic axis. In other words, the inclination angle tends to decrease (Philippov et al. 2014; Tong & Kou 2017). Tong & Xu (2012) applied the spin-down model of Contopoulos & Spitkovsky (2006) to simulate the rotational evolution of SGR 0418+5729, which is a so-called low magnetic field magnetar. They pointed out that SGR 0418+5729 may be a normal magnetar but with a small inclination angle. From the point of view of inclination angle evolution, this is possible. PSRs J1119–6127 and PSR J1846–0258 are two young pulsars with high magnetic fields about 8 × 10¹³ G and 9.7 × 10¹³ G at the pole. Both of them are observed with magnetar-like bursts (Gavriil et al. 2008; Göğüş et al. 2016). Besides, the radio emission of PSR J1119–6127 disappeared during the bursts, but reappeared about two weeks after the bursts (Majid et al. 2017). Connections between high magnetic field pulsars and magnetars are discussed (Archibald et al. 2017). The wind braking model is applied to explain their observations and to simulate their long-term rotational evolution (Ou et al. 2016). From the point of view of magnetosphere evolution, PSR J0250+5854 may be a high magnetic field pulsar or a magnetar on the death edge; it has a small inclination angle but a strong magnetic field. Its true age may be much smaller than the characteristic age. Due to its relatively high magnetic field, it may be burst active but with a relatively longer waiting time (Perna & Pons 2011).

As the Hall and ohmic diffusion are both temperature dependent, the magnetothermal evolution is studied (Pons et al. 2009; Viganò et al. 2013). For PSR J0250+5854, an age of 1.4 × 10⁵ yr is calculated with an initial magnetic field of 4 × 10¹⁵ G. The expected thermal luminosity of PSR J0250+5854 is about 10³⁴ erg s⁻¹ via the magnetothermal evolution models of Viganò et al. (2013). However, the Swift/X-Ray Telescope nondetection places an upper limit on the luminosity of 1.5 × 10³² erg s⁻¹ on kT = 85 eV for ${N}_{H}=9\times {10}^{21}\,{\mathrm{cm}}^{2}$ (Tan et al. 2018). One possibility is that the evolution of PSR J0250+5854 is dominated by magnetospheric evolution. However, because the inner super-dense matter of neutron star is still a mystery, it is also possible that PSR J0250+5854 has undergone a rapid cooling process (Yakovlev et al. 2008).

In fact, both the effect of magnetospheric evolution and magnetic field decay may work on the pulsar; the dominant mechanism could also be identified by the period derivative and the distributions of long-period pulsars. In the case of magnetospheric evolution, the spin-down torque continues to decrease in the later stage, which will result in a continuous decrease in the period derivative (Figure 2). Hence, a relatively low period derivative will be measured—for example, the radio pulsar PSR J2144–3933, with a period of 8.5 s and a relatively low period derivative of 4.96 × 10⁻¹⁶ s/s (Young et al. 1999). However, in the case of magnetic field decay, the magnetic field tends to be constant at its late stage, and so does the spin-down torque. Pulsars will have a relatively high period derivative in this case. Clustering distribution of more similar long-period pulsars will help to distinguish these two cases.

A fallback disk model is also commonly applied to explain the observations of pulsars and pulsar-like objects. Due to the significant increase in the period during the propeller phase, the evolution of a long-period neutron star could be understood by considering fallback disk torque (Menou et al. 2001; Fu & Li 2013; Tong et al. 2016). After the propeller phase, the neutron star tends to be in rotational equilibrium with the fallback disk at the equilibrium period of ${P}_{\mathrm{eq}}=2\pi {\left(\tfrac{{R}_{{\rm{m}}}^{3}}{{GM}}\right)}^{1/2}$, where G is the gravitational constant and M is the mass of the neutron star.

We also calculated the possible evolution of PSR J0250+5854 under the fallback disk model of Tong et al. (2016). Assuming that PSR J0250+5854 is spinning at the rotational equilibrium phase, a relatively large age about 10⁷ yr could be estimated according to its period derivative: $\dot{P}\approx \tfrac{3\alpha P}{7t}$, where α is the power-law index of the mass accretion rate, and α = 1.25 is taken here. The model implies that the magnetic field of PSR J0250+5854 is about 10¹² G, which is much smaller than its dipole-magnetic field. Similarly, a lower magnetic field of about 10¹¹ G and a larger age of about 10⁸ yr are expected for PSR J2144–3933, because of its low period derivative. However, the typical lifetime of a fallback disk is generally taken as about a few thousand years (Menou et al. 2001). The model-expected period derivatives would be much larger than the observed values for these two long-period pulsars. One possibility is that the fallback disk around them is not active, and the pulsar is braked by magnetic dipole radiation at the moment; their ages and magnetic fields are much closer to their characteristic ages and characteristic magnetic fields. The pulsar will spin down along the isodynamic lines. Future measurements of the pulsar age and magnetic field would help to probe the possibility of braking under torques due to a fallback disk for these long-period pulsars.

In conclusion, PSR J0250+5854 may be an old, high magnetic field pulsar or a magnetar. Considering the magnetospheric evolution, PSR J0250+5854 may be on the death edge, it may have a small inclination angle and a relatively high magnetic field, and magnetar-like emissions may be observed. However, in the case of magnetic field decay, the magnetic-field-powered magnetar-like emission may be difficult to observe. Besides, from the point of view of evolution studied, pulsars will have low period derivatives in the case of magnetospheric evolution but have relatively higher period derivatives in the case of magnetic field decay. Clustering distribution of more similar long-period pulsars will also help to distinguish these two cases.

This work is supported by the National Key R&D Program of China under grant No. 2018YFA0404703 and the Open Project Program of the Key Laboratory of FAST, NAOC, Chinese Academy of Sciences, and the West Light Foundation of CAS (2018-XBQNXZ-B-023). H.T. is supported by the NSFC (11773008). R.X.X. is supported by the National Key R&D Program of China (grant No. 2017YFA0402602), the National Natural Science Foundation of China (grant Nos. 11673002 and U1531243), and the Strategic Priority Research Program of Chinese Academy Sciences (grant No. XDB23010200). X.Z. is supported by the National Natural Science Foundation of China (Nos. 11873040, 11373006, and U1838108). The FAST FELLOW-SHIP is supported by Special Funding for Advanced Users, budgeted and administrated by Center for Astronomical Mega-Science, Chinese Academy of Sciences (CAMS).