Grouping bi-chi-squared method for pulsar navigation experiment using observations of Rossi X-ray timing explorer

Hin SUN, Jinyu SU, Zhonwn DENG, Liron SHEN,*,Wimin BAO,c, Xiopin LI, Linshn LI, Zh SU, Wncon WANG

a School of Aerospace Science and Technology, Xidian University, Xi’an 710126, China

b Shaanxi Key Laboratory of Space Extreme Detection, Xi’an 710126, China

c Peng Cheng Laboratory, Shenzhen 518000, China

d Academy of Advanced Interdisciplinary Research, Xidian University, Xi’an 710126, China

e Beijing Institute of Control Engineering, Beijing 100190, China

f Institute of Space Electronic Information Technology, Xi’an 710100, China

g Shandong Institute of Aerospace Electronic Technology, Yantai 264032, China

KEYWORDS Model buildings;Navigation analysis;Navigation testing;Pulsars;RXTE observations

Abstract X-ray pulsar-based navigation is a revolutionary technology which is capable of providing the required navigation information in the solar system. Performing as an important pulsar-based navigation technique, the Significance Enhancement of Pulse-profile with Orbit-dynamics (SEPO)can estimate orbital elements by using the distortion of pulse profile.Based on the SEPO technique,we propose a grouping bi-chi-squared technique and adopt the observations of Rossi X-ray Timing Explorer(RXTE)to carry out experimental verification.The pulsar timing is refined to determine spin parameters of the Crab pulsar (PSR B0531+21) during the observation periods, and a short-term dynamic model for RXTE satellite is established by considering the effects of diverse perturbations.Experimental results suggest that the position and velocity errors are 16.3 km (3σ) and 13.3 m/s(3σ)with two adjacent observations split by one day(exposure time of 1.5 ks),outperforming those of the POLAR experiment whose results are 19.2 km(3σ)and 432 m/s(3σ).The approach provided is particularly applicable to the estimation of orbital elements via using high-flux observations.

1. Introduction

Serving as a complementary and an alternative technology to Deep Space Network (DSN), X-ray Pulsar Navigation(XPNAV) is a promising navigation method that enables advanced navigation capability for spacecrafts to travel through the solar system and beyond.1In recent decades,with the continuous advancement of X-ray detector technology,2,3an ever-growing interest has been focusing on XPNAV, most of which concentrates on perspectives of near-Earth and deep space applications. This trend becomes more apparent ever since the following X-ray observatories and XPNAV demonstration satellites have been launched. These satellites include X-ray pulsar navigation-I (XPNAV-1, November 10,2016),4,5Station Explorer X-ray Timing and Navigation Technology (SEXTANT, June 3, 2017),6as well as Hard X-ray Modulation Telescope (HXMT, June 15, 2017).7,8

With respect to XPNAV demonstration, existing experimental frameworks can be categorized into two types: realtime and non-real-time. The former depends upon measuring the time differences of pulse Time of Arrivals(TOAs)between the predicted and the observed TOAs at a given epoch and reference location.The latter relies on the analysis of pulse profile deformation. The SEXTANT mission is the world’s first realtime and autonomous XPNAV demonstration which was initialized by intentionally degrading orbit position and velocity.By processing the observations of multiple millisecond pulsars,a navigation accuracy within 10 km was realized compared with the position information supplied by on-board Global Positioning System (GPS) receiver. The SEXTANT experiment is a milestone in the XPNAV history, which lays a solid foundation for implementing autonomous XPNAV.9,10In 2017, Runnels and Gebre proposed a recursive range estimation using the Crab pulsar (PSR B0531+21) observations, in which an extended Kalman filter was adopted to estimate the time difference of pulse TOAs at the spacecraft and the origin.The estimation performance was validated through the flight data obtained from the X-ray observatory Suzaku, by which the estimated standard deviations of position reached 54 km(1σ).11The research group of the XPNAV-1 conducted a ground pulsar-based experiment using 85-day X-ray observations from the Time-resolved Soft X-ray Spectrometer(TSXS),and verified the feasibility of XPNAV using PSR B0531+21.Despite the small effective area of TSXS (merely 2.4 cm2@1.5 keV), it accumulated a large amount of continuous observations on the Crab pulsar. The ranging measurements relative to the Earth along the pulsar line-of-sight was introduced as control points into the orbit propagation of the XPNAV-1 for the states of the spacecraft cannot be completely observed by relying on only a single pulsar’s observations.The results have proven that the XPNAV-1 is capable of conducting self-location with an average navigation error of 38.4 km at the control points by observing the Crab pulsar for a period of nearly three months.6

