Stringent Tests of Gravity with Highly Relativistic Binary Pulsars in the Era of LISA and SKA

Currently, general relativity (GR) is the best tested theory of gravity and has passed all of the tests with flying colors (Will 2014). But there are tentative alternative theories that could possibly explain the "missing mass" problem in galaxies and the accelerating expansion of our universe without introducing dark matter or dark energy (Jain & Khoury 2010; Clifton et al. 2012). In addition, there are theoretical arguments that GR may not be the exactly correct gravity theory in the ultraviolet end, and it is necessary to test GR and alternative theories more precisely. There are many gravity tests in the solar system and these tests provide a very tight restriction on the parameter space of theories beyond GR. Great achievements were made in bounding the parameterized post-Newtonian (PPN) parameters (Will 2018). Because different theories of gravity correspond to different values of PPN parameters in the weak field, bounding PPN parameters can help us test theories of gravity in a systematic way. The solar system experiments probe gravitation in the weak field regime. However, for some gravity theories, such as the scalar−tensor gravity of Damour & Esposito-Farèse (1993), gravity can be reduced to GR in the weak field, while it is very different from GR in the strong field (Damour & Esposito-Farèse 1992, 1996; Shao & Wex 2016). So testing theories of gravity in the strong field can help us understand gravity theory more comprehensively.

Radio pulsars, rotating neutron stars (NSs) with radio emission along their magnetic poles, are ideal astrophysical laboratories for strong field gravity tests, because of the strong self-gravity of NSs. Binary pulsar systems, in which one or more bodies are NSs, provide us an opportunity to test gravity in the strong field. We can measure the parameters of binary pulsar systems to a very high precision, which benefits from an extremely accurate measurement technology, the so-called pulsar timing (Lorimer & Kramer 2005; Wex 2014). Pulsar timing precisely measures the times of arrival (TOA) of radio signals at a telescope on the Earth and fits an appropriate timing model to these TOAs to get a phase-connected solution. In 1974, Hulse & Taylor (1975) discovered the first binary pulsar system, PSR B1913+16, which is a double NS (DNS) system. This system proved the existence of gravitational wave (GW) radiation for the first time (Taylor et al. 1979), using the fact that the observed orbital period decay rate, ${\dot{P}}_{b}^{\mathrm{obs}}$, is consistent with the predicted value from GR, ${\dot{P}}_{b}^{\mathrm{GR}}$. With more and more binary systems discovered, we can use them to test a variety of theories of gravity (see, e.g., Shao & Wex 2012, 2013; Shao et al. 2013; Shao 2014a, 2014b; Manchester 2015; Kramer 2016; Miao et al. 2020; Xu et al. 2020, 2021). For example, ${\dot{P}}_{b}$ of binary pulsar systems can be used to restrict the time variation of the gravitational constant G (Wex 2014; Zhu et al. 2019); it can also be used to bound the mass of the graviton in massive graviton theories (Finn & Sutton 2002; Miao et al. 2019; Shao 2021).

At present, more than 300 binary pulsar systems have been discovered ⁷ (Manchester et al. 2005), including pulsar-NS (PSR-NS) systems, pulsar−white dwarf (PSR-WD) systems, pulsar−main-sequence star (PSR-MS) systems, and so on. Among them, PSR-NS systems are ideal to test gravity theory, because most pulsars in these systems are millisecond pulsars (MSPs) which come from the "reversal mechanism" (namely, the so-called "recycled pulsars"; Tauris et al. 2017). MSPs are generally more precise timers than normal pulsars, because the timing precision of normal pulsars is usually dominated by red noise, spin irregularities, glitches, etc. (Lyne et al. 2010). Therefore MSPs can provide more precise gravity tests.

Up to now, 19 DNS systems have been discovered from radio telescopes ⁸ (Manchester et al. 2005). One of them is a double pulsar system, PSR J0737−3039A/B (Kramer et al. 2006), which provides a wealth of relativistic tests. The DNS system with the shortest orbital period is PSR J1946+2052 with P_b = 0.078 days, and this system has the largest periastron advance, $\dot{\omega }=25.6\pm 0.3\,\deg \,{\mathrm{yr}}^{-1}$ (Stovall et al. 2018). Binary pulsar systems with shorter orbital periods can provide more relativistic gravity tests because of the higher relative velocity v. They also can provide smaller fractional uncertainties for some binary orbital parameters. For example, the expected dependence of fractional uncertainty of ${\dot{P}}_{b}$ is proportional to P_b ³ (Damour & Taylor 1992; Lorimer & Kramer 2005). So in most cases shorter orbital period PSR-NS systems can help us test gravity better. But we have not detected DNS systems with orbital periods of several minutes by radio observation yet, because the fast-changing Doppler shift caused by orbital acceleration of the pulsar can smear the pulsar signal in the Fourier domain. Therefore, this selection effect makes the detection of pulsars in very tight binary systems difficult (Tauris et al. 2017).

In addition to radio observations, GW observations from LIGO/Virgo have also detected two DNS systems, GW170817 in the second observing run (O2; Abbott et al. 2017) and GW190425 in the third observing run (O3; Abbott et al. 2020). ⁹ Because of the high frequency range of LIGO and Virgo, we can only detect GWs from DNS systems near their phase of merger. Therefore, the orbital period of the discovered DNSs has a gap between 1.88 hr (PSR J1946+2052) and ∼1 ms at the phase of merger. If we can detect a DNS system with an orbital period in this gap, it can help us better understand the population and formation mechanism of DNS systems. That said, these systems with orbital periods of several minutes are highly relativistic, so they can provide complementary gravity tests to those done with current pulsar and GW observations.

In 2018, Kyutoku et al. (2019) proposed a multimessenger strategy, combining the Laser Interferometer Space Antenna (LISA) and the Square Kilometre Array (SKA), to detect radio pulses from Galactic DNSs with a period shorter than 10 minutes. LISA is a future space-based GW detector, and it is sensitive at mHz bands (Amaro-Seoane et al. 2017), which makes DNS systems with ${P}_{b}\sim \min$ a target for detection. SKA (Phase 2), under construction, will greatly improve the detection sensitivity of radio observations in the future. ¹⁰ According to the schedule, LISA and SKA will perform their missions simultaneously for a duration of time. Kyutoku et al. (2019) considered that we can utilize LISA to discover DNS systems, and provide their sky location and binary parameters with a high accuracy. Then we can use SKA to probe whether or not a visible pulsar exists in the DNS system. The sky location information can help SKA locate the DNS, and the information of binary parameters can help recover pulse signals modulated by orbital motion. With the help of LISA, SKA can detect PSR-NS systems with ${P}_{b}\sim \min$ more easily with an affordable computational cost. So to some extent, the multimessenger strategy can help us reduce the difficulty in searches of short-orbital-period PSR-NS systems.

As a future high-sensitivity radio instrument, SKA can improve the measurement precision of TOAs of pulses to an unprecedented level, which can help improve the measured accuracy of binary pulsar parameters. The high-precision binary parameters can further help us improve gravity tests. The theoretically projected relation between P_b and the fractional error of ${\dot{P}}_{b}$ in a measurement is $\delta ({\dot{P}}_{b})\propto {P}_{b}^{3}$ (Lorimer & Kramer 2005). So for the short-orbital-period PSR-NS systems, $\delta ({\dot{P}}_{b})$ can be well measured. Then we can use the high-precision ${\dot{P}}_{b}$ to provide tighter bounds on the parameter space of gravity theories beyond GR.

In general, short-orbital-period binary pulsars are short lived. Recent literature have estimated the number of DNS systems that can be detected by LISA. Based on different Galaxy DNS merger rates, Andrews et al. (2020) considered that LISA will detect on average 240 (330) DNS systems in the Galaxy for a 4 yr (8 yr) mission with signal-to-noise ratio (S/N) above 7. Lau et al. (2020) estimated that around 35 DNSs will accumulate S/N > 8 over a 4 yr LISA mission. Although the estimated number for DNSs to be detected by LISA are different in the above two papers, they confirm that a certain number of DNS systems are likely to be detected by LISA. Therefore, we expect that the multimessenger strategy can detect short-orbital-period PSR-NS systems in the future.

In this work, we assume that in the future some short-orbital-period PSR-NS systems are detected by the multimessenger strategy. Based on this assumption, we simulate the fractional error of parameters of these short-orbital-period PSR-NS systems using the standard pulsar-timing software TEMPO2, ¹¹ and use the simulated precision of the parameter ${\dot{P}}_{b}$ to constrain specific theories of gravity, including varying-G theory, massive gravity, and scalar−tensor theories. Results of our simulation show significant improvement in the constraints of specific parameters of gravity theories in the strong field. We also provide a simple estimate of orbital periastron advance rate, ${\dot{\omega }}_{\mathrm{LT}}$, due to the Lense−Thirring effect. Our result suggests that we are able to measure the value of ${\dot{\omega }}_{\mathrm{LT}}$ in the short-orbital-period PSR-NS systems, which can help us restrict the equations of state (EOSs) of NSs. We need to underline that we have ignored higher-order contributions in our simulation; thus our results are indicative, but they provide a reasonable estimate of the magnitude of the precision of parameters.

The paper is organized as follows. In the next section, we overview the multimessenger strategy and list the estimated number of DNSs that could be detected by LISA from several works. In Section 3, we simulate the parameters of PSR-NS systems to be detected by the joint observations of LISA and SKA. We give the simulated fractional uncertainty of parameters. In Section 4, we use our simulated fractional uncertainty of the parameter ${\dot{P}}_{b}$ to constrain the parameter spaces of some specific theories of gravity. In Section 5, we provide an estimate on the order of magnitude of orbital periastron advance rate, ${\dot{\omega }}_{\mathrm{LT}}$. Finally we give a discussion in Section 6.

At present, the detected DNS system of the shortest orbital period is PSR J1946+2052 with P_b = 0.078 days (Stovall et al. 2018). The merging timescale is provided by Peters (1964) and Liu et al. (2014),

