© The Authors 2025

1. Introduction

The discovery of the first binary pulsar, PSR B1913+16 (Hulse & Taylor 1975), opened the era of gravitational wave (GW) astronomy. Importantly, the detailed timing of the pulsar enabled the detection of three relativistic effects in its orbital motion: (1) an increase in the longitude of periastron with time (parametrised by its rate, $\dot{\omega }$), (2) the combination of the varying relativistic time dilation and gravitational redshift with orbital phase (parametrised by a physical quantity known as the amplitude of the Einstein delay, γ), and (3) a decrease of the orbital period with time (again, parametrised by its rate, Ṗ_b, Taylor & Weisberg 1982, 1989). These ‘Post-Keplerian’ (PK) parameters quantify relativistic effects in the orbital motion and in the propagation of light for all fully conservative, boost-invariant gravity theories (Damour & Taylor 1992). In general relativity (GR), they depend only on the two masses of the components of the system, at least to the leading order. Thus, if we are sure the measured orbital effects are relativistic¹ and assume GR as the correct theory of gravity, then with two measured PK parameters we can determine the two masses in the binary. Measuring a third PK effect, the Ṗ_b, provided a test of the self-consistency of the theory, and indeed the observed Ṗ_b matched the GR prediction for the orbital decay of the system caused by GW damping. This was the first test of the radiative properties of gravity and the first test of the gravitational properties of strongly self-gravitating objects. Notably, these two aspects of gravity cannot be studied in the laboratory of the Solar System (for a review, see Freire & Wex 2024).

The observation of the orbital decay in this system represented the first evidence of GWs. This was especially important, as the reality of GWs had not yet been fully established theoretically. Its many astrophysical implications (e.g. the inevitability of double neutron star mergers) paved the path for the development of ground-based GW detectors.

In the following decades, a few hundred binary radio pulsars were discovered, but of these, only 24 have been confirmed as double neutron star systems (DNSs)². Of special importance for the discussion in this paper is the Double Pulsar system PSR J0737−3039A/B. With an orbital period of 2h 27 m, it was at the time of its discovery the DNS with the shortest coalescence time known, about 86 Myr, a fact that significantly increased the expected rate of such mergers (Burgay et al. 2003). Apart from that, this is the only DNS known where both neutron stars are radio pulsars (Lyne et al. 2004). Furthermore, its orbital inclination is the highest known for any binary pulsar, yielding extremely precise measurements of the Shapiro delay.

Continued radio timing of this system has yielded not only the most precise neutron star masses but a total of several independent tests of GR, all of which the theory passes (Table V in Kramer et al. 2021). These include the most precise test of the quadrupole formula for GW damping, which is 25 times more precise than the second best system (the aforementioned test with the Hulse-Taylor pulsar, Weisberg & Huang 2016); the test of the propagation of light in a spacetime with a curvature that is 10³ times larger than near Sagittarius A*; and the detection, for the first time, of next-to-leading order effects (see the improvement of the test in the signal propagation with the MeerKat data in Hu et al. 2022). Also for the first time, the system has also allowed the derivation of a useful upper limit of the moment of inertia (I < 3 × 10⁴⁵ g cm²; 90% C. L.) of a radio pulsar, PSR J0737−3039A. The results from this system are directly relevant to this work and are discussed in more detail below.

The focus of this work, PSR J1946+2052, was discovered in 2017 by the Arecibo L-Band Feed Array pulsar (PALFA) survey (Cordes et al. 2006; Lazarus et al. 2015) and reported by Stovall et al. (2018), henceforth Paper 2018. Of all known Galactic DNSs, PSR J1946+2052 has the shortest orbital period, 1h 53 m, and the shortest coalescence time, 46 Myr. The pulsar’s spin period is the smallest for any known member of a DNS, ∼17 ms. The known spin-down means that the pulsar has the largest spin parameter at merger among all Galactic DNSs ( Paper 2018). Its orbital eccentricity (e  =  0.063848(9)) is the smallest for any Galactic DNS, even lower than that of PSR J1325−6253 (e = 0.0640091(7), Sengar et al. 2022). Interestingly, the orbital period and eccentricity of PSR J1946+2052 are very similar to those that the Double Pulsar system will have in about 40 Myr, when it will be almost halfway through its merger. Given the characteristic age of PSR J1946+2052, the upper limits in the orbital eccentricity and period at birth are e < 0.14 and P_b < 0.17 days (Stovall et al. 2018), which makes it very likely that they have decreased significantly since the formation of the system, owing to GW emission.

The published results on the PSR J1946+2052 system were based on a very small timing baseline of 71 days. Despite that short baseline, a remarkably accurate measurement of the spin period and spin period derivative could be established, thanks to the precise position obtained with the Very Large Array. The spin parameters imply a characteristic age of about 290 Myr and a surface magnetic field of about 4  ×  10⁹ G. Also remarkable was by far the largest periastron advance for any known DNS system in the Galaxy, $\dot{\omega }=25.6\pm 0.3deg{\mathrm{yr}}^{-1}$. This implies (assuming the validity of GR) a total system mass of 2.50 ± 0.04 M_⊙, making this one of the least massive DNSs known (see also Martinez et al. 2017).

Since the publication of Paper 2018, we carried out regular radio timing observations of PSR J1946+2052 by using the Arecibo 305-m radio telescope until the end of its operational time in 2020. Since 2019, we have also been observing the system with the Five-hundred-meter Aperture Spherical radio Telescope (FAST; Jiang et al. 2020). Based on the FAST data, Meng et al. (2024) reported the significant profile evolution of PSR J1946+2052, which indicates that its spin axis is misaligned with the orbital angular momentum and is undergoing relativistic spin precession. The first report of the relativistic spin precession was in PSR B1913+16, where Weisberg et al. (1989) and Kramer (1998) observed a secular variation of the relative amplitude of the two prominent components of PSR B1913+16’s profile. Such a pulse evolution will introduce time offsets when using one standard template to generate times-of-arrival (ToAs), which should be compensated correctly.

This paper discusses some of the results of these radio observations, especially the timing. Its structure is as follows: In Sect. 2, we describe the observations themselves, the data processing and ToA generation and analysis. In Sect. 3, we present the new timing results, in particular the proper motion of the system and PK parameters. In Sect. 4, we discuss the masses of the two neutron stars in this system and the three GR tests that are now possible with this system. The higher-order contributions in $\dot{\omega }$ are also considered, leading to the correction to the total mass. Finally, we summarise our results in Sect. 5.

2. Observations and data analysis

Most Arecibo observations in this work used the L-band receiver; for details, see Paper 2018. This set of observations ended on April 4, 2020. In addition, we have made a set of six observations with the Arecibo P-band receiver, which has a centre frequency of 327 MHz. All of these used the Puerto Rican Ultimate Pulsar Processing Instrument (PUPPI) back-end, which was based on its Green Bank predecessor (GUPPI, DuPlain et al. 2008). We have also obtained one full-orbit observation with the 800 MHz receiver of the Green Bank Telescope, with GUPPI as a back-end.

The FAST observations were carried out with the centre beam of FAST’s 19-beam receiver at a central frequency of 1250 MHz with a bandwidth of 500 MHz. The frequency resolution is 0.122 MHz. The signal is digitised in 8-bit, converted into 4-polarisations with a sample time of 49.152 μs and de-dispersed incoherently. Most of the FAST observations are 2 hours long. This yields high signal-to-noise ratio pulse profiles and, by covering full orbits, results in precise measurements of the orbital parameters. In Fig. 1, we displayed several examples of the observed pulse profiles from FAST observations and the template we used to derive ToAs.

Fig. 1. Frequency-averaged template used to fit ToAs and some examples of the observed pulse profile (Stokes I) displayed in red and black solid lines, respectively. We generated the template with the pulse profile of 2019-03-29 and aligned each profile with the centre of the two Gaussian functions that we used to fit the pulse, indicated by the blue dashed lines. The less severe profile evolution in the main pulse can be seen in this figure compared to that in the interpulse, which makes it reasonable to fit ToAs only with the main pulse. One can also notice that the separation between the main pulse and the interpulse is increasing over time.

The initial timing analysis for PSR J1946+2052 was described in Paper 2018. As in that analysis, we generated multi-frequency ToAs and used the TEMPO timing software³ to analyse them. First, the ToAs measured at the telescope are converted to the Bureau International des Poids et Mesures (BIPM) 2023 timescale. To subtract the telescope’s motion relative to the Solar System barycenter, we have used the International Earth Rotation Service data and the Jet Propulsion Laboratory’s DE 440 Solar System ephemeris (Park et al. 2021).

As reported in Meng et al. (2024), in the FAST observations, the pulse profile of PSR J1946+2052 is seen to change with time due to the relativistic spin precession, while in Arecibo observations, there’s no sign of profile evolution. Dealing with the time offsets introduced by profile evolution is necessary when using a single template to generate ToAs. The same situation can be found in several other relativistic binaries, such as PSR B1913+16 (Weisberg & Huang 2016), PSR B1534+12 (Fonseca et al. 2014) and PSR J1906+0746 (Desvignes et al. 2019). Unlike PSR B1913+16, but similarly to what happens in PSR J1906+0746, the profile of PSR J1946+2052 has a main pulse and an interpulse. The separation between these two main components is increasing, as shown in Fig. 1, and it is hard to identify the absolute movement of each pulse. We chose to make a standard template with only the main pulse, owing to its greater observed stability.