Some similar ideas that have potential to perform flight experiments have also been developed for XPNAV, which include but are not limited to introducing pulsar phase and Doppler frequency estimation, or realizing an integrated navigation system by utilizing external observation information.Pulsar phase and Doppler frequency estimation method was introduced for XPNAV by using on-orbit epoch folding, and a linearized pulsar phase model was creatively put forward for a deep space explorer mission by dividing the received phase into predicted and estimated parts. This method solves the initial phase problem greatly well and provides theoretical foundation for practical applications.12,13Two-stage estimation method of Doppler frequency and phase delay was presented to overcome the Doppler effects,and the results turn out that within an observation interval, a constant Doppler frequency can be estimated to overcome the effect of the initial velocity errors.14Some advanced filtering techniques were applied to reduce the large initial position uncertainties and to estimate the spacecraft position and velocity for steady state operations.15In terms of integrated navigation based on pulsars’ observations, some observational information such as starlight angular distance, solar Doppler16,17and inertial integrated navigation are adopted as the navigation measurement to aid the XPNAV method.18,19

Distinguishing from the first framework, the HXMT team reported a non-real-time method being referred to as ‘‘Significance Enhancement of Pulse-profile with Orbit-dynamics(SEPO)” to estimate the orbital elements using the distortion of the pulse profile. The 31-day observations of the Crab pulsar from the POLAR detector (about 200 cm2) on-board China’s TianGong-2 space station were used to test pulsar navigation. The orbital elements were determined by the orbit deviation within 20 km (1σ). Analogous experimental demonstration was also performed using the observations of the Crab pulsar from Insight-HXMT satellite. By combining a five-day period of observations from High Energy X-ray Telescope,Medium Energy X-ray Telescope and Low Energy X-ray Telescope,the position and velocity were pinpointed within 10 km and 10 m/s,respectively.16,17The SEPO is an off-line technique to determine the orbital elements. Instead of generating the standard pulse profile,the advantage of SEPO is to use the distortion information of the observed pulse profile.

Despite the above achievements,certain challenges must be overcome for SEPO. First, hardly can a long-term orbitdynamics model be built accurately. As the model error accumulates, the inversion method for deducing orbital elements may fail due to the long-term error averaging effect. Second,the rationale of the SEPO technique is to estimate orbital elements by examining the variation of pulse profile, during which one of the elements is to be changed while the rest five elements are assumed to remain accurate and unchanged.Intriguingly, these six orbital elements are closely combined to determine the position of the spacecraft, for which reason the assumption for fixing the parameters seem unreasonable in the SEPO technique.In fact,the accuracy of the SEPO technique depends heavily on the X-ray flux rather than the time duration of observations. This paper proposes a decoupling method for six-dimension parameter estimation via statistical analysis of grouping bi-chi-squared technique.The key contributions of this paper are summarized as follows.

First,based on the principle of SEPO,we propose an inversion method of grouping bi-chi-squared and a processing framework, which can reduce the orbital elements from the distortion information of the X-ray pulse profiles.

Second,distinguishing from the existing experiments whose results are deduced from long observations (ten days or even longer),we use the RXTE observations to realize higher accuracy within a shorter time.

Third, to minimize the calculation complexity,we establish a short-term orbit dynamics model of the RXTE satellite and rearrange a barycentric TOA transformation.

The remainder of this paper is organized as follows. Section 2 formulates the framework of conducting the inversion method of deducing orbital elements. Section 3 presents the screening criterion for the RXTE observations and determines the spin parameters of the Crab pulsar. Section 4 analyzes the experimental results. Section 5 summarizes conclusions of this paper.

2. Technical scheme

In this section, we propose a technical scheme of orbital elements inversion as shown in Fig.1.Large assumed orbit errors lead to distorted epoch-folded pulse profiles, suggesting that inaccurate orbit parameters result in the deviation of X-ray TOAs and in the blurry pulse profile. On the contrary, the pulse profile is closely adjacent to the real one when the given orbital elements are accurate enough. Similar to the SEPO method, the significance analysis of the pulse profile and orbit dynamics are combined to estimate the orbital elements by using the Crab pulsar observations.20,21

