Prospects of constraining the dense matter equation of state from the timing analysis of pulsars in double neutron star binaries: the cases of PSR J0737-3039A and PSR J1757-1854

The Lense-Thirring effect from spinning neutron stars in double neutron star binaries contribute to the periastron advance of the orbit. This extra term involves the moment of inertia of the neutron stars. Moment of inertia, on the other hand, depends on the mass and spin of the neutron star as well as the equation of state of the matter. If at least one member of the double neutron star binary (better the faster one) is a radio pulsar, then accurate timing analysis might lead to the estimation of the contribution of the Lense-Thirring effect to the periastron advance, which will lead to the measurement of the moment of inertia of the pulsar. Combination of the knowledge on the values of the moment of inertia, the mass, and the spin of the pulsar, will give a new constraint on the equation of state. Pulsars in double neutron star binaries are the best for this purpose as short orbits and moderately high eccentricities make the Lense-Thirring effect substantial, whereas tidal effects are negligible (unlike pulsars with main sequence or white-dwarf binaries). The most promising pulsars are PSR J0737-3039A and PSR J1757-1854. The spin-precession of pulsars due to the misalignment between the spin and the orbital angular momentum vectors affect the contribution of the Lense-Thirring effect to the periastron advance. This effect has been explored for both PSR J0737-3039A and PSR J1757-1854, and as the misalignment angles for both of these pulsars are small, the variation in the Lense-Thirring term is not much. However, to extract the Lense-Thirring effect from the observed rate of the periastron advance, more accurate timing solutions including precise proper motion and distance measurements are essential.


Introduction
Timing analysis of binary pulsars lead to the measurement of pulsar's spin, Keplerian orbital parameters (the orbital period P b , the orbital eccentricity e, the longitude of periastron ω, the projected semi-major axis of the orbit x p = a p sin i where a p is the semi-major axis of the pulsar orbit and i is the angle between the orbit and the sky-plane, and the epoch of the periastron passage 1 ), as well as post-Keplerian (PK) parameters like the Einstein parameter (γ), Shapiro range (r) and shape (s) parameters, the rate of the periastron advanceω, the rate of change of the orbital periodṖ b , the relativistic deformation of the orbit δ θ , etc. [4]. Sometimes, the Shapiro delay is parametrized 1 The remaining Keplerian parameter, the longitude of the ascending node ϕ does not come in the standard pulsar timing algorithm, it can be measured via proper motion only in very special cases, see [1][2][3] for details.
Submitted to Universe, pages 1 - 15 www.mdpi.com/journal/universe differently, with parameters h 3 and ς [5], or with r and z s [6], and these parameters can easily be expressed in terms of conventional parameters s and r. Measurement of PK parameters leads estimation of masses of the pulsars and their companions, as well as tests of various theories of gravity [7,8]. Note that, in principle, measurements of only two PK parameters are enough to extract two unknowns, i.e., the masses of the pulsar and the companion while measurements of more than two PK parameters lead tests of gravity theory through consistency. However, the uncertainty in measurement is not equal for every PK parameter and usually the two most accurate parameters are used to obtain the best mass estimates (e.g., Fig 1 of [9]). Additionally, relativistic binary pulsars have the potential to constrain the Equation of State (EoS) of matter at extreme densities. Measurements of masses of two pulsars being around 2 M ⊙ have already ruled out many soft EsoS [10,11]. On the other hand, the recent analysis of the gravitational wave event GW170817 from the merger of two neutron stars ruled out some extremely stiff EsoS [12,13]. Still a large number of EsoS are allowed, including several hybrid [14] and strange quark matter [15] EsoS in addition to standard hadronic ones. So, further progress in this issue is essential, and that is expected in the near future as in one hand more neutron star−neutron star mergers are expected to be detected in future runs of the advanced LIGO, and on the other hand, many binary pulsars are being timed regularly and accurately. Besides, upcoming radio telescopes like MeerKAT and SKA will lead significant improvement in pulsar timing.
Among all of the PK parameters, in the present article, I concentrate mainly on the periastron advance, which has a potential to constrain the dense matter equation of state due to the Lense-Thirring effect. Neutron star−neutron star binaries or 'double neutron star' (DNS) systems are the best for this purpose as these systems are compact enough to display effects of strong field gravity yet wide enough to have negligible tidal effects. Moreover, due to high compactness of neutron stars, spin induced quadrupole moments are also negligible. On the other hand, as these neutron stars are rapidly spinning (spin periods of most of the pulsars in DNSs are less than 100 milliseconds), the Lense-Thirring effect is significant. That is why I concentrate on DNSs in this article, although the mathematical formulations are valid for any kind of general relativistic binaries with a possibility of additional effects depending on the nature of the objects in the binaries. There are about sixteen DNSs known at the present time, including the first discovered binary pulsar -the Hulse-Taylor binary, PSR B1913+16. For one DNS, PSR J0737−3039A/B, both members are radio pulsars and the system is known as the double pulsar 2 .
For a long period of time, the double pulsar was the most relativistic binary. However, the recently discovered DNS, PSR J1757−1854 shows even stronger effects of general relativity in some aspects, i.e.,Ṗ and γ [17] due to its high eccentricity combined with a small orbital period. In fact, its eccentricity is 6.9 times and orbital period is 1.8 times larger than those of PSR J0737−3039A/B. Noteworthily,ω is larger for PSR J0737−3039A/B than PSR J1757−1854 where the contribution of eccentricity is less dominant, see Lorimer & Kramer [4] for expressions for PK parameters.
PSR J1906+0746 is another DNS having the second smallest value of P b , just after PSR J0737−3039A/B. However, its P b is slightly (1.1 times) and e is significantly (7.1 times) smaller than those of PSR J1757−1854. That is why although PSR J1906+0746 has the third largest value ofω, other PK parameters are small, even smaller than most of the other DNSs. Finally, the latest discovered DNS, PSR J1946+2052, has broken all records by having the smallest value of P b and the largest value ofω [18].