$$\begin{eqnarray}&&{T}_{\mathrm{merge}}=\displaystyle \frac{5}{256}{\left(\displaystyle \frac{{P}_{b}}{2\pi }\right)}^{8/3}{T}_{\odot }^{-5/3}\displaystyle \frac{{\left({m}_{p}+{m}_{c}\right)}^{1/3}}{{m}_{p}{m}_{c}}G(e),\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}G(e) & = & 1-3.6481{e}^{2}+5.1237{e}^{4}\\ & & -\ 3.5427{e}^{6}+1.3124{e}^{8}+\cdots .\end{array}\end{eqnarray} \tag{ 2 }$$

In Equation (1), T_⊙ ≡ GM_⊙ c⁻³ ≈ 4.9254909 μs where M_⊙ is the solar mass, m_p and m_c are the masses of the pulsar and the companion in solar units, respectively, and e is the eccentricity of the binary system. By calculating Equation (1), we get the merging timescale of PSR J1946+2052, T_merge = 46 Myr, which means that this system will take 46 Myr to merge. The GW events from DNS mergers suggest that the Galactic merger rate of DNSs is about one in 10⁴ yr (Abbott et al. 2017, 2020). So there must be undetected DNS systems with P_b < 0.078 days in the Galaxy. The reason why we have not observed DNS systems with an orbital period less than one hour is that we are constrained by our present searching methods. Radio observations have contributed the most to DNS detection, as 19 DNS systems have been detected by radio telescopes. Orbital motion of the short-orbital-period binary systems can lead to a severe Doppler smearing that makes it difficult for radio observations to recover the pulse signals of such systems. In GW observations, LIGO/Virgo detected two DNS systems, GW170817 in O2 and GW190425 in O3 (Abbott et al. 2017, 2020). However, LIGO/Virgo are only sensitive to the merger phase. Therefore, constrained by the existing technology and observations, currently it is difficult to detect DNS systems with ${P}_{b}\sim \min$.

LISA, a future space-based GW detector, which is expected to be operational in the 2030s, is sensitive at mHz bands, so the DNS systems with P_b ∼ minutes are potential targets for LISA. Recently, several works have provided estimates of the number of DNS systems that could be detected by LISA. Andrews et al. (2020) used the Galaxy DNS merger rate of 210 Myr⁻¹, which is derived from Abbott et al. (2017, 2019), to estimate the number of DNSs that LISA could detect. They showed that LISA will detect on average 240 (330) DNS systems in the Galaxy with S/N > 7 for a 4 yr (8 yr) mission. They also used a conservative DNS merger rate, ${42}_{-14}^{+30}\,{\mathrm{Myr}}^{-1}$ (Pol et al. 2019), and estimated that LISA will detect on average 46 (65) DNS systems in the Galaxy for a 4 yr (8 yr) mission with S/N > 7. The smaller merger rate could suffer from a selection effect (Tauris et al. 2017). At the same time, Lau et al. (2020) provided an estimate that around 35 DNSs will accumulate S/N above 8 over a 4 yr LISA mission. They gave a smaller detectable number of DNSs than that in Andrews et al. (2020). It is because Lau et al. (2020) used a smaller DNS merger rate, 33 Myr⁻¹. Overall, whether we use a conservative or an optimistic DNS merger rate, LISA has a high probability of detecting DNS systems with ${P}_{b}\sim \min$, and there may be PSR-NS systems among them.

Kyutoku et al. (2019) provided a multimessenger strategy, using LISA and SKA, to detect PSR-NS systems with P_b ∼ 10 minutes in the Galaxy. LISA can detect DNSs with P_b ∼ minutes and provide information on the location and orbital parameters of DNSs. This information can help SKA locate DNS systems and recover the radio pulse signal that could exist in the system (Kyutoku et al. 2019). So, the joint observation strategy can significantly reduce the number of trials used to search for DNS systems with radio telescopes alone. Kyutoku et al. (2019) estimated the 1σ statistical errors of parameters when S/N ρ = 200 for a 2 yr observation mission, and the uncertainties of parameters can be found in their Equations (9)–(11). With the estimated errors of parameters, Kyutoku et al. (2019) showed that LISA can locate a DNS within an error range ΔΩ ≈ 0.036 deg² (see their Equation (11)) and then SKA only needs to take 240 pointings to locate the DNS system precisely by a tied-array beam. So the multimessenger strategy can greatly reduce the number of points that an all-sky survey needs. After locating the DNS system, we can use the orbital frequency and binary parameters provided by LISA to recover a possible pulse signal of the pulsar, which also can greatly reduce the number of required trials when using a coherent integration. ¹² We should notice that an all-sky survey is not confined to detect PSR-NS systems with P_b ∼ 10 minutes. It can be used for other radio source searching as well. The multimessenger observation can specifically and efficiently help us detect PSR-NS systems with P_b ∼ 10 minutes.

Based on Kyutoku et al. (2019), Thrane et al. (2020) offered an extended analysis for the joint detection of LISA and SKA. They considered a 4 yr observation mission of LISA where the S/N can be improved to ρ = 360, so the fractional errors of parameters can be further reduced. With the improved precision of parameters, we can use the joint detection to detect short-orbital-period PSR-NS systems more efficiently.

In conclusion, with a reasonable merger rate for Galactic DNSs, the joint detection of SKA and LISA can provide an efficient and more targeted way for detecting PSR-NS systems with very short orbital periods.

The multimessenger strategy can offer a specific approach to search for the short-orbital-period PSR-NS systems. These systems can provide a more relativistic situation in which to test gravity. In our paper, we assume that we will have detected PSR-NS systems by this way in the future, and we simulate the precision of orbital parameters of the short-orbital-period PSR-NS systems for SKA.

SKA as the next generation radio telescope will enable the precision of TOAs to reach an unprecedented level. Liu et al. (2014) provided an approximate measurement precision on TOAs at 1.4 GHz with 10 minutes integration, which is listed in Table 1 of their work. These results are based on the presumed instrumental sensitivities. In practice, the precision of TOAs depends not only on instrumental sensitivities, but also on interstellar medium evolution, intrinsic pulse shape changes, etc. (Liu et al. 2011). Based on various effects, Liu et al. (2011) predicted that a TOA precision of a normal-brightness MSP can achieve 80–230 ns at 1.4 GHz with 10 minutes integration by the SKA.

The timing precision of young pulsars is normally dominated by red noise, spin irregularities and glitches, etc., so the timing precision could not be significantly improved by increasing the system sensitivity of telescopes. Hence, we assume that the pulsar in a PSR-NS system is recycled, for which TOAs have a stable precision in general, and we use a spin period P = 20 ms. For simplification, we only consider white noise. In our work, considering the possible contributions of various effects, we use two different measurement precisions of TOAs to perform the simulations, σ_TOA = 100 ns and σ_TOA = 1 μs. LISA is sensitive to DNS systems with P_b ∼ 10 minutes, so the target of the multimessenger observation strategy is mainly PSR-NS systems with P_b ∼ 10 minutes. Meanwhile, radio observation can also provide us with an opportunity to detect systems with a period around 100 min by an "acceleration search" (Johnston & Kulkarni 1991) and a "jerk search" (Bagchi et al. 2013; Andersen & Ransom 2018). So we simulate PSR-NS systems with P_b ∈ [8, 120] minutes. For the masses of the NSs, we choose two sets of parameters, m_p = 1.3 M_⊙, m_c = 1.7 M_⊙ and m_p = 1.35 M_⊙, m_c = 1.44 M_⊙. Andrews et al. (2020) used the equations from Peters (1964) to show how orbits of the 20 known PSR-NSs ¹³ will evolve over the next 10 Gyr as gravitational radiation causes them to circularize and decay. Figure 1 in Andrews et al. (2020) shows that the orbital eccentricity goes smaller than 0.1 when the orbital period P_b of the 20 PSR-NS systems evolves to 10 minutes. So in our simulation, we set eccentricity e = 0.1 for the simulation of all systems. ¹⁴

We use the FAKE plug-in in TEMPO2 (Hobbs 2012) to simulate pulsar-timing data. We need to provide a parameter file that contains parameters of the pulsar and binary orbit as input. In this work, we are interested in how the short orbit PSR-NS systems combined with the high-precision measurement of SKA can improve the measurement precision of orbital parameters. Hence, we only consider simulating the time delay of TOAs related to orbital motion. These time delay terms of orbital motion are the Römer delay caused by the orbital motion of the pulsar, the Einstein delay caused by the gravitational redshift and the second-order Doppler effect, and the Shapiro delay, which is due to the gravitational field of the companion star. ¹⁵ For the timing model of orbital motion, different models have different orbital parameters to fit. Damour & Deruelle (1986) introduced a phenomenological parameterization, which is called the parameterized post-Keplerian (PPK) formalism, to describe the orbital motion of pulsars. The post-Keplerian (PK) parameters are solutions of the first-order post-Newtonian (PN) equations of orbital motion and they are generic for fully conservative gravity theories (Damour & Deruelle 1985, 1986). In our simulation, we choose the Damour−Deruelle (DD) timing model.

In the DD timing model, we are mainly able to measure 5 PK parameters—$\dot{\omega },\gamma ,r,s,{\dot{P}}_{b}$—in binary pulsar systems. In GR, they can be expressed as functions of two component masses (Lorimer & Kramer 2005),

