Constraining the dense matter equation-of-state with radio pulsars

Radio pulsars provide some of the most important constraints for our understanding of matter at supranuclear densities. So far, these constraints are mostly given by precision mass measurements of neutron stars (NS). By combining single measurements of the two most massive pulsars, J0348$+$0432 and J0740$+$6620, the resulting lower limit of 1.98 $M_\odot$ (99% confidence) of the maximum NS mass, excludes a large number of equations of state (EOSs). Further EOS constraints, complementary to other methods, are likely to come from the measurement of the moment of inertia (MOI) of binary pulsars in relativistic orbits. The Double Pulsar, PSR J0737$-$3039A/B, is the most promising system for the first measurement of the MOI via pulsar timing. Reviewing this method, based in particular on the first MeerKAT observations of the Double Pulsar, we provide well-founded projections into the future by simulating timing observations with MeerKAT and the SKA. For the first time, we account for the spin-down mass loss in the analysis. Our results suggest that an MOI measurement with 11% accuracy (68% confidence) is possible by 2030. If by 2030 the EOS is sufficiently well known, however, we find that the Double Pulsar will allow for a 7% test of Lense-Thirring precession, or alternatively provide a $\sim3\sigma$-measurement of the next-to-leading order gravitational wave damping in GR. Finally, we demonstrate that potential new discoveries of double NS systems with orbital periods shorter than that of the Double Pulsar promise significant improvements in these measurements and the constraints on NS matter.