Precession in double neutron star binaries
Precessions (of both the spin and the orbit) of neutron stars in DNSs are very important. Following Barker & O'Connell [19], the rate of the change of the unit spin vector (s a ) of a spinning neutron star (a) in a binary can be written as: where the angular spin-precession frequency − → Ω sa can be written as: (2) − → Ω sPN a is the angular spin-precession frequency of the neutron star due to the space-time curvature around its companion and can be calculated within the 'post-Newtonian' (PN) formalism.
− → Ω sLT a is the angular spin-precession frequency of that neutron star due to the Lense-Thirring effect of its spinning companion. We can further write: where k = − → L /| − → L | is the unit vector along the orbital angular momentum, − → L is the orbital angular momentum, and Here if a is the neutron star under consideration, a + 1 is its companion. The amplitudes in Eqns. 3 and 4 are given by: where G is the gravitational constant, c is the speed of light in vacuum, and n = 2π/P b is the angular orbital frequency. M a is the mass, I a is the moment of inertia, P sa is the spin period of the a-th neutron star, and For pulsars in DNSs, the Lense-Thirring term ( − → Ω sLT a ) can be ignored. Even for the case of the slow pulsar (B) of the double pulsar, where the Lense-Thirring effect due to the spin of the fast pulsar (A) contributes to − → Ω s B , I find A PN B = 5.07481 deg yr −1 and A LT B = 3.48865 × 10 −5 deg yr −1 (using the values of the parameters given in Table 1). This leads to Ω sLT B ∼ −6.97731 × 10 −5 deg yr −1 using the fact that s A is almost parallel to k [20,21], and hence Ω s B = 5.07474 deg yr −1 , which is close to the observed median value of Ω s B = 4.77 deg yr −1 [22]. It is unlikely that companions of other pulsars would be much faster than PSR A, and hence − → Ω sLT a can be neglected. So, Eqn. 1 becomeṡ where χ a is the angle between k and s a and u is a unit vector perpendicular to the plane containing k and s a in the direction given by the right-hand rule. So, the measurement of the spin-precession of a pulsar helps estimate χ a when other parameters involved in Eqn. 5 are already known.
Similarly, neglecting the tidal and the spin-quadrupole effects, as well as the spin-spin interaction, the orbital angular precession frequency can be written as: where − → Ω bPN and Ω bLT are contributions from the space-time curvature and the Lense-Thirring effect respectively. Note that, both members of the binary contribute to each term.
− → Ω b leads to the precession of both the Laplace-Runge-Lenz vector A as well as L, but none is directly observable. Damour & Schäfer [23] first studied the manifestation of − → Ω b in terms of observable parameters by decomposing − → Ω b as: where h a is the unit vector along the line-of-sight, i.e., from the earth to the a-th neutron star (the pulsar), Υ a = h a ×k |h a ×k| is the unit vector along the line of the ascending node, ϕ a is the longitude of the ascending node of the a-th neutron star and i is the angle between the orbit and the sky-plane. Damour & Schäfer [23] has also given: As already mentioned,ω a is the parameter of interest as it has the potential to put constraints on the dense matter EoS. So, I explore properties ofω a in the next and subsequent sections.

