In-orbit demonstration of X-ray pulsar navigation with the Insight-HXMT satellite

TL;DR

This paper demonstrates the successful in-orbit use of X-ray pulsar navigation with the Insight-HXMT satellite, achieving precise orbit determination using a novel method and data from a single pulsar observation.

Contribution

It introduces the SEPO method for pulsar navigation and validates it through in-orbit demonstration with Insight-HXMT using only one pulsar.

Findings

Orbit determined within 10 km (3 sigma)

Velocity pinpointed within 10 m/s (3 sigma)

Successful in-orbit demonstration of X-ray pulsar navigation

Abstract

In this work, we report the in-orbit demonstration of X-ray pulsar navigation with Insight-Hard X-ray Modulation Telescope (Insight-HXMT), which was launched on Jun. 15th, 2017. The new pulsar navigation method 'Significance Enhancement of Pulse-profile with Orbit-dynamics' (SEPO) is adopted to determine the orbit with observations of only one pulsar. In this test, the Crab pulsar is chosen and observed by Insight-HXMT from Aug. 31th to Sept. 5th in 2017. Using the 5-day-long observation data, the orbit of Insight-HXMT is determined successfully with the three telescopes onboard - High Energy X-ray Telescope (HE), Medium Energy X-ray Telescope (ME) and Low Energy X-ray Telescope (LE) - respectively. Combining all the data, the position and velocity of the Insight-HXMT are pinpointed to within 10 km (3 sigma) and 10 m/s (3 sigma), respectively.

Tables

Table 1. Table 1: The pulsar parameters of the Crab pulsar

RA	${05}^h{34}^m{31}^s.973$
DEC	${22}^{\circ }{00}^{\prime }{52.06}^{\prime \prime }$
$\nu$	29.639022542326
$\dot{\nu }$	$-3.6867\times {10}^{-10}$
$\ddot{\nu }$	0
EPOCH	57992.16548081821

Table 2. Table 2: The orbit parameter spaces. ’ Δ ​ a Δ 𝑎 \Delta a ’, ’ Δ ​ e Δ 𝑒 \Delta e ’, ’ Δ ​ i Δ 𝑖 \Delta i ’, ’ Δ ​ ω Δ 𝜔 \Delta\omega ’, ’ Δ ​ Ω Δ Ω \Delta\Omega ’ and ’ Δ ​ w Δ 𝑤 \Delta w ’ are the deviations for semi-major axis, eccentricity, inclination angle, RAAN, argument of perigee, and mean anomaly, respectively.

Orbital element deviation	Min	Max	Step
$\Delta$ a (m)	-200	200	20
$\Delta$ e	-0.0018	0.0182	0.001
$\Delta$ i (∘)	-0.5	0.5	0.05
$\Delta$ $\omega$ (∘)	-1.0	1.0	0.1
$\Delta$ $\Omega$ (∘)	-1.0	1.0	0.1
$\Delta$ w (∘)	-1.0	1.0	0.1

Table 3. Table 3: Best estimated value of the elements and errors (3 σ 𝜎 \sigma )

Payload	Orbital element deviation	$\Delta$a	$\Delta$e	$\Delta$i	$\Delta \omega$	$\Delta \Omega$	$\Delta$w
		m	${10}^{-3}$	∘	∘	∘	∘
LE	deviation	7.02	-0.38	0.01	0.01	0.003	0.003
	error ($3\sigma$)	11.75	1.03	0.02	0.06	0.09	0.09
ME	deviation	3.07	-0.11	0.005	0.01	0.007	0.006
	error ($3\sigma$)	16.44	1.28	0.03	0.064	0.11	0.11
HE	deviation	4.89	-0.58	0.006	0.009	0.002	0.002
	error ($3\sigma$)	7.58	0.59	0.01	0.03	0.05	0.05
all detectors	deviation	5.25	-0.50	0.007	0.01	0.002	0.002
	error ($3\sigma$)	6.40	0.52	0.01	0.03	0.04	0.04

Table 4. Table 4: Best estimated value of Δ Δ \Delta a and errors (3 σ 𝜎 \sigma )

enlargement of $\delta$		0%	1%	5%	10%
$\Delta$a (m)	deviation	-3.15	-6.12	-5.4	-6.8
	error ($3\sigma$)	24.96	25.12	26.12	27.15

Equations

$${\chi }^2=\sum _{i=1}^n\frac{(P({\varphi }_i)-\bar{P}{)}^2}{\bar{P}},$$