In order to measure the time offsets caused by profile evolution, we introduced three assumptions: (1) the shape of the emission beam is a circular hollow cone (Rankin 1983; Kramer 1998), (2) the centre of the main pulse is stable in the spin phase because it is the region with the fastest change in the position angle of the linear polarisation, which according to the rotating vector model (Radhakrishnan & Cooke 1969) makes it the point in the rotation when one of the magnetic poles is closest to our line of sight; furthermore this region showed less significant evolution in Meng et al. (2024) and (3) we made a simplifying assumption, which is that, relative to this fiducial point, any time offsets in different frequency channels of each epoch are caused by changes in the dispersion measure (DM) only, not due to the pulse profile evolution with time⁴.

Firstly, we scrunched the time, polarisation and frequency information to get the high signal-to-noise ratio pulse profile of each observation. We chose the first FAST observation (2019-03-29) as the template after making the template frequency-resolved by dividing it into 4 sub-bands and smoothing the pulse profile with paas in PSRCHIVE (Hotan et al. 2004). The frequency-averaged template is displayed in the top panel of Fig. 1. Secondly, we used two Gaussian functions to fit every main pulse and aligned each main pulse with the centre of the two Gaussian functions. This procedure ensures that the time offset would be 0 if the pulse profile of PSR J1946+2052 were stable. Finally, we used the Fourier domain algorithm of Taylor (1992) with Markov chain Monte Carlo (FDM, implemented in the pat program of PSRCHIVE) to fit the phase offsets in each main pulse, and convert them to ToA offsets.

The result of the ToA offsets is displayed in Fig. 2, wherein we also present the evolution of the half-width of the main pulse, denoted by width offsets. The significant correlation between the ToA offsets and width offsets indicates how profile evolution will affect the ToA measurement. Then we applied the linear relation between MJD and ToA offsets to all the ToAs we derived from FAST observations.

Fig. 2. Time-of-arrival offsets derived from the standard template of the first observation and the integrated pulse profile displayed with blue points. The pulse width (the unit is the same as the ToA offset) of each integrated pulse profile is plotted with red points. The correlation between the ToA offsets and pulse widths indicates the strong influence introduced by the profile evolution on measuring the ToAs. The blue solid line is the linear fit between MJD and ToA offsets, which is used to correct ToAs.

The pulse profile template from FAST observation uses a two-dimensional template that tracks the pulse profile as a function of both spin phase and radio frequency, whereas the Arecibo and GBT datasets rely on simpler one-dimensional templates. Therefore, these template definitions do not share the exact same fiducial pulse longitude as a function of radio frequency; combining the datasets leaves small, systematic offsets in both pulse phase and DM.

Regarding the definition of the reference longitude on the NS, this problem is difficult to solve, given that the radio pulses at low frequency (like L-band) are significantly scatter-broadened relative to high frequency (S-band). If we use TEMPO with more than one iteration, the phase offsets are translated to time offsets, which can bias the estimates of times of passage through periastron and thus of estimates of P_b and Ṗ_b (for a detailed discussion, see Guo et al. 2021). This is not a problem in this work because we always ran TEMPO with a single iteration. In this case, any offsets between ToA datasets are treated as phase offsets, originating from different definitions of the reference longitude of the NS. This can be done under the assumption that there are no significant offsets between the absolute timing of the different timing systems, which is known from other works on pulsars with better timing precision.

Regarding the DM offset between the FAST and Arecibo/GBT data, its exact value depends on the exact alignment of the pulse profiles at different frequencies in the FAST 2D template, which is affected by the DM that was used to create that profile. In Fig. 3, we display the DM measurements for a set of different epochs using the DMX model, a piecewise-constant function to describe the DM variation (Demorest et al. 2013). The uncorrected measurements are shown in panel (a). In this panel, the DM offset is very clear, and has a value of about 0.04 pc cm⁻³. After applying this DM offset onto the Arecibo and GBT ToAs (panel b), the DM variation during the entire observation period is still significant. This variation is mostly caused by the motion of the Earth and the pulsar: the radio waves travelling through the different paths from the pulsar to the Earth probe slightly different regions of the ionised interstellar medium (IISM), which will add to slightly different electron column densities. If the varying DM is not taken into account, it will introduce systematic uncertainties to the measurement of the pulsar properties. These DM variations can be used to probe the property of the IISM, which is shown in Appendix A.

Fig. 3. Dispersion measure variations derived from the DMX model with a time bin of 1 day. Panels (a) and (b) represent the DM variation before and after the DM correction, respectively. DM measurements from Arecibo and GBT are represented by black points, and those from FAST are represented by red points. The red dashed line indicates the final measurement of the DM, which is 93.9281 pc cm−3. We display the 10-order DM derivative fit in panel (b) with the blue dashed line.

We employed DM derivatives to describe these variations. The reason we chose DM derivatives rather than the DMX values in Fig. 3 is to reduce the number of parameters we fit in the model. To determine the number of DM derivatives to use, we compared the χ² and the reduced χ² from the fit of ToAs with different numbers of DM derivatives. The result is displayed in Fig. 4, which reveals that after the tenth DM derivative, the quality of the fit starts degrading, indicating over-fitting. We also display the DM variation predicted by our best 10-DM derivative model in panel (b) of Fig. 3, which provides a good description for the measured DM variation.

Fig. 4. Reduced χ2 and χ2 with different numbers of DM derivatives shown with black solid circles and red solid stars, respectively. The ten-order DM derivative fit generates the lowest χ2 and reduced χ2. Using a higher-order DM derivative will overfit, as indicated by the larger χ2 and reduced χ2 after ten DM derivatives.

Most of the results below are very robust, in the sense that they depend only weakly on the particular DM model used, with one exception, which is the measurement of the Shapiro delay parameters: given its small amplitude, the Shapiro delay is easily contaminated by systematic errors and thus depends significantly on the DM model. To deal with this issue, we decided in advance to report the solution with the number of derivatives that produces the lowest reduced χ², without knowing in advance whether the values of the Shapiro delay parameters in that solution would match the predictions of GR or not⁵. All results reported below are based on that model.

We used two binary models to examine the results for timing analysis, both based on the binary model of Damour & Deruelle (Damour & Deruelle 1986). The first is the ‘DDFWHE’ model, which, like the standard Damour & Deruelle (DD) model, is theory-independent. The only difference is that the Shapiro delay is parametrised differently; in this model the ‘orthometric parametrisation’ is used instead (Freire & Wex 2010); this is useful since it minimises the correlation of Shapiro delay parameters, providing a better description of the possible location of the system in the mass-mass (or mass – inclination) planes, and providing, as we show, improved tests of gravity theories. This model was implemented in TEMPO by Weisberg & Huang (2016). These theory-independent models are necessary not only for understanding which relativistic effects are detected in the data but also to quantify them precisely and test gravity theories.

The second binary model, known as the “DDGR” model, is a theory-dependent model that assumes that GR provides a correct description of all relativistic effects in the system (Taylor & Weisberg 1989). In this model, the only two unknowns are the two masses, or in the specific formalism of the DDGR model, the total system mass (M⁶) and the companion mass (m_c). Apart from immediate estimates of the masses, this model has the advantage of more easily detecting effects caused by, for instance, the acceleration of the system in the Galactic field, or alternatively the effects of spin-orbit coupling on the orbit of the system, because the difference between the observed value and the prediction from GR, such as in Ṗ_b (XPBDOT) can be measured directly.

3. Results

The parameters of our TEMPO fits are presented in Table 1. The ToA residuals (ToA minus the prediction of the DDFWHE timing solution in that table) are shown in Fig. 5. They display no observable trends, which implies that the DDFWHE solution provides, within their measurement accuracy, an adequate description of the observed ToAs. We now discuss some of the parameters in this timing solution that will be especially relevant for the following discussions.

Fig. 5. Residuals obtained using the DDFWHE timing solution in Table 1. The residuals in blue, orange, green and red are derived from L-band/PUPPI data, P-band/PUPPI data, single GBT observation and FAST data. Top: Residuals as a function of time. Bottom: Residuals as a function of orbital phase. The post-fit residuals’ root mean square is consistent with the ToA uncertainties, and the RMS of the residuals from Arecibo L-band, Arecibo P-band, GBT and FAST are 66.088 μs, 87.709 μs, 138.936 μs and 13.581 μs, respectively. No unmodelled trends are seen in the ToA residuals, indicating that, within measurement uncertainty, the DDFWHE timing solution provides an adequate description of the timing of the system.

3.1. Position and proper motion

