Since A. Hewish and his graduate student J. Bell discovered the radio pulse signal with a strict cycle in 1967^[1], the research on pulsars has never been stopped. The pulsar signals have extraordinarily high long-term stability and can be used as a natural clock signal. The idea of utilizing pulsars for navigation was proposed in the 1970s^[2]. High-precision measurement of pulsar angle information has the potential to establish a stellar navigation network with meter-level orbiting accuracy^[3], which may fundamentally improve the capability of mankind to pilot into the solar system and beyond. Very long baseline interferometry is the primary method used to determine the angular position of pulsars^[4]. Due to the constraints of baseline length and signal wavelength, the measurement accuracy of pulsar angular positions remains at the millisecond level. Visible light can also be implemented for pulsar angle measurement based on interferometry approaches, yet the resolution is restricted by the observation wavelength and the telescope aperture^[5-8]. The wavelength of X-ray radiation is several orders of magnitude shorter than that of visible light and radio waves, which provides the possibility to achieve higher resolution via X-ray measurement. Researchers have been investigating X-ray pulsar navigation since the 1980s^[9-13]. NASA conducted the first conceptual space demonstration of X-ray navigation in 2018^[14]. However, conventional interferometry methods based on first-order interference extract the phase difference from the directly recorded interference fringes, so the fabricating accuracy of the optical devices needs to be smaller than the wavelength^[15], which is a huge challenge in the X-ray region. Due to the very short wavelength of X-ray radiation, it is difficult to manufacture high-precision X-ray equipment for space interferometry^[16].

Different from first-order interferometry, the methods based on second-order interference detect the intensity fluctuations of optical fields and can indirectly obtain the interference fringes by intensity correlation calculation^[15]. Thus, the requirement for the fabricating accuracy of optical devices is not so high^[17], and they are less affected by environmental noises^[18-20]. The concept of intensity correlation originated from the Hanbury Brown-Twiss (HBT) experiment^[21,22], which was proposed to measure the angular diameter of Sirius^[21]. It relies on coincidence counting measurements, imposing high temporal resolution requirements on the detector. Subsequently, Glauber clarified the nature of higher-order coherence of optical fields^[23], laying the theoretical foundation for intensity correlation measurement. In the past decades, the application of visible light intensity correlation has developed rapidly in remote sensing, super-resolution imaging, spectral imaging, three-dimensional imaging, etc.^[24-31]. Now it has been expanded to the X-ray region^[32-38], and some progress in pulsar observation has been achieved^[39,40]. Nevertheless, these X-ray techniques require a detecting time that is shorter than the coherence time of the radiation field. The coherence time of pulsar X-ray radiation is much less than picoseconds, while existing X-ray space detectors, such as silicon drift detectors (SDDs), have a temporal resolution of only microseconds. This means that X-ray intensity correlation measurement cannot be easily achieved.

In this paper, we propose a method to extract the pulsar angle information based on spatially modulated X-ray intensity correlation (SMXIC). By inserting a modulator in front of the detector, the radiation field is spatially modulated. Then the observing angle of the target pulsar can be obtained by measuring the intensity correlation between two points in the detected radiation field. This approach utilizes the spatial average instead of the temporal ensemble average adopted in previous X-ray correlation techniques. It can significantly reduce the temporal resolution requirements for X-ray detectors, allowing the utilization of detectors with microsecond temporal resolution. The theoretical model of this method is established, and a proof-of-concept X-ray experiment is demonstrated successfully. Related simulation is performed to discuss the potential of this method, which may be helpful in the future application of pulsar navigation.

Figure 1 illustrates our method of measuring the observing angle of a target pulsar based on SMXIC. The X-ray radiation emitted from the pulsar is detected in two directions, namely the reference beam and the test beam. The reference beam is recorded by a single-photon detector $D_r$, and a modulator $M_r$ is positioned in front of the detector. The distance between $D_r$ and $M_r$ is $z_0$. The arrangement in the test beam is the same as that in the reference beam. The difference lies in that the detector $D_t$ can be moved along the $X_t$-direction. ${\beta }_t$ and ${\beta }_r$ represent the incident angles of the X-rays illuminating the modulator $M_r$ and $M_t$, respectively. It should be noted that the two modulators are identical and positioned in the same plane. The angle $\gamma$ between the two beams is the observing angle of the target pulsar to be measured. According to the geometric relationship, the angle $\gamma$ can be determined by $\gamma ={\beta }_t-{\beta }_r$.

Due to the extremely long distance from the pulsar and the limited size of the modulators, the optical field on each modulator can be regarded as spatially coherent and described by a plane wave as $$E_0(x,{\beta }_k)\propto \mathrm{rect}\left(\frac{x}{a}\right)\exp \left(j\frac{2\pi }{\lambda }\sin {\beta }_{k}x\right),k=r,t.$$(1)Here, $x$ is the coordinate on the modulator plane, $a$ is the size of the modulators, and $\lambda$ is the X-ray wavelength.