$$\begin{eqnarray}&&\dot{\omega }=3{T}_{\odot }^{2/3}{\left(\displaystyle \frac{{P}_{b}}{2\pi }\right)}^{-5/3}\displaystyle \frac{1}{1-{e}^{2}}{\left({m}_{p}+{m}_{c}\right)}^{2/3},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\gamma ={T}_{\odot }^{2/3}{\left(\displaystyle \frac{{P}_{b}}{2\pi }\right)}^{1/3}e\displaystyle \frac{{m}_{c}\left({m}_{p}+2{m}_{c}\right)}{{\left({m}_{p}+{m}_{c}\right)}^{4/3}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&r={T}_{\odot }{m}_{c},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&s=\sin i={T}_{\odot }^{-1/3}{\left(\displaystyle \frac{{P}_{b}}{2\pi }\right)}^{-2/3}x\displaystyle \frac{{\left({m}_{p}+{m}_{c}\right)}^{2/3}}{{m}_{c}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\dot{P}}_{b}=-\displaystyle \frac{192\pi }{5}{T}_{\odot }^{5/3}{\left(\displaystyle \frac{{P}_{b}}{2\pi }\right)}^{-5/3}f(e)\displaystyle \frac{{m}_{p}{m}_{c}}{{\left({m}_{p}+{m}_{c}\right)}^{1/3}},\end{eqnarray} \tag{ 7 }$$

with

$$\begin{eqnarray}&&f(e)\equiv \displaystyle \frac{1+(73/24){e}^{2}+(37/96){e}^{4}}{{\left(1-{e}^{2}\right)}^{7/2}}.\end{eqnarray} \tag{ 8 }$$

In the above equations, $x\equiv (a/c)\sin i$ is the projected semimajor axis of the pulsar orbit and i is the orbital inclination; $\dot{\omega }$ is the relativistic precession of periastron of the pulsar; γ is a PK parameter that is related to the Einstein delay; r and s are PK parameters that are related to the Shapiro delay; and ${\dot{P}}_{b}$ is the orbital period decay rate caused by GW quadrupole radiation in GR. For the five PK parameters, except the parameters r and s that do not depend on P_b , the expected fractional uncertainty of ${\dot{P}}_{b}$ is proportional to P_b ³, the expected fractional uncertainty of $\dot{\omega }$ is proportional to P_b , and the expected fractional uncertainty of γ is proportional to ${P}_{b}^{4/3}$; these dependent relationships are listed in Table 2 in Damour & Taylor (1992). So PSR-NS systems with P_b ∼ 10 minutes can provide better measurement precision for PK parameters ${\dot{P}}_{b}$, $\dot{\omega }$, and γ, which can help us improve the ability of gravity tests. In principle, there are extra relativistic PK parameters, e.g., δ_θ . We have omitted them in this work.

As mentioned previously, we should provide the orbital parameters of PSR-NS systems for the parameter file. We use the parameters of the example PSR-NSs listed in Table 1 and Equations (3)–(7) to provide the values of the PK parameters. In fact, the PK parameters should be generic; however, current measurements of PK parameters in binary pulsars show no significant deviation from the predicted values in GR. In general, we use the error of the observations of PK parameters as a deviation from GR to limit the parameter spaces of alternative theories of gravity. Here, we assume that GR is correct and predict the bounds on other theories of gravity constrained by these PSR-NS systems.

We assume a 4 yr observation and ten TOAs per week for SKA. Ten TOAs per week means that we need to observe a PSR-NS system for 100 minutes per week, when each TOA is obtained with a 10 minutes integration time. ¹⁶ With the parameter file and fitting conditions determined, TEMPO2 can generate simulated TOAs and give a phase-connected solution using the FAKE plug-in by a weighted least-squares algorithm (for details, see Hobbs et al. 2006). The simulation can provide values and measurement uncertainties of parameters, and we can utilize the simulated results to predict how the short-period PSR-NS systems can limit specific theories of gravity.

In our simulation, we need to give a value of the longitude of periastron ω₀ at the epoch of periastron passage T₀ in the parameter file. Figure 3 of Damour & Taylor (1992) shows that the fractional uncertainties of five PK parameters vary with the value of ω₀. The reason is that there is some degeneracy between different time delay terms, e.g., the Römer delay and Einstein delay, and different values of ω₀ can lead to different degrees of degeneracy. Therefore we simulate the fractional uncertainties of five PK parameters for different values of ω₀. We exhibit results for the "PSR-NS I" system with δ_TOA = 100 ns (see Table 1) for different P_b 's in Figure 1. As P_b increases, different values of ω₀ can result in a significant difference in the fractional errors of 4 PK parameters, except ${\dot{P}}_{b}$. So we take the effect of ω₀ into account in the simulation and calculate the fractional error with ${\omega }_{0}\in \left[0^\circ ,360^\circ \right)$.

Zoom In Zoom Out Reset image size

Figure 2 shows the fractional uncertainties of the five PK parameters as functions of P_b for a 4 yr observation. The band of the expected fractional uncertainty is caused by simulating with different ω₀'s. In the upper panels of Figure 2, we present the expected measurement precisions of the five PK parameters with σ_TOA = 1 μs. In the lower panels of Figure 2, we present the results with σ_TOA = 100 ns. As shown in Figure 2, the simulated fractional uncertainties of PK parameters with σ_TOA = 100 ns are an order of magnitude smaller than those with σ_TOA = 1 μs. SKA can provide a high-precision measurement for TOAs of pulsars. The expected sensitivity could be even smaller than 100 ns (Liu et al. 2014), so the precision of PK parameters might be even higher for real detection in the future. At the same TOA precision, the simulated results of the PSR-NS I system are slightly better than those of the PSR-NS II system, which could be caused by a more significant mass difference and thus some reduced degeneracy among parameters.

Zoom In Zoom Out Reset image size

Based on a 4 yr observation, we notice that the simulated fractional uncertainty of ${\dot{P}}_{b}$ as a function of P_b is largely consistent with what theory predicts, namely $\propto {P}_{b}^{3}$ (Damour & Taylor 1992; Lorimer & Kramer 2005). ¹⁷ So the predicted fractional uncertainty of ${\dot{P}}_{b}$ can offer reasonable predictions of gravity tests. The simulated fractional uncertainties of PK parameters r and s are mainly consistent with the theoretical prediction, and are not dependent on P_b , except that there is a slight rise in P_b ∈ [40, 60] minutes. We notice that the slope of PK parameter $\dot{\omega }$ changes significantly in P_b ∈ [30, 60] minutes for a 4 yr observation, which could be caused by the degeneracy between different time delay terms. For parameter γ in a 4 yr observation, our simulated results decrease with P_b when P_b < 60 minutes, and then increase with P_b . So the changing trend of $\delta \left(\gamma \right)$ is not strictly consistent with the theoretical prediction, $\propto {P}_{b}^{4/3}$, when P_b < 60 minutes. For DD model, according to Equations (2.2b) and (2.2c) in Damour & Taylor (1992), the Einstein delay is always degenerate with the Römer delay, and the degeneracy depends on the change of relativistic advance of periastron (Liu et al. 2014). So we think that the trend of $\delta \left(\gamma \right)$ could come from the degeneracy between the Römer delay and the Einstein delay. For PSR-NS systems with P_b ∼ 10 minutes, they have a large $\dot{\omega }$, and considering a 4 yr observation, their periastrons have precessed several rounds. Hence the periodic change of ω cannot eliminate the degeneracy effectively enough and leads to the trend of $\delta \left(\gamma \right)$ that is different from the theoretical prediction. For binary systems with P_b > 60 minutes, the change of ω is less than 2π, so the trend of $\delta \left(\gamma \right)$ varying with P_b is roughly consistent with the theoretical prediction when considering different values of ω₀. ¹⁸ So the DD model cannot avoid the degeneracy because of the specific form of the time delay terms, which makes γ less well measured when Δω > 2π in observations.

We have showed the expected fractional uncertainties of PK parameters of short-orbital-period PSR-NS systems that could be detected by LISA and SKA in joint observations in the future. But we need to emphasize that, strictly speaking, the DD timing model is not sufficient to describe the orbital motion of PSR-NS systems with a very short orbital period ${P}_{b}\sim \min$, because it only takes into account the leading PN order terms. A PSR-NS system with P_b = 8 minutes is an extremely relativistic system, so, in principle, we need to take higher-order effects into account in our simulation. Since we mainly use the fractional error of ${\dot{P}}_{b}$ to test gravity in the following, we provide an extended analysis for ${\dot{P}}_{b}$. Equation(7) gives the 2.5 PN formula for ${\dot{P}}_{b}$, ${\dot{P}}_{b}^{2.5\,\mathrm{PN}}$, and we actually have appended it with the 3.5 PN contribution (Blanchet & Schaefer 1989),

$$\begin{eqnarray}&&\delta {\dot{P}}_{b}^{3.5\,\mathrm{PN}}={\dot{P}}_{b}^{2.5\,\mathrm{PN}}{X}^{3.5\,\mathrm{PN}},\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&\begin{array}{l}{X}^{3.5\,\mathrm{PN}}=\displaystyle \frac{{\left({{GmM}}_{\odot }\right)}^{2/3}}{336\left(1-{e}^{2}\right){c}^{2}f(e)}{\left(\displaystyle \frac{{P}_{b}}{2\pi }\right)}^{-2/3}\\ \quad \times \ \left[1273+\displaystyle \frac{16495}{2}{e}^{2}+\displaystyle \frac{42231}{8}{e}^{4}+\displaystyle \frac{3947}{16}{e}^{6}\right.\\ \quad -\ \left(924+3381{e}^{2}+\displaystyle \frac{1659}{4}{e}^{4}-\displaystyle \frac{259}{4}{e}^{6}\right)\displaystyle \frac{{m}_{p}{m}_{c}}{{m}^{2}}\\ \quad +\ \left.\left(3297{e}^{2}+4221{e}^{4}+\displaystyle \frac{2331}{8}{e}^{6}\right)\displaystyle \frac{{m}_{p}-{m}_{c}}{m}\right].\end{array}\end{eqnarray} \tag{ 10 }$$

For a PSR-NS system with P_b = 8 minutes and e = 0.1, one has X^{3.5 PN} ∼ 1 × 10⁻⁴. Our simulated fractional error of ${\dot{P}}_{b}$ can reach 10⁻⁷ for σ_TOA = 1 μs. So, comparing with the simulated fractional error, the contribution from $\delta {\dot{P}}_{b}^{3.5\,\mathrm{PN}}$ can be significant for an 8 minutes system. However $\delta {\dot{P}}_{b}^{3.5\,\mathrm{PN}}$ is still a small value relative to ${\dot{P}}_{b}^{2.5\,\mathrm{PN}}$, so we use the DD timing model to simulate this situation, as it can give a correct estimate of the magnitude of the expected fractional uncertainties of the PK parameters, and it can predict the ability of a highly relativistic binary pulsar to test gravity. We hope that we can provide a more suitable timing model to simulate or process TOAs for short-orbital-period PSR-NS systems in future work.

In this section, we utilize the simulated results of ${\dot{P}}_{b}$ to do some specific gravity tests. Before we use the fractional error of ${\dot{P}}_{b}$ to test gravity, we need to emphasize that we should subtract the nonintrinsic contribution in ${\dot{P}}_{b}$ and use the intrinsic value ${\dot{P}}_{b}^{\mathrm{intr}}$ in a real gravity test. The value of ${\dot{P}}_{b}$ from actual pulsar timing is ${\dot{P}}_{b}^{\mathrm{obs}}$, which contains the dynamical contributions from the Shklovskii effect (Shklovskii 1970) and the differential Galactic acceleration between the solar system barycenter and the pulsar system (Damour & Taylor 1991; Lazaridis et al. 2009). So, in general, we need to subtract the contribution of the Shklovskii effect and the Galactic acceleration difference in ${\dot{P}}_{b}^{\mathrm{obs}}$ first via

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\dot{P}}_{b}}{{P}_{b}}\right)}^{\mathrm{intr}}={\left(\displaystyle \frac{{\dot{P}}_{b}}{{P}_{b}}\right)}^{\mathrm{obs}}-{\left(\displaystyle \frac{{\dot{P}}_{b}}{{P}_{b}}\right)}^{\mathrm{Shk}}-{\left(\displaystyle \frac{{\dot{P}}_{b}}{{P}_{b}}\right)}^{\mathrm{Gal}}.\end{eqnarray} \tag{ 11 }$$