We have improved the precision of the position measurement by combining TOAs from FAST, Arecibo and GBT. The timing solution provides a position with R.A. (J2000) 19:46:14.13475(4) and DEC. (J2000) +20:52:24.829(1) on 2017 Aug. 17. Combining the position with the DM, we are able to estimate the distance to the pulsar as d = 4.2 kpc based on the NE2001 model (Cordes & Lazio 2002) and 3.5 kpc based on the YMW16 model (Yao et al. 2017). The result indicates a relatively low galactic height of less than 0.2 kpc, suggesting that the system is not likely to have a large vertical velocity.

After 7 years of timing, we are able to constrain the proper motion with high precision: −1.5(1) mas yr⁻¹ in R.A. and −4.0(2) mas yr⁻¹ in Dec., leading to a total proper motion of μ = 4.2(2) mas yr⁻¹. With this proper motion and the distance from the YMW16 model, we derive the heliocentric transverse velocity of the system of v_T = 70(3) km s⁻¹. The position angle of the proper motion, in Galactic coordinates, is Θ_μ  =  260 ± 2 deg, which is 10 degrees South from the Western direction of the Galactic plane. The implications of this measurement for the evolution of the system will be discussed in a future publication. In what follows, the value of μ is used the most.

3.2. Spin period derivative

The measurement of the spin period derivative has been improved by four orders of magnitude. Its value is 1-σ consistent with the measurement in Paper 2018, 0.9(2)×10⁻¹⁸. This and the measurement of μ provide us with a chance to estimate the intrinsic spin-down rate. The observed spin-down rate is modified from the intrinsic one by the unknown radial velocity of the pulsar system in the form of the Doppler factor (D), as given by (Phinney 1993):

$$\begin{array}{c}\hfill P_{\mathrm{obs}}=D^{-1}{P}_{\mathrm{int}}={[1+({\mathbf{V}}_{\mathrm{PSR}}-{\mathbf{V}}_{\mathrm{SSB}})·\mathbf{n}/c]}^{-1}{P}_{\mathrm{int}},\end{array}$$(1)

where V_PSR is the pulsar system velocity, V_SSB is the velocity of the Solar System Barycentre (SSB), n is the unit vector pointing from the SSB to the pulsar system, P_obs is the observed spin period and P_int is the intrinsic spin period. By differentiating the equation above, one gets (Phinney 1993)

$$\begin{array}{c}\hfill {\dot{P}}_{\mathrm{obs}}={\dot{P}}_{\mathrm{int}}+{\dot{P}}_{\mathrm{gal}}+{\dot{P}}_{\mathrm{shk}},\end{array}$$(2)

where the second term results from the differential Galactic acceleration between the Solar and the pulsar system, and the third term is the Shklovskii term (Shklovskii 1970). The last two terms can be calculated according to

$$\begin{array}{c}\hfill \frac{{\dot{P}}_{\mathrm{gal}}}{P}=-\frac{({\mathbf{a}}_{\mathrm{PSR}}-{\mathbf{a}}_{\mathrm{SSB}})·\mathbf{n}}{c}\text{and}\frac{{\dot{P}}_{\mathrm{shk}}}{P}=-\frac{{\mu }^{2}d}{c},\end{array}$$(3)

where a_PSR and a_SSB are the accelerations of the pulsar system and the SSB in the gravitational potential of the Milky Way (see also Damour & Taylor 1991).

To calculate the accelerations of the pulsar system and the SSB, we adopted the Galactic potential model MWPotential2014 deployed in the Python package galpy⁷ (Bovy 2015), scaled such that the distance from the Sun to the Galactic centre is R₀ = 8.275(34) kpc (GRAVITY Collaboration 2021), and the circular velocity of the Sun’s local standard of rest is Θ₀ = 240.5(41) km s⁻¹ (Guo et al. 2021). For our calculations of the differential Galactic acceleration, we further assumed the Sun’s height above the local Galactic mid-plane as Z_⊙ = 20.8(3) pc (Bennett & Bovy 2019).

Then we used the two distances determined by YMW16 and NE2001 as Gaussian distributions with 20% uncertainty (Cordes & Lazio 2002) and derived the contribution from the Galactic acceleration and the Shklovskii effect. The results are presented in Table 2. The resulting Ṗ_int values from the two DM models are consistent with each other as 1.116(2)×10⁻¹⁸ s s⁻¹, and are larger than Ṗ_obs by a factor of 1.002(1). Using the expressions in Lorimer & Kramer (2005), this results in a characteristic age of 0.24 Gyr, a B-field of 4.4  ×  10⁹ G and, by using a moment of inertia as I = 1.31 × 10⁴⁵ g cm² (see in Appendix B), a spin-down energy of 1.183(2)×10³⁴ erg s⁻¹.

3.3. Rate of advance of periastron

As mentioned earlier, the defining feature of PSR J1946+2052 is its extremely short orbital period, P_b = 0.0785 d. This parameter primarily accounts for the fact that, as outlined in Paper 2018, this binary pulsar exhibits the highest recorded rate of periastron advance, $\dot{\omega }=25.6(3)deg{\mathrm{yr}}^{-1}$.

This parameter has been massively improved in this work: our current measurement, $\dot{\omega }=25.79205(40)deg{\mathrm{yr}}^{-1}$. This is 1-σ compatible with the measurement presented in Paper 2018, but three orders of magnitude more precise. Assuming that this is as predicted by GR, then the total mass of the system in solar mass parameters is, to the leading post-Newtonian order, given by Robertson (1938)

$$\begin{array}{c}\hfill M=\frac{1}{T_{\odot }}{\left[\frac{\dot{\omega }}{3}(1-e^2)\right]}^{\frac{3}{2}}{n}_{\mathrm{b}}^{-\frac{5}{2}},n_{\mathrm{b}}\equiv \frac{2\pi }{P_{\mathrm{b}}},\end{array}$$(4)

where $T_{\odot }\equiv {(\mathcal{G}\mathcal{M})}_{\odot }^{\mathrm{N}}/c^3=4.925490947641266978...\mathrm{\mu }\mathrm{s}$ is an exact quantity, the nominal solar mass parameter ${(\mathcal{G}\mathcal{M})}_{\odot }^{\mathrm{N}}$ in time units (Prša et al. 2016)⁸ and e is the orbital eccentricity. Using the values for these quantities from Table 1, we obtain M  =  2.531858(60) M_⊙. This is the lowest total mass measured among known DNS systems in the Galactic field, with the possible exception of PSR J1411+2551 (M  =  2.538(22) M_⊙, Martinez et al. 2017). Also, as noted in the Introduction, this DNS has the lowest eccentricity for any known DNS in the Galaxy.

3.4. Einstein delay

Whereas the timing solution in Paper 2018 is sensitive to only one PK parameter ($\dot{\omega }$), our enhanced timing precision and extended timeline of the dataset enabled the detection of four additional PK parameters. The first is the Einstein delay amplitude, $\gamma =0.2591\pm 0.0006\mathrm{ms}$, which is ∼430-σ significant. In GR, this is given by

$$\begin{array}{c}\hfill \gamma =g\frac{m_{\mathrm{c}}(M+m_{\mathrm{c}})}{M^{\frac{4}{3}}},\end{array}$$(5)

where

$$\begin{array}{c}\hfill g=e{n}_{\mathrm{b}}^{-\frac{1}{3}}{T}_{\odot }^{\frac{2}{3}},\end{array}$$(6)

and where m_c is the companion mass. Although the small values of e and $n_{\mathrm{b}}^{-\frac{1}{3}}$ make our value of g the smallest for any DNS system, which could result in a large relative uncertainty, our measurement of γ is still extremely precise with an uncertainty of 600 ns, owing to the large precession angle ($\mathrm{\Delta }\omega =T\dot{\omega }=184.09deg$) that has been covered over the length of our timing baseline (T).

With a measurement of M and γ, we can make a first estimate of the individual masses of the components. The mass of the companion is given by (Ridolfi et al. 2019)

$$\begin{array}{c}\hfill m_{\mathrm{c}}=\frac{1}{2}(\sqrt{M^2+4M^{\frac{4}{3}}\frac{\gamma }{g}}-M),\end{array}$$(7)

from which we get m_c = 1.2480(21) M_⊙. The mass of the pulsar is given by m_p = M − m_c = 1.2838(21) M_⊙. This indicates that PSR J1946+2052 is a slightly asymmetric system. Rewriting the mass function equation, we obtain

$$\begin{array}{c}\hfill sini=\frac{x}{m_{\mathrm{c}}}{T}_{\odot }^{-\frac{1}{3}}{\left(n_{\mathrm{b}}M\right)}^{\frac{2}{3}},\end{array}$$(8)

where i is the orbital inclination. From this we obtain sin i = 0.9599(16), which corresponds to i = 73.71(33) deg or i = 106.29(33) deg. This is 2-σ consistent with i = 63°_−3°⁺⁵° derived by Meng et al. (2024), which includes the small misalignment between the spin axis of the pulsar and the orbital angular momentum that they estimated.

3.5. Shapiro delay

In our timing, we also detect, for the first time, the Shapiro delay. To quantify this detection, we adopted the orthometric parametrisation introduced by Freire & Wex (2010). The orthometric amplitude h₃ and the orthometric ratio ς of Shapiro delay are given in GR by