Specifically, the X-ray observations are first screened according to the given conditions, presenting errors of the six orbital elements (semi-major axis, eccentricity, inclination,Right Ascension of Ascending Node (RAAN), argument of Perigee and true anomaly) to calculate the position and velocity of the satellite at a fixed initial epoch.

Fig. 1 Orbital elements inversion using X-ray observations.

Then the position and velocity of the satellite are predicted using the orbit dynamic model that is established by considering the perturbation terms. According to the assumed orbit parameters, the photon TOAs at the spacecraft are converted to the corresponding Barycentric Dynamical Time (TDB) at the Solar System Barycenter (SSB), after which the X-ray phases are computed to derive a pulse profile with the help of the refined timing model, and the significance of the pulse profile is then to be calculated using the chi-squared test with the assumed parameters. Eventually, the parameterdecoupling technique and least squares fitting are performed to obtain the orbital elements.

2.1. Short-term orbit dynamics of RXTE satellite

Owing to the influence of various perturbations, the orbit of RXTE satellite changes continuously, under which circumstance the trajectory will not locate in a closed plane.22To avoid excessive computation during orbit determination, we establish a short-term prediction model for RXTE satellite,which is composed of the following perturbation accelerations

where aEis the non-spherical perturbation of the Earth,aSand aMare the gravitational perturbation from the Sun and the Moon respectively, aAis the atmospheric drag perturbation,aDis the solar radiation pressure perturbation.

(1) Earth Gravity Harmonics

The equatorial bulge exerts a pulling force on the Low-Earth Orbit (LEO) spacecraft that moves toward the equatorial plane. The pulling force exerted is three orders of magnitude lower than the central gravity, causing precessional motion of the orbital plane and shifting of the line of node by several degrees per day. The asphericity of the Earth also induces further perturbations, affecting all the six orbital elements of LEO satellites in particular.22,23To ensure the prediction accuracy of satellite orbit while minimizing the amount of calculation, we adopted JGM-3 model for orbit prediction,and set the degree and order as 20 × 20 in the JGM-3 model.

(2) Atmospheric Drag Effects

The atmospheric drag force mainly affects the semi-major axis and the eccentricity, causing slow attenuation of the RXTE’s orbit. This paper utilizes the Harris-Priester (H-P)atmospheric density model to estimate a rough atmospheric density ρ, which is a relatively simple model based on properties of the Earth’s upper atmosphere under quasi-hydrostatic conditions.

(3) Solar and Lunar perturbations

Since the Sun and the Moon are far away from the Earth,they are generally regarded as particles. The position coordinates of the Sun, the Moon and the Earth relative to SSB are obtained from the Planetary Ephemeris DE405 which refers to the International Celestial Reference System (ICRF)realized by a series of radio sources.

(4) Solar Radiation Pressure

As a perturbation of surface force,the solar radiation pressure typically affects all orbital elements of low-altitude satellite(especially eccentricity and argument of perigee).22,23Here,we adopt the following expression to describe the perturbation acceleration of solar pressure acting upon satellite

where P1AU≈4.56×10-6N/m2is the solar radiation pressure at low-altitude satellite, nSis the unit vector pointing to the direction of solar radiation,μ is the exposure coefficient(which ranges from complete absorption 0 to complete reflection 1),S is the cross-section area of spacecraft perpendicular to the direction of solar radiation. Similar to the atmospheric drag coefficient, the radiation pressure coefficient CRis estimated as a free parameter in orbit determination program.

2.2. Barycentric TOA transformation

In this subsection, we use a version of ’Barycorr’ adopted by the observations of Neutron star Interior Composition Explorer (NICER) project.24,25The adopted time transformation is expressed as