$$\Phi -{\Phi }_0=\nu (t-t_e)+\frac{1}{2}\dot{\nu }{(t-t_e)}^2+\frac{1}{6}\ddot{\nu }{(t-t_e)}^2,$$

$${\sigma }_p=\frac{\frac{1}{2}W}{SNR},$$

$$SN{R}_{\mathrm{filtered}}=\frac{1000SNR}{1000+SNR}.$$

$${\chi }^2$$

$$=\sum _{i=1}^n\sum _{k=1}^m\sum _{k^{\prime }=1}^{m}p({\varphi }_i-\delta {\varphi }_k)p({\varphi }_i-\delta {\varphi }_{k^{\prime }})/\bar{P}-n\bar{P}.$$

$$\delta =\frac{w}{d},$$

Full text

stars: pulsar; techniques: miscellaneous

1 Introduction

Space navigation plays an important role for spacecraft launched to the earth orbit, lunar and Mars vicinities, and deep space. At present, the navigation of the spacecraft in low earth orbits mainly relies on global navigation satellite systems (GNSS) and the ground-based tracking systems, while the navigation in deep-space is mainly based on radio technologies (e.g. the US deep-space network). With the increasing number of space tasks, the ground stations are becoming overloaded and the operating costs are getting higher. Navigation also can not be handled timely, due to the communication delay or unexpected malfunction of ground stations. Pulsar navigation, an autonomous navigation technology, is receiving more and more attention as it is less dependent on the support of ground equipments and meets the continuous navigation requirements for space missions in different orbits.

The first pulsar was discovered in 1967 (Hewish et al., 1968), and up until today more than 2000 pulsars have been discovered with frequency covering from radio, infrared, optical, ultraviolet, X-ray to gamma-ray. They are called ‘celestial GPS satellites’ because of their long-term timing stability which is comparable to atomic clocks on the earth (Taylor, 1991). The pulsar navigation using the radio and X-ray pulsars was proposed in 1974 (Downs, 1974) and in 1981 (Chester & Butman, 1981), respectively. A large number of researches on the theory, algorithms and simulations have been carried out for spacecrafts in earth orbits, Mars orbits and deep space orbits (e.g. Sheikh et al (2007); Hanson et al (2008); Emadzadeh & Speyer (2011); Becker et al. (2013); Wei et al. (2013); Wang et al. (2014); Shemar et al. (2016); Cui et al. (2016)), and the tracking phase or positioning errors for different pulsars are predicted in detail (Sheikh et al, 2007; Hanson et al, 2008). In 1999, the first in-orbit test was carried out with the Unconventional Stellar Aspect (USA) on the ARGOS satellite (Wood, 1993; Sheikh, 2005). In 2016, the results of testing pulsar navigation with Gamma-Ray Bursts Polarimeter (POLAR), which was launched onboard the Chinese space laboratory TG-2 on Sept. 15th, 2016, were reported (Zheng et al., 2017). With the effective area of $\sim$ 200 cm2 and wide field of view of more than 2$\pi$ Sr (Produit et al., 2005; Li et al., 2018), POLAR has monitored the Crab pulsar for a long time. With 31-day-long observations of the Crab pulsar, the TG-2 orbit was determined successfully. The orbital elements were determined with the orbit deviation within 20 km (1 $\sigma$). On Jun. 3th, 2017, the Neutron Star Interior Composition Explorer (NICER) was launched on the International Space Station (ISS). It comprises 56 identical X-ray telescopes with effective area of 2000  cm2 in total, each of the telescopes consisting of a concentrating X-ray optic and a single-pixel silicon drift detector (SDD) (Paul et al., 2017). By measuring the Time of Arrival (TOAs) of five millisecond pulsars, NICER pinpointed its location within 5 km (Alexandra, 2017).