$$\begin{array}{c}\hfill \varsigma =\frac{sini}{1+\sqrt{1-{sin}^{2}i}},h_3=m_{\mathrm{c}}{T}_{\odot }{\varsigma }^3.\end{array}$$(9)

The advantage of using such a parametrisation relative to the r, s parametrisation in the DD model is that it reduces the correlation between the two parameters of the Shapiro delay. Furthermore, as we show below, the h₃ test also provides a more precise test of GR compared to the r parameter used in the DD model.

This parametrisation has been implemented in the DDFWHE timing model in TEMPO. The values measured from the timing of PSR J1946+2052 derived the values of h₃ = 2.58(14) μs and ς = 0.760(36), i.e. this is a highly significant detection of the Shapiro delay. From these, we derive sin i = 0.963(13) and m_c = 1.21(19) M_⊙. These values are consistent with the values derived above from $\dot{\omega }$ and γ, but less precise. The larger deviation from i = 90° implies that the Shapiro delay signal is much weaker than in the Double Pulsar (Kramer et al. 2021), hence the much larger relative uncertainties in the measurement of the Shapiro delay parameters.

Using the DD model, we obtain s = 0.962(13) and m_c, S = r/T_⊙ = 1.21(21), which are again consistent with the values derived from the DDFWHE solution. The values of i derived from the ς in the DDFWHE solution and s in the DD solution are both i = 74 ± 2 deg. This is not surprising given the exclusive dependence of ς and s on i, but is important for the interpretation of the GR tests made with ς and s.

3.6. Orbital decay

As a result of our timing, we now have a 1200-σ significant measurement or the observed variation of the orbital period: Ṗ_b,obs = −1829.6 ± 1.5 fs/s. Most of this is due to orbital decay induced by GW damping. In GR, this orbital decay (Ṗ_b,GR) is given by (Peters 1964):

$$\begin{array}{c}\hfill {\dot{P}}_{\mathrm{b},\mathrm{GR}}=-\frac{192\pi }{5}{\left(T_{\odot }{n}_b\right)}^{\frac{5}{3}}f(e)\frac{m_{\mathrm{c}}(M-m_{\mathrm{c}})}{M^{\frac{1}{3}}},\end{array}$$(10)

where

$$\begin{array}{c}\hfill f(e)=\frac{1+(73/24)e^2+(37/96)e^4}{{(1-e^2)}^{\frac{7}{2}}};\end{array}$$(11)

these expressions were re-written as a function of M because of its high precision but also because m_c and m_p are not determined independently.

Similar to Ṗ_obs, the value of Ṗ_b,obs is affected by the variation of the Doppler factor exactly the same way as $\dot{P}$, i.e. with contributions from the relative Galactic accelerations (Ṗ_b,gal) and the Shklovskii effect (Ṗ_b,shk). Therefore, we can modify Eq. (2) as

$$\begin{array}{c}\hfill {\dot{P}}_{\mathrm{b},\mathrm{obs}}={\dot{P}}_{\mathrm{b},\mathrm{int}}+{\dot{P}}_{\mathrm{b},\mathrm{ext}}={\dot{P}}_{\mathrm{b},\mathrm{int}}+{\dot{P}}_{\mathrm{b},\mathrm{gal}}+{\dot{P}}_{\mathrm{b},\mathrm{shk}},\end{array}$$(12)

where Ṗ_b,int is the intrinsic variation of the orbital period, which is dominated by the emission of GW; Ṗ_b,ext is the external contribution of the orbital decay; and

$$\begin{array}{c}\hfill {\dot{P}}_{\mathrm{b},\mathrm{gal}}=\frac{{\dot{P}}_{\mathrm{gal}}}{P}{P}_{\mathrm{b}},{\dot{P}}_{\mathrm{b},\mathrm{shk}}=\frac{{\dot{P}}_{\mathrm{shk}}}{P}{P}_{\mathrm{b}}.\end{array}$$(13)

As in Sect. 3.2, we used the distances from YMW16 and NE2001 and assumed a 20% uncertainty. Considering the uncertainty of the proper motion and using the distance from YMW16, Ṗ_b,gal and Ṗ_b,shk can be determined as −1.93(62) fs s⁻¹⁹ and .05(23) fs s⁻¹, respectively.

The sum of the external effects (for this distance and a 20% uncertainty) is thus given by −0.88(66) fs s⁻¹, which is consistent with ΔṖ_b = Ṗ_b − Ṗ_b,GR calculated in Sect. 4.1, the latter is indicated by the gray bar in Fig. 6. The precision of this sum is mostly limited by the uncertainty in the distance, although the uncertainty of the proper motion remains important.

Fig. 6. Contributions of Ṗb from differential Galactic acceleration and the Shklovskii effect as a function of the distance to the pulsar (see text for details). The orange and blue areas indicate the 1-σ confidence of these two effects. The solid lines represent the nominal values, and the dotted curves are the ±1-σ uncertainties. The total external contribution of Ṗb is shown by the corresponding red area and curves. The uncertainties come from the uncertainty of R0, Θ0 and the proper motion. The dashed red and black lines are the distances derived from YMW16 and NE2001. The grey area displays ΔṖb = Ṗb,obs − Ṗb,GR. This quantity, estimated in detail in Sect. 4.1, represents a measurement of the total external contribution that assumes the validity of GR; its width represents its uncertainty, which is dominated by the error of the observed Ṗb.

This implies that Ṗ_b,int = −1828.8(16) fs s⁻¹, which is consistent with Ṗ_b,GR. The precision of this parameter is mostly limited by the precision of Ṗ_b,obs, but this will improve quickly in the near future. All Ṗ_b values derived assuming the NE2001 distance are 1-σ consistent, and they are presented in Table 3.

Another contribution that could change the orbital period is the mass loss due to the spin-down of the pulsar. Damour & Taylor (1991) estimated it as

$$\begin{array}{c}\hfill \frac{{\dot{P}}_{\mathrm{b}}^{{\dot{m}}_{\mathrm{p}}}}{P_{\mathrm{b}}}=8{\pi }^2\frac{I_{\mathrm{p}}{\dot{P}}_{\mathrm{int}}}{M{c}^2{P}^3},\end{array}$$(14)

where I_p is the moment of inertia of the pulsar. By using the total mass derived in Sect. 3.3 and the timing parameters in Table 1, ${\dot{P}}_{\mathrm{b}}^{{\dot{m}}_{\mathrm{p}}}$ can be calculated as 3.5(2)×10⁻¹⁷ s s⁻¹ by using I_p = (1.31 ± 0.08)×10⁴⁵ g cm² from Appendix B, which is completely ignorable in the current analysis in the orbital decay.

Owing to the short spin period (17 ms) and the small magnetic field strength at the surface (4.4 × 10⁹ G), it is clear that the pulsar is the first-born neutron star in this binary system and the companion is the second-born NS. Calculating the orbit into the past, we can see there were a few times when it crossed the Galactic plane within the characteristic age of the pulsar. By assuming the system formed at the time of the last crossing, we can estimate a conservative lower limit of the characteristic age of the companion with tan b/μ_b ≈ 10 Myr, where b is the Galactic latitude and μ_b is the Galactic latitude component of the proper motion.

Since the companion is the second-born (hence non-recycled) NS, the location of the companion in the P−Ṗ diagram will be in the region of normal pulsars. Therefore, we can simply estimate the lower limit of the companion’s spin period (P_c) and the first derivative of the spin period (Ṗ_c) by using the characteristic age calculated above in the P−Ṗ diagram (Lorimer & Kramer 2005), resulting in P_c ≥ 180 ms and Ṗ_c ≈ 6 × 10⁻¹⁵ s s⁻¹. Utilising Eq. (14), the upper limit on the change of the orbital period caused by the mass loss of the companion could be estimated as 1.8 × 10⁻²⁰ s s⁻¹, which can be completely neglected in the current analysis of the orbital decay.

4. Masses and tests of general relativity

The measurement of 5 PK parameters enables us to determine the masses of the pulsar and the companion, and perform three tests of GR by testing the self-consistency of Eqs. (4), (5), (9) and (10) for the Keplerian and PK parameters of the system.

This consistency is depicted graphically in the mass-mass diagram in Fig. 7. In this diagram, we highlight the regions where the masses are 1-σ consistent with the measurements of $\dot{\omega }$, γ, h₃, ς and Ṗ_b using the aforementioned equations; these form bands in the diagram. As we can see, all bands meet at a small common region perfectly even in the enlarged panel, indicating GR passed the three tests. Next, we elaborate on these two aspects.

Fig. 7. Mass-mass diagram for PSR J1946+2052. The regions consistent with the measured ${\dot{P}}_{\mathrm{b}},\dot{\omega },\gamma ,h_3,\varsigma$, and their 1-σ uncertainties are displayed in red, green, blue, purple, and orange. Note that Ṗb in this diagram has been corrected by removing the external contributions. The grey area is excluded by the mass function and sin i ≤ 1. The inset is an expanded view, showing in more detail the intersection of all PK parameters and the self-consistency of GR. The $\dot{\omega }$-line is based on the leading-order Equation (4). The (small) modifications related to the next-to-leading order contributions to the periastron advance are discussed in detail in Sect. 4.4.