The second term of Equation (11) is the contribution from the Shklovskii effect, which is due to the transverse motion of the pulsar system (Shklovskii 1970), and the specific form is shown by Equation (9) in Lazaridis et al. (2009). The third term of Equation (11) is the contribution from the Galactic acceleration difference (Damour & Taylor 1991; Lazaridis et al. 2009), and the specific form is given by Equation (16) in Lazaridis et al. (2009).

The actual fractional error of ${\dot{P}}_{b}^{\mathrm{intr}}$ depends on the errors of ${\dot{P}}_{b}^{\mathrm{Shk}}$ and ${\dot{P}}_{b}^{\mathrm{Gal}}$. So, the fractional errors of ${\dot{P}}_{b}^{\mathrm{Shk}}$ and ${\dot{P}}_{b}^{\mathrm{Gal}}$ will affect gravity tests. In other words, their fractional errors can determine how tight the PSR-NS system can bound the parameter space of non-GR parameters. For ${\dot{P}}_{b}^{\mathrm{Shk}}$ and ${\dot{P}}_{b}^{\mathrm{Gal}}$, their fractional errors all rely on the precision of distance D. At present, we use very long baseline interferometry (VLBI) or the dispersion measure (DM) to measure the distance of binary pulsars. For the nearby binary pulsar systems, VLBI can provide an accurate measurement for D. Once we have a measured distance with a small error, the errors of ${\dot{P}}_{b}^{\mathrm{Shk}}$ and ${\dot{P}}_{b}^{\mathrm{Gal}}$ could contribute little to the total error, and the precision of ${\dot{P}}_{b}^{\mathrm{intr}}$ would largely depend on ${\dot{P}}_{b}^{\mathrm{obs}}$, such as in the case of PSRs J2222−0137 (Cognard et al. 2017) and J0737−3039 (Kramer et al. 2006). In the opposite case, for remote binary pulsar systems, we cannot provide an accurate measurement of D by VLBI or DM. The large uncertainty of D results in large errors for ${\dot{P}}_{b}^{\mathrm{Shk}}$ and ${\dot{P}}_{b}^{\mathrm{Gal}}$, which restrict the precision of ${\dot{P}}_{b}^{{\rm{intr}}}$; for example, see PSR B1913+16 (Weisberg & Huang 2016; Deller et al. 2018).

So the accurate measurement of the distance of a binary pulsar system helps us reduce the error contribution from astrophysical terms in ${\dot{P}}_{b}$. Based on a specific DNS system, Kyutoku et al. (2019) provided an estimated fractional uncertainty of distance for a 2 yr mission of LISA,

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}D}{D}=0.01{\left(\displaystyle \frac{\rho }{200}\right)}^{-1}.\end{eqnarray} \tag{ 12 }$$

In Thrane et al. (2020), considering a 4 yr mission of LISA, the fractional uncertainty of distance can reach 0.006 for a specific DNS system with ρ = 360. So LISA could improve the accuracy of pulsar distance measurement significantly and then maybe avoid introducing large errors caused by ${\dot{P}}_{b}^{\mathrm{Shk}}$ and ${\dot{P}}_{b}^{\mathrm{Gal}}$ (with a well modeled Galactic potential). Therefore, if LISA can detect a PSR-NS system with a high S/N, it can improve the measure of the distance of a PSR-NS system.

To simplify the simulation, we did not include the possible error contribution from ${\dot{P}}_{b}^{\mathrm{Gal}}$ and ${\dot{P}}_{b}^{\mathrm{Shk}}$, so we may derive optimistic results when testing gravity. In Section 4.2, to be more complete, we provide a simple calculation for the possible error contribution to ${\dot{P}}_{b}^{\mathrm{intr}}$ from ${\dot{P}}_{b}^{\mathrm{Shk}}$ and ${\dot{P}}_{b}^{\mathrm{Gal}}$, and show the possible impact on gravity tests. We assume LISA can work well in distance measurements, so we use a fractional error of ΔD/D = 0.01. Actually, D may not be measured with such a high precision, so even if the nonintrinsic effects are considered, the results may still be optimistic because D is assumed to be measured to a high precision.

4.1. Theoretical Framework

Based on the expected measurement precision of the Five-hundred-meter Aperture Spherical radio Telescope and SKA, Liu et al. (2014) simulated the timing precision for black hole−pulsar (BH-PSR) binary systems. Seymour & Yagi (2018) used the simulated parameter precision of ${\dot{P}}_{b}$ from Liu et al. (2014) to constrain specific gravity theories. In this section, following Seymour & Yagi (2018), we use the simulated precision of ${\dot{P}}_{b}$ of the short-orbital-period PSR-NS systems to constrain alternative theories of gravity.

In Seymour & Yagi (2018), the authors gave a generic formalism that provides a mapping between generic non-GR parameters in the orbital decay rate ${\dot{P}}_{b}$ and theoretical constants in various alternative theories of gravity. Details can be found in Table 1 of Seymour & Yagi (2018). A parameterized non-GR modification to ${\dot{P}}_{b}$ reads,

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{P}}_{b}}{{P}_{b}}=\displaystyle \frac{{\dot{P}}_{b}}{{P}_{b}}\left|{}_{\mathrm{GR}}\left(1+{\rm{\Gamma }}{v}^{2n}\right),\right.\end{eqnarray} \tag{ 13 }$$

where v is the relative velocity of two compact objects in a binary system, $v={\left(\tfrac{2\pi M}{{P}_{b}}\tfrac{G}{{c}^{3}}\right)}^{1/3}$, and M is the total mass of the binary system. The term $\left({\dot{P}}_{b}/{P}_{b}\right)\left|{}_{\mathrm{GR}}\right.$ is the quantity evaluated in GR (see Equation (7)). The non-GR parameters are Γ and n. Term Γ gives the overall magnitude of the correction, and n shows how the correction depends on v. In Equation (13), a correction term proportional to v²ⁿ corresponds to an n-PN order correction to GR. For different theories of gravity, the non-GR parameter Γ corresponds to different forms, which are listed in Table 1 of Seymour & Yagi (2018), where a specific non-GR theory can be mapped to a set of (Γ, n). In our work, we use the same mapping relation to test gravity. We want to estimate how tightly the short-orbital-period PSR-NS systems can constrain the parameters of non-GR theories in the future.

Following Equation (4) in Seymour & Yagi (2018), we define a fractional error δ for the observed ${\dot{P}}_{b}/{P}_{b}$, and we can use δ to bound Γ by the following equation,

$$\begin{eqnarray}&&\left|{\rm{\Gamma }}\right|\lt \displaystyle \frac{\delta }{{v}^{2n}}.\end{eqnarray} \tag{ 14 }$$

In Section 3, we have simulated a series of fractional errors of ${\dot{P}}_{b}$ for different P_b 's, and these fractional errors of ${\dot{P}}_{b}$ correspond to the parameter δ here. Hence we can use the data to predict the constraints on Γ for different PN orders.

We choose a PSR-NS system with an 8 minute orbital period. Figure 3 presents the predicted upper bounds on Γ as a function of PN orders. For comparison, we also present the result from a 3 day orbital period PSR-BH system (Seymour & Yagi 2018). For a PSR-NS system with an 8 minute orbital period, its relative velocity v approximately equals 5.7 × 10⁻³, which is 5 times larger than the velocity of a 3 day orbital period PSR-BH system, for which v ∼ 1 × 10⁻³. For negative PN, Γ ∝ v^∣2n∣, so for a larger value of n, a binary system with a larger v is more disadvantageous for bounding Γ. For example, if n = −4, for an 8 minute PSR-NS system, its value of v^∣2n∣ is 4 × 10⁵ times larger than that of a 3 day PSR-BH system, which greatly reduces the restriction on Γ. However, the 8 minute orbital period PSR-NS system has a smaller fractional error δ, which makes it possible to provide a stronger constraint on Γ. Combining those two factors, the constraint from an 8 minute PSR-NS system is tighter than the constraint from a 3 day PSR-BH system at some negative PN orders, as shown in Figure 3. At positive PN orders, a larger v is more competitive, so an 8 minute PSR-NS system can naturally provide a tighter constraint than a 3 day PSR-BH system. Like Seymour & Yagi (2018), we also show the limits from GW observations in Figure 3, using GW150914 and GW151226 as examples (Yunes et al. 2016). Because GWs detected by LIGO/Virgo come from the phase of merger, which indicates a larger relative velocity, it leads to a tighter restriction on Γ at positive PN orders. But due to the small δ of the 8 minute PSR-NS system, it can give a slightly tighter restriction than GW observations at positive PN order when n ≲ 2.