Periastron advance
The observed rate of the periastron advance (ω a,obs ) of the a-th member of a relativistic binary can be written as:ω a,obs =ω a +ω Kop a =ω PN a +ω LT a +ω LT a+1 +ω Kop a (12) whereω PN a is due to the space-time curvature caused by both members of the binary,ω LT a is due to the Lense-Thirring effect of the star a,ω LT a+1 is due to the Lense-Thirring effect of the star a + 1. These three terms together comes into the expressions ofω a (Eqn. 11b) caused by the precession of the orbit.ω Kop a is a secular variation due to the gradual change of the apparent orientation of the orbit with respect to the line-of-sight due to the proper motion of the barycenter of the binary [24]. The expressions for different terms are: where only the first and the second order terms are retained in the expression ofω PN a . µ α and µ δ are the proper motion of the barycenter of the binary (which is measured as the proper motion of the visible object) in the right ascension and the declination respectively, both expressed in the units of milliarcseconds per year. All other parameters are in the SI units. The parameters introduced in Eqns. 13a,b,c are as follow: where X a = M a /M, M = M a + M a+1 is the total mass of the system. β sa is defined in Eqn. 7, and where λ a is the angle between h a and s a (see Fig. 1). It is obvious from Eqn. 16 that the maximum value of g sa occurs when s a is parallel to the vector (3 sin 2 i − 1) k + cos i h a , giving If s a k, i.e., χ a = 0 and λ a = i: Using the expression of g sa, in Eqn. 13b, I getω LT a , , and using the expression of g sa, max , I getω LT a , max . The negative sign in the Lense-Thirring term (Eqn. 13b) implies the fact that actuallẏ ω LT a , max is the minimum value. However, depending on the values of i, χ p , λ p ; g sa (Eqn 16) can be negative makingω LT a (Eqn 13b) positive. Note that, the Lense-Thirring effects from both the pulsar and the companion come in the total periastron advance rate (Eqns. 12, 13b). But if the companion is much slower than the pulsar, as in the case of the double pulsar, then the contribution from the companion can be neglected. Due to the non-detection of any pulsation, we do not know the values of the spin periods of the companions for other DNSs, and cannot rule out the possibility of significant Lense-Thirring effect from those companions. However, most of the pulsars in DNSs are recycled (at least mildly), suggesting the fact that pulsars were born as neutron stars earlier than their companions. The companions, i.e., the second born neutron stars in DNSs, are expected to be slow; because even if neutron stars are born with spin periods in the range of a few tens of milliseconds up to a few hundred milliseconds, they quickly spin down to periods of a few seconds unless they gain angular momentum via mass accretion from their Roche Lobe filling giant companions and become fast rotators. Such spin-up is possible only for the first born neutron star in a DNS. That is why in the present article, I ignore the Lense-Thirring effect from the companions of pulsars in DNSs.
In such a case, i.e. whenω LT a+1 is negligible,ω LT a can be extracted by subtractingω PN a +ω Kop a fromω a,obs (Eqn. 12). One can estimate the value of I a fromω LT a if all other relevant parameters like P b , e, M a , M a+1 , sin i, χ a , λ a are known (Eqns 7, 14, 15, and 16). Then, from the knowledge of I a , M a , and P sa , a new constraint on the EoS can be placed as it is a well-known fact that the moment of inertia of a neutron star depends on its mass, spin-period and the EoS [25, and references therein]. For this purpose, DNSs with large values ofω LT a (larger than the measurement uncertainty inω a,obs ) are preferable. Large values of the orbital eccentricity e, and small values of P b and P s makeω LT a larger.
The procedure is actually a bit trickier than it seems first. Usually, timing analysis of a binary pulsar reports the most accurate measurement ofω a,obs out of all PK parameters, and the masses of the pulsar and the companion are estimated by equatingω a,obs withω PN a and using the second most accurately measured PK parameter. Ifω LT a is large, then this procedure would become erroneous. One should instead estimate masses using two PK parameters other thanω a,obs , and these mass values should also be accurate enough to give a precise estimate ofω PN a . Accurate measurements of proper motion is also needed to evaluateω Kop a . Only then,ω PN a +ω Kop a can be subtracted from the observedω a,obs to get the value ofω LT a . In short, one needs very precise measurements of at least three PK parameters,ω a,obs and two more, as well as the distance and the proper motion of the pulsar (the a-th neutron star).
Because of its small values of P b and P s , PSR J0737−3039A is expected to have a large value oḟ ω LT a in spite of its low eccentricity. Similarly, as PSR J1757−1854 has a slightly larger P b , slightly smaller P s , and much larger e, it is also expected have a large value ofω LT a . For this reason, in the next section, I exploreω LT a for these two DNSs in detail. Note that, although, PSR J1946+2052 might be a useful system for this purpose, it is not possible presently to extend my study for this system due to the lack of knowledge of masses of the pulsar and the companion.
One should also remember the dependence ofω LT a (via g sa ) on the orientation of s a with respect to k and h a (Eqn 16). If I assume that the a-th neutron star is the pulsar, then all 'a's in the subscripts of above equations will be replaced by 'p's, e.g., s p , h p , etc. Most of the time, the value of g sp, is estimated in the literature. However, at least for some pulsars χ p = 0. In such cases, the values of χ p and λ p can be measured by analyzing the change of the pulse profile shape and the polarization. Conventionally, χ p is denoted by δ and λ p is denoted by ζ. Note that, even if χ p is a constant over time, λ p would vary due to the spin precession and can be written as [26]: where Ω sp is the angular spin precession frequency (amplitude of − → Ω sp in Eqn. 2), t is the epoch of the observation, T p0 is the epoch of the precession phase being zero. Note that, i itself might change with time as given in Eqn. (11c).
In Table 1, with other relevant parameters, I compile values of χ p wherever available. I do not compile values of λ p although available in some cases as it is a time-dependent parameter. From Fig.  1, it is clear that the spin would precess in such a way that it would lie on the surface of a fiducial cone with the vertex at the pulsar and having the half-opening angle as χ p . The constraint on λ p is: (i) if χ p < i, then λ p,max = i + χ p and λ p,min = i − χ p (Fig. 1a), (ii) if χ p > i, then λ p,max = χ p + i and λ p,min = χ p − i (Fig. 1b).