where ΔCrefers to the clock error correction caused by the atomic clock. The intrinsic delay of the RXTE’s instruments needs to be included in high-precision clock corrections. For Proportional Counter Array (PCA), the time delay is about 16 to 20 μs. However, the above step is unnecessary for NICER observations,because the satellite clock of the NICER has been synchronized with GPS clock. The Einstein delay ΔEis a part of the relativistic effect,which can be attributed to the gravitational redshift and time expansion caused by the motion of the Earth and other celestial bodies. The observation clock is affected by the time-varying gravitational potential and Doppler frequency shift. After correcting the term of Einstein delay, the Terrestrial Time (TT) can be converted to the TDB. The Roemer delay ΔRrepresents the vacuum propagation time from RXTE satellite to SBB. The Shapiro delay ΔSmeans the time delay where the light bends due to the gravitational force of celestial bodies (mainly the Sun, Saturn and Jupiter).

The Roemer delay ΔRcan be regarded as the time of photons propagating from spacecraft to SSB, which can be calculated by

where c is the speed of light, the vector rscis the spacecraft detection position relative to SSB,npdenotes the unit direction vector pointing to pulsar,which depends primarily on the position of pulsars in space. For the Earth-orbiting satellite, the vector rsccan be further decomposed into two parts: the inertial position of the Earth (rE) and the relative position (rsc/E),and rsc=rE+rsc/E. Generally, rEis specified by planetary ephemeris,whereas rsc/Eis determined by the orbital dynamics model.The two vectors should be expressed in inertial coordinate system.

The Einstein delay ΔEis a modified transformation from TT to TDB, which can be further divided into two terms where vEis the Earth’s velocity vector, and ΔEGis the analytical expression for the time transformation TDB-TT which containing 127 sinusoidal periodic terms with an accuracy at the 100 ns level.26

The first-order solar gravitational delay is expressed as

where the universal gravitational constant G and the mass of the Sun MSare given by the Planetary Ephemeris DE405. The variable ? is the angle that is jointly formed by the Crab pulsar, the Sun and the RXTE satellite. The above time transformation equation is derived under the assumption that the pulsar is stationary and that the X-ray photons emitted by the pulsar are parallel during the whole observation period. Comparing with the complete time transformation equation, the accuracy of the simplified version can reach a level of microsecond, which is adequate and acceptable for orbital elements determination.

2.3. Orbital elements extracting

To identify RXTE satellite’s position, we divide parametric grids according to the range of the given six orbital elements(i.e. semi-major axis a, eccentricity e, inclination i, RAAN Ω,argument of perigee ω,and true anomaly θ).The grid numbers of the corresponding six orbital elements are Na，Ne，Ni，NΩ，Nωand Nθ. The position and velocity vectors in the orbital coordinate system can be calculated by using the orbital elements,which can be expressed as follows23

Based on the calculated initial position, the satellite’s position at any time can be predicted in inertial coordinate system according to the aforementioned orbital dynamics model in Subsection 2.1.Then,the detected photon TOAs are converted to the SSB via using the time transformation equation presented in Subsection 2.2. Subsequently, the photon TOAs are converted into the phase of photons φ（t ） using the timing model, which can be expressed as

where f （k-1） represents the （k-1）th derivative of the pulsar spin frequency,the time difference t-t0represents the change between the requested epoch t and the initial reference epoch t0, φ（t ）-φ0represents the total accumulated pulsar rotation phase from the initial reference phase φ0. This phaseconversion process is equivalent to the normalization of the pulse period, and the average pulse profile is obtained by

where φiis the i th phase bin(one phase is evenly divided into Nbbins), Ciis the number of accumulated photons in the i th phase bin, Npis the number of the phase cycles for a single observation period, cj（φi） is the number of photons in the i th phase bin and the j th phase cycle. We adopt chi-squared test to calculate the statistical characteristics of the observed profile as follows

The largest chi-squared statistics corresponds to a optimal set of six orbital elements.However,these parameters are coupled together,to which we use a grouping bi-chi-squared technique to extract the orbital elements from the dataset of the chi-squared statistics. The semi-major axis is first obtained by the following equation

where the set of the remaining chi-squared statistics is divided into Negroups and N=NiNΩNωNθ. The notation χekrepresents result of bi-chi-squared statistics for the eccentricity.The superscripts ’res’ marks the remaining chi-squared statistics corresponding to the group where the optimal semi major axis locates. By analogy, the first five orbital elements can be estimated one by one, and the chi-squared test quantity involved in the calculation gradually becomes less. The last orbital element of the true anomaly is estimated by