4.1. Masses of the two neutron stars

Using the DDGR model, which assumes the validity of GR (and furthermore, that the observed PK parameters are relativistic), we can fit directly for the masses. By doing so, we obtained M = 2.531837(24) M_⊙ and m_c = 1.2476(20) M_⊙, and thus m_p = 1.2842(21) M_⊙. Unsurprisingly, these values coincide exactly with our estimates above based on $\dot{\omega }$ and γ; this indicates that it is the presence of the relativistic effects associated with $\dot{\omega }$ and γ in the timing of the pulsar that is chiefly responsible for the precise DDGR mass estimates. It corresponds to sin i = 0.960(2) and the orbital inclination angle i = 73.7(3)° (or 106.3(3)°).

In order to better estimate the uncertainties in these parameters, we made a χ² map in a three-dimensional space where the coordinates are M, cos i and ΔṖ_b. This used the same basic Bayesian methodology as Splaver et al. (2002), where the χ² obtained by a DDGR model at each point of that space (with the corresponding values of MTOT, M2 and XPBDOT fixed to their values at that point) is used to derive a 3D probability distribution function (pdf). In Fig. 8, this is projected into four 2D spaces that have cosi, ΔṖ_b, m_c and m_c as axes, and finally the pdfs are marginalised along those axes. From this analysis, we obtain: M = 2.53186(4) M_⊙, m_c = 1.248(2) M_⊙, m_p = 1.284(2) M_⊙, cosi = 0.280(5) and $\mathrm{\Delta }{\dot{P}}_{\mathrm{b}}=-0.9_{-1.5}^{+1.4}\mathrm{fs}{\mathrm{s}}^{-1}$. These estimates and their uncertainties agree with the simpler estimates above, indicating that the uncertainties obtained by the simple DDGR fit are essentially accurate. Apart from measuring the masses and their uncertainties, the detailed mapping of the allowed masses enables a precise prediction of Ṗ_b,GR and an estimate of ΔṖ_b, where the correlation between m_c and m_p is taken into account.

Fig. 8. Constraints on the orbital inclination angle, masses of the pulsar and the companion, and the Ṗb deviation from the measurement and the prediction of GR (ΔṖb). The black, grey, and light grey contours represent the 1-σ, 2-σ, and 3-σ confidence levels of these parameters. We show the probability density of each parameter on the upper and right edges, where the grey area and solid grey line represent the 1-σ confidence area and the nominal values; the dashed lines at the edges represent the ±1-σ limits.

4.2. Testing general relativity with the orbital decay

The masses we used in the previous section are calculated from $\dot{\omega }$ and γ that are derived from the DDFWHE model. Then, we use these masses to calculate the predicted orbital decay under GR (Ṗ_b,GR). The difference between Ṗ_b,int and Ṗ_b,GR (Ṗ_b,xs) is at most 0.3 ± 1.6 fs s⁻¹, which indicates a non-detectable difference. The ratio of ${\dot{P_{\mathrm{b}}}}^{\mathrm{int}}$ to Ṗ_b,GR is given by

$$\begin{array}{c}\hfill \frac{{\dot{P}}_{\mathrm{b},\mathrm{int}}}{{\dot{P}}_{\mathrm{b},\mathrm{GR}}}=\frac{-1.8288(16)\times {10}^{-12}}{-1.8287(1)\times {10}^{-12}}=1.00005(91),\end{array}$$(15)

and thus Ṗ_b,int is fully consistent with the predicted value, Ṗ_b,GR.

This radiative test is about twice as precise as the test with the Hulse-Taylor pulsar (Weisberg & Huang 2016; Deller et al. 2018), making it the second most precise test of the quadrupole formula; however, it is still one order of magnitude behind the precision achieved with the Double Pulsar (Kramer et al. 2021).

In Fig. 9, we see the cumulative shift of periastron time. This is the difference between the time of passage through the periastron measured at any particular time and what its value would be if the orbital period had stayed constant at its value at the start of the timing. The linear change in the orbital period predicted by GR results in a quadratic shift in the time of periastron that is indicated by the red solid line. In black, we measured the times of passage of periastron with subsets of TOAs, where all other timing parameters are fixed. As we can see, there is an excellent match between the GR prediction and the measurements of T₀. This plot does not only tell us that the magnitude of the orbital decay is as expected by GR, but also shows that this decay proceeds linearly with time as expected from a constant rate of loss of orbital energy.

Fig. 9. Cumulative shift of periastron time relative to the predictions of a constant orbit caused by the decay of the orbital period. In the upper panel, the black points are the measurements of the periastron time shift, and the red solid line is the GR prediction, which is the sum of the GW damping calculated from GR and the calculated contribution from the Galactic acceleration and the Shklovskii effect. We display the residual between the measurements and the GR prediction in the lower panel.

Note that the result is derived from 7 years of observations with Arecibo, GBT and FAST. Considering the observation duration of 31 years for the Hulse-Taylor pulsar and 16 years for the Double Pulsar, with the follow-up observations of FAST, PSR J1946+2052 has the potential to increase the precision of its test significantly as the measurement for Ṗ_obs improves according to $T_{\mathrm{obs}}^{5/2}$ where T_obs is the observed time span (Damour & Taylor 1992).

4.3. Testing general relativity with the Shapiro delay

The $\dot{\omega }$–γ–Ṗ_b test of the previous section is a mixed test which combines quasi-stationary and radiative strong field effects. The additional detection of the Shapiro delay in the PSR J1946+2052 system allows for a different type of test, a purely quasi-stationary test that combines orbital and signal propagation effects. Therefore, such a test probes different aspects of gravity compared to the $\dot{\omega }$–γ–Ṗ_b test, in particular the propagation of electromagnetic signals in the spacetime of a strongly self-gravitating body (cf. Taylor et al. 1992; Damour 2009; Wex 2014). Although the measurements of the Shapiro delay parameters h₃ and ς are much less significant than the measurement of Ṗ_b, they still provide two tests for GR with good precision. The GR predictions are calculated from Eq. (7), (8) and (9), where M and m_c were derived from $\dot{\omega }$ and γ using the GR equations. Dividing the measured parameters by their GR predictions, we obtained

$$\begin{array}{cc}& \frac{h_3}{h_{3,\mathrm{GR}}}=\frac{2.58(14)\times {10}^{-6}}{2.59(4)\times {10}^{-6}}=0.997(56),\hfill \end{array}$$(16)

$$\begin{array}{cc}& \frac{\varsigma }{{\varsigma }_{\mathrm{GR}}}=\frac{0.760(36)}{0.750(4)}=1.01(5),\hfill \end{array}$$(17)

indicating a good consistency.

We can compare these with the GR tests provided by the r,s parametrisation of the DD model:

$$\begin{array}{cc}& \frac{m_{\mathrm{c},\mathrm{S}}}{m_{\mathrm{c},\mathrm{GR}}}=\frac{1.21(21){\mathrm{M}}_{\odot }}{1.248(2){\mathrm{M}}_{\odot }}=0.97(15),\hfill \end{array}$$(18)

$$\begin{array}{cc}& \frac{sin{i}_{\mathrm{S}}}{sin{i}_{\mathrm{GR}}}=\frac{0.962(13)}{0.9599(16)}=1.002(13).\hfill \end{array}$$(19)

The results are more precise than those in the Hulse-Taylor pulsar (Weisberg & Huang 2016) but less precise than PSR B1534+12 (Fonseca et al. 2014) and the Double Pulsar (Kramer et al. 2021), owing to the intermediate orbital inclination angle of 73.7(3)°.

These measurements provide a direct comparison of the h₃ and r tests, with the former being 2.5 times more precise, apart from being less correlated with ς than r is to s. This illustrates the fact highlighted by Freire & Wex (2010) that, by decreasing the correlation between the parameters that quantify the Shapiro delay, the orthometric parametrisation provides an improved GR test relative to the “classical” r-s parametrisation. The difference between these two tests can be illustrated in Fig. 7: in this diagram, the r constraint from Shapiro delay (m_c = 1.21(21) M_⊙) would occupy a much wider region of the diagram compared to the h₃ band.

Note, on the other hand, that the s test seems to be more precise than the ς test. This is a purely numerical artefact, caused by the fact that as angles come closer to 90°, the variations of sin i become very small. As discussed in Sect. 3.5, the ς and s constraints result in very similar orbital inclinations, which implies similar bands for these two parameters in Fig. 7.

4.4. Effect of higher-order contributions on the measurement of the total mass

Given the precise measurement of $\dot{\omega }$, we can explore the higher-order corrections on this parameter. The $\dot{\omega }$ can be written as

$$\begin{array}{c}\hfill \dot{\omega }={\dot{\omega }}^{1\mathrm{PN}}+{\dot{\omega }}^{2\mathrm{PN}}+{\dot{\omega }}^{\mathrm{LT}},\end{array}$$(20)