Results
For the purpose of demonstration, I compute the value of the moment of inertia (I p ) for the two most interesting pulsars, PSR J0737−3039A and PSR J1757−1854 using the APR equation of state [27] and the RNS code 3 [28]. I find that, both of these pulsars have I p = 1.26 × 10 45 gm cm 2 . I use this value in my calculations, remembering the fact that the true values of I p would be different depending on the true EoS and we are actually seeking an answer whether it would be possible to know the value of I p by singling outω LT p from the total observedω p,obs .
I tabulate values ofω PN p ,ω LT p , , andω LT p , max for all pulsars in DNSs in Table 2. I exclude candidate DNSs, i.e., PSR B1820−11 and PSR J1753−2240 from this calculation due to poor mass constraints. For other cases, I use the limiting mass values where actual values are unavailable. Onlẏ ω LT p , can be calculated for the systems with undetermined/unpublished values of sin i. As expected, PSR J0737−3039A shows the largest Lense-Thirring effect, followed by PSR J1757−1854. In fact, |ω LT p , max |/|ω PN p | is slightly larger for PSR J1757−1854 (2.92 × 10 −5 ) than that for PSR J0737−3039A (2.81 × 10 −5 ). This fact makes it a very interesting system, and I investigate this system in detail. The discovery paper [17] also mentioned large Lense-Thirring effect causing significant amount of di/dt (Eqn. 11c). Orientation of different vectors relevant for the estimation ofω LTa where the subscript a refers to the object being observed, i.e., the pulsar (so the companion would be represented by the subscript a + 1). The vectors are as follow: k is the unit vector along the orbital angular momentum, h a is the unit vector along the line-of-sight, and s a is the unit spin vector. The angles in the figures are: i is the inclination angle between the orbit of the pulsar and the sky-plane as well as the angle between k and h a , χ a is the angle between k and s a , and λ a is the angle between h a and s a . The left and right panels are for χ a < i and χ a > i, respectively. In both of the panels, two positions of s a are shown corresponding to maximum and minimum values of λ a (see text). The pulsar, i.e., the object a is located at the vertex of the fudicial cone (gray in color) made by the precessing s a around k.
I find that, for PSR J1757−1854, |ω LT p , max | = 3.03170 × 10 −4 deg yr −1 is achieved for χ p = 3 • , λ p = 99 • . But as the value of χ p for this system is not yet known (although Cameron et al. [17] obtained a constraint χ p ∼ 25 • based on their simulation of the kick velocity), I vary χ p in the range of 0 • − 60 • . It is obvious thatω LT p depends more on χ p than on λ p as the absolute value of the factor (see Eqn. 16) with cos χ p , i.e., 3 sin 2 i − 1 is larger than that of the factor with cos λ p , i.e., cos i. Also, the measured value of sin i gives two values of i (as sin θ = sin(π − θ)), I use the larger one to get a bigger range of λ p . The minimum value of |ω LT p | = 1.37597 × 10 −4 deg yr −1 is achieved for χ p = 60 • , λ p = 36 • . If I fix χ p strictly as 25 • , then the maximum and minimum values of |ω LT p | are 2.81108 × 10 −4 deg yr −1 and 2.67663 × 10 −4 deg yr −1 , for λ p = 121 • and λ p = 71 • , respectively. In Fig. 2 As Table 2 shows, |ω LT p , max | for PSR J0737−3039A is around 1.4 times smaller than the presently published accuracy,ω p,obs = 16.89947 ± 0.00068 deg yr −1 where the uncertainty is twice the formal 1σ value obtained in the timing solution [9]. The uncertainties in other PK parameters are even poorer (Table 1 of [9]), so if we excludeω p,obs , the mass estimates would be less accurate and should not be used to calculateω PN p . So, to achieve the goal of estimatingω LT p , lowering only the uncertainty iṅ ω p,obs would not be enough, one needs to improve the accuracy of at least two other parameters that can be used to get masses at least as precise as the ones already published. It is impossible to improve the accuracy of the Keplerian parameter R used by Kramer et al. [9] in combination withω p,obs to report the values of the masses, because R = M p /M c = x c /x p involves x c , the projected semi-major axis of the companion (PSR B), which is not visible presently and even when it was visible, it was not a good timer as PSR A. Solving the expressions for PK parameters, I find that the combination of s anḋ P b is the best choice left whenω p,obs and R are excluded, i.e., gives the narrowest ranges of masses M p = 1.34 ± 0.02 M ⊙ , M c = 1.251 ± 0.007 M ⊙ using the published values of the uncertainties. I also find the fact that, to get masses as accurate as the ones reported by Kramer et al. [9], the uncertainty inṖ b should be reduced at least to 8.0 × 10 −16 assuming that s and relevant Keplerian parameters (P b , e, x p ) would also be improved by at least an order of magnitude. The published uncertainty iṅ P b is 1.7 × 10 −14 . However, the question arises whether for such an improved measurement ofṖ b , the dynamical contribution coming from the proper motion of the pulsar (Shklovskii effect) and the relative acceleration between the pulsar and the solar system barycenter, could be ignored. Using the proper motion reported by Kramer et al. [9] and the distance (1.1 kpc) estimated by Deller et al. [29] using VLBI parallax measurement, I find thatṖ b, dyn = 4.91 × 10 −17 where I have used a realistic model of the Galactic potential [30]. The improvement required in the timing solution of PSR J0737−3039A is expected by 2030 [31,32] accompanied by an improvement inω p,obs . Fortunately, ω Kop p is smaller than |ω LT p , | -using the proper motion reported in Kramer et al. [9], I find that it can be at most 1.18 × 10 −6 deg yr −1 (I could not calculate the actual value, as ϕ p is unknown). However, one will still need better estimates of χ p and λ p to be able to extract the value of the moment of inertia of the pulsar from theω LT p .
On the other hand, for PSR J1757−1854, both |ω LT p , max | and |ω LT p , | are about 1.5 times larger than the presently published accuracy,ω p,obs = 10.3651 ± 0.0002 deg yr −1 . However, even this system cannot be used at present to estimateω LT p , not only because of the unknown values of χ p and λ p , but also because of the less accurate values of other PK parameters. If one uses γ andṖ b , then the uncertainties in these parameters should be improved by one and two orders of magnitudes respectively to get the uncertainty in masses of about 0.0008 M ⊙ . This seems plausible, as the present accuracy inṖ b is only 0.2 × 10 −12 [17]. However, due to the lack of measurements of the proper motion and the parallax distance, I could not calculateṖ b, dyn for this pulsar, only can give rough estimates based on the distance guessed from the dispersion measure using two different models of the Galactic electron density, e.g., NE2001 model [33,34] and YMW16 model [35]. The contribution due to the relative acceleration is −2.11 × 10 −15 for the NE2001 distance (7.4 kpc), and −1.65 × 10 −14 for the YMW16 distance (19.6 kpc). Both of the values seem to be too large to be canceled out by the Shlovskii term. Measurements of the proper motion and the distance (using parallax) are essential to estimateṖ b, dyn . The measurement of the proper motion will also help estimateω Kop p .