Zoom In Zoom Out Reset image size

Overall, the results in Figure 3 show that the 8 minute PSR-NS system can provide a tight bound for the non-GR theories of negative PN orders when n ≳ −4. However, when comparing these results in a strict way, one should be aware that the prediction may depend on the details of specific theories. NSs and BHs can behave quite differently in alternative gravity theories. We leave this to a future detailed study. In any case, the system that we considered here will provide very useful and complementary constraints. We underline that here the corrections of different PN orders are relative to GR's quadrupole radiation, not to the equations of motion.

4.2. Projected Constraints

In Figure 3, we show the constraint from an 8 minute PSR-NS system on parameter Γ as a function of PN orders. In Section 4.1, the generic formalism (Equation (13)) offers a mapping between generic non-GR parameters in the orbital decay rate and theoretical parameters in various alternative theories of gravity (Seymour & Yagi 2018). We can use this mapping to restrict specific theories of gravity, and then predict how the short-orbital-period PSR-NS systems constrain these theories.

4.2.1. Varying-G Theory

In GR, the Newtonian gravitational constant G does not change with time. But in some alternative theories, the locally measured Newtonian gravitational constant G may vary with time, because of the evolution of the universe (Will 2014). The gravitational constant that varies with time violates local position invariance, which leads to a violation of the strong equivalence principle (SEP). Hence proving whether the gravitational constant varies with time helps us test the SEP. If the gravitational constant is changing with time, it can cause changes in the compactness and Arnowitt–Deser–Misner mass of a compact star, which would result in an additional contribution to the orbital period decay rate ${\dot{P}}_{b}$ (Nordtvedt 1990).

The relation between ${\dot{P}}_{b}$ and $\dot{G}$ is given by Nordtvedt (1990),

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{P}}_{b}^{\dot{G}}}{{P}_{b}}=-2\displaystyle \frac{\dot{G}}{G}\left[1-\displaystyle \frac{2m+{m}_{c}}{2m}{s}_{p}-\displaystyle \frac{2m+{m}_{p}}{2m}{s}_{c}\right],\end{eqnarray} \tag{ 15 }$$

where s_i denotes the "sensitivity" of the NS, which depends on the mass and the EOS of the NS,

$$\begin{eqnarray}&&{s}_{i}\equiv -\displaystyle \frac{\partial \mathrm{ln}{M}_{i}}{\partial \mathrm{ln}G},\quad \quad (i=p,c),\end{eqnarray} \tag{ 16 }$$

where M_i = m_i M_⊙ (notice that in this paper m_i is the mass in units of solar masses). In GR, s_i = 0. We use the value of s_i as per Damour & Esposito-Farèse (1992), ¹⁹

$$\begin{eqnarray}&&{s}_{i}=0.16\left(\displaystyle \frac{{M}_{i}}{1.33{M}_{\odot }}\right).\end{eqnarray} \tag{ 17 }$$

With Equations (14)–(15), we get the corresponding expression of Γ,

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{5}{48}\displaystyle \frac{\dot{G}}{{c}^{3}}\displaystyle \frac{{m}^{3}{M}_{\odot }}{{m}_{p}{m}_{c}}\left[1-\displaystyle \frac{2m+{m}_{c}}{2m}{s}_{p}-\displaystyle \frac{2m+{m}_{p}}{2m}{s}_{c}\right]\displaystyle \frac{1}{f(e)},\end{eqnarray} \tag{ 18 }$$

and it is a −4 PN order correction relative to GR. So combining Equations (13) and (18), we get the upper bound on $\dot{G}$ by

$$\begin{eqnarray}&&\displaystyle \frac{| \dot{G}| }{G}\lt -\displaystyle \frac{\delta }{2}{\left(\displaystyle \frac{\dot{{P}_{b}}}{{P}_{b}}\right)}_{\mathrm{GR}}{\left(1-\displaystyle \frac{2m+{m}_{c}}{2m}{s}_{p}-\displaystyle \frac{2m+{m}_{p}}{2m}{s}_{c}\right)}^{-1}.\end{eqnarray} \tag{ 19 }$$

There is an associated dipole radiation of GWs for most gravitational theories that violate SEP. The dipole radiation of GWs can contribute an extra ${\dot{P}}_{b}^{\mathrm{dipole}}$ that is degenerate with ${\dot{P}}_{b}^{\dot{G}}$. The degeneracy indicates that we need to use multiple binary pulsar systems to break the degeneracy. In our work, for the sake of simplicity, we ignore the possible degeneracy and consider only one violation of GR at a time.

Figure 4(a) shows the predicted upper limit on $\dot{G}$ from short-orbital-period PSR-NS systems with P_b ∈ [8, 120] minutes to be detected by SKA. The best bound is $| \dot{G}| /G\lt 4.8\times {10}^{-14}\,{\mathrm{yr}}^{-1}$, which is from the 8 minute PSR-NS I system with σ_TOA = 100 ns.

Zoom In Zoom Out Reset image size

The simulated fractional uncertainty of ${\dot{P}}_{b}$ varies with P_b , which is largely consistent with what is theoretically predicted (namely $\delta \propto {P}_{b}^{3}$; Damour & Taylor 1992; Lorimer & Kramer 2005). Because of the degeneracy of some delay terms, there are still slight differences between theoretic prediction and our simulation. In Figure 4(a), the black dotted line is the trend from theoretical prediction, $\delta \propto {P}_{b}^{3}$, and we use the PSR-BH I system with σ_TOA = 100 ns to illustrate it. Because we show the constraint on $| \dot{G}| /G$, which is a cumulative effect over one year, there is a noticeable deviation in P_b ∈ [40, 60] minutes.

We do not plot the bounds of a PSR-BH system from Seymour & Yagi (2018), and their results are similar to ours. Although the 8 minute orbital period system has a very high precision for ${\dot{P}}_{b}$, as mentioned before, the varying-G theory is a −4 PN order correction relative to GR. A system with a large v is disadvantageous to provide a tighter bound on $\dot{G}$, so the two systems (PSR-NS and PSR-BH) have similar results.

To provide a more complete analysis, we also consider the fractional error that contains error contributions from ${\dot{P}}_{b}^{\mathrm{Shk}}$ and ${\dot{P}}_{b}^{\mathrm{Gal}}$ to restrict $\dot{G}$ (see Equation (11)). We assume that LISA can help us determine the distance of the PSR-NS system very well. So for the fractional error of D, we adopt 0.01. ²⁰ Considering the error contributions from ${\dot{P}}_{b}^{\mathrm{Shk}}$ and ${\dot{P}}_{b}^{\mathrm{Gal}}$, we show in Figure 4(a) the bound from PSR-NS I systems with δ_TOA = 100 ns on $\dot{G}$ with the red dotted line. We find that by introducing the errors of ${\dot{P}}_{b}^{\mathrm{Shk}}$ and ${\dot{P}}_{b}^{\mathrm{Gal}}$, our results become worse. Although our distance accuracy is high, the precision of other parameters that we adopt from PSR J0737−3039 are low, which prevents us from limiting the non-GR parameters to a high precision. If we, likely with SKA, further improve the measurement accuracy of these parameters in the future, we will constrain $\dot{G}$ better.

For comparison, we show the current bounds from the solar system experiments in Figure 4(a). We notice that the estimated bound of an 8 minute PSR-NS system is slightly weaker than the solar system result from NASA Messenger (red dashed line; Genova et al. 2018). For now, the best constraint comes from the solar system experiment, which is from a weak field. However, the sensitivity s_i is a body-dependent quantity, so it could lead to an enhancement of $\dot{G}$ for strongly self-gravitating bodies. NSs are strongly self-gravitating bodies, so it provides the constraint on $\dot{G}$ from a strong field. The constraint on $\dot{G}$ from PSR-NS systems can provide a complementary constraint compared with the constraint from the solar system experiment. One point that needs to be emphasized is that we can improve the precision of the parameters through long-term observations. For ${\dot{P}}_{b}$, the expected dependence of fractional uncertainty on observing time, T_obs, is proportional to ${T}_{\mathrm{obs}}^{-5/2}$. Therefore, according to the theoretical prediction, compared with a 4 yr observation, the measurement precision of ${\dot{P}}_{b}$ will increase by about 10 times for a 10 yr timing. We simulate a 10 yr timing for the PSR-NS I system with δ_TOA = 100 ns by using the same cadence, and the result is consistent with the theoretical prediction; that is, the measurement precision is improved by an order of magnitude. Hence, a 10 yr observation could increase the projected limit to $| \dot{G}| /G\lt 5\times {10}^{-15}\,{\mathrm{yr}}^{-1}$ for an 8 minute PSR-NS I system with δ_TOA = 100 ns. Thus, an 8 minute PSR-NS system has the potential to give a tighter bound than the present result from the solar system for $\dot{G}$.

We also provide the previous upper limit on $\dot{G}$ from binary pulsar systems that combined the results from PSRs J0437−4715, J1738+0333, and J1713+0747 (blue dashed line; Zhu et al. 2019). Our simulated results show that the constraint on $\dot{G}$ in a strong field can be improved significantly.

4.2.2. Fierz−Pauli-like Massive Gravity