Intriguingly, the grouping bi-chi-squared technique is not necessarily required in estimating the last orbital element θ.Instead,it needs merely to utilize the residual chi-squared data,thereby estimating the orbital elements one by one.

3. Experimental results

3.1. RXTE observations reduction

The Crab pulsar data were obtained using the PCA mounted on RXTE satellite and is sensitive to the X-ray energy band 2-60 keV with a total collecting area of 6500 cm2and 1° field of view. We investigated the RXTE observations of the Crab pulsar from the data archive of the High Energy Astrophysics Science Archive Research Center (HEASARC) for about 16 years from February 1996 to January 2012.To get available observations for navigation test,we set two important rules as follows.

(1) The interval separation of every two adjacent observations does not exceed two days.

(2) No timing glitches exist in the observation during the whole period of observations.

Most of the Crab pulsar observations obtained by the PCA were separated by one or two weeks, for which reason it was rendered difficult to acquire sufficient data to conduct navigation test. Taking into account the above rules, we selected six eligible data subsets with the event mode between MJD 55463 and MJD 55469, involving about 2.5×107events with a total exposure time of about 5000 s, as listed out in Table 1.

The radio timing model for the Crab pulsar is usually updated once in a month. Unfortunately, the X-ray observation times were just between the two radio timing models,thereby introducing large errors into timing residuals. To obtain a more accurate timing model, we refined a pulsar timing model using the Crab pulsar observations of the RXTE.

Table 1 RXTE observations of the Crab pulsar.

The background flux is counted due to its significant effect on the statistic analysis of the observed pulse profile. According to the pulsar parameters given in Table 2,we obtained the pulse profile by epoch-folding all the RXTE observations.The ratio of the background flux and the pulsar flux is calculated in the range 2-60 keV, and the result is 8.8 which means that background flux accounts for ～90% of the detected X-ray photons. According to this ratio, we can get the pulse flux and background flux with the events data listed in Table 1.The corresponding timing analysis of the Crab pulsar was accomplished by performing not only common phasecoherent analysis of the data, but also partial phase-coherent analysis that is less sensitive to the timing noise superposed on the deterministic spin-down.27By using the known position of the Crab pulsar(RA=05h34m31s.972,Dec=22°00＇52.07＇＇in J2000.0 coordinates and the solar system ephemeris DE405),the time series were reduced to TDB at SSB.The unit direction vector of the Crab pulsar can be calculated by the given RA and DEC.28The corresponding spin parameters were fitted,as listed out in Table 2. Through adopting the fitting techniques,27-29we obtained a refined timing model whose Root Mean Square (RMS) of timing residuals is 1.5 μs.

3.2. Results

We used the short-term dynamic model to calculate the position and velocity vectors of the RXTE satellite, to which the prediction results are shown in Fig.2.It can be seen that small orbit error exists compared with the orbit file from RXTE archive, and the orbit prediction errors gradually diverge as the observation time increases. The position errors are generally less than 400 m, and the velocity errors are less than 0.45 m/s within 38 hours. The increasing deviation is mainly caused by the geopotential model. It is known that highaccuracy orbit prediction is of decisive significance for precise inversion of the long-term RXTE satellite’s orbital elements.However,with respect to the short-term orbital elements inversion of the RXTE satellite, high-accuracy orbit prediction is not necessary as long as the observations can noticeably reflect the changes of orbit parameters.

Fig. 2 Orbital prediction error of the RXTE satellite within about two days.

Based on the six sets of X-ray observations listed in Table 1,the orbit determination experiments were conducted with every two adjacent groups. We statistically averaged the RMS of the six orbital elements, to which the results are presented in Table 3.Although the samples of X-ray observations are limited, it can be apparently observed that the proposed method is effective to obtain an excellent parameter within a period of less than two days. In the proposed grouping bichi-squared technique, we assumed that all the six orbital elements are changeable, and we only restrict them to change within a reasonable range as listed in Table 4. At each grid point of the orbital elements, the observed pulse profile was obtained and calculated by chi-squared test. The calculation results of all the parameters are shown in Fig.3.It can be seen that these orbital elements are coupled together into chisquared dataset, implying that each chi-squared value corresponds to the evaluation of the combination of all the elements. In reality, the periodic fluctuations in chi-squared dataset were caused by the repetitive grids of orbital elements.The local optimal solutions of orbital elements can be obtained by the peaks in the dataset. The middle part of the dataset is relatively flat, indicating that the values of some orbital elements are accurate enough, to which chi-squared test is not applicable to observably reflect the changes of parameters.