On Jun. 15th, 2017, Insight-Hard X-ray Modulation Telescope (Insight-HXMT) was launched in China (Zhang et al., 2018; Li et al., 2018). It consists of three X-ray slat-collimated telescopes: High Energy X-ray Telescope (HE), Medium Energy X-ray Telescope (ME) and Low Energy X-ray Telescope (LE). HE contains 18 cylindrical NaI(Tl)/CsI(Na) Phoswich detectors with the energy band 20-250 keV and a total geometrical area of about 5000 cm2. ME is composed of 1728 Si-PIN detectors with the energy band 5-30 keV and a total geometrical area of about 952 cm2. LE is composed of 96 Swept Charge Device (SCD) detectors which are sensitive in the 1-15 keV with a total geometrical area of about 384 cm2. Due to the wide energy band and large detection area, Insight-HXMT can also be used to make in-orbit test of pulsar navigation. In this paper, the new pulsar navigation method ‘Significance Enhancement of Pulse-profile with Orbit-dynamics’ (SEPO), which has been tested in POLAR (Zheng et al., 2017), is adopted to determine the orbit of Insight-HXMT autonomously. The Kalman Filter method can also be used to make in-orbit demonstration of X-ray pulsar navigation (Wang et al., 2018).

2 Navigation method

Conventionally, at least three pulsars are needed in pulsar navigation models to estimate the position of a spacecraft. However, it results in complexity of the payload and poses risks for control system. The one-pulsar navigation method SEPO has been proposed by combining the significance analysis of pulse profile and orbit dynamics. The orbital dynamics model can produce a high-precision orbit forecast in a short period of time but with a long-term drift. The pulse profile will be ‘deformed’ due to the orbit deviation as explained below, resulting in the decrease of the significance of the profile signal, which is used to update the orbital parameters continuously.

2.1 The orbit dynamics

The orbit is described by six orbital elements: semi-major axis ($a$), eccentricity ($e$), inclination angle ($i$), right ascension of the ascending node (RAAN, $\omega$), argument of perigee ($\Omega$), and mean anomaly ($w$). Considering the altitude of the Insight-HXMT orbit (about 550 km) and the complexity of the earth’s perturbation force, the orbit forecast model is constructed with the 70th order gravity field, the atmospheric drag, and the solar radiation pressure, and the track integrator is Runge-Kutta-Fehlberg 7(8) with an integration step of 60 seconds (Wang & Gao, 2016). In such conditions, the position error is less than 3.5 km for the forecast orbit in five days. It could be further simplified according to the different conditions of a spacecraft.

2.2 The significance of the pulse profiles

To get the ‘standard’ pulse profile, the arrival time of each photon on the local site is corrected to the Solar System Barycenter (SSB) and folded with the pulsar ephemeris with certain bins (Ge et al., 2012). In this process, the orbit of the spacecraft plays an important role for the time correction. If the orbit deviates from the correct position, the calculated profile will be deformed from the ‘standard’ profile due to the wrong phase.

The significance of the pulse profile is defined as follows;

[TABLE]

where $P(\phi_{i})$ is the counts of the profile at $\phi_{i}$ , $\bar{P}$ is the mean counts of the profile and $n$ is the total bin number of the profile.

It is expected that the significance of the calculated profile varies with spacecraft orbit deviation. The more the orbit deviates from the real one, the more the profile will be deformed, and the less the significance. $\chi^{2}$ will reach the maximum for the zero deviation, i.e., the true orbit.

We make a simulation for five-day observations of the Crab pulsar with HE. The true orbit ‘ob0’ is given by the GPS on the satellite, and the orbit ‘ob1’ is calculated with the dynamics of orbits by changing the ‘semi-major axis’ with 200 m. As shown in Fig. 1, the deviation will be larger, to more than 50 km after 2 days. Using the predicted orbit ‘ob1’, the profile is deformed and the significance is decreased from $1.237\times 10^{6}$ to $1.226\times 10^{6}$, i.e. by about $1\%$.

2.3 The orbit fitting with grid search

Considering the correlations between the significance of the profile and the orbit deviation, the orbit determination algorithm is described by the processing flow in Fig. 2. In this algorithm, the orbit forecast and the profile’s significance analysis are combined. The predicted orbit is obtained with orbital dynamics model to complete the profile folding. At the same time, the long-term stability characteristics of the profile can be used to continuously update the orbital parameters, thus the long-term drift of the orbit model can be corrected.

On one hand, the orbit forecast model is constructed by considering the orbit dynamics as described in Section 2.1. Then, with the estimated orbit parameters $p_{0}$ ($a$, $e$, $i$, $\omega$, $\Omega$, $w$) at the current time ($t_{0}$) and the certain parameter spaces $\delta p_{k}\,(k=1,2,..N)$, the $k^{\rm th}$ predicted orbit is generated for the time span [$t_{s}$, $t_{0}$] ($t_{s}$ is the beginning time for the used historic observation data). On the other hand, the used data with the time span [$t_{s}$, $t_{0}$] are chosen and filtered according to the selection criteria as described in Section 3.