In GR, gravitation is mediated by a massless spin-2 graviton, and the graviton is massless in most theories of gravity. The existence of the graviton is not experimentally confirmed yet. There are some theories of gravity in which the graviton has mass and they cannot be completely ruled out yet. So we can use different theories to bound the graviton mass (de Rham et al. 2017), and we show some results in Figure 4(b). If the graviton is massive, the gravitational potential possesses an exponential Yukawa suppression. In the solar system, which can be viewed as a static field, the Yukawa potential can cause planets to behave differently from the Newtonian potential, which is ∝ 1/r. Bernus et al. (2020) used the planetary ephemeris INPOP19a to constrain the graviton mass and provided m_g < 3.62 × 10 eV/c², shown with a red dashed line in Figure 4(b). The constraint on m_g from the solar system is a bound in a weak field; we also provide two constraints from a strong field. The first bound comes from GW data. If gravity is propagated by a massive field, the massive graviton can lead to a modified dispersion relation. In other words, the velocity of the graviton depends on the GW frequency. Using the modified dispersion relation, the LIGO/Virgo Collaboration bounds the graviton mass to be m_g < 1.8 × 10⁻²³ eV c⁻² by combining the data of 31 GW events from the catalog GWTC-2 (Abbott et al. 2021). This result is represented by a blue dashed line in Figure 4(b). This bound coming from the direct detection of GWs is an effect of GWs' phase accumulation. The second bound is from binary pulsar systems, which is based on a Fierz−Pauli-like massive gravity. Finn & Sutton (2002) presented a phenomenological term in the Lagrangian for linearized gravity with a mass term $\sim {m}_{g}^{2}\left({h}_{\mu \nu }^{2}-\tfrac{1}{2}{h}^{2}\right)$, which is similar to the Fierz−Pauli mass term $\sim {m}_{g}^{2}\left({h}_{\mu \nu }^{2}-{h}^{2}\right)$; (Fierz & Pauli 1939), but with a different coefficient for h². The choice of this linear mass term ensures that the theory reduces to GR when m_g → 0 (avoiding the van Dam−Veltman−Zakharov (vDVZ) discontinuity) and the derived wave equation takes the standard form (Finn & Sutton 2002)

$$\begin{eqnarray}&&\left(\square -{m}_{g}^{2}\right){\bar{h}}_{\mu \nu }+16\pi {T}_{\mu \nu }=0,\end{eqnarray} \tag{ 20 }$$

with an h-independent source T_{μ ν} . We have defined ${\bar{h}}_{\mu \nu }\equiv {h}_{\mu \nu }-\tfrac{1}{2}{\eta }_{\mu \nu }h$. Finn & Sutton (2002) worked out the relation between m_g and δ,

$$\begin{eqnarray}&&{m}_{g}^{2}\lt \displaystyle \frac{24}{5}F(e){\left(\displaystyle \frac{2\pi {\hslash }}{{c}^{2}{P}_{b}}\right)}^{2}\delta ,\end{eqnarray} \tag{ 21 }$$

where $F(e)=f(e){\left(1-{e}^{2}\right)}^{1/2}$ (see Equation (8) for f(e)), and ℏ is the reduced Planck constant. With Equation (21), we can use binary pulsars to restrict m_g . Finn & Sutton (2002) provided a method to constrain the graviton mass in a dynamic regime. The Fierz−Pauli-like massive gravity gives a −3 PN order correction relative to GR. By comparison with Equation (14), we can obtain the form of Γ,

$$\begin{eqnarray}&&{\rm{\Gamma }}\lt \displaystyle \frac{5}{24}{m}_{g}^{2}{m}^{2}\displaystyle \frac{{M}_{\odot }^{2}{G}^{2}}{{c}^{2}{{\hslash }}^{2}}\displaystyle \frac{1}{F(e)}.\end{eqnarray} \tag{ 22 }$$

Finn & Sutton (2002) used PSRs B1913+16 and B1534+12 to bound m_g . Recently, using nine well-timed binary pulsars and the Bayesian theorem, Miao et al. (2019) provided an improved bound, m_g < 5.2 × 10⁻²¹ eV c⁻², which is shown with a light blue dashed line in Figure 4(b).

Here, we use the simulated δ of ${\dot{P}}_{b}$ of the short-orbital-period PSR-NS systems by SKA, to predict the constraints on the graviton mass. The results of the bounds are shown in Figure 4(b). The results show that the PSR-NS system with an 8 minute short orbital period can give a tighter limit than previous results in Miao et al. (2019) by roughly a factor of 4: m_g < 1.2 × 10⁻²¹ eV c⁻². Although the projected bounds from PSR-NS systems are looser than the bounds from GWs and the solar system, they reflect a dynamical process of two compact objects and a strong field effect. The constraint from GWs is based on a modified dispersion relation that is an effect of the accumulation in GW phases, so it is a kinematic effect from propagation. The constraint from the solar system reflects an effect from a static weak field. Hence the restrictions from the binary pulsars are complementary to these constraints. For an 8 minute PSR-NS I system with δ_TOA = 100 ns, simulated results for a 10 yr observation span improve it to m_g < 4.0 × 10⁻²² eV c⁻². Because m_g ∝ δ^1/2, the increase in the precision of ${\dot{P}}_{b}$ cannot improve the limit on m_g significantly in this case. We also show the trend of the theoretical prediction, $\delta \propto {P}_{b}^{3}$, with the black dotted line in Figure 4(b). Our results do not significantly deviate from the theoretical trend. In addition, we give the restriction on m_g when there is an error contribution from ${\dot{P}}_{b}^{\mathrm{Shk}}$ and ${\dot{P}}_{b}^{\mathrm{Gal}}$, and it is shown with a red dotted line in Figure 4(b).

4.2.3. Cubic Galileon Massive Gravity

If the graviton is massive, there may be extra scalar and vector degrees of freedom. The vector degrees of freedom can be ignored, because they usually do not couple to matter in the decoupling limit (de Rham et al. 2017). The extra scalar degree of freedom can couple to matter, which leads to a fifth force, and the fifth force would make the theory of massive gravity fail to recover GR when m_g → 0 with the so-called vDVZ discontinuity. However the vDVZ discontinuity can be solved by introducing screening mechanisms, e.g., the Vainshtein mechanism (Vainshtein 1972). The scalar degree of freedom becomes strongly coupled and is strongly suppressed, because of the Vainshtein mechanism, so the deviation from GR is suppressed. In other words, the fifth force is suppressed and the metric field propagates the gravitational force.

The simplest model used to display the Vainshtein mechanism is called the cubic Galileon (Luty et al. 2003), which is a Lorentz invariant massive gravity model. The massive graviton causes extra GW radiation, and there is not only quadrupole radiation but also monopole and dipole radiation that need to be considered. de Rham et al. (2013) calculated the contributions of these three radiation types, and found that the dipole radiation is orders of magnitude smaller than the other two types and the quadrupole radiation has the largest contribution in the vast majority of cases. In fact, we calculated the contribution of monopole radiation and found that it can be ignored, indeed. So we just consider quadrupole radiation in our work.

The extra contribution to ${\dot{P}}_{b}$ from the quadrupole radiation power in the cubic Galileon model is (de Rham et al. 2013)

$$\begin{eqnarray}&&{\dot{P}}_{b}^{\mathrm{quad}}=-\displaystyle \frac{15}{\sqrt{2}}{\lambda }^{2}{\pi }^{3/2}\displaystyle \frac{{\left(\mathrm{GM}\right)}^{5/6}}{{n}_{b}^{1/6}{\hslash }{c}^{1/2}}\displaystyle \frac{{m}_{p}{m}_{c}{{ \mathcal Y }}^{2}}{{m}^{2}}{m}_{g}{f}_{\mathrm{cg}}(e),\end{eqnarray} \tag{ 23 }$$

where ${ \mathcal Y }\equiv \left({m}_{p}^{1/2}+{m}_{c}^{1/2}\right)/{m}^{1/2}$, the constant λ = 0.2125 (Shao et al. 2020), and ${f}_{\mathrm{cg}}(e)={\sum }_{n=0}^{\infty }{\left|{I}_{n}^{\mathrm{quad}}(e)\right|}^{2}$,

$$\begin{eqnarray}&&{I}_{n}^{\mathrm{quad}}(e)=\displaystyle \frac{{\left(1-{e}^{2}\right)}^{3/2}}{2\pi }{n}^{7/4}{\int }_{0}^{2\pi }\displaystyle \frac{\exp \left[i(2-n)x\right]}{{\left(1+e\cos x\right)}^{3/2}}{\rm{d}}x.\end{eqnarray} \tag{ 24 }$$

We choose n = 30 in our calculation, which is enough for e = 0.1 (Shao et al. 2020). With Equation (13) and Equation (23), we can provide the constraint on m_g ,

$$\begin{eqnarray}&&{m}_{g}\lt \displaystyle \frac{8}{5{\lambda }^{2}}\displaystyle \frac{1}{{f}_{\mathrm{cg}}(e)}\displaystyle \frac{{m}^{1/2}{M}_{\mathrm{Pl}}}{{M}_{\mathrm{quad}}^{2}}\displaystyle \frac{{{ \mathcal P }}_{\mathrm{GR}}}{{n}_{b}^{1/2}{\left({n}_{b}a\right)}^{3}}{\left({\hslash }{M}_{\odot }\right)}^{1/2}\delta ,\end{eqnarray} \tag{ 25 }$$

where n_b = 2π/P_b is the orbital angular frequency, a is the semimajor axis, ${M}_{\mathrm{Pl}}\equiv {\left({\hslash }c/8\pi G\right)}^{1/2}$ is the reduced Planck mass, ${{ \mathcal P }}_{\mathrm{GR}}$ is the radiation power of GR,

$$\begin{eqnarray}&&{{ \mathcal P }}_{\mathrm{GR}}=\displaystyle \frac{32m}{5}\displaystyle \frac{{m}_{p}^{2}{m}_{c}^{2}}{{a}^{5}}\displaystyle \frac{{G}^{4}{M}_{\odot }^{5}}{{c}^{5}}f(e),\end{eqnarray} \tag{ 26 }$$

and

$$\begin{eqnarray}&&{M}_{\mathrm{quad}}\equiv \displaystyle \frac{{m}_{p}{m}_{c}}{m}{ \mathcal Y }{M}_{\odot }.\end{eqnarray} \tag{ 27 }$$