INTRODUCTION
Neutron stars (NSs) are among the most compact and exotic objects in nature, comprised of extraordinarily dense matter that is not accessible in laboratory experiments. Determining the properties and structure of the cold dense matter inside NSs is therefore a tremendous challenge in nuclear physics. Thus far, a variety of equations of state (EOSs) have been proposed to describe the pressure -density relation inside NSs (see e.g., Lattimer & Prakash 2001. Constraining the EOS is crucial for understanding aspects of fundamental physics, such as the internal structure of NSs, the dynamics of binary mergers, and r-process nucleosynthesis (for a recent review see Özel & Freire 2016).
Various observational methods have emerged to measure the macroscopic properties of NSs, which promise to increase our knowledge of the EOS. The gravitational wave (GW) observation of a binary NS merger with LIGO/Virgo offers the possibility of measuring the tidal deformability (Abbott et al. 2017(Abbott et al. , 2018. X-ray observations of emissions from the hot regions on NS surface with E-mail: huhu@mpifr-bonn.mpg.de NICER (Watts et al. 2016) allows a joint mass-radius estimation (Riley et al. 2019;Miller et al. 2019).
The largest number of known NSs, however, can be observed as radio pulsars. Currently about 3000 pulsars are known, and the ability of radio astronomers to measure pulsar properties precisely via a technique known as "pulsar timing", suggests that important information about the EOS of NSs can also be derived from such measurements. This is indeed the case. The most direct and best known route is to measure the masses of NSs precisely. This is possible in binary pulsars using relativistic orbital effects, potentially combined with other information. The mass range, especially the maximum mass observed, must obviously be consistent with the range of masses supported by a proposed EOS. In addition, there are other orbital effects that also offer the possibility of measuring the moment of inertia (MOI) in binary pulsars via relativistic spinorbit coupling, as was first suggested by Damour & Schäfer (1988). The MOI of a NS depends crucially on the EOS and hence allows us to constrain or even identify it (Morrison et al. 2004;Lattimer & Schutz 2005;Greif et al. 2020). Accessing the MOI of isolated NSs, in contrast, may be possible if one can reliably derive or measure the total loss in rotational energy, E, which relates the MOI with the observed period and period derivative.
In this work, we provide insight into the various methods using binary pulsars and their current status in Section 2, before we focus specifically on the possibility of using the Double Pulsar (Burgay et al. 2003;Lyne et al. 2004) for MOI measurements. We will provide an in-depth study of the relevant factors in Section 3, where we explain how Lense-Thirring (LT) precession affects the periastron advance. Section 4 describes the intrinsic and extrinsic contributions to the orbital period decay. We describe how we simulate future timing observations in Section 5 and evaluate how the Double Pulsar can measure the MOI and constrain the EOS in Section 6. Prospects of testing LT precession and constraining theories of gravity is discussed in Section 7 by assuming the EOS is known. We investigate potential constraints on next-to-leading order GW damping in Section 8, and potential constraints from future discoveries of more relativistic binary pulsars in Section 9. Finally, we conclude in Section 10.

Mass measurements
A given EOS i can only sustain a NS up to a certain maximum mass, M max i . Finding a massive NS of mass M j , consequently excludes all EOS with M max i < M j . This was, for instance possible, by using a Shapiro delay measurement in PSR J1614−2230, where Demorest et al. (2010) determined a mass M = 1.97 ± 0.04 M . We note that a recent update on continued timing observations (Arzoumanian et al. 2018), implies a significantly lower mass of 1.908 ± 0.016 M for this pulsar. As pointed out by Cromartie et al. (2019), a Shapiro delay measurement and the determined uncertainty can be affected by the exact orbital sampling (see also Hu et al. in prep.).
In 2013, Antoniadis et al. (2013) could determine the mass of PSR J0348+0432 without using a Shapiro delay measurement. They combined radio timing measurements of the orbit of the pulsar with precise spectroscopy data of the white dwarf companion in the optical regime to derive a mass of 2.01 ± 0.04 M , confirming the existence of 2-M NSs via a complementary method.
Recently, Cromartie et al. (2019) used a Shapiro delay measurement in PSR J0740+6620 to determine a pulsar mass of 2.14 +0.10 −0.09 M . We can use the masses of these latter two most massive pulsars, J0348+0432 (fully accounting for the rather asymmetric probability density distribution found by Antoniadis et al. (2013)) and J0740+6620, to obtain a 99% confidence lower limit for the maximum mass of a NS, 1.98 M < M max . Such a constraint already rules out a number of soft EOSs as shown in Figure 1. 1 We can compare this lower limit derived from pulsar timing with an upper limit placed by the NS-NS merger GW170817 observed by LIGO (Abbott et al. 2017). Assuming that the NS-NS merger resulted in the formation of a black hole, one finds an upper limit of about 2.3 M for the maximum mass of a NS (Rezzolla et al. 2018;Shibata et al. 2019). Figure 1. The mass of a NS as function of its radius for different EOS (Lattimer & Prakash 2001). The horizontal bands indicate the 2-σ range for the masses of the two most massive radio pulsars known to date, PSR J0348+0432 (Antoniadis et al. 2013) in blue and PSR J0740+6620 (Cromartie et al. 2019) in yellow.

Relativistic spin-orbit coupling
Unlike in Newtonian gravity, the gravitational field of a body in general relativity (GR) has contributions from the mass currents related to the body's proper rotation. Lense & Thirring (1918) with substantial help from Albert Einstein (see Pfister 2007) -have shown that the rotation of the Sun has, in principle, an effect on the planetary orbits. This relativistic spin-orbit coupling, also known as LT precession, has since been well tested in the gravitational field of the rotating Earth with the help of satellite laser ranging (Ciufolini & Pavlis 2004;Ciufolini et al. 2019). Similarly, in relativistic binaries, the spin of a compact rotating body is expected to couple gravitationally with the orbital motion of the system (Barker & O'Connell 1975), leading to a precession of the orbit, while the total angular momentum is conserved. 2 This LT precession of the orbit is potentially observable, hence providing a route to access the MOI of the pulsar (Damour & Schäfer 1988). An MOI measurement, even with an accuracy of ∼10%, would offer important constraints of the EOS (Morrison et al. 2004;Lattimer & Schutz 2005).
The LT precession of the orbit may be detected via the variation in the orbital inclination angle, i, and hence in the (observable) projected semi-major axis of the pulsar obit, x = a p sin i/c (a p is the semi-major axis, and c the speed of light). However, for this to be detectable, the misalignment angle between pulsar spin and angular momentum vector must be sufficiently large. Also, the orbital inclination angle must not be too close to 90 degrees ("edgeon" geometry), since the precession leads to a contribution to the rate of change of the projected semi-major axis given by where (di/dt) LT is given by Eq.(3.27) in Damour & Taylor (1992). For nearly edge-on systems, i.e. i ≈ 90 deg, this contribution becomes small and most likely undetectable since cot i 1. However, in relativistic binary systems with smaller inclination angles, such as PSR J1757−1854, this measurement appears to be possible by Cameron et al. (2018). To achieve this, two challenges will have to be addressed successfully. Firstly, the precession is expected to cause a variation in the pulse profile with time due to a change in the viewing geometry (e.g. Kramer 1998). Special care in the timing procedure is then needed to obtain sufficient precision and to properly account for possible systematic errors (e.g. Stairs et al. 2002;Bhat et al. 2008;van Leeuwen et al. 2015). Moreover, since we require access to the direction of the pulsar spin vector (Damour & Schäfer 1988;Damour & Taylor 1992), the geometry of the binary system and the pulsar needs also to be measured. This is, however, possible via polarisation measurements, as shown previously (e.g. Kramer 1998;Stairs et al. 2004;Desvignes et al. 2019;Venkatraman Krishnan et al. 2020).
Alternatively, rather than using a contribution to x, one can exploit LT precession also via its contribution to the advance of periastron. If it is possible to isolate the contribution of ω LT from the total periastron advance, then the MOI can be determined (Damour & Schäfer 1988). This method was suggested for the Double Pulsar, PSR J0737−3039A/B (Lyne et al. 2004;Kramer et al. 2006). Kramer & Wex (2009) concluded that a MOI measurement of ∼10% accuracy is possible by ∼2030, with the timing accuracy achievable at the time. Later, Kehl et al. (2017) simulated timing data from emerging telescopes, i.e., the Square Kilometer Array (SKA; e.g. Kramer & Stappers 2015) and its precursor MeerKAT (Camilo 2018;Bailes et al. 2018), which greatly improve the timing precision, and predict a MOI measurement with an accuracy well below 10% by 2030. However, the timeline of the nominal operation of the SKA assumed by Kehl et al. (2017) was optimistic compared to the current estimates. With MeerKAT in operation since about 2018 (albeit initially with limited capability), operations of the first phase of the SKA (SKA 1, initially expected to have about 10% of the full SKA's sensitivity) are not expected before 2027. However, first useful data from commissioning observations may be already available in 2025. 3 In addition, compared to Kehl et al. (2017), we now already have about two years of Double Pulsar timing observations with MeerKAT, and therefore have more realistic numbers for the timing precision and cadence of observations, not only for the current MeerKAT configuration but also for future extensions. Moreover, Kehl et al. (2017) did not incorporate the contribution of spin-down mass loss of pulsar A to the orbital period derivative into the simulations. As we will show below, considering this effect is important, and its impact on our ability to measure the MOI needs to be studied in a fully consistent analysis. Hence, more complete simulations of the MOI measurement in the Double Pulsar should give us a more realistic estimate of the system's (near) future capability to constrain the EOS of ultra-dense matter inside a NS.
Consequently, in what follows, we present new and important details of how to measure the MOI of radio pulsars using the method of isolating the LT contribution to the advance of periastron. Using the Double Pulsar as the most promising system for this kind of experiment, we simulate timing data of PSR J0737−3039A that can be expected from MeerKAT and future extensions, to assess our ability to measure its MOI in the next 10 years.

LENSE-THIRRING EFFECT IN THE DOUBLE PULSAR
The Double Pulsar is the only system to-date where both NSs have been observed as pulsars (Burgay et al. 2003;Lyne et al. 2004), with an orbital period of only 2.4 h. Breton et al. (2008) used the 3 See skatelescope.org for updates. system to provide a 13%-test of spin-orbit interaction of strongly self-gravitating bodies using the relativistic spin precession in pulsar B. The compact, relativistic nature of the system also allows the measurement of several post-Keplerian (PK) parameters to an unparalleled level of accuracy. This not only enables some of the most stringent tests of GR related to strong-field gravity Kramer & Wex 2009;, but it is also crucial for the efforts to measure the MOI and to constrain the EOS of a NS.

Spin-orbit coupling contribution to the periastron advance
To simplify the problem, we neglect the LT contribution of pulsar B, since it spins about 122 times slower than pulsar A. Such a simplification is well justified, as will become clear below. In addition, the long term observations of the pulse profile of PSR J0737−3039A shows that the misalignment angle between the spin axis of pulsar A and the orbital angular momentum has an upper limit of 3.2° (Ferdman et al. 2008(Ferdman et al. , 2013. Therefore, for all practical purposes, we can assume that the spin of pulsar A is parallel to the orbital angular momentum, which is consistent with evolutionary considerations for the Double Pulsar system and a low-kick supernova formation (e.g. Stairs et al. 2006;Tauris et al. 2017). Pol et al. (2018) confirmed that pulsar A is indeed rotating prograde in its orbit, using the modulation of pulsar B's radio emission by the interaction with the wind of pulsar A. Consequently, the spin of pulsar A only induces a change to the advance of periastron, and does not lead to a change in the orbital inclination, more specifically, the projected semi-major axis. Following Damour & Schäfer (1988), the total intrinsic contribution to the periastron advance in the Double Pulsar system can be written, with sufficient precision, as where n b is the orbital frequency, and e T is the proper-time eccentricity used as the observed eccentricity in the standard timing model (Manchester et al. 2015) and defined in Damour & Deruelle (1986). The factor in front of the right-hand side of Eq. (2) is the first post-Newtonian (1PN) contribution; the higher order corrections due to 2PN effects and LT precession caused by pulsar A are indicated by the second and third term in the square brackets respectively. The following notations are used to simplify Eq. (2), (1 − e 2 T ) 1/2 The subscript A stands for pulsar A. G is the Newtonian gravitational constant, and M = m A + m B is the total mass defined as the sum  Kramer et al. (2006). I (45) A = I A /(10 45 g cm 2 ). The current measurement precision for ω is already ∼10 −5 deg yr −1 (Kramer et al. in prep.), which is about 40 times smaller than ω LT,A .  Table 1 lists the values of each term contributing to ω intr , using the Keplerian parameters and masses (m A = 1.3381 M , m B = 1.2489 M ) measured in Kramer et al. (2006). We note that the contribution due to the LT precession ω LT,A depends on the MOI, whereby is written as a function of I are around unity for realistic EOSs. It is evident that the contribution from the LT effect is comparable to that of 2PN, but with opposite signs.
The analysis of timing data from relativistic binary pulsars is based on a particularly simple and elegant solution of the post-Newtonian equations of motion, the so called Damour-Deruelle (DD) model (Damour & Deruelle 1985;Damour & Deruelle 1986;Manchester et al. 2015). In the quasi-Keplerian parametrization of the DD model one can see that the advance of periastron is proportional to the true anomaly. This behaviour is modified by two periodic terms as part of the generalised quasi-Keplerian parametrization, which is a natural extension of the DD model when including 2PN and spin-orbit terms (Damour & Schäfer 1988;Schäfer & Wex 1993;Wex 1995). However, these periodic terms will remain well below measurability for the foreseeable future, for any of the known binary pulsars. For that reason, we will ignore them in our analysis.
Besides the coupling to the orbital angular momentum (spinorbit coupling), the spin of pulsar A also couples to the spin of pulsar B (spin-spin coupling) (Barker & O'Connell 1975). However, the spin of pulsar B is about a factor of 3 × 10 6 smaller than the orbital angular momentum. Hence, spin-spin coupling is totally irrelevant here.
Finally there are, at least in principle, also contributions from the rotationally induced mass quadrupole moments of pulsars A and B to the orbital dynamics (Barker & O'Connell 1975). These spin-squared contributions give rise to an additional change in the advance of periastron (Smarr & Blandford 1976;Wex 1998). The contribution from the quadrupole moment of pulsar A is estimated to be ∼ 3 × 10 −8 deg yr −1 , where we have used the relations in Bauböck et al. (2013) to calculate the mass quadrupole. This is four orders of magnitude smaller than the second order effects. The contribution from pulsar B is even smaller (about 10 4 times) due to its slower rotation. Hence we can totally ignore such contributions in this study.

The proper motion contribution to the observed periastron precession
Apart from the intrinsic contributions to the periastron advance, the proper motion of a binary system also can change the apparent geometrical orientation of the orbit, and hence the observed periastron advance (Kopeikin 1996). As a consequence, the observed value of periastron advance is shifted from its intrinsic value by Here, ω K is the Kopeikin term that satisfies where i is the orbital inclination as defined in Damour & Taylor (1992), µ α and µ δ are the proper motion in right ascension and declination, and Ω is the longitude of the ascending node (measured from East, in the sense of rotation towards North). Using the parameters measured by Kramer et al. (2006) and the estimated Ω = 25(2)°b y Rickett et al. (2014) Given the current measurement precision ∆ ω ∼ 10 −5 deg yr −1 (Kramer et al. in prep.), the Kopeikin term is a small correction to the intrinsic periastron advance that we use in this study. However, since it is three orders of magnitude smaller than ω LT,A (see Table 1), it does not have a significant influence on the LT measurement.

Challenges on extracting the Lense-Thirring contribution and measuring the MOI
Although the current measurement precision ∆ ω is already ∼ 40 times smaller than ω LT,A , it is not that straightforward to extract the LT contribution from ω obs with Eqs.
(2) and (7), as the two masses (m A , m B ) are needed to calculate ω 1PN and ω 2PN . The masses need to be obtained from any other two PK parameters, where the best two here are the Shapiro delay shape parameter s and the orbital period derivative P b (see Figure 2). For the Double Pulsar, we already have sufficient precision for s, so the limitation is mainly from P b (Kramer et al. in prep.). The measurement precision of P b will improve over time, especially with the addition of MeerKAT and the SKA. However, the observed value of P b is influenced by extrinsic acceleration effects, which depend on the distance of the pulsar and the Galactic gravitational potential. Moreover, the spin-down mass loss of the pulsars also have an impact on P b , which itself depends on the MOI, meaning the masses can not be determined independently from I A . These are the challenges for measuring the MOI. An alternative option to P b could be the time dilation amplitude γ, whose fractional error is about one order of magnitude larger than P b (see Figure 2). However, based on the assumptions of observing plan in Section 5, it would take at least two decades from now to obtain a 1σ-measurement of ω LT,A using only γ, s and ω, a precision that is already reached now with P b (Kramer et al. in prep.). Hence, a comprehensive understanding of the individual contributions to P b is needed, which will be discussed in detail in the following section.

THE INTRINSIC AND EXTRINSIC CONTRIBUTIONS TO THE ORBITAL PERIOD DECAY
The observed value of the orbital period decay comprises several effects (Damour & Taylor 1991). For the purpose of this study, we only consider the dominant terms where gravitational wave damping (GR) and mass loss of pulsar A ( m A ) are intrinsic contributions, and Galactic acceleration (Gal) and Shklovskii effect (Shk) are extrinsic contributions. Thereby, the intrinsic orbital period decay can be extracted from the observed value using Consequently, the uncertainty in the intrinsic orbital period decay also depends on the error in the pulsar distance and the uncertainty in the Galactic gravitational potential at the location of the pulsar and the Earth.

Gravitational wave damping
The binary system loses energy in the form of GW emission, which shrinks the orbit of the system, and in turn gradually reduces the orbital period. The post-Newtonian approximation is employed to describe the orbital dynamics of the binary system (see e.g. Damour 1987; Blanchet 2014), i.e. the equations of motion are expanded with respect to v/c, where v denotes a typical orbital velocity. The change of the orbital period due to GW damping enters at order (v/c) 5 , i.e. the 2.5PN approximation. The corresponding change in the orbital period is given by (Peters & Mathews 1963;Esposito & Harrison 1975;Wagoner 1975) where η = m A m B /M 2 is the symmetric mass ratio. Later, Blanchet & Schäfer (1989) extended the expression to the next-to-leading order (3.5PN), where δm denotes the mass difference of the timed pulsar and its companion, in our case, δm = m A − m B . Eq. (12) can be written in a simplified form as where the relative correction of the 3.5PN order, X 3.5PN , is 1.40 × 10 −5 for the Double Pulsar. To date, only the leading order contribution to the orbital period decay is considered in the analysis and interpretation of any of the known binary pulsars. The higher order correction, however, will need to be included in the future, when we reach the necessary timing precision with emerging powerful radio telescopes such as the SKA. We will evaluate future measurability of the 3.5PN contribution to P b in Section 8. Besides the damping of the binary period, the emission of GWs in principle has an additional effect on the observed P b . Junker & Schäfer (1992) have shown that a double NS system with asymmetric masses in an eccentric orbit becomes accelerated due to the GW recoil. Since any acceleration along the line of sight leads to an apparent change in the orbital period (Damour & Taylor 1991), the GW recoil at 3.5PN order will also affect the observed orbital period at some level. As Junker & Schäfer (1992) have pointed out, the recoil acceleration changes its direction with the advance of periastron, in our case on a timescale of about 21 years. However, using Eq. (103) in Junker & Schäfer (1992) we find a maximum shift in the observed P b due to GW recoil of 4.6 × 10 −24 , which is seven orders of magnitude below the current measurement precision.

Galactic acceleration and Shklovskii effect
The contribution of Galactic acceleration can be calculated with (Damour & Taylor 1991;Nice & Taylor 1995;Lazaridis et al. 2009) where  (Cordes & Lazio 2002). We note that new, preliminary timing and VLBI measurements indicate a distance closer to the DM distance (Kramer et al. in prep.). Hence, for our simulation, we consider an intermediate distance of 0.8 kpc with a 10% error. We will see in Section 6, using a different distance does not have a big influence on our results. The vertical contribution of the Galactic acceleration K z for Galactic height z ≡ |d sin b| 1.5 kpc can be approximated with the expression (Holmberg & Flynn 2004;Lazaridis et al. 2009) where z kpc ≡ z[kpc]. For K z , we consider a typical error of about 10% (Holmberg & Flynn 2004;Zhang et al. 2013). The Galactic parameters R 0 is the distance from the Sun to the Galactic center, and Θ 0 is the Galactic circular velocity at the location of the Sun. In our calculation, we adopt the recent result in Gravity Collaboration et al. (2019), where R 0 = 8.178 ± 0.026 kpc 6 and Θ 0 = 236.9 ± 4.2 km s −1 . The slope parameter at the radius of the Sun is defined as (Damour & Taylor 1991 We note this term is often ignored in other studies, as the rotation curve is nearly flat in the vicinity of the Sun. Its uncertainty, however, could be relevant for measuring the MOI, and as such is included in our study. The slope of the Galactic rotation curve at the location of the Sun estimated by Reid et al. (2014) is −0.2 ± 0.4 km s −1 kpc −1 , corresponding to b 0 = 0.007 ± 0.014. Lately, Eilers et al. (2019) found a slope significantly different from zero, i.e. −1.7 ± 0.1 km s −1 kpc −1 (b 0 = 0.0603 ± 0.0035), with a systematic uncertainty of 0.46 km s −1 kpc −1 . Both results will be employed later in our simulation, but for Eilers et al. (2019) we only consider the statistical error and assume that the systematic error can be well understood in the future. Besides the Galactic acceleration, there are additional accelerations due to masses in the vicinity of the Sun or the pulsar, primarily giant molecular clouds (GMCs), but also stars, black holes, and other external masses (Damour & Taylor 1991;Kehl 2015). These have most likely, if at all, only a small influence on the result, which should not limit our ability of measuring the MOI with a precision lower than 10% (Kehl 2015). The influence of these masses mostly depends on the distance to the Double Pulsar, with which we can, for instance, trace and restrict the presence and influence of GMCs using Galactic carbon monoxide (CO) surveys (Neininger et al. 1998;Glover & Mac Low 2011). We expect that a more firmly established distance measurement in the future will allow a refined analysis to confirm our conclusions.
Finally, the transverse motion of a pulsar leads to an apparent change in the orbital period. This is known as the Shklovskii effect (Shklovskii 1970), and is given as

Mass loss
A pulsar looses mass due to its energy emission, which changes the orbital period by (Jeans 1924(Jeans , 1925) Although the emission process of pulsars is not fully understood, the mass-energy loss can be calculated (with sufficient precision) from the loss in rotational kinetic energy, i.e., E rot j m j c 2 (Damour & Taylor 1991), where E rot j = I j Ω j Ω j , with Ω j the angular velocity of the (proper) rotation of body j ( j = A or B), given in terms of the spin period by Ω j = 2π/P j . Hence, Clearly, the mass-loss correction to the rate of orbital period decay also depends on the MOI, and therefore on the EOS. Table 2 lists the predicted value of each contribution to P obs b , where the mass-loss contributions are written as a function of I (45) j . The contribution due to the mass loss of pulsar A is one order of magnitude smaller than that of the Galactic acceleration and the Shklovskii effect, and of the same order of magnitude as the current measurement precision (Kramer et al. in prep.), hence must be considered. The mass-loss contribution of pulsar B, however, is nearly four orders of magnitude smaller than that of pulsar A and thus can be safely ignored.

SIMULATIONS
In order to investigate the capability of measuring the MOI and testing GR with radio pulsars, we developed a simulation framework to generate and analyse time-of-arrivals (TOAs) for binary pulsars. In this section, we will describe how we simulate TOAs from emerging Table 2. Contributions to the rate of orbital period decay in the Double Pulsar, calculated with Keplerian parameters and masses measured in Kramer et al. (2006). The Galactic acceleration is computed using Galactic measurements by Gravity Collaboration et al. (2019)  A . The current measurement precision for P b is below 0.1 fs/s (Kramer et al. in prep.).

Contribution
[fs/s] telescopes for PSR J0737−3039A based on realistic assumptions, and how to measure PK parameters and timing parallax.
To simulate TOAs of PSR J0737−3039A from current and future telescopes, knowledge of the sensitivity of the telescopes, as well as (realistic) assumptions about a future observing plans are needed. We consider the best telescopes for observing this pulsar, i.e., MeerKAT and its future arrays. Unfortunately, this pulsar is not in the field of view of the Five-hundred-meter Aperture Spherical radio Telescope (FAST; Nan et al. 2011), the largest radio telescope today and in the near future.
MeerKAT is a precursor for the mid-frequency array of the SKA, which comprises 64 dishes, each with a diameter of 13.5 m. This corresponds to an effective diameter ( / eff ) of 108 m. Regular timing observations for the Double Pulsar started in 2019 as a part of the MeerTIME project (Bailes et al. 2018). The MeerKAT extension, hereafter MeerKAT+, is a joint collaboration of the South African Radio Astronomy Observatory (SARAO) and the Max-Planck-Society (MPG) to extend MeerKAT by the addition of 20 SKA-type dishes, each 15 m in diameter, to MeerKAT. MeerKAT+ is expected to operate from 2022, providing an increase in sensitivity by 50% (Kramer, priv. comm.) The first phase of the SKA midfrequency array, SKA 1-mid, is planned to build 112 additional dishes with 15 m diameter, extending MeerKAT+ further, with first data from 2025 and full operation after 2027. We summarise the observing plans and the effective diameters of these telescopes in Table 3.
In order to estimate the TOA uncertainty of each observing phase, we need to consider noise contributions for pulsar A. The TOA uncertainty of pulsar A with real MeerKAT observations at L-band is about 1.06 µs for a 5 minutes integration over the full bandwidth (Bailes et al. 2020). Since the system performance of MeerKAT+ and SKA 1-mid are expected to be similar to that of MeerKAT, and the radiometer noise σ rn reduces in reverse proportional to the effective collection area of the telescope A eff , we can therefore calculate the radiometer noise using the relation where the superscript "MK" stands for MeerKAT. We are not considering noise budgets other than the radiometer noise, because: 1) The phase jitter has not been detected in the current MeerKAT observations and must be rather small (Bailes et al. 2020, Hu et al. in prep.). It may become important in the future observing phase as the radiometer noise reduces, but the influence of jitter can potentially be reduced using Bayesian methods (Imgrund et al. 2015) or binning and combining the data in orbital phase.
2) The contributions from scintillation and other effects are expected to be one or more orders of magnitude smaller than the radiometer noise of SKA 1-mid, hence are neglected. As a result, in our simulation, we adopt the TOA uncertainties solely based on the radiometer noise estimation for each observing phase, which can be found in Table 3.
Based on the above assumptions, we generate TOAs of PSR J0737−3039A that mimic observations with MeerKAT, MeerKAT+, and SKA 1-mid from 2019 to 2030 covering two full orbits per month (∼5 h), and combine them with the existing TOAs from multiple telescopes (Kramer et al. in prep.) to form a longrange dataset . Technically speaking, we only use the observing cadence and TOA uncertainties from the existing TOAs, since the data analysis by Kramer et. al (in prep.) is still ongoing, and in the next steps all TOAs will be simulated to fit our model, under the assumption of Gaussian white noise.
The first step is to create a parameter file (model) for pulsar A. For this, we take precisely measured masses from Kramer et al. (2006), m A = 1.3381 M , m B = 1.2489 M , and assume EOS AP4 (see Lattimer & Prakash 2001). This particular choice of EOS satisfies the current lower limit of 1.98 M (99% confidence level, hereafter C.L.) for the maximum mass of a NS (see Section 1 for details), and also lies in the MOI ranges obtained for pulsar A by Gorda (2016); Lim et al. (2019); Greif et al. (2020). The MOI of pulsar A, under this assumption, is therefore I AP4 A = 1.24 × 10 45 g cm 2 . We create a parameter file by taking the well measured Keplerian parameters of the Double Pulsar  and the PK parameters computed from m A , m B , and I A . For the advance of periastron ω, we consider first and second order PN terms and the LT contribution. As for the orbital period decay P b , we consider leading order (2.5PN) GW emission, Galactic acceleration, Shklovskii effect and mass loss in pulsar A. The 3.5PN GW term is only considered in Section 8.
We then adjust the TOAs to perfectly match with our model, and add a Gaussian white noise to each TOA, according to its σ TOA . The red noise from DM variations is not considered in our simulation, since it can be in principle corrected for with multi-frequency data. In a final step, we use the pulsar timing software TEMPO 7 to 7 http://tempo.sourceforge.net/ fit for the timing parameters and obtain their uncertainties, including the PK parameters, which are of particular importance here. From 2018 to 2030, the dataset is split with a step size of 6 months, so as to demonstrate how the measurements improve with time.
The predicted fractional errors of the PK parameters are shown in Figure 2. As part of the simulation, we also measure the timing parallax π x , which gives an idea of the precision of future distance measurement from timing parallax. The predicted uncertainty of π x is shown in Figure 3. For the uncertainty of pulsar distance, which enters the Galactic acceleration and the Shklovskii effect, we adopt the value calculated from timing parallax when its uncertainty is smaller than what we assumed in Section 4.2, which is from mid-2021. Aside from timing parallax measurement, in the future, the VLBI parallax measurements with the SKA can potentially provide an accurate distance measurement (Smits et al. 2011).

MEASURING THE MOI AND CONSTRAINING THE EOS
Based on our TOA simulation, we predict the future timing measurement of PK parameters ( Figure 2). The three best measured parameters, P obs b , ω obs and s, are promising for the determination of I A . With Eqs. (7) and (10), we obtain the intrinsic periastron advance ω intr (m A , m B , I A ) and the intrinsic orbital period decay P intr b (m A , m B , I A ). Since both now ω intr and P intr b depend on the MOI, we can not directly use P intr b and s to determine the masses and hence measure I A from ω intr as in Kehl et al. (2017). Instead, a self-consistent method is employed to solve for the masses (m A , m B ) and I A jointly from P intr b (m A , m B , I A ), s (m A , m B ) and ω intr (m A , m B , I A ). To estimate the probability distribution function for I A , we perform a Monte Carlo simulation to randomise the observed parameters according to their uncertainties. This process is repeated for the measurements from 2018 to 2030. Figure 4 shows the predicted measurements of I A with time, where the new telescopes clearly help to narrow down the uncertainty of I A . Here we adopt the Galactic measurements (R 0 , Θ 0 ) by Gravity Collaboration et al. (2019) and the slope measurement by Reid et al. (2014). The predicted uncertainty of I A with time is also illustrated as the blue line in Figure 5. In this case, we expect to achieve an MOI measurement with 25% precision at 68% C.L. by the year 2030. Our simulation shows that, although the uncertainty of P obs b is initially higher than the Galactic acceleration, it decreases with additional years of precise timing observations (see Figure 2), and by 2030, the error in the Galactic acceleration is three times higher than the error in P obs b , which becomes the limiting factor for measuring the MOI.
However, the measurements of the Galactic potential is expected to improve through various observational methods, such as  Eilers et al. (2019) are similar to the measurements used in the previous case (blue line), but here we assume the systematic errors can be well understood in the near future, and only consider the statistical errors. With this assumption, we expect to measure the MOI with 11% precision at 68% C.L. in 2030. This is nearly the same as using an error-free Galactic model, which is indicated by the red line in Figure 5. Therefore, with future measurements of the Galactic potential and a better understanding of the systematic errors, a MOI measurement with 11% precision from the Double Pulsar seems realistic.
One important factor for the result is the influence of the mass loss in pulsar A, which was neglected in the previous study by Kehl et al. (2017). Without considering this contribution, the uncertainty of I A significantly reduces and reaches 7% by 2030 (see the grey line in Figure 5), in contrast to the red line. In addition, we find that increasing the observing cadence does not significantly improve the precision of MOI measurements.
As mentioned in Section 4.2, different approaches provide very different measurement of the distance of the Double Pulsar, and a compromise distance of 0.8 kpc is thereby employed in our study. To investigate how distance influences the MOI measurement, we consider two extreme cases, d = 0.4 kpc and d = 1.6 kpc, with the same setups as in the d = 0.8 kpc simulations. Using the current Galactic measurements, we find that the uncertainty of the MOI measurement reaches 17% by 2030 when d = 0.4 kpc, and has a much higher uncertainty (43%) when d = 1.6 kpc. However, with negligible error in the Galactic potential, both predict ∼11% measurements by 2030, same as for the case of d = 0.8 kpc. Since an improved Galactic model is expected in the near future, the value we employ for the distance should not have a significant impact on the prediction of the MOI uncertainty.
An 11% precision measurement of the MOI would further improve the constraints of the EOS of NSs (Lattimer & Schutz 2005;Greif et al. 2020). Figure 6 shows the MOIs of a number of EOSs, which are scaled by a factor of M 3/2 in order to reduce the range of the ordinate (cf. Lattimer & Schutz 2005). The 11% measurement predicted from our simulation is illustrated by the red bar centered at the assumed EOS AP4, and located at the precisely measured mass of pulsar A. To compare with the constraints from other methods, we mark the curves in different styles. The observations of the binary neutron-star merger event GW170817 by LIGO/Virgo (Abbott et al. 2018) placed a constraint for the radii of both NSs, 11.9 ± 1.4 km (90% C.L.), which excludes the EOSs in grey dashed curves. Recently, a more stringent constraint combining GW170817 with nuclear theory was obtained by Capano et al. (2020), where they found the radius for a 1.4M NS is 11.0 +0.9 −0.6 km (90% C.L.). This further excludes the EOSs in blue dashed curves. The remaining promising EOSs from this constraint are marked in blue solid curves, which is already very close to our 11% prediction from the MOI measurement in 2030. With more and more binary NS mergers expected to be detected in the coming years, tighter constraints on the EOS are likely to be achieved. Meanwhile, recent NICER observation delivered a joint mass-radius measurement for PSR J0030+0451 from two independent analyses. Riley et al. −0.14 M and 13.02 +1.24 −1.06 km. This is a weak constraint on the EOS, but is expected to improve with more observations in the near future. The upcoming X-ray missions, such as eXTP (Zhang et al. 2016) and ATHENA (Barret et al. 2013), are also promising to improve our understanding of the mass-radius relation for NSs.
Therefore, it is fair to assume that the GWs and X-ray observations will place a more stringent constraint on the EOS within the next 10 years, and if the EOS can be known with sufficient precision, we can in turn use this information as an input to our analysis, test the LT precession and constrain theories of gravity with the Double Pulsar. We will discuss this scenario in detail in the next section.

TESTING LENSE-THIRRING PRECESSION
As discussed in the previous section, the MOI measurement of PSR J0737−3039A is expected to reach 11% accuracy by 2030, whereas GWs and X-ray observations are likely to give a better constraint on the EOS. In this section, we discuss the prospects of testing LT precession and constraining theories of gravity using the Double Pulsar, if the EOS is known.
We again adopt EOS AP4 and this time assume that a precision of 5% could be achieved when calculating the MOI of pulsar A, based on a (hypothesized) future improvement in our understanding of super-dense matter. Given I A as an input to our simulations, only the masses are unknown for the intrinsic orbital period decay P intr b and the Shapiro shape parameter s. With the masses measured from ( P intr b , s) and the given I A , we can directly test the LT contribution to the periastron advance ω LT,A . To discuss the physical meaning of such a test, we use the generic framework for relativistic gravity theories introduced by Damour & Taylor (1992), which is fully conservative and based on a Lagrangian that includes a generic term L SO for spin-orbit interaction. As in Damour & Taylor (1992), we will make no assumption about the (strong-field) coupling function Figure 6. Constraints of EOSs from an 11% measurement of the MOI of PSR J0737−3039A (red). EOS AP4 was assumed in the simulation (curve through red dot). The grey dashed curves indicate EOSs that are disfavoured by the LIGO/Virgo observations of the GW170817 binary neutron-star merger (Abbott et al. 2018). The blue dashed curves are additionally excluded by the refined (combined with nuclear theory) GW170817 analysis by Capano et al. (2020). The following EOSs have been plotted (ascending in their intersection with the left border): WFF1, WFF2, AP4, BSk20, AP3, SLy4, BSk25, MPA1, BSk21, SLy9, BL, BSk22, H4, PAL1, MS2, MS0 (https://compose.obspm.fr). All these EOSs are able to support a NS of 1.98 M , the current lower limit for the maximum mass (see Section 1 for details). Γ B A , which enters L SO . Since the spin axis of pulsar A has been found to be practically parallel to the orbital angular momentum, the general form of the LT contribution to the periastron advance can be written as where σ A is a generic strong-field spin-orbit coupling constant, defined by In GR, the generalised gravitational constant G equals G, and the coupling function Γ B A equals 2G (Damour & Taylor 1992), so that But in other theories, Γ B A is expected to deviate from 2G, including modifications by self-gravity contributions from the strongly selfgravitating masses in the system.
We define a parameter δ LT to measure the relative deviation of the theory-independent parameter σ A /G from its GR prediction, By inserting Eq. (22) into the above definition, one obtains for the spin-orbit coupling function To assess potential constraints on a non-GR spin-orbit coupling, we multiply the expression of ω LT,A in GR (last term in Eq. (2)) by (1 + δ LT ), and solve for the parameter δ LT using the three PK parameters and ω intr (m A , m B , δ LT ). One has to keep in mind that, for simplicity, we make here the assumption that the non-spin related parts of the orbital dynamics and signal propagation are (to sufficient approximation) given by their GR expressions. It goes without saying, that in practice one has to conduct a fully self-consistent analysis within a given class of alternative gravity theories. For a discussion that purely focuses on the measurability of a potential deviation in the LT contribution, our approach is sufficient. Figure 7 shows the expected decrease in the uncertainty of δ LT with future observations. With R 0 and Θ 0 measurements from Gravity Collaboration et al. (2019) and the slope measurement from Reid et al. (2014), we expect to measure δ LT with 18% precision at 68% C.L. by 2030, which is indicated by the blue line. The red line adopts the Galactic measurements from Eilers et al. (2019), where we expect to achieve a 9% precision by 2030. In the ideal case, we assume that the Galactic potential, the distance to the Double Pulsar, and the MOI can be precisely measured in the future, so that we could leave out the errors. In this scenario, we expect to measure δ LT with 7% precision by 2030 (green line). We have seen in Section 6 that change from the Galactic measurements by Eilers et al. (2019) to an error-free Galactic model has little enhancement on the measurements of the MOI, and the uncertainty of the timing parallax is relatively small, therefore, the improvement from 9% (red line) to 7% (green line) is to a fair fraction (nearly half) related to the uncertainty of the MOI. Breton et al. (2008) have conducted a different experiment for spin-orbit coupling in the Double Pulsar system. Studying the geodetic precession of pulsar B, they were able to show that σ B /G is in agreement with GR, with a precision of about 13%. Analogously to Eq. (22), σ B is related to Γ A B . A priori there is no reason to assume that Γ A B and Γ B A are equal (see discussion in Damour & Taylor 1992). Consequently, a LT test with pulsar A would nicely complement the geodetic precession test of Breton et al. (2008), when investigating the relativistic interaction between the proper rotation of the two NSs and their orbital motion.
Finally, short range modifications of gravity, related to the strong gravitational field of a NS, could significantly change the structure of the star and therefore its MOI, without any "direct" impact on the orbital dynamics or the signal propagation in a binary pulsar system. Examples of such theories are scalar-tensor theories with a massive scalar field having a sufficiently short Compton wavelength (see e.g. Ramazanoǧlu & Pretorius 2016;Yazadjiev et al. 2016). While in such a scenario, PK parameters related to time dilation, GW damping, and Shapiro delay remain (practically) unaffected (see e.g. Alsing et al. 2012), one could still expect a deviation in the precession of periastron of the Double Pulsar. The reason is that due to the modification of the MOI the spin of pulsar A and therefore the spin-orbit coupling is modified. Testing the LT precession in the Double Pulsar can therefore be used to constrain such deviations from GR. It is important to note, that P m b would also be modified accordingly, and therefore has to be accounted for. Hence, limits on δ LT would consequently be somewhat weaker than given above (cf. Section 6). In such a scenario it could generally be difficult to disentangle uncertainties in the EOS and deviations from GR by astronomical observations. For this, a combination of various experiments, like GWs from binary neutron-star mergers, X-ray observations, and radio pulsar timing might turn out to be necessary. Nonetheless, the future measurement of the LT precession in the Double Pulsar is expected to provide important contributions when constraining such deviations from GR.

NEXT-TO-LEADING ORDER GRAVITATIONAL WAVE DAMPING
In GR, the loss of energy of a material system due to GWs is to leading order sourced by a time-dependent mass quadrupole (Einstein 1918;Eddington 1922). This also holds for binary systems where a change in the mass quadrupole is driven by gravity itself. It enters the two-body equations of motion at the 2.5PN order (see e.g. Damour 1987). When computing the next-to-leading order contribution to GW damping, one also has to account for the mass-octupole and the current quadrupole moments (Thorne 1980). Next-to-leading order contributions enter the equations of motion at 3.5PN (O(c −7 )), and therefore correspond to the 1PN corrections in the radiation reaction force (Iyer & Will 1995;Pati & Will 2002;Königsdörffer et al. 2003;Nissanke & Blanchet 2005). The corresponding change in the orbital period of a binary system has been determined out by Blanchet & Schäfer (1989) and is given by Eq. (12). In this section we will investigate if next-to-leading order corrections to the GW damping are expected to become important in the near future for the timing observation of the Double Pulsar. Again we assume EOS AP4 and a 5% error in the knowledge of the MOI I A . We implement the 3.5PN contribution into our model by using Eq. (12), and adjust the TOAs accordingly. After running simulations as described in Section 5, we obtain the measured PK parameters. We use Eq. (13) to solve for the relative correction of the 3.5PN order X 3.5PN using the three PK parameters P intr b (m A , m B , X 3.5PN ), s (m A , m B ), and ω intr (m A , m B ). Figure 8 illustrates the predicted uncertainty of X 3.5PN with observing phase, which will fall below its theoretical value X theo which will reach a precision of 85% at 68% C.L. by 2030. Adopting the Galactic measurements (statistical errors) by Eilers et al. (2019), the red line shows that X 3.5PN can be constrained with a precision of 42% by 2030. By contrast, in the ideal case where there are no errors in the Galactic model, the distance and the MOI, X 3.5PN can be constrained with a precision of 33% by 2030, where nearly half of the improvement is contributed from the MOI.

POTENTIAL NEW DISCOVERIES
Large pulsar surveys with MeerKAT, FAST and the forthcoming SKA, such as TRAPUM (Stappers & Kramer 2016) and CRAFTS (Li et al. 2018), can potentially discover more relativistic double neutron star (DNS) systems, preferably with a more compact orbit than PSR J0737−3039. An example of such a system, PSR J1946+2052, with a more relativistic orbit than the Double Pulsar (P b 1.88 h) and larger periastron advance ( ω ≈ 26 deg yr −1 ) and LT precession ( ω LT ≈ 0.001 deg yr −1 ), was recently discovered in the PALFA survey (Stovall et al. 2018). In its orbital parameters, the PSR J1946+2052 system resembles a system similar to the Double Pulsar, but that has evolved further due to GW damping, by about 40 Myr. While it is still unclear, if for PSR J1946+2052 the necessary precision in the mass determination can be reached to rival the Double Pulsar in the tests proposed here 8 , it certainly adds confidence to the hope of finding more relativistic "cousins" of the Double Pulsar in the coming years. Such binary pulsars would quite likely enable MOI measurements with superior precision within a comparably short period of time, and improve the constraints of the EOS.
Here we consider two scenarios, one with an orbital period of 100 minutes and one with 50 minutes, which are within the expected acceleration searches by MeerKAT. Assuming such systems can be found in 2020 and we start timing them regularly from 2021, with two orbits per month, we run our simulation again to predict the measurements of the MOI. To simplify the simulation, we assume these systems satisfy the conditions of the Double Pulsar (inclination i close to 90 degrees, similar distance and brightness, etc.) but with modified orbital parameters, assuming that these systems had an orbit like the Double Pulsar some time in the past, and then evolved by GW damping to an orbital period of 100 or 50 minutes. In reality, these systems are likely to be further away. Nonetheless, it is also possible that such systems are bright and nearby, but were missed in the past surveys due to their high acceleration (see Johnston & Kulkarni 1991;Ransom 2001;Jouteux et al. 2002;Ng et al. 2014;Cameron et al. 2018).
We calculate the evolved semi-major axis using Kepler's third law and the evolved eccentricity using the a − e relation in Peters (1964), for the orbital period of 100 minutes and 50 minutes, respectively. Then we calculate the PK parameters and run simulations as described in Section 5 and 6. Assuming the same distance as the Double Pulsar, we convert the uncertainty of timing parallax into an uncertainty for the distance. The Galactic measurement by Eilers et al. (2019) is adopted in the simulation and, as before, we assume the systematic uncertainties can be well understood in the future.
Our results show that, for the DNS system with an orbital period of 100 minutes, we could measure the MOI with 12% precision by 2030 and with 4.5% by 2035 at 68% C.L. As for an orbital period of 50 minutes, we expect an MOI measurement with 1.5% precision by 2030 and with 0.5% by 2035 at 68% C.L. Such measurements would probably be comparable to the by then available constraints from other methods (GWs and X-ray observations, nuclear physics, etc.) and help for determining the EOS of NSs.
Furthermore, LISA has the potential to discover ultra relativistic DNS systems with a characteristic orbital frequency of 0.8 mHz (Lau et al. 2020). Thrane et al. (2020) suggested that following up such systems with SKA for 10 years could potentially measure the mass-radius relation with a precision <1%. To this end, we perform a simulation for a DNS system with 20 minute orbital period, and find an MOI precision of ∼0.2% (68% C.L.) may be possible with 10 years of timing with SKA 1-mid.
However, there is a low chance that the new discovered DNS systems will be edge-on to our line-of-sight, as is the case for PSR J0737−3039, hence a precise measurement of s might not be possible. Instead, we may need to use γ to constrain the masses and MOI, whose fractional error is usually a few orders of magnitude larger than s (see Figure 2). This is indeed the case for PSR J1946+2052, despite its relativistic nature, determining the masses with sufficient precision will be challenging.
Moreover, not all DNS systems are ideal to test the Lense-Thirring precession in terms of periastron advance ω LT . Systems like the aforementioned PSR J1757−1854 have a large eccentricity most likely caused by a large kick (Tauris et al. 2017) causing a significant misalignment between the spin of pulsar and the orbital angular momentum, and hence ω LT can not be determined as straightforwardly as in the Double Pulsar. However, as pointed out in Section 2.2, this allows an alternate test using the contribution of LT precession to the rate of change of the projected semi-major axis x LT (Cameron et al. 2018) if profile changes due to geodetic precession can be accounted for in the timing process and the spin orientation can be determined with sufficient precision.

CONCLUSION
In this paper, we have developed a consistent method to measure the MOI of radio pulsars, which has been applied to mock data for the Double Pulsar. We simulated TOAs of PSR J0737−3039A assuming future observations with MeerKAT, MeerKAT+ and SKA 1-mid which cover two orbits per month. We found a MOI measurement with 11% accuracy (68% C.L.) could be achievable by the end of this decade, if we have sufficient knowledge of the Galactic gravitational potential (e.g., from Gaia mission (Gaia Collaboration 2016)). We also found that the mass loss of pulsar A has a considerable impact on the measurement of the MOI. Neglecting this contribution to the orbital period change leads to an overoptimistic prediction. This is the main reason why, even with the better timing precision used in this paper as compared to Kramer & Wex (2009), by ∼2030 we would still only reach the same accuracy level as predicted by Kramer & Wex (2009). Additionally, the assumptions made in this paper are more realistic compared to Kehl et al. (2017), with timing precision from MeerKAT observation, as well as updated timeline and size of upcoming telescopes.
In the second part of the paper, Section 7 and 8, we have assumed that a better constraint on the EOS might be achieved with GWs and X-ray observations in the future, so as to investigate the capability of testing LT precession and 3.5PN order contributions to the GW damping. This assumption coincides with Landry et al. (2020) where they found that constraints from GWs and X-ray observations are likely to have larger contributions in constraining the EOS than the MOI measurement of J0737−3039A. Assuming a 5% error in the determination of the MOI, we simulated measurements of the relative deviation of the theory-independent spin-orbit coupling parameter σ A /G from GR's prediction. We found a 9% precision measurement is possible by 2030 with an improved Galactic model, whereas a 7% precision measurement in the ideal caseno errors in the Galactic model, the distance, and the MOI. This test is a complement to Breton et al. (2008), where they found a 13% constraint on σ B /G. This measurement would enable a constraint for the coupling function Γ B A that enters the spin-orbit Lagrangian of the two-body interaction for strongly self-gravitating masses. Such a measurement could be sensitive to short range deviations from GR, which otherwise would not show up in the orbital dynamics of such systems.
We have also studied the measurability of GR's next-to-leading (3.5PN) order GW-damping contribution. We predicted that the uncertainty of the 3.5PN order correction X 3.5PN will fall below its theoretical value at the beginning of SKA 1-mid (∼2026) and a measurement of X 3.5PN with 3σ-significance is possible in ∼10 years, if by then we have sufficient knowledge of the Galactic gravitational potential, pulsar distance, and the EOS. This means that from the SKA 1-mid era, we will have to include the 3.5PN term in our analysis in order to avoid any bias. Binary mergers detected by LIGO/Virgo do provide constraints on post-Newtonian (PN) terms (Abbott et al. 2016). Their way of counting the PN contributions is relative to the Einstein quadrupole formula, i.e. the order they enter the radiation reaction force (Blanchet 2014). Their 1PN term therefore contains 3.5PN contributions from the equations of motion. As a comparison to our 3.5PN 3-σ result, (Abbott et al. 2019) provide a ∼10% measurement (90% C.L.) of the (radiative) 1PN coefficient with GW170817. Future merger events will most likely lead to even more precise measurements of this term. While at the 2.5PN (0PN radiative) level, the Double Pulsar is still many orders of magnitude more precise than LIGO/Virgo mergers (Kramer 2016, Kramer et al. in prep.). When it comes to higher order PN contributions, we conclude that binary pulsars are not expected to be competitive, simply because of the much smaller orbital velocity.
Finally, we discussed potential new discoveries of DNS systems with radio telescopes like MeerKAT, FAST, and SKA, as well as the space-based future GW observatory LISA. We demonstrated that for a DNS system which mimics the evolved PSR J0737−3039 with an orbital period of 50 minutes, the MOI measurement is expected to reach 1.5% precision (68% C.L.) after 10 years observation with MeerKAT, MeerKAT+ and SKA 1-mid, and 0.5% precision after 15 years. Moreover, LISA is expected to find DNS systems with a characteristic orbital period of 20 minutes in the near future (Lau et al. 2020). Such discoveries can significantly tighten the constraints for the EOS.
To conclude, although the EOS constraints resulting from a future MOI measurement with the Double Pulsar are not likely to exceed those with LIGO/Virgo mergers and X-ray observations in the coming years, we still anticipating other aspects of science coming from this unique gravity laboratory in future studies based on an improved understanding of the NS EOS as an input. Furthermore, the discovery of more relativistic binary pulsars, possible with the unprecedented surveying capabilities of new and upcoming radio telescopes and advances in data analysis (e.g. Lentati et al. 2018), could ultimately lead to EOS constraints quite competitive with other methods.