where ${\dot{\omega }}^{1\mathrm{PN}}$ and ${\dot{\omega }}^{2\mathrm{PN}}$ are the first and second post-Newtonian contributions, ignoring spin contributions, and ${\dot{\omega }}^{\mathrm{LT}}$ is the Lense-Thirring (LT) precession contribution caused by the coupling between the spins of the binary components and the orbital motion. Note that the LT term in Eq. (20) is completely dominated by the spin of the (observed) pulsar due to the expected large spin period of the companion (see details in Sect. 3.6).

The 1PN and 2PN contributions of $\dot{\omega }$ are given by (Robertson 1938; Damour & Deruelle 1985, 1986; Damour & Schäfer 1987, 1988)

$$\begin{array}{cc}& {\dot{\omega }}^{1\mathrm{PN}}=3n_{\mathrm{b}}^{\frac{5}{3}}{(T_{\odot }M)}^{\frac{2}{3}}{(1-e^2)}^{-1},\hfill \end{array}$$(21)

$$\begin{array}{cc}& {\dot{\omega }}^{2\mathrm{PN}}=3n_{\mathrm{b}}^{\frac{7}{3}}{(T_{\odot }M)}^{\frac{4}{3}}{(1-e^2)}^{-1}{f}_{\mathrm{O}},\hfill \end{array}$$(22)

and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill f_{\mathrm{O}}=& \frac{1}{1-e^2}(\frac{39}{4}{X}_{\mathrm{p}}^2+\frac{27}{4}{X}_{\mathrm{c}}^2+15X_{\mathrm{p}}{X}_{\mathrm{c}})\hfill \\ \hfill & -(\frac{13}{4}{X}_{\mathrm{p}}^2+\frac{1}{4}{X}_{\mathrm{c}}^2+\frac{13}{3}{X}_{\mathrm{p}}{X}_{\mathrm{c}}),\hfill \end{array}\end{array}$$(23)

where X_p ≡ m_p/M and X_c ≡ m_c/M.

If the misalignment angle between the spin axis and the orbital angular momentum is small (in the case of PSR J1946+2052 Meng et al. 2024 estimate a misalignment angle of only ∼0.2deg), then ${\dot{\omega }}^{\mathrm{LT}}$ is given by (Barker & O’Connell 1975; Damour & Schäfer 1988; Freire & Wex 2024)

$$\begin{array}{c}\hfill {\dot{\omega }}^{\mathrm{LT}}=-n_{\mathrm{b}}^2(T_{\odot }M){(1-e^2)}^{-3/2}(3+X_{\mathrm{p}})X_{\mathrm{p}}{\chi }_{\mathrm{p}}.\end{array}$$(24)

Here, χ_p is the dimensionless spin of the pulsar, given by

$$\begin{array}{c}\hfill {\chi }_{\mathrm{p}}=\frac{G}{c^5}{(T_{\odot }{m}_{\mathrm{p}})}^{-2}2\pi {\nu }_{\mathrm{p}}{I}_{\mathrm{p}},\end{array}$$(25)

and ν_p is the spin frequency of the pulsar.

From these equations, and the masses of the components of the system, we derive ${\dot{\omega }}^{2\mathrm{PN}}=0.00078422(2){\mathrm{deg}\mathrm{yr}}^{-1}$ and ${\dot{\omega }}^{\mathrm{LT}}=-0.000884(1)\times I_{\mathrm{p}}^{(45)}{\mathrm{deg}\mathrm{yr}}^{-1}$. One thing that is important to emphasise is that both values are significantly larger than the corresponding values for PSR J0737−3039A. In particular, for the latter pulsar, Kramer et al. (2021) estimate ${\dot{\omega }}^{\mathrm{LT}}=-0.000377\times I_{\mathrm{A}}^{(45)}{\mathrm{deg}\mathrm{yr}}^{-1}$, a factor of 2.3 smaller than for PSR J1946+2052, which can only be in a small part compensated by the slightly larger value of I_A compared to I_p. The reason is that, as shown by Eq. (24), the effect increases with the square of the orbital frequency and thus (according to Kepler’s law) with the inverse cube of the orbital separation of the masses. Furthermore, as shown by Eq. (25), the spin frequency is important, and PSR J1946+2052 has the fastest spin of any known pulsar in a DNS system.

A consequence of the magnitude of these effects is that they are already larger than the uncertainty of $\dot{\omega }$, $\delta \dot{\omega }=0.00040{\mathrm{deg}\mathrm{yr}}^{-1}$. This implies that to measure M and its uncertainty accurately (in GR), we have to take the effect of ${\dot{\omega }}^{2\mathrm{PN}}$ and ${\dot{\omega }}^{\mathrm{LT}}$ into account, as in the case of the PSR J0737−3039 system (Kramer et al. 2021). The value of $I_{\mathrm{p}}^{(45)}$ is determined as 1.31(8) by employing different equations of state (see Appendix B for details), indicating that ${\dot{\omega }}^{\mathrm{LT}}=-0.00117(7){\mathrm{deg}\mathrm{yr}}^{-1}$, which is 3 times larger than $\delta \dot{\omega }$.

The total mass of the system in Sects. 3 and 4 are derived from $\dot{\omega }$ with only 1PN contribution. Now we start to consider the 2PN and LT contributions. One way to distinguish the significance of the LT contribution in $\dot{\omega }$ is to ignore the LT term and assume the observed $\dot{\omega }$ consists of 1PN and 2PN contributions. By doing this, we can plot the curve of $\dot{\omega }$ in the mass-mass diagram, as indicated by ${\dot{\omega }}_0$. Then we plot the observed $\dot{\omega }$ with the full expression, including ${\dot{\omega }}^{2\mathrm{PN}}$ and ${\dot{\omega }}^{\mathrm{LT}}$. From Fig. 10, one notes that the 1-σ areas of ${\dot{\omega }}_0$ and $\dot{\omega }$ are already separated, which indicates different values of the total mass of the system when we consider the LT contribution or not. In Table 4, we present the calculations of the masses with different contributions in $\dot{\omega }$. The masses of the pulsar and the companion are all consistent at -1σ confidence level, indicating that the analyses in the previous sections will not be affected by the higher-order contributions of $\dot{\omega }$. However, the total mass shows a significant difference, which means the 2PN and LT contributions have already had impact on the measurement of the total mass.

Fig. 10. Enlarged mass-mass diagram for PSR J1946+2052. To demonstrate the influence introduced by the LT effect, we display the total mass constraint for ${\dot{\omega }}_0$ with black dashed lines, where we ignore the LT contribution. The black solid lines represent $\dot{\omega }$ with the full contribution, including the 2PN and LT effect. The significant deviation between ${\dot{\omega }}_0$ and $\dot{\omega }$ indicates that the LT contribution significantly affects our measurement of the total mass. The magnitude of this shift depends on the MOI of the pulsar. Owing to the much larger uncertainty of Ṗb, the total mass cannot be measured independently with sufficient precision to determine the MOI.

Another possible source of uncertainty is the contribution from the proper motion. Rewriting the equations of Kopeikin (1996), one obtains:

$$\begin{array}{c}\hfill {\dot{\omega }}_{\mathrm{K}}=\frac{\mu }{sini}cos({\mathrm{\Theta }}_{\mu }-\mathrm{\Omega }),\end{array}$$(26)

where Θ_μ and Ω are the position angles of the proper motion and the line of nodes, respectively. Using the values above, we obtain at most a contribution of 1.3 × 10⁻⁶ deg yr⁻¹, which is two orders of magnitude smaller than the measurement uncertainty.

5. Conclusions

In this work, we have presented a comprehensive timing analysis of the DNS system PSR J1946+2052 using the observations from Arecibo, GBT, and FAST. After dealing with the profile evolution caused by the relativistic spin precession of the pulsar, we derived ToAs from the 7-year period of observations. These were used to derive precise timing solutions using the DDFWHE and DDGR binary models.

These solutions led to precise measurements of the proper motion and of five PK parameters, enabling a robust determination of precise masses for the pulsar and the companion and three tests of GR. The observed orbital decay deviates from the GR prediction by less than 9 × 10⁻⁴ (68% C. L.), making this the second-most precise test of the quadrupole formula to date. The Shapiro delay measurement in this pulsar produced two additional GR tests. Overall, GR passed all three tests.

Owing to the short orbital period, the higher-order contributions (2PN and LT effect) to $\dot{\omega }$ are significantly larger than those for the Double Pulsar system. By adopting multiple EOSs and multi-messenger constraints on the radius of the neutron star, we derived the moment of inertia and the LT effect on $\dot{\omega }$, which is three times larger than the measurement uncertainty for $\dot{\omega }$. This means that these higher-order contributions must be taken into account to derive accurate estimates of the total mass of the system.

Future observations from FAST will continue to rapidly improve the precision of Ṗ_b,obs. However, the lack of precise distance measurements that are independent from the Galactic electron density models will preclude a significant reduction of the uncertainty of Ṗ_b,ext, which is now of the order of 0.7 fs/s, or about half the current uncertainty of Ṗ_b,obs. Thus, to improve the precision of the radiative test (i.e. of Ṗ_b,ext and Ṗ_b,int) by more than a factor of two, a precise independent measurement of the distance will be necessary. This uncertainty of the measurement of Ṗ_b,int precludes a precise independent determination of the total mass of the system, which would be necessary to extract ${\dot{\omega }}^{\mathrm{LT}}$ from the total observed $\dot{\omega }$ and determine the MOI of the pulsar. The future VLBI observations operated by SKA will greatly improve the measurement of the distance and improve the precision of Ṗ_b,int, which will then provide a more precise radiative test and possibly allow a measurement of the MOI for this pulsar.