The cubic Galileon massive gravity model provides a $-\tfrac{11}{4}$ PN order correction relative to GR, and Γ is

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{25{\lambda }^{2}}{{2}^{8}}\displaystyle \frac{{M}_{\mathrm{quad}}^{2}}{{M}_{\mathrm{Pl}}}\displaystyle \frac{{m}^{3}{m}_{g}}{{m}_{p}^{2}{m}_{c}^{2}}{\left(\displaystyle \frac{G}{{\hslash }c}\right)}^{1/2}\displaystyle \frac{{f}_{\mathrm{cg}}(e)}{f(e)}\displaystyle \frac{1}{{M}_{\odot }}.\end{eqnarray} \tag{ 28 }$$

With Equation (25), we show constraints from our simulated results in Figure 4(c). The light blue dashed line is the bound from a combination of binary pulsar systems, which gives m_g < 2.0 × 10⁻²⁸ eV c⁻² (Shao et al. 2020). Compared with it, the short-orbital-period PSR-NS systems can greatly improve the constraint on m_g , and a projected bound from an 8 minute PSR-NS system is m_g < 1.9 × 10⁻³¹ eV c⁻². Because m_g ∝ δ, the deviation with respect to the trend predicted by theory (black dotted line) is more significant than that of the Fierz−Pauli-like massive gravity. The red dotted line in Figure 4(c) is the restriction that considers the possible error contributions from ${\dot{P}}_{b}^{\mathrm{Shk}}$ and ${\dot{P}}_{b}^{\mathrm{Gal}}$. In Figure 4(c), the best restriction is from the red dashed line, m_g < 10⁻³² eV c⁻² (Dvali et al. 2003). It is a bound from the lunar laser ranging experiments in the weak field. For a 10 yr observation, an 8 minute PSR-NS I system with δ_TOA = 100 ns can provide a bound of m_g < 2.0 × 10⁻³² eV c⁻², which is comparable to that of the lunar laser ranging experiments. So if we can observe these systems for decades, we could provide a tighter constraint than the current result. So the short-orbital-period PSR-NS systems have the potential to significantly improve the restriction on m_g in the strong field.

4.2.4. Scalar−tensor Theory

In GR, gravity is mediated only by a long-range, spin-2 tensor field, g_{μ ν} , but in alternative theories there can be scalar or vector fields involved in the gravitational interaction. The scalar−tensor gravity, which adds a spin-0 scalar field, φ, is the most natural alternative theory. In our work, we take the Damour and Esposito-Farése (DEF) theory as an example (Damour & Esposito-Farèse 1992, 1993, 1996). In the weak field, the DEF gravity theory can be reduced to GR, which means that it can pass the solar system tests. However, in the strong field, DEF theory can be completely different from GR, due to a strong field spontaneous scalarization phenomenon (Damour & Esposito-Farèse 1993, 1996), so binary pulsar systems provide an ideal laboratory to study the DEF theory.

In the DEF theory, the GW damping is different from that in GR. There is an additional contribution of GW radiation from the scalar field, which leads to an extra orbital decay rate ${\dot{P}}_{b}^{\mathrm{scalar}}$, and it contains monopole, dipole, and quadrupole contributions. We only consider the dipole contribution which is a −1 PN correction relative to GR, and the expression of ${\dot{P}}_{b}^{\mathrm{dipole}}$ is (Damour & Esposito-Farèse 1993, 1996)

$$\begin{eqnarray}&&{\dot{P}}_{b}^{\mathrm{dipole}\ }=-2\pi {T}_{\odot }^{* }d(e)\left(\displaystyle \frac{2\pi }{{P}_{b}}\right)\displaystyle \frac{{m}_{p}{m}_{c}}{{m}_{p}+{m}_{c}}{\left({\alpha }_{p}-{\alpha }_{c}\right)}^{2},\end{eqnarray} \tag{ 29 }$$

where α_i is the coupling strength between the scalar field and matter, $d(e)\equiv \left(1+\tfrac{1}{2}{e}^{2}\right){\left(1-{e}^{2}\right)}^{-5/2}$, ${T}_{\odot }^{* }={G}^{* }{M}_{\odot }/{c}^{3}$, ${G}^{* }=G/\left(1+{\alpha }_{0}^{2}\right)$ ²¹ , and α₀ is the linear coupling of the scalar field to matter fields. For simplification, we ignore the difference between G* and G. Because the current weak field constraint on ${\alpha }_{0}^{2}$ has reached 10⁻⁵ (Bertotti et al. 2003), we set G* = G. The effective scalar couplings α_p and α_c depend on the specific EOS of NSs.

In our work, we only use the simulated fractional error to constrain $\left|{\alpha }_{p}-{\alpha }_{c}\right|$, and we do not further analyze α_p by using specific EOSs. ²² The spontaneous scalarization of the DEF theory happens in a strong field, so it is conducive to test the DEF theory in PSR-NS systems. With Equation (13) and Equation (29), we can provide the constraint of $\left|{\alpha }_{p}-{\alpha }_{c}\right|$,

$$\begin{eqnarray}&&{\left({\alpha }_{p}-{\alpha }_{c}\right)}^{2}=\displaystyle \frac{96}{5}{T}_{\odot }^{2/3}{\left(\displaystyle \frac{{P}_{b}}{2\pi }\right)}^{-2/3}D(e){\left({m}_{p}+{m}_{c}\right)}^{2/3}\delta ,\end{eqnarray} \tag{ 30 }$$

where D(e) ≡ f(e)/d(e). Using Equation (30), we provide the results in Figure 4(d). We notice that the best constraint on $\left|{\alpha }_{p}-{\alpha }_{c}\right|$ can reach 1.6 × 10⁻⁶. In the work of Shao & Wex (2016), they showed constraints of ∼10⁻³ for several best measured PSR-WD systems where $\left|{\alpha }_{p}-{\alpha }_{c}\right|=\left|{\alpha }_{p}-{\alpha }_{0}\right|$ due to the companion being a WD. ²³ The simulated 8 minute PSR-NS system seems to provide a tighter bound on $\left|{\alpha }_{p}-{\alpha }_{c}\right|$, but the companion in the system is an NS, so α_c could be very close to α_p , which makes $\left|{\alpha }_{p}-{\alpha }_{c}\right|$ small for most cases. Our predicted constraints on $\left|{\alpha }_{p}-{\alpha }_{c}\right|$ come from a simplified analysis, and the real situation needs to be discussed with the EOS dependence included. In the figure, we also provide the trend of the prediction from theory (black dotted line) and the results containing the possible error contributions from ${\dot{P}}_{b}^{\mathrm{Shk}}$ and ${\dot{P}}_{b}^{\mathrm{Gal}}$ (red dotted line).

In Section 4, we have only used the PK parameter ${\dot{P}}_{b}$ to test gravity theories, but in our simulated results, the fractional error of $\dot{\omega }$ could reach the highest precision compared with the other four PK parameters. So, using $\dot{\omega }$ to test gravity is also powerful. There is another advantage of using $\dot{\omega }$ to test gravity, in that the observed $\dot{\omega }$ of an isolated binary pulsar system is relatively "clean," because, in general, the nonintrinsic contribution to $\dot{\omega }$ can be ignored. So our simulated fractional error of $\dot{\omega }$ is probably more realistic. In this section we provide a possible application of the high-precision $\dot{\omega }$ measurement, just for an illustrative purpose.

In relativistic gravity theories, which are different from Newtonian gravity, a rotating body in a binary system can have an effect on spin and orbital dynamics. For PSR-NS systems, the two components have proper rotation, and in leading order contributions we only consider the effect from spin−orbital (SO) coupling (Wex 2014). The relativistic SO coupling in PSR-NS systems can lead the orbital and the two spins to precess. For the precession of spin, the main effect is a secular change in the orientation of spin and this effect is due to the spacetime curvature and is spin-independent. We generally call this effect geodetic precession and it has been detected in some binary pulsar systems, such as PSR J0737−3039 (Breton et al. 2008). The precession of orbit contains two effects that need to be considered. The first effect is due to the mass−mass interaction. The second effect is due to the Lense−Thirring effect related to the spins of the two stars.

Damour & Schaefer (1988) decomposed the effect of precession of the orbit in terms of observable parameters of pulsars and calculated the contribution to periastron advance ${\dot{\omega }}_{p}$,

$$\begin{eqnarray}&&{\dot{\omega }}_{p}={\dot{\omega }}_{\mathrm{PN}}+{\dot{\omega }}_{{\mathrm{LT}}_{p}}+{\dot{\omega }}_{{\mathrm{LT}}_{c}},\end{eqnarray} \tag{ 31 }$$

where ${\dot{\omega }}_{\mathrm{PN}}$ is caused by the mass−mass interaction, which results from the two bodies of the binary system, and ${\dot{\omega }}_{{\mathrm{LT}}_{p}}$ and ${\dot{\omega }}_{{\mathrm{LT}}_{c}}$ are due to the Lense−Thirring effect of the pulsar and the companion, respectively. For the actual observed value of ${\dot{\omega }}_{p}^{\mathrm{obs}}$, it should include the Kopeikin term, which is caused by the proper motion of a binary system (Kopeikin 1996), but in general, the contribution of the Kopeikin term is small and can be ignored or subtracted (Bagchi 2018; Hu et al. 2020). The expression of the first term in Equation (31) is

$$\begin{eqnarray}&&{\dot{\omega }}_{\mathrm{PN}}=3{T}_{\odot }^{2/3}{\left(\displaystyle \frac{{P}_{b}}{2\pi }\right)}^{-5/3}\displaystyle \frac{{\left({m}_{p}+{m}_{c}\right)}^{2/3}}{1-{e}^{2}}\left(1+{f}_{{\rm{O}}}{v}^{2}\right),\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}\begin{array}{rcl}{f}_{{\rm{O}}} & = & \displaystyle \frac{1}{1-{e}^{2}}\left(\displaystyle \frac{3}{2}{X}_{p}^{2}+\displaystyle \frac{3}{2}{X}_{p}+\displaystyle \frac{27}{4}\right)\\ & & +\left(\displaystyle \frac{5}{6}{X}_{p}^{2}-\displaystyle \frac{23}{6}{X}_{p}-\displaystyle \frac{1}{4}\right),\end{array}\end{eqnarray} \tag{ 33 }$$