The intensity distribution detected in the two beams can be written as $$I_k(x_k,{\beta }_k)={\left|\int \mathrm{d}x{E}_0(x,{\beta }_k)t_m(x)h_k(x,x_k)\right|}^2,k=r,t.$$(2)Here, $x_r$ and $x_t$ are the position coordinates of $D_r$ and $D_t$, respectively. $t_m(x)$ denotes the transmittance of the two identical modulators, and $h_k(x,x_k)$ corresponds to the free-space transfer function from the modulator to the detector plane. Under the paraxial approximation, we have $$h_k(x,x_k)\propto \exp \left[\frac{j\pi }{\lambda z_0}{(x-x_k)}^2\right],k=r,t.$$(3)

Then, the correlation between the intensity fluctuations at $D_r$ and $D_t$ can be calculated as $$\begin{gathered} \mathrm{\Delta }{G}_s^{(2)}(x_r,x_t)={⟨I_r(x_r)I_t(x_t)⟩}_s-{⟨I_r(x_r)⟩}_s{⟨I_t(x_t)⟩}_s \\ =\int \mathrm{d}{x}_1\mathrm{d}{x}_1^{\prime }\mathrm{d}{x}_2\mathrm{d}{x}_2^{\prime }{E}_0(x_1,{\beta }_r)E_0^{*}(x_1^{\prime },{\beta }_r)E_0(x_2,{\beta }_t)E_0^{*}(x_2^{\prime },{\beta }_t) \\ \times {⟨t_m(x_1)t_m^{*}(x_1^{\prime })t_m(x_2)t_m^{*}(x_2^{\prime })⟩}_s \\ \times h_r(x_1,x_r)h_r^{*}(x_1^{\prime },x_r)h_t(x_2,x_t)h_t^{*}(x_2^{\prime },x_t) \\ -\int \mathrm{d}{x}_1\mathrm{d}{x}_1^{\prime }{E}_0(x_1,{\beta }_r)E_0^{*}(x_1^{\prime },{\beta }_r) \\ \times {⟨t_m(x_1)t_m^{*}(x_1^{\prime })⟩}_s{h}_r(x_1,x_r)h_r^{*}(x_1^{\prime },x_r) \\ \times \int \mathrm{d}{x}_2\mathrm{d}{x}_2^{\prime }{E}_0(x_2,{\beta }_t)E_0^{*}(x_2^{\prime },{\beta }_t) \\ \times {⟨t_m(x_2)t_m^{*}(x_2^{\prime })⟩}_s{h}_t(x_2,x_t)h_t^{*}(x_2^{\prime },x_t), \end{gathered}$$(4)where ${⟨·⟩}_s$ means the spatial average^[38] and $*$ represents the complex conjugate.

Assuming that the modulators will introduce field fluctuations that can be described by a complex circular Gaussian ergodic random process with zero mean^[41], then we have $$\begin{gathered} {⟨t_m(x_1)t_m^{*}(x_1^{\prime })t_m(x_2)t_m^{*}(x_2^{\prime })⟩}_s \\ =G_s^{(1)}(x_1,x_1^{\prime })G_s^{(1)}(x_2,x_2^{\prime })+G_s^{(1)}(x_1,x_2^{\prime })G_s^{(1)}(x_2,x_1^{\prime }), \end{gathered}$$(5)where $G_s^{(1)}$ is the first-order spatial intensity correlation function. Supposing the feature size of the modulator is sufficiently small, $G_s^{(1)}$ can be regarded as a delta function.

Substituting Eqs. (1)–(3), and (5) into Eq. (4), after some calculations, we obtain $$\begin{gathered} \mathrm{\Delta }{G}_s^{(2)}(x_r,x_t) \\ \propto \int \mathrm{d}{x}_1\mathrm{d}{x}_2\mathrm{rect}\left(\frac{x_1}{a}\right)\mathrm{rect}\left(\frac{x_2}{a}\right) \\ \times \exp [\frac{j2\pi }{\lambda }(\sin {\beta }_r-\sin {\beta }_t+\frac{x_t-x_r}{z_0})(x_1-x_2)] \\ \propto {\mathrm{sinc}}^2\left[\frac{a}{\lambda }(\sin {\beta }_t-\sin {\beta }_r-\frac{x_t-x_r}{z_0})\right]. \end{gathered}$$(6)