Summary and Conclusions
Measurement of the moment of inertia of a pulsar will provide yet another constraint on the dense matter equation of state. This issue has been discussed in the literature in the past [32,51,52, and references therein], but unfortunately the effect of the misalignment between the spin and the orbital angular momentum vectors has not been explored in detail mainly because these two vectors are almost parallel for PSR J0737−3039, the only suitable system existed until recently.
However, even a small misalignment angle would have significant consequences, because for a particular neutron star, the theoretical values of the moment of inertia using different EsoS are sometimes very close (depending on the stiffness of the EoS). One can see   whereas for PSR J1757−1854ω LT p lies in the range of −1.376 × 10 −4 to −3.032 × 10 −4 deg yr −1 . Future measurement of χ p for PSR J1757−1854 will narrow down the range ofω LT p . For both of the systems, significant improvements in timing solution are needed. Moreover, the proper motion and the parallax measurements for PSR J1757−1854 is crucial. Additionally, it is essential to measure χ p and λ p for PSR J1757−1854, hence long-term studies of polarization properties and profile variation would be very useful.
Note that, although all of the calculations presented in this article is within the framework devised by Barker & O'Connell [19] and Damour & Schäfer [23], alternative formalisms to incorporate the Lense-Thirring effect inω p,obs exist. One example is Iorio [54], but for systems with i ∼ 90 • , χ p ∼ 0 • these two formalisms take similar forms. Precisely, Iorio [54] reportsω LT p , = −3.7 × 10 −4 deg yr −1 for PSR J0737−3039A using I p = 1.00 × 10 45 gm cm 2 , which is close to what I get, i.e. ω LT p , = −3.766 × 10 −4 deg yr −1 if I use the same value of I p instead of 1.26 × 10 45 gm cm 2 used in this article. A fundamentally different idea, i.e., to incorporate the Lense-Thirring effect into the delays in the pulse arrival times has been proposed recently [55]. If this can be implemented in future algorithms for pulsar timing, the effect of the Lense-Thirring effect could be measured directly.
Finally, if the measurement accuracy ofω p,obs improves by a few orders of magnitudes, one will need to subtract the Kopeikin term (the secular variation due to the proper motion). The measurement of the proper motion of PSR J1757−1854 will help us calculate this term. Moreover, in order to useṖ b to estimate masses accurately, one needs to subtract the contributions from the proper motion and the relative acceleration between the pulsar and the solar system barycenter. A package to perform these tasks based on a realistic potential of the Milky Way has been recently developed [30], and gradual improvement in the model of the Galactic potential is expected with the help of Gaia data.

Acknowledgement
The author thanks the organizers of CSQCD VI (especially Prof. David Blashcke) for an excellent hospitality, all three reviewers for extensive suggestions on the earlier version of the manuscript, and ANI Technologies Pvt. Ltd. and Uber Technologies Inc. for helping survive in Chennai.