Another method to determine the MOI is measuring the LT contribution in $\dot{x}({\dot{x}}^{\mathrm{LT}})$, a result of the LT-induced secular variation of the orbital inclination i (see e.g. Damour & Schäfer 1988; Damour & Taylor 1992). If future observations confirm the result in Meng et al. (2024), at least in its order of magnitude, then the maximum of ${\dot{x}}^{\mathrm{LT}}$ is of the order of 3 × 10⁻¹⁶ s s⁻¹, which is considerably smaller than the current constraint on ẋ with a level of 3 × 10⁻¹⁴ s s⁻¹.

The ongoing observations of PSR J1946+2052 will continue to monitor the observed changes of the pulse profile and polarisation (Meng et al. 2024), which not only will improve the constraints on the orbital geometrical parameters, including solid measurement of the misalignment angle, but might also result in an additional test of GR in this system via the measurement of the geodetic precession rate. This test might be particularly favoured since for PSR J1946+2052 we detected a strong interpulse (see Fig. 1). In most cases similar to this, one can achieve an improved determination of the spin geometry from the rotating vector model, which results in improved measurements of the geodetic precession rate (Stairs et al. 2004; Desvignes et al. 2019).

Continued observations will also greatly improve the precision of the proper motion, which will improve the estimate of the pulsar’s transverse velocity relative to that of its Galactic co-rotating frame. The combination of the spatial velocity and the misalignment angle will yield estimates of the kick associated with the second supernova in this system and provide an improved understanding of the evolution of the PSR J1946+2052 system.

Acknowledgments

We thank Cijie Zhang and Yukai Zhou for deploying and maintaining the server used throughout this research. This work was supported by the National Natural Science Foundation of China (12041303, 12421003, 12203072, 12203070), the National SKA Program of China (2020SKA0120200, 2020SKA0120300), Beijing Nova Program (No. 20250484786) and the CAS-MPG LEGACY project. Pulsar research at UBC is supported by an NSERC Discovery Grant and by the Canadian Institute for Advanced Research. LM gratefully acknowledges the support of the China Scholarship Council (CSC) during the visit to the Max-Planck Institute for Radio Astronomy in Bonn, Germany. PCCF gratefully acknowledges continuing support from the Max-Planck-Gesellschaft and the hospitality of the Academia Sinica Institute of Astronomy and Astrophysics in Taipei, where part of this work was conducted while he was a Visiting Scholar. EP is supported by a Juan de la Cierva fellowship (JDC2022-049957-I). LS was supported by the Max Planck Partner Group Program funded by the Max Planck Society. MAM is supported by the NANOGrav Physics Frontiers Center (NSF Award Number 2020265). JN is supported by the National Natural Science Foundation of China (NSFC, No. 12503055) and the Postdoctoral Fellowship Program of CPSF under Grant Number GZB20250737. JY was supported by the National Science Foundation of Xinjiang Uygur Autonomous Region (2022D01D85), the Major Science and Technology Program of Xinjiang Uygur Autonomous Region (2022A03013-2), the Tianchi Talent project, and the CAS Project for Young Scientists in Basic Research (YSBR-063), the Tianshan talents program (2023TSYCTD0013), and the Chinese Academy of Sciences (CAS) “Light of West China” Program (No. xbzg-zdsys-202410 and No. 2022-XBQNXZ-015). This work made use of the data from FAST (Five-hundred-meter Aperture Spherical radio Telescope). FAST is a Chinese national mega-science facility, operated by National Astronomical Observatories, Chinese Academy of Sciences. While taking data for this work, the Arecibo Observatory was operated by SRI International under a cooperative agreement with the National Science Foundation (NSF; AST-1100968), and in alliance with Ana G. Méndez-Universidad Metropolitana, and the Universities Space Research Association. The National Radio Astronomy Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities. We thank the staff of the Arecibo Observatory for all the help during this project, especially Arun Venkataraman.

References

Appendix A: Observed dispersion measure variations and their power spectrum

Dispersion measure variation is a common phenomenon in pulsar astronomy, starting from the first detection of such epoch-dependent variations in the Crab Pulsar (Rankin & Roberts 1971). The solar wind can contribute to DM variation, which can be corrected, to first order, by TEMPO by assuming a constant electron density of the Solar wind at 1 au from the Sun; this is by default set at 10 cm⁻³. Due to the elliptical latitude of 41.25°, the impact of the solar wind could be negligible. The more significant contribution could be the variation of the large DM of 93.9281(9) pc cm⁻³ of PSR J1946+2052.

Fig. A.1. Measurements of DM structure functions as a function of time lag are displayed in black circles. The black diamond represents the DM structure function when the time lag equals the time scale of DISS (τd ≈ 3 minutes), which is included in the fitting. The red and blue dashed lines are the fitting result with Dϕ(τ)∝τβ − 2 and fixed β of 11/3, respectively. Only points to the left of the grey line are involved in the fitting due to the possible bias introduced by the sudden decrease.

As mentioned previously, these variations are due to the fact that the motion of the pulsar and the Earth changes the path followed by the radio waves from the pulsar to the Earth; these changing paths imply that, at different times, the radio waves go through slightly different regions of the IISM. The resulting stochastic variation of the DM is an interesting probe of the spatial structure of the IISM.

One way to characterise this stochastic variation is by using the spatial power spectrum, which is a function of power-law wavenumber(Rickett 1990):

$$\begin{array}{c}\hfill P(q)=C_n^2{q}^{-\beta };2\pi /l_{\mathrm{o}}<q<2\pi /l_{\mathrm{i}},\end{array}$$(A.1)

where $C_n^2$ is the scaling factor of the power spectrum, β is the power-law index, l_o and l_i are the outer and inner scales of the turbulence. This equation assumes that the electron-density fluctuations within the interstellar medium are isotropic and the magnitude of the wavenumber determines the spectrum (Lam et al. 2016). Then the power spectrum can be estimated by the phase structure function (SF) with

$$\begin{array}{c}\hfill D_{\varphi }(\tau )=⟨[\varphi (t+\tau )-\varphi {(t)}^2]⟩,\end{array}$$(A.2)

where ϕ(t) is the geometrical phase delay that indicates the spatial variations in the electron density along the line of sight between the pulsar and the observer in a specific time, τ is the time lag between phases, and the angle brackets denote an ensemble average. If the spectrum index satisfies 2 < β < 4 and the spatial scale is between outer and inner of the IISM that governs the power spectrum of electron density fluctuations, the phase SF can be given by D_ϕ(τ) = (τ/τ_d)^β − 2 (Rickett 1977), where τ_d is the time scale of the diffractive interstellar scintillation (DISS). Based on this, we can also define DM SF to describe DM variations:

$$\begin{array}{c}\hfill D_{\mathrm{DM}}(\tau )=⟨[\mathrm{DM}(t+\tau )-\mathrm{DM}{(t)}^2]⟩,\end{array}$$(A.3)

where the τ here is the time lag between the DM measurements. The DM SF can be related to the phase SF with

$$\begin{array}{c}\hfill D_{\varphi }(\tau )={\left(\frac{2\pi C}{f}\right)}^2{D}_{\mathrm{DM}}(\tau ),\end{array}$$(A.4)

where C = 4.148808 × 10³ MHz² pc⁻¹ cm³ and f is the frequency of the observation in MHz.

The observed DM variations are shown in Fig. 3. Given the systematic DM offset between DM measurements from Arecibo and FAST, we only use FAST measurements to analyse the DM variation. We display the result of the DM SF in Fig. A.1, where we included in the fit the DM SF when the time lag equals the time scale of DISS. This timescale (τ_d) is estimated to be ∼3 minutes by using the dynamic spectrum of the observation with the highest signal-to-noise ratio. The value of the DM SF for τ = τ_d is then calculated with D_ϕ(τ_d) = 1. In addition, Fig. A.1 shows the power-law fit to the DM SF. The resulting spectral index of the SF is 3.69 ± 0.05, which is consistent with the typical Kolmogorov spectrum with an index of 11/3 (∼3.667), indicating that the interstellar medium between the Earth and the pulsar is incompressible, homogeneous and isotropic.

Appendix B: The moment of inertia of the pulsar

To calculate the Lense-Thirring contribution to the periastron advance caused by the fast-rotating pulsar, one needs to know the MOI of PSR J1946+2052, I_p. Given the rather precise value for the pulsar mass m_p (see Tab. 1), the MOI is known with high precision once an EOS has been chosen to describe the density-pressure relation for neutron star matter. Unfortunately, our knowledge of the properties of matter at supranuclear densities is still limited, leading to a considerable uncertainty in our knowledge of the true EOS for neutron stars. This is reflected in a corresponding inaccuracy in our knowledge of the MOI of a neutron star with given mass. The purpose of this section is to obtain a probability distribution for I_p that takes into account this imperfect knowledge of the density structure of PSR J1946+2052.