Normally, ${\beta }_r$ and ${\beta }_t$ are so small that we have $\sin {\beta }_t-\sin {\beta }_r={\beta }_t-{\beta }_r=\gamma$. The initial position of $D_r$ relative to $M_r$ and the initial position of $D_t$ relative to $M_t$ are the same. Fix the position of $D_r$, and move $D_t$ along the $X_t$-direction to form a shift $\mathrm{\Delta }x=x_t-x_r$ relative to the initial position. The normalized second-order intensity correlation between the two optical fields at $D_r$ and $D_t$ can be calculated to yield $$\begin{gathered} g_s^{(2)}(\mathrm{\Delta }x)=\frac{{⟨I_r(x_r)I_t(x_t)⟩}_s}{{⟨I_r(x_r)⟩}_s{⟨I_t(x_t)⟩}_s} \\ =1+{\mathrm{sinc}}^2\left[\frac{a}{\lambda }(\gamma -\frac{\mathrm{\Delta }x}{z_0})\right]. \end{gathered}$$(7)

As the shift $\mathrm{\Delta }x$ changes, $g_s^{(2)}(\mathrm{\Delta }x)$ changes and reaches its maximum under the condition of $\gamma -\mathrm{\Delta }x/z_0=0$. Therefore, the angle $\gamma$ can be obtained by calculating the intensity correlation between $D_r$ and $D_t$. It should be mentioned that $\gamma$ is only related to the peak position of the intensity correlation results according to Eq. (7).

A proof-of-concept experiment was demonstrated on a tabletop X-ray system operating in vacuum. Figure 2(a) shows the schematic of our experimental layout. A magnesium anode X-ray source was adopted to produce X-rays with a photon energy of 1.25 keV (1 nm wavelength). In the experiment, the tube voltage was set to 5.4 kV, and the current was 0.48 mA. A porous golden film serving as the modulator was positioned at a distance of 1.37 m from the source. The thickness of the film was 280 nm, and the size of the randomly distributed holes on the film was about 100 nm. This film was mounted on a motorized translation stage and can be moved between two positions. To simplify the experiment, we used a CCD camera instead of the single-photon detector for scanning detection. The CCD camera with a pixel size of 15 µm was placed 77 cm downstream of the film to capture the X-ray signals. Initially, the film was regarded as the modulator $M_r$. The X-rays passed through the film, and the intensity distribution of the reference beam $I_r(x_r)$ was recorded by the CCD camera. After moving in the direction indicated by the red arrow in Fig. 2(a), the film was treated as the modulator $M_t$. Then the intensity distribution of the test beam $I_t(x_t)$ was acquired by the CCD camera.

In the experiment, the position of the CCD camera was fixed, and we utilized a specific area ($100\times 100$ pixels) on the CCD plane as the signal detection area, which is denoted by the blue rectangle in Fig. 2(c). In practical scenarios, this specific area can also be acquired by scanning with a single-photon detector. For the reference beam, only one frame of intensity distribution $I_r(x_r)$ was recorded. For the test beam, by scanning the detection area pixel by pixel along the direction indicated by the blue arrow, a series of $I_t(x_t)$ at different $\mathrm{\Delta }x$ was acquired. An example of the intensity distribution pairs of $I_r(x_r)$ and $I_t(x_t)$ is displayed in Fig. 2(b). The spatial intensity correlation between each pair of $I_r(x_r)$ and $I_t(x_t)$ was calculated. After normalization, the values of $g_s^{(2)}(\mathrm{\Delta }x)$ at different $\mathrm{\Delta }x$ were obtained. By fitting these values according to Eq. (7), the observing angle can be extracted from the peak position of the fitting curve.

First, the experiment was performed at $\gamma =0$. The experimental intensity correlation data of $g_s^{(2)}(\mathrm{\Delta }x)$ and the corresponding fitting curve are plotted in Fig. 3. The value of the peak position $\mathrm{\Delta }{x}_p$ can be obtained when $g_s^{(2)}(\mathrm{\Delta }x)$ reaches its maximum value. For comparison, the theoretical curve calculated according to Eq. (7) is also plotted in Fig. 3 with the dashed line. It can be observed that the peak positions of the experimental fitting curve and the theoretical curve are consistent. Compared with the theoretical curve, the experimental curve shows slight widening, which may be attributed to some experimental factors, such as limited pixel size. Fortunately, the peak position of the curve remains unchanged, and extra errors will not be introduced in the angle measurement. It should be mentioned that the correlation is weak at positions with large $\mathrm{\Delta }x$, and the experimental results may be affected by noise.

Then, the experiment of measuring different observing angles was carried out. Six different angles ($\gamma =1.5$, 3.0, 4.5, 6.0, 7.5, and 9.0 arcsec) were measured. For each angle, the intensity correlation values are plotted in Fig. 4. By fitting these data according to Eq. (7) and then calculating the position of the maximum value, we have the value of $\mathrm{\Delta }{x}_p$ for each angle, and the experimental result of $\gamma$ can be calculated from $\mathrm{\Delta }{x}_p/z_0$. The solid curves with different colors in Fig. 4 give the fitting results for the six angles. The peak positions of the six curves are $\mathrm{\Delta }{x}_p=5.29$, 10.84, 16.34, 21.94, 27.09, and 31.74 µm, respectively. Hence, the measured values for the six angles are 1.42, 2.90, 4.38, 5.88, 7.26, and 8.50 arcsec, respectively. These values are close to the corresponding theoretical values.