With the predicted orbit at different grid point, the arrival times of photons are corrected to SSB and folded to generate the ‘predicted’ profiles using phases calculated with the following equation (Ge et al., 2012);

[TABLE]

where $\Phi_{0}$, $\nu$, $\dot{\nu}$ and $\ddot{\nu}$ are the phase, frequency, frequency derivative and the second order derivative of frequency at the reference epoch $t_{e}$, respectively. Then, all the significances for different profiles are calculated. Finally, the ‘significance analysis’ is carried out, i.e. the trend between the significance of profiles and the orbit parameters is fitted with the Gaussian function and the optimal orbit parameters are obtained.

3 Data reduction

Almost five-day observation data of the Crab pulsar from 2017-08-31 10:00:00 to 2017-09-05 08:00:00 (UTC) are chosen, including all the data collected by HE, ME, and LE. The data are analyzed with the Insight-HXMT Data Analysis software (HXMTDAS) software v2.0111http://www.uu-world.cn/hxmt/hxmtsoft/index2.html. The following data selection criteria as recommended in the software are used:

(1) Pointing offset angles $\leq 0.05^{\circ}$;

(2) Elevation angles $\geq 12^{\circ}$;

(3) Geomagnetic cutoff rigidity $\geq 6$;

(4) Only the detectors with small FoVs selected.

In addition to the above criteria, intervals with high background induced by the charged particles or bright earth are excluded manually. Then we get the filtered data with the total exposure duration of 180 ks, 190 ks, and 120 ks for HE, ME, LE, respectively. The arrival times of the photons are corrected to SSB, then the pulsed profiles are folded for HE, ME, and LE. The pulsar ephemeris, shown in Table 1, is obtained from Insight-HXMT data with the timing residual of about 10 $\mu$s (1$\sigma$). As the $\ddot{\nu}$ parameter of the Crab pulsar is about $10^{-20}\,{\rm s}^{-3}$, the affected pulsar phase offset after 5 days is negligible ($\frac{1}{6}\cdot\ddot{\nu}\cdot dt^{3}=1.3\times 10^{-4}$), thus we set $\ddot{\nu}=0$.

4 Results

The orbital elements of Insight-HXMT are given at 07:59:00 (UTC) on Sept. 5th, 2017, and the space grids of parameters are chosen as listed in Table 2. For an easy display, only the deviations of the parameters ($a$, $e$, $i$, $\omega$, $\Omega$, $w$) from the true values ($6922.8781$ km, $0.00181017$, $42.9715^{\circ}$, $207.0229^{\circ}$, $116.9049^{\circ}$, $22.5215^{\circ}$) are listed.

First, the orbit determination method is demonstrated with LE. At each grid point of the parameter space, the ‘predicted’ profile is folded and its significance is obtained. As shown in Fig. 3, the significance of the profile varies significantly with the deviation of the orbital parameters and has a maximum around the true value (zero deviation). Then the maxima obtained by fitting a Gaussian function curve give the optimal values for the six orbital parameters separately, i.e. each time one parameter is free with the other five fixed, and the results are shown in Table 3. The process obtaining the best parameters is very similar to that used by Sheikh et al to get the TOAs from the observed pulse profiles (Sheikh, 2005). Therefore, the uncertainties of these parameter values can be obtained in a similar way as follows;

[TABLE]

where W is Full Width Half Maximum (FWHM), and SNR is the signal to noise ratio obtained by the bootstrap method (Diaconis & Efron, 1983). As mentioned by Sheikh (Sheikh, 2005), the SNR is limited to a maximum of 1000 by:

[TABLE]

Then similar demonstrations are also done for HE and ME, shown in Table 3. It shows that the errors for HE are less than for ME and LE, because it has collected more pulsed photons. Combining all the data from HE, ME and LE, the errors become smaller. Transfer the orbit elements to the cartesian parameters, we obtain the 3$\sigma$ errors of the current position and velocity, which are 3.15 km, 6.61 km, and 4.46 km for x, y, z and 0.0073 km s*-1*, 0.0033 km s*-1*, and 0.0042 km s*-1* for $v_{x}$, $v_{y}$, $v_{z}$, respectively. In short, the position and velocity are pinpointed within 10 km (3$\sigma$) and 10 m s*-1* (3$\sigma$).

5 Discussions