In a first step, we use the piecewise polytropic approach of Read et al. (2009) and calculate the MOI of PSR J1946+2052 for all EOSs in Table III of Read et al. (2009) and Table V of Kumar & Landry (2019) that have a maximum mass of at least 1.96 M_⊙. A maximum mass below that value is excluded with 95% confidence by the mass measurement of Fonseca et al. (2021).¹⁰ Figure B.1 shows this result for neutron stars with radii in the range of 10 km to 14 km, where I_p can be approximated with sufficient precision as a linear function of the radius of PSR J1946+2052 R_p. Note that on both ends the R_p-range of Fig. B.1 stretches well beyond the range of EOSs allowed by a combination of a large set of constraints by nuclear physics and multi-messenger astrophysics, as given in Koehn et al. (2024).

An important quantity, when presenting constraints on EOSs is the radius of a “canonical” neutron star of 1.4 M_⊙. Koehn et al. (2024) find $R_{1.4}=12.{27}_{-0.94}^{+0.83}\mathrm{km}$ with 95% credibility as their conservative result. To a good approximation, we can adopt these limits for PSR J1946+2052, since for most EOSs capable of supporting a ∼2 M_⊙ neutron star the radius changes very little when the neutron star mass is lowered from 1.4 M_⊙ to 1.284 M_⊙ (see e.g. Fig. 2 in Lattimer & Prakash (2001)). For simplicity, we now assume a Gaussian probability distribution for R_p with the same 2-σ bounds as R_1.4 in Koehn et al. (2024). Using the linear relation between R_p and I_p of Fig. B.1 (red line), this probability distribution can be converted to a Gaussian probability distribution for the MOI of PSR J1946+2052, which is (1-σ error)

$$\begin{array}{c}\hfill I_{\mathrm{p}}=(1.31\pm 0.08)\times {10}^{45}\mathrm{g}c{m}^2.\end{array}$$(B.1)

The ∼6% uncertainty for the MOI of the pulsar leads to the same fractional error for the calculated ${\dot{\omega }}_{\mathrm{LT}}$. Given that ${\dot{\omega }}_{\mathrm{LT}}$ is only three times the 1−σ error of the observed periastron advance ($\delta \dot{\omega }$), this error is practically negligible for the considerations in the main text, in particular Sec. 4.4. In retrospect, this also justifies some of the idealised approaches in this section.

Fig. B.1. Moment of inertia of PSR J1946+2052 (mass = 1.284 M⊙) as a function of the pulsar’s radius, for 50 EOSs, which all have a maximum neutron-star mass exceeding 1.96 M⊙. The red line shows the best linear fit to the data points. The grey range indicates the conservative 95% credibility interval of R1.4 given by Koehn et al. (2024). See text for more details.

All Tables

All Figures

Fig. 1. Frequency-averaged template used to fit ToAs and some examples of the observed pulse profile (Stokes I) displayed in red and black solid lines, respectively. We generated the template with the pulse profile of 2019-03-29 and aligned each profile with the centre of the two Gaussian functions that we used to fit the pulse, indicated by the blue dashed lines. The less severe profile evolution in the main pulse can be seen in this figure compared to that in the interpulse, which makes it reasonable to fit ToAs only with the main pulse. One can also notice that the separation between the main pulse and the interpulse is increasing over time.

In the text

Fig. 2. Time-of-arrival offsets derived from the standard template of the first observation and the integrated pulse profile displayed with blue points. The pulse width (the unit is the same as the ToA offset) of each integrated pulse profile is plotted with red points. The correlation between the ToA offsets and pulse widths indicates the strong influence introduced by the profile evolution on measuring the ToAs. The blue solid line is the linear fit between MJD and ToA offsets, which is used to correct ToAs.

In the text

Fig. 3. Dispersion measure variations derived from the DMX model with a time bin of 1 day. Panels (a) and (b) represent the DM variation before and after the DM correction, respectively. DM measurements from Arecibo and GBT are represented by black points, and those from FAST are represented by red points. The red dashed line indicates the final measurement of the DM, which is 93.9281 pc cm−3. We display the 10-order DM derivative fit in panel (b) with the blue dashed line.

In the text

	Fig. 4. Reduced χ2 and χ2 with different numbers of DM derivatives shown with black solid circles and red solid stars, respectively. The ten-order DM derivative fit generates the lowest χ2 and reduced χ2. Using a higher-order DM derivative will overfit, as indicated by the larger χ2 and reduced χ2 after ten DM derivatives.
In the text

Fig. 5. Residuals obtained using the DDFWHE timing solution in Table 1. The residuals in blue, orange, green and red are derived from L-band/PUPPI data, P-band/PUPPI data, single GBT observation and FAST data. Top: Residuals as a function of time. Bottom: Residuals as a function of orbital phase. The post-fit residuals’ root mean square is consistent with the ToA uncertainties, and the RMS of the residuals from Arecibo L-band, Arecibo P-band, GBT and FAST are 66.088 μs, 87.709 μs, 138.936 μs and 13.581 μs, respectively. No unmodelled trends are seen in the ToA residuals, indicating that, within measurement uncertainty, the DDFWHE timing solution provides an adequate description of the timing of the system.

In the text

Fig. 6. Contributions of Ṗb from differential Galactic acceleration and the Shklovskii effect as a function of the distance to the pulsar (see text for details). The orange and blue areas indicate the 1-σ confidence of these two effects. The solid lines represent the nominal values, and the dotted curves are the ±1-σ uncertainties. The total external contribution of Ṗb is shown by the corresponding red area and curves. The uncertainties come from the uncertainty of R0, Θ0 and the proper motion. The dashed red and black lines are the distances derived from YMW16 and NE2001. The grey area displays ΔṖb = Ṗb,obs − Ṗb,GR. This quantity, estimated in detail in Sect. 4.1, represents a measurement of the total external contribution that assumes the validity of GR; its width represents its uncertainty, which is dominated by the error of the observed Ṗb.

In the text

Fig. 7. Mass-mass diagram for PSR J1946+2052. The regions consistent with the measured ${\dot{P}}_{\mathrm{b}},\dot{\omega },\gamma ,h_3,\varsigma$, and their 1-σ uncertainties are displayed in red, green, blue, purple, and orange. Note that Ṗb in this diagram has been corrected by removing the external contributions. The grey area is excluded by the mass function and sin i ≤ 1. The inset is an expanded view, showing in more detail the intersection of all PK parameters and the self-consistency of GR. The $\dot{\omega }$-line is based on the leading-order Equation (4). The (small) modifications related to the next-to-leading order contributions to the periastron advance are discussed in detail in Sect. 4.4.

In the text

Fig. 8. Constraints on the orbital inclination angle, masses of the pulsar and the companion, and the Ṗb deviation from the measurement and the prediction of GR (ΔṖb). The black, grey, and light grey contours represent the 1-σ, 2-σ, and 3-σ confidence levels of these parameters. We show the probability density of each parameter on the upper and right edges, where the grey area and solid grey line represent the 1-σ confidence area and the nominal values; the dashed lines at the edges represent the ±1-σ limits.

In the text

Fig. 9. Cumulative shift of periastron time relative to the predictions of a constant orbit caused by the decay of the orbital period. In the upper panel, the black points are the measurements of the periastron time shift, and the red solid line is the GR prediction, which is the sum of the GW damping calculated from GR and the calculated contribution from the Galactic acceleration and the Shklovskii effect. We display the residual between the measurements and the GR prediction in the lower panel.

In the text

Fig. 10. Enlarged mass-mass diagram for PSR J1946+2052. To demonstrate the influence introduced by the LT effect, we display the total mass constraint for ${\dot{\omega }}_0$ with black dashed lines, where we ignore the LT contribution. The black solid lines represent $\dot{\omega }$ with the full contribution, including the 2PN and LT effect. The significant deviation between ${\dot{\omega }}_0$ and $\dot{\omega }$ indicates that the LT contribution significantly affects our measurement of the total mass. The magnitude of this shift depends on the MOI of the pulsar. Owing to the much larger uncertainty of Ṗb, the total mass cannot be measured independently with sufficient precision to determine the MOI.

In the text

Fig. A.1. Measurements of DM structure functions as a function of time lag are displayed in black circles. The black diamond represents the DM structure function when the time lag equals the time scale of DISS (τd ≈ 3 minutes), which is included in the fitting. The red and blue dashed lines are the fitting result with Dϕ(τ)∝τβ − 2 and fixed β of 11/3, respectively. Only points to the left of the grey line are involved in the fitting due to the possible bias introduced by the sudden decrease.

In the text

	Fig. B.1. Moment of inertia of PSR J1946+2052 (mass = 1.284 M⊙) as a function of the pulsar’s radius, for 50 EOSs, which all have a maximum neutron-star mass exceeding 1.96 M⊙. The red line shows the best linear fit to the data points. The grey range indicates the conservative 95% credibility interval of R1.4 given by Koehn et al. (2024). See text for more details.
In the text