In order to evaluate the measurement accuracy and reduce random errors, we conducted five repeated measurements for each angle. The results are summarized in Table 1. The first column in the table is the theoretical values for the six observing angles. Each row presents the five measuring results of an angle. By calculating the mean and standard deviation of each angle, the results with error bars are provided in Fig. 5. It is clear that the experimental results are in agreement with the theoretical values. The relative error of the measurements at different observing angles is about 5%. It is mainly caused by the relative errors associated with the distance $z_0$ and the position shift $\mathrm{\Delta }x$. If the relative errors remain stable, a decrease in absolute error can be expected when measuring smaller angles.

To explore the potential for future applications, we also made some numerical simulations based on the theoretical model. According to Eq. (7), the observing angle is obtained at the maximum value of $g_s^{(2)}(\mathrm{\Delta }x)$, where $\gamma =\mathrm{\Delta }x/z_0$. It means that the control error of the position shift $\mathrm{\Delta }x$ and the measurement error of the distance $z_0$ will directly influence the measurement accuracy of the observing angle.

The error equation of the angle measurement process can be derived through the total differential of $\gamma =\mathrm{\Delta }x/z_0$^[42]: $$\mathrm{\Delta }\gamma =|\frac{1}{z_0}\delta (\mathrm{\Delta }x)|+|-\frac{\mathrm{\Delta }x}{z_0^2}\delta (z_0)|.$$(8)

Here, $\delta (\mathrm{\Delta }x)$ represents the control error of the position shift $\mathrm{\Delta }x$, and $\delta (z_0)$ represents the measurement error of the distance $z_0$. For simplicity, $\delta (\mathrm{\Delta }x)$ and $\delta (z_0)$ are referred to as position error and distance error, respectively, in the following text. Their influence on the angle measurement can be investigated through simulation based on Eq. (8).

Typically, X-ray pulsar observations are conducted in space utilizing X-ray apparatus installed on satellite platforms^[13]. At present, the distance $z_0$ can be measured via satellite laser ranging techniques with centimeter precision^[43], and the position control accuracy of satellite formation flight is expected to reach the micrometer level soon^[44-48]. In our simulation, taking the Crab Pulsar as an example, it is situated more than 1000 light years away from the Earth^[49]. The simulation used an X-ray wavelength of 1 nm (1.2 keV), which falls within the spectral range of the pulsar. The distance $z_0$ was set to 800 km, and the position shift $\mathrm{\Delta }x$ was 100 m. The boundary of the two errors adopted in the simulation was $\delta (z_0)=1–10\mathrm{cm}$ and $\delta (\mathrm{\Delta }x)=1–10\mathrm{\mu m}$. The results are shown in Fig. 6. It can be observed that reducing the position error $\delta (\mathrm{\Delta }x)$ and the distance error $\delta (z_0)$ can significantly improve the angle measurement accuracy. From Fig. 6, it can be found that, when the distance error and position error decrease to 3 cm and 3 µm, respectively, the angle measurement accuracy at the level of micro-arcsecond can be expected. Precise position control and ranging of satellites during the SMXIC measurement process are crucial for achieving accurate pulsar angle information.

We have proposed a pulsar angle measurement approach based on spatially modulated X-ray intensity correlation (SMXIC), which overcomes the obstacle of high temporal resolution requirements for X-ray detectors in previous X-ray correlation techniques. It is theoretically proved that the observing angle of a target pulsar can be obtained by measuring the intensity correlation between two points in the detected radiation field. A proof-of-concept experiment based on SMXIC measurement was carried out using a tabletop X-ray system. Six different observing angles have been measured, and the experimental results are consistent with the theoretical values. The impact of satellite position control and ranging errors in future applications has been investigated through simulation, and it may be possible to achieve micro-arcsecond precision angle measurement of pulsars just using available X-ray apparatus.

At this short distance of our proof-of-concept experiment, some imaging methods may also give quite good results. However, they are difficult to apply in space due to the inherent constraints. The accuracy of current pulsar measurement methods is theoretically restricted by the receiving aperture of the detecting system. Our approach provides the feasibility of achieving high-precision measurement of pulsar angle information based on multiple satellites. In practical scenarios, the typical distance of pulsars will be much longer, and the satellites may need to separate far enough to achieve the desired measurement accuracy. From an engineering perspective, it may bring additional challenges, such as the fabrication of modulators and the requirement for high-precision device positioning.