Conventionally, a minimum of three pulsars are required to get the absolute navigation. Considering the clock time-offset on the spacecraft, another pulsar is needed. Thus at least four pulsars should be observed simultaneously, which would require at least four detectors, increasing the technical complexity and the cost of mass, energy and money. An alternate choice is to observe different pulsars sequentially using a single detector, which would increase the risk of the control system. NICER has adopted the latter strategy. With the large effective areas and low background rates, it could get the TOAs of some millisecond pulsars in several kiloseconds. By observing five millisecond pulsars sequentially, the spacecraft position was determined under real flight conditions successfully (Alexandra, 2017).

The method SEPO, determining the orbit by using only one pulsar, has been proposed and demonstrated with POLAR (Zheng et al., 2017) and Insight-HXMT in this paper, respectively. As POLAR has an effective area of about 200 cm2 and high background rates (about 4000 counts s*-1*), the orbit was determined within 20 km (1 $\sigma$) by monitoring the Crab pulsar for 1 month. However, it should be pointed out that the POLAR’s response for photons varies with their incident angles, thus the obtained different profiles are mixed and become ‘broader’, decreasing the navigation precision which is shown in appendix A. For Insight-HXMT, with the pointed observation mode, the response to the incident pulsar emission is almost the same. With the observations of the Crab pulsar for 5 days, the position is determined within 10 km. So it is feasible to use one detector for navigation by monitoring one pulsar for a long time, which would result in low request of mass, energy, and less spacecraft control. It is favorable for the navigation applications, particular for the deep space navigation during the cruise phase of flight.

In addition, the standard pulse profile and parameters are important external inputs in X-ray pulsar navigation. However, it is hard to achieve the fixed ‘standard’ profile in space X-ray detection. For example, if the energy band or the response of the instrument is different, the ‘standard’ profiles become different. Furthermore, the response of the instrument may vary in the long-term operation in space, thus the profile varies too. If we still use the previous ‘standard’ profile, it would inevitably affect the navigation accuracy and even result in unpredictable systematic deviation. For SEPO, neither the standard pulse profile nor the continuous update of the pulsar parameters is needed. Therefore, it is no longer dependent on the standard pulse profile, nor affected by the pulsar observation conditions and the detector performance changes. It can be applied to different detectors and different mission scenarios.

For the in-orbit detections of a pulsar, the obtained profile could be regarded as the summation of all the profiles at different point of the orbit. That is, $P(\phi_{i})=\sum_{k=1}^{m}p_{k}(\phi_{i})=\sum_{k=1}^{m}p(\phi_{i}-\delta{\phi_{k}})$, where $\delta{\phi_{k}}$ represents the phase offset due to the position deviation. Thus, considering the significance of the profile defined in Sec. 2.2,

[TABLE]

As shown in equation 5, $\chi^{2}$ of the profile will reach the maximum with $\delta\phi_{k}=\delta\phi_{k^{\prime}}$, i.e. $\phi_{k}=const$. $\phi_{k}=0$ represents the true orbit with zero deviation. If $\phi_{k}=const$ (nonzero), it means that the orbit has fixed deviations from the real orbit. However, constrained by the orbital dynamics, there is few such completely ‘parallel’ tracks, thus ensuring to find the true orbit. For the cruising orbit in deep space, there indeed exist approximately parallel tracks. In such conditions, more measurements besides the significance of the profiles, TOAs for example, are needed to discriminate the true orbit from the others.

The proposed SEPO method, however, has not been proven mathematically. Thus some simulations have been performed and the detailed results are shown in the appendix, which shows that the method works with very different pulse profiles. In this in-orbit demonstration with the SEPO method, the position is pinpointed to within 10 km (3$\sigma$), while the theoretical bound is $\sim$ 250 meters (3$\sigma$) for 1920 m${}^{2}\cdot$s of LE (the effective area of 160 cm2 for the small FOV with the exposure duration of 120 ks) (Hanson et al, 2008). Possible reason are the effects of orbit forecast error ($<3.5$ km), absolute time accuracy of Insight-HXMT (Li et al., 2018), or the intrinsic timing noise of the Crab pulsar. The ephemeris used in this demonstration has an timing residual of $\sim$10 $\mu$s (1$\sigma$), resulting in position error of $\sim$9.0 km (3$\sigma$). See the detailed discussions on the effects of timing noise at the end of this section.