where X_p ≡ m_p /m. The first term on the right side of Equation (32) is the 1 PN term of ${\dot{\omega }}_{\mathrm{PN}}$ and is given by Equation (3); the second term is the 2 PN term of ${\dot{\omega }}_{\mathrm{PN}}$. We calculate the fractional value of the 2 PN term for the 8 minute PSR-NS system II and get f_O v² ≈ 2 × 10⁻⁴, so the contribution from ${\dot{\omega }}_{2\mathrm{PN}}$ is four orders of magnitude smaller than ${\dot{\omega }}_{1\mathrm{PN}}$. Our simulated fractional error of ${\dot{\omega }}_{1\mathrm{PN}}$ can reach 10⁻⁸ for an 8 minute PSR-NS system with σ_TOA = 1 μs, so the contribution from ${\dot{\omega }}_{2\mathrm{PN}}$ should be detectable for an 8 minute binary.

The expressions of the second and third terms in Equation (31) are (Bagchi 2018; Hu et al. 2020)

$$\begin{eqnarray}\begin{array}{rcl}{\dot{\omega }}_{{\mathrm{LT}}_{p}}+{\dot{\omega }}_{{\mathrm{LT}}_{c}} & = & -3{T}_{\odot }^{2/3}{\left(\displaystyle \frac{{P}_{b}}{2\pi }\right)}^{-5/3}\displaystyle \frac{{\left({m}_{p}+{m}_{c}\right)}^{2/3}}{1-{e}^{2}}v\\ & & \times \left({g}_{p}^{{\rm{s}}}{\beta }_{p}^{{\rm{s}}}+{g}_{c}^{{\rm{s}}}{\beta }_{c}^{{\rm{s}}}\right),\end{array}\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{\beta }_{p}^{{\rm{s}}}=\displaystyle \frac{{{cI}}_{p}}{{{GM}}_{p}^{2}}\displaystyle \frac{2\pi }{{P}_{p}},\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}\begin{array}{rcl}{g}_{p}^{{\rm{s}}} & = & \displaystyle \frac{{X}_{p}\left(4{X}_{p}+3{X}_{c}\right)}{6{\left(1-{e}^{2}\right)}^{1/2}{\sin }^{2}i}\left[\left(3{\sin }^{2}i-1\right)\cos {\chi }_{p}\right.\\ & & \left.+\cos i\cos {\lambda }_{p}\right],\end{array}\end{eqnarray} \tag{ 36 }$$

where I_p is the moment of inertia of the pulsar and depends on the EOS of NSs, and P_p is the spin period of the pulsar. χ_p is the angle between s _p and k , where k is the unit vector along the orbital angular momentum and s _p is the unit spin vector of the pulsar; λ_p is the angle between s _p and h _p where h _p is the unit vector along the line of sight. To get ${\beta }_{c}^{{\rm{s}}}$ and ${g}_{c}^{{\rm{s}}}$, we just need to interchange the subscript "c" and "p." Generally, the pulsar in a PSR-NS system is a recycled pulsar and the companion is a normal pulsar, so the rotation of the companion is much slower than the pulsar, and we can ignore the contribution of ${\dot{\omega }}_{{\mathrm{LT}}_{c}}$; see, e.g., PSR J0737−3039 (Kramer & Wex 2009).

To estimate the magnitude of ${\dot{\omega }}_{{\mathrm{LT}}_{p}}$ contributed by the Lense−Thirring effect, we simplify the calculation a little bit. If s _p is parallel to k , namely χ_p = 0 and λ_p = i, Equation (36) can be simplified to

$$\begin{eqnarray}&&{g}_{p\,\parallel }^{{\rm{s}}}=\displaystyle \frac{\left(\tfrac{1}{3}{X}_{p}^{2}+{X}_{p}\right)}{{\left(1-{e}^{2}\right)}^{1/2}}.\end{eqnarray} \tag{ 37 }$$

With these equations, we can estimate the contribution from ${\dot{\omega }}_{{\mathrm{LT}}_{p}}$ for periastron advance ${\dot{\omega }}_{p}$. For the 8 minute PSR-NS system II,

$$\begin{eqnarray}&&{\dot{\omega }}_{{\mathrm{LT}}_{p}}=-0.14\,\deg \,{\mathrm{yr}}^{-1}\left(\displaystyle \frac{{I}_{p}}{{10}^{45}\,{\rm{g}}\,{\mathrm{cm}}^{2}}\right)\left(\displaystyle \frac{20\,\mathrm{ms}}{{P}_{p}}\right).\end{eqnarray} \tag{ 38 }$$

So if I_p = 10⁴⁵ g cm² and P_p = 20 ms, we have ${\dot{\omega }}_{{\mathrm{LT}}_{p}}=-0.14\,\deg \,{\mathrm{yr}}^{-1}$. For the 8 minute PSR-NS system II, $\left|{\dot{\omega }}_{{\mathrm{LT}}_{p}}/{\dot{\omega }}_{1\mathrm{PN}}\right|\approx 6.1\times {10}^{-5}$. Remember that the fractional error of $\dot{\omega }$ can reach 10⁻⁸ with σ_TOA = 1 μs in our simulation. So the contribution of ${\dot{\omega }}_{{\mathrm{LT}}_{p}}$ can be significant compared with the error of ${\dot{\omega }}^{\mathrm{obs}}$, and it means that we are very likely to detect the value of ${\dot{\omega }}_{{\mathrm{LT}}_{p}}$ contributed by the Lense−Thirring effect. Notice that the value of ${\dot{\omega }}_{{\mathrm{LT}}_{p}}$ depends on the moment of inertia I_p , as shown in Equation (35). When m_p and spin are determined, I_p only depends on the EOS. So if we can separate ${\dot{\omega }}_{{\mathrm{LT}}_{p}}$, it would help us restrict the EOS. This is another way to restrict the EOS in addition to using the maximum mass of NSs (Lattimer & Schutz 2005; Hu et al. 2020).

In our work, we assume that one can use the LISA and SKA multimessenger strategy (Kyutoku et al. 2019; Andrews et al. 2020; Lau et al. 2020) to detect short-orbital-period PSR-NS systems with ${P}_{b}\sim \min$ in the future. The short-orbital-period PSR-NS systems can provide extreme gravity tests that may show new information about gravity. We use the plug-in FAKE of TEMPO2 to simulate the TOAs of short-orbital-period PSR-NS systems and fit TOAs to get the fractional errors of five PK parameters. Although we do not consider the higher PN order contribution in our simulation, considering that v ≪ c, the main contribution from the lowest PN order can provide reliable errors for parameters that can help us constrain gravity theories. As our results show, compared with the previous bounds from binary pulsar systems, short-orbital-period systems can significantly improve the constraints on parameters of some specific gravity theories. The constraints of parameters we predicted are looser than the bound from the solar system, but highly relativistic PSR-NS systems give constraints from a strong field, which is a complementary gravity regime. Our simulation is simplified and we expect to provide a more complete simulation in the future, such as to study a more complete timing model.

For our simulation, we choose two different groups of masses (see Table 1). We found that the PSR-NS I system (m_p = 1.3 M_⊙, m_p = 1.7 M_⊙) can provide a slightly tighter bound than the PSR-NS II system (m_p = 1.35 M_⊙, m_p = 1.44 M_⊙) with the same σ_TOA. It indicates that PSR-NS systems with larger mass differences are more suitable for limiting these theories of gravity that we discussed, in particular for the dipole radiation in the DEF gravity, where the asymmetry of masses will play an important role (see, e.g., Zhao et al. 2021 in the context of GWs).

We simulate five PK parameters, and among them, r and s are parameters related to the Shapiro delay. The precision of r and s depends on the orbital inclination i, which we have set to i = 60° in our simulation. If the PSR-NS system is nearly edge-on (i ∼ 90°), the effect of the Shapiro delay can be more significant, such as in PSR J0737−3039 (Kramer et al. 2006), and it can make r and s more precise in an actual measurement. The higher precision of r and s can help us measure the masses of PSR-NS systems more precisely and test GR more stringently.

In the future, TianQin (Luo et al. 2016) and Taiji (Luo et al. 2020) could also join the multimessenger observation strategies. The sensitive frequency of Taiji is slightly lower than that of LISA, which allows the joint observation of Taiji and SKA to detect more PSR-NS systems with a bit larger P_b . In addition, the mission times of TianQin or Taiji might be different from that of LISA, so in that case we can have a longer total time to detect short-orbital-period PSR-NS systems by different combinations of joint observation strategies. This may increase the detection numbers.

On the astrophysical side, if we can detect short-orbital-period PSR-NS systems in the future, it can fill in the gap in current observations and allow us to better model the DNS population. Detecting such systems can help us understand the formation of DNSs better, and possibly verify new formation channels of DNSs, such as a fast-merging channel (Andrews et al. 2020).

Furthermore, if LISA and SKA can detect a PSR-NS system, as we discussed earlier, LISA can help SKA to locate and detect PSR-NS systems. Meanwhile, SKA can provide a more accurate measurement of parameters, which in turn, could help us analyze the GW waveforms of PSR-NS systems in the LISA data. We leave these studies to future publications.

We thank the anonymous referee for constructive comments that improved the work. This work was supported by the National SKA Program of China (2020SKA0120300), the National Natural Science Foundation of China (11975027, 11991053, 11721303), the Young Elite Scientists Sponsorship Program by the China Association for Science and Technology (2018QNRC001), and the Max Planck Partner Group Program funded by the Max Planck Society. It was partially supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB23010200), and the High-Performance Computing Platform of Peking University.

Facilities: LISA - , SKA. -

Software: TEMPO2 (Hobbs et al. 2006; Edwards et al. 2006).