Table 2 Spin parameters of Crab pulsar during the selected observation period.

Table 3 Estimated orbit errors of RXTE satellite compared with POLAR experiment (3σ).

Table 4 Errors of initial orbital elements.

Fig. 3 Chi-squared test by traversing all the orbital elements.

To obtain the optimal semi-major axis, we firstly divided the chi-squared dataset into Na(Na=9) groups according to the parameter range. Then we obtained nine groups of distribution curves by histogram statistics of the chi-squared values in each group, as shown in Fig. 4(a). Each distribution curve illustrates the variation of the pulse profile with different errors of the semi-major axis, and the changes of the distribution curve also reflect the influence of the other five orbital elements. It can be observed that smaller semi-major axis error implies more concentrated distribution curve, otherwise the distribution curve is more scattered, indicating that a better statistical significance of pulse profile is obtained under a small semi-major axis error. It is also appropriate to select a curve with the most concentrated distribution and take its corresponding semi-major axis error as the optimal semi-major axis error.The grouping bi-chi-squared technique was used to evaluate the neutralization and dispersion of the distribution curves. The least square fitting for parabola function is used to avoid the calculation errors of the optimal semi-major axis caused by limited sampling values,as shown in Fig.4(b).It can be seen that the semi-major axis parameter can further be refined to obtain better results.However,the amount of calculation will increase significantly with the increase of grid point data.In response to the solution of RAAN,the process is similar, and the results are presented in Fig. 5.

The chi-squared data in Fig. 3 are actually the result of pulse profile significance analysis after iterating over all the orbital elements. In the previous processing, each parameter calculation will strip off part of the chi-square data, which can be regarded as the process of decoupling and there is a sharp decline in the number of data involved in calculations.When the parameter of the true anomaly is to be estimated,only nine remaining chi-square values are left as shown in Fig. 6. We can directly estimate the optimal true anomaly by using parabolic curve fitting based on the least squares technique.Note that there are several orders of magnitude for the bichi-squared statistics among the Figs.4-6.It is mainly because that the bi-chi-square values given by Eq. (14) depend on the number of the remaining chi-squared statistics involving in parameters estimation. When estimating the parameter semimajor axis, a total number of N=NeNiNΩNωNθchi-square values are computed and the corresponding residuals are summed. However, only Nθchi-square values are utilized to estimate the last orbital element.

In order to be able to compare these two methods to some extent, we listed the experiment results of POLAR (by SEPO method) and RXTE (by our method) in Table 3. In these two experiments, the POLAR experiment utilized 4.8 billion X-ray events with total exposure time 782 ks, and RXTE experiment utilized 25 million X-ray events with exposure time of 1.5 ks.The POLAR experimental results were carried out by SEPO method, and the position error and velocity error are 19.2 km and 432 m/s. The RXTE experimental results were carried out by our method, and the position error of 16.3 km (3σ) and a velocity error of 13.3 m/s (3σ) with only two-orbit observations for each experiment.It is seen from the experimental results, our method has a higher accuracy,although the number of observed photons between these two experiments are widening by about 200 times.

Fig. 4 Grouping bi-chi-squared technique results for the grouped semi-major axis.

Fig. 5 Grouping bi-chi-squared technique results for the grouping RAAN.

Fig. 6 Decoupled chi-square values of the pulse profiles for different true anomaly errors.

Furthermore, we calculated the chi-squared values of the pulse profile with respect to each orbital element by SEPO method. When analyzing an orbital element, we assume that the other five contain random errors and the errors’ amplitudes are the half interval of each parameter multiplied by the random number between (-1,1). This assumption is reasonable for the parameter ranges are the only priori information about the number of orbit elements. The SEPO method calculates the chi-squared value of the pulse profile with respect to the errors of each orbital element within the parameter ranges listed in Table 4. The corresponding results obtained by directly using SEPO method are plotted in Fig.7.It can be seen that the chi-square values of pulse profile with respect to each orbital element significantly deviate from the centre point (theoretically, when the error is 0, the corresponding chi-square value is the largest). It is because that the inaccuracy of the other five orbital elements dominates the distortion of the pulse profile and the influence of a single orbital element on the pulse profile is insignificant, greatly reducing the estimation accuracy of orbital elements.