By now, the clock time-offset on the spacecraft has not been considered in the SEPO method. The short-term drift of clock within the integration time would deform the calculated profile, change the pulse significance and result in additional systematic errors. While the long-term drift just changes the absolute TOA, which has little effects on the pulse significance. Thus the clock time-offset can be taken into account and corrected with the significance analysis in future. In addition, the pointing error also decreases the significance of the profiles by changing the effective areas. However, the pointing error has no effects on the pulse phase of the pulsar, thus has little effect on the solution in the SEPO method, even if the pointing error has a dependence on orbit position. In this work, the observation of Insight-HXMT was performed in the ‘pointed mode’ with the pointing error $\leq 0.028^{\circ}$ (3$\sigma$). In this situation, the effective area changes by only about 0.5%, which can be ignored completely.

The timing noise of pulsars (Hobbs et al, 2010) has a great influence on the accuracy of pulsar navigation as shown above. The timing noise of a pulsar includes white noise component and red noise component. White noise can be suppressed by accumulating observation data over a long period of time. However, the red noise is difficult to suppress or eliminate. Actually, the effects of the red noise is quite similar to the clock time-offset and can be corrected as mentioned above. On the other hand, combining with the study of the mechanism of noise (such as magnetic field evolution), it is possible to give a more accurate prediction of the spin evolutions of pulsars, reducing the impact of timing noise (especially red noise) on the accuracy of pulsar navigation (Yi & Zhang, 2015; Gao et al, 2016; Zhang & Xie, 2012a, b). In addition, the glitch, that is, a sudden change in the rotation period, will also have to be considered in future.

6 Conclusion

A new navigation method SEPO has been proposed that combines the observed pulse profile with orbit dynamics. It has been demonstrated with the Insight-HXMT observations of only one pulsar (Crab) and the orbit has been determined successfully. Combining all the data, we obtain the best position estimations within 10 km (3$\sigma$) and velocity within 10 m s*-1* (3$\sigma$), respectively.

This work is supported by the National Key R&D Program of China (2016YFA0400800) and the National Natural Science Foundation of China under grants U1838101 and 11503027. This work made use of the data from the Insight-HXMT mission, a project funded by China National Space Administration (CNSA) and the Chinese Academy of Sciences (CAS). Appendix

Appendix A Simulations on the SEPO method

We perform some simulations to evaluate the feasibility of the new navigation method SEPO, especially for different pulse-profiles. The duty cycle of the profile $\delta$ is defined as follows;

[TABLE]

where $w$ is the width of the pulse (above 10% maximum of the peak), $d$ is the pulsar period. Thus $\delta$ describes the sharpness of a pulse profile.

First, three different pulse profiles from three pulsars, i.e., PSRs B1821-24, B1509-58, J1811-1925 and a supposed sinusoidal-like pulse shape are chosen for test, as shown in Fig. 4 (the left panel). To exclude the effects of other factors, the parameters such as position and period are fixed to be the same as the Crab pulsar. The pulse count rate is assumed to be 130 cnts$\cdot$ s*-1* and the background count rate to be 1200 cnts$\cdot$ s*-1*, as obtained from the LE data. Then the Insight-HXMT/LE observation data of the four pulsars are simulated with the true orbit of Insight-HXMT from Aug. 30th to Sept. 5th, 2019.

Then the method SEPO is carried out and the significance of the profiles with the deviations of the elements (semi-major axis for example) are shown in Fig. 4 (the right panel). It shows that the significance of all the four different pulse-profiles varies with $\Delta$a, and the optimal value can be found by Gaussian-fitting, which is near the true value (zero). The pulse-profiles do have effects on the position error. The sharper is the pulse-profile, the larger is the calculated significance, and the smaller is the position error. The duty cycle $\delta$ is 0.16, 0.55, 0.72, 0.80 for PSRs B1821-24, B1509-58, J1811-1925 and sinusoidal-like pulsar. The profiles like PSR B1821-24 are the best candidate due to its the navigation for the small duty cycle. For PSR J1811-1925, the significance trend is not smooth, however, the optimal value can still be obtained by fitting.

To estimate the effects of pulse shapes in details, a double-peak profile is generated with $\delta=0.16$, and $\delta$ is increased by 1%, 5%, and 10%, respectively. Then the observation data are simulated and the method SEPO is carried out as described above. The results are shown in Fig. 5, and the optimal semi-major axis and the errors ($3\sigma$) are shown in Table 4. It also shows that the significance decreases as the profile gets broader. With $\delta$ becoming larger by 10%, the resulted error of the semi-major axis is increased by 8.8%.