Fig. 7 Chi-squared values of the pulse profile with respect to each orbital element calculated by SEPO method (the other five orbital elements have random errors).

Fig.8 Chi-squared values of the pulse profile with respect to each orbital element calculated by proposed method (the centre values of orbital elements are given by the proposed method).

We also calculated the chi-squared values of the pulse profile with respect to each orbital element by our method.The six plots of Fig. 8 show the relation between the pulse profile distortion and the errors of orbital elements whose centre values are estimated by the proposed method.It can be seen that with the increase of the orbital element errors,the distortions of the pulse profiles are large, resulting in small chi-square values.With the decrease of orbital element errors, the pulse profiles are close to the standard template, and a larger chi-square value is obtained,so as to form a‘‘bell”curve.The results indicate that there is a close relation between the orbital elements and pulse profile distortion in the selected parameter ranges.By calculating the location of the the maximum chi-square value in the ranges, we can obtain the orbital elements.

Both the proposed method and the SEPO method involve complex operations such as barycentric TOA transformation,orbit dynamics and parameter iteration which makes it difficult to directly evaluate the number of floating-point operations. At present, both methods require to estimate the orbital parameters with the help of computer clusters. For more intuitive comparison, the computational complexity of the two methods is compared by the time consumed by performing calculations of one million of photons with a single CPU.The time consumption of SEPO method is 1.25 h/Mcnts,and that of our method is 64 h/Mcnts. Although the barycentric TOA transformation and orbit dynamics have been simplified in this paper without affecting the calculation accuracy, the computational complexity of our method is significantly higher than that of SEPO method due to multidimensional cyclic iteration.It is a key issue to reduce the iteration and the computational complexity of our method in the future.

In fact,due to the short time duration in each observation,the accuracy of orbit estimation is still low.Our future research interest will be focusing on achieving more accurate orbital parameter estimation of 1 km by using large-area detectors with more intensive observation and shorter observation period. What needs to be reminded is that more attention must be paid to the following situations when implementing our proposed method.First,longer duration of the whole observation interval implies stronger average effect of orbit parameters inversion,suggesting a degrading in the accuracy of parameter estimation. This also explains why the experimental results of this study outperform those of the POLAR whose X-ray observation lasts for a time duration of one month. Second,the faster the observation changes with the orbit, the higher the inversion accuracy increases with the orbital elements.Imaging an extreme situation in which the orbit plane is parallel to the radiation direction of the pulsar, our proposed method may exhibit its best performance. However, when the orbit plane is perpendicular to the radiation direction of the pulsar, the orbit motion cannot be reflected in the X-ray observations and the orbit parameters are not observable any more.

4. Conclusions

(1) Based on the idea of SEPO technique, we propose a grouping bi-chi-squared method and a processing framework, which can reduce the orbital elements from the distortion information of the pulse profiles. This technique is applicable to the estimation of orbital elements via using high-flux observation on the ground,nevertheless, it has certain difficulties of being applied in orbit due to the low efficiency of grid search.

(2) The pulsar timing was refined to determine spin parameters of the Crab pulsar (PSR B0531+21) during these observational periods. The RMS of timing residuals for the refined timing model performed better than 1.5 μs.

(3) We establish a short-term orbit dynamics model of the RXTE satellite by considering the effects of diverse perturbations and rearrange a barycentric TOA transformation to minimize the calculation complexity.

(4) Distinguishing from the existing experiments whose results are deduced from a long observations (ten days or even longer), we take use of the RXTE observations to realize higher accuracy within 2-day observations.Experimental results suggest that the position and velocity errors are 16.3 km (3σ) and 13.3 m/s (3σ) with two adjacent observations split by one day (exposure time of 1500 s), outperforming those of the POLAR experiment whose results are 19.2 km (3σ) and 432 m/s (3σ).

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This research adopted PCA observations from RXTE missions published on the website. Funding support for this work was provided by the National Natural Science Foundation of China(No.62103313);the Space Optoelectronic Measurement& Perception Laboratory of BICE, China (No. LabSOMP-2020-06); and the Central Universities, China (No. JC2007).