© The Authors 2024

1 Introduction

General relativity (GR) proved to be the most successful theory of gravity for more than a century. Up to now, it has passed all experimental tests exceptionally well. The weak field regime has been verified by the Solar System experiments (Will 2014), whereas strong field effects have been tested especially well by the timing of binary pulsars (Wex & Kramer 2020). Among other things, pulsar tests have enabled precise tests of the radiative properties of gravity and the strong equivalence principle (SEP). Furthermore, the large-scale behaviour of gravity (low spacetime curvature and temporal variation) has been tested in cosmological observations (Clifton et al. 2012). Finally, in the last few years, ground-based gravitational wave (GW) detectors observed GWs from coalescing binary black holes (BHs) (Abbott et al. 2016, 2017c, 2019a, 2021b), binary neutron stars (NSs) (Abbott et al. 2017a,b, 2019b), and also BH-NS binaries (Abbott et al. 2021a), probing the hitherto unexplored highly dynamic, strong-field regime of gravity. Thus far, GR has accounted for all observed effects both in weak and strong gravitational fields, in the quasi-static and the highly dynamical regime, and on small as well as on large scales¹.

Despite such successes, there are still convincing reasons to investigate modified theories of gravity (Berti et al. 2015). GR describes a gravitational interaction by a single massless, spin-two tensor field 𝑔_µν, known as the metric of spacetime. Some of the best motivated alternative gravity theories, the scalar-tensor gravity (STG) theories, incorporate additional scalar degrees of freedom to mediate gravity. Such scalar fields arise naturally in higher dimensional theories in their low energy limits, for example in the old Kaluza-Klein theory (Kaluza 1921; Klein 1926) and string theories (Fujii & Maeda 2007). Scalar-tensor gravity theories are also motivated by cosmological questions of inflation, dark energy, and even the thus far undiscovered theory of quantum gravity (Clifton et al. 2012).

A scalar-tensor gravity theory which for many years was considered as the only natural competitor of GR is known as Jordan–Fierz–Brans–Dicke (JFBD) gravity (Jordan 1955, 1959; Fierz 1956; Brans & Dicke 1961)². JFBD gravity is a metric theory (matter and non-gravitational fields couple only to one specific spacetime metric) and it therefore fulfils the Einstein equivalence principle by design (Will 1993; Damour 2012). However, an additional scalar field is non-minimally coupled to matter via the choice of a special conformal coupling function; thus, SEP is violated (Will 2014). This coupling function depends on a single parameter ω_BD, which by now is tightly constrained, in particular with the pulsar in the stellar triple system (Voisin et al. 2020). Bergmann (1968) and Wagoner (1970) presented the most general STG theory with one scalar field which, while in action, is quadratic at most. The most general mono-scalar–tensor theory with second-order field equations is Horndeski gravity (Horndeski 1974). In 1992, Damour & Esposito-Farèse (1992b) presented a generic class of STG theories with an arbitrary number of scalar fields. In the scope of this paper, we mainly focus on Damour–Esposito–Farèse (DEF) gravity (Damour & Esposito-Farèse 1993). Damour–Esposito–Farèse gravity is a Bergman-Wagoner theory which has a massless scalar field and a quadratic coupling function with two arbitrary parameters.

Apart from being well-motivated, DEF gravity predicts interesting effects which are not present in weak fields of the Solar System but could be tested by employing pulsar astronomy. Neutron stars have strong gravitational fields and large gravitational binding energies due to their high compactness. Because of this, in STG theories they can have large gravitational form factors (also known as scalar charges), unlike weakly self-gravitating objects such as normal stars and white dwarfs (WDs). Consequently, while STG theories generally predict the emission of dipolar and monopolar GWs, which increase the rate of orbital period decay relative to the GR expectation, this is expected to be especially strong in asymmetric binary systems containing a NS and a WD, due to the significant differences in their compactness and scalar charges.

Furthermore, Damour & Esposito-Farèse (1993, 1996) found a specific phenomenon in NSs called ‘spontaneous scalarisation’. This fully non-perturbative effect excites the scalar field above the cosmological background value in NSs, allowing scalar charges of order unity, even if there is no deviation from GR in the weak-field regime. This could produce, for binary systems containing NSs with specific masses, a highly enhanced rate of dipolar GW emission.

Some of the tightest tests of gravity with strongly self-gravitating objects are provided by the high precision timing of binary pulsars (Stairs 2003; Wex 2014; Shao & Wex 2016; Wex & Kramer 2020; Kramer et al. 2021). These tests are performed in a quasi-stationary strong-field gravity regime, where the gravitational field is strong near and inside NSs and produces high curvature, whereas the velocities are small compared to the speed of light υ/c ~ 10⁻³. The first timing model (BT) was introduced by Blandford & Teukolsky (1976) and it allowed for information from the timing of binary pulsars to be extracted. An extended model covering all relativistic effects in the dynamics of a binary system up to the first post-Newtonian level was proposed later by Damour & Deruelle (1985, 1986) (DD).

The DD model made it possible to perform self-consistent tests of gravity by means of the parametrised post-Keplerian (PPK) formalism (see, for instance, applications by Taylor & Weisberg 1989). In this generic framework, theory-independent Keplerian and post-Keplerian (PK) timing parameters, which quantify the relativistic motion of the pulsar and the propagation of its radio signals, are fitted to the observational timing data. Apart from the well-measured Keplerian parameters of the orbit, each PK parameter only depends on the masses of the pulsar and its companion. As soon as the expressions for the PK parameters in a particular gravity theory are known, these two masses can be derived from the PK parameters once at least two of them are known. If more than two PK parameters become available, then a test of the consistency of that gravity theory can be made. Later, Damour & Taylor (1992) expanded the DD model and introduced pulse structure parameters. They also showed that the DD model can be applied to a large set of fully conservative theories of gravity and presented phenomenological PK expressions in a framework of generic boost-invariant theories (the modified Einstein-Infeld-Hoffmann formalism).

However, such an independent measurement of PK parameters means that any correlations between them inadvertently reduce the sensitivity of their measurements. Taylor & Weisberg (1989) solved this problem for GR with the introduction of the DDGR model, where the GR expressions for PK parameters are incorporated in the model, and two masses are directly fitted to the observational timing data. This mitigates or breaks any correlations between the PK parameters.

In this work we present a new pulsar timing model – the ‘DDSTG’ model which is called after Damour, Deruelle and STG. The idea is very similar to that of the DDGR model: fitting directly for the masses of the components from pulsar timing data, but instead of GR using a particular member of the two-parametric DEF gravity. Apart from the theoretical predictions for all PK parameters within that theory, the model uses pre-calculated grids of physical parameters (gravitational form factors) of NSs in that particular gravity theory. The direct fit of two masses of the companions reduces the model to a minimal set of parameters and naturally solves all the issues of possible correlations between the observed parameters. This feature, achieved by the construction, makes the DDSTG model superior to the traditional methods, for example the ‘PK method’ based on the measurement of PK parameters.

As an application of the new DDSTG model, we use the most recently published timing data for PSR J2222–0137 from Cognard et al. (2017) and Guo et al. (2021b). This binary pulsar is interesting for several reasons: It is one of the closest binary pulsars known and has a good timing precision. A precise measurement of the variation of the orbital period (Ṗ_b) and the high asymmetry in compactness between the pulsar and its companion allows for tight constraints on the emission of dipolar GWs. Moreover, the pulsar’s mass m_p = 1.831(10) M_⊙ (Guo et al. 2021b) lies in the range (m_p ≳ 1.5 M_⊙) subject to the spontaneous scalarisation effect; the non-detection of dipolar GW emission in this system has introduced strong constraints on the existence of that phenomenon, at least within the DEF framework (Zhao et al. 2022).

We organise the paper as follows: in Sect. 2, we introduce the basics of DEF gravity and present our machinery preparations necessary for the implementation of the new model. In Sect. 3, we discuss the standard timing models used thus far and introduce our new DDSTG model for testing STG theories. The subsequent sections are devoted to the applications of our new timing model. In Sect. 4, we provide a brief description of the Guo et al. (2021b) dataset on PSR J2222–0137 and apply the DDSTG model to it, discuss the limits on DEF gravity, and demonstrate that our new method does indeed provide tighter and more reliable constraints than previous methods. In Sect. 5, based on simulated datasets for the largest radio observatories, we discuss what limits on these theories might be achieved in the near future with continued timing of PSR J2222–0137. In Sect. 6, we show possible future constraints from reasonable pulsar-black hole (PSR-BH) binary systems. We also briefly discuss the possible origin of PSR-BH systems and how it affects the test if such a system is located in a globular cluster (GC). Section 7 is devoted to a discussion of other alternative gravity theories and possibilities of the DDSTG model extensions beyond DEF gravity. In Sect. 8 we discuss the effect of the time-varying gravitational constant as a perspective future extension of the DDSTG model. In Sect. 9 we summarise our results.

2 Damour–Esposito–Farèse gravity

This paper investigates STG theories, natural alternatives to GR. Specifically, we mainly focus on DEF gravity. In the following we summarise aspects of that class of gravity theories that are relevant for this paper.

2.1 Theoretical aspects

Damour–Esposito–Farèse gravity is a STG theory with one long range, massless scalar field ϕ non-minimally coupled to the curvature scalar. The theory is fully described by the action, which is most simply presented in the Einstein frame (Damour & Esposito-Farèse 1993, 1996) $$S=\frac{c^4}{16\pi G_{*}}{\displaystyle \int \frac{{\text{d}}^{4}x}{c}}{g}_{*}^{1/2}\left(R_{*}-2g_{*}^{\mu v}{\partial }_{\mu }\phi {\partial }_v\phi \right)+S_m\left[{\psi }_m;A^2(\phi )g_{\mu v}^{*}\right],$$(1)

where G_* is the bare gravitational coupling constant, 𝑔_* and R_* are the determinant and the Ricci scalar associated with the Einstein metric ($g_{\mu v}^{*}$). The last term in Eq. (1) describes an action associated with any matter fields (ψ_m). A(ϕ) is the coupling function and takes a specific exponential form in DEF gravity: $A(\phi )=\exp \left[{\alpha }_0\left(\phi -{\phi }_0\right)+\frac{1}{2}{\beta }_0{\left(\phi -{\phi }_0\right)}^2\right]$. Introduced quantities {α₀,β₀} are two arbitrary parameters of the theory and ϕ₀ is the scalar field at spatial infinity. The coupling function A(ϕ) plays the role of a conformal factor connecting the ‘physical metric’ ${\tilde{g}}_{\mu \nu }=A^2(\phi )g_{\mu \nu }^{*}$ with the Einstein metric. The metric ${\tilde{g}}_{\mu v}$ is the one measured by laboratory clocks and rods due to the universal coupling of matter to this metric. The bare gravitational constant is simply related to the one measured in a Cavendish experiment by $G_{\text{Cav}}=G_{*}\left(1+{\alpha }_0^2\right)$. GR corresponds to vanishing coupling parameters, that is α₀ = β₀ = 0, and JFBD gravity is recovered by setting β₀ = 0 while keeping α₀ ≠ 0.

The field equations are derived by variation of the action in Eq. (1) and are most simply formulated in terms of the pure-spin field variables $\left(g_{\mu v}^{*},\phi \right)$ in the Einstein frame: $$R_{\mu \nu }^{*}=2{\partial }_{\mu }\phi {\partial }_{\nu }\phi +\frac{8\pi G_{*}}{c^4}\left(T_{\mu \nu }^{*}-\frac{1}{2}{T}^{*}{g}_{\mu \nu }^{*}\right),$$(2) $${\square }_{g^{*}}\phi =-\frac{8\pi G_{*}}{c^4}\alpha (\phi )T_{*},$$(3)

with a material stress-energy tensor $T_{*}^{\mu \nu }=2c{g}_{*}^{-1/2}\delta S_m/\delta g_{\mu v}^{*}$ and α(ϕ) = ∂ In A(ϕ)/∂ϕ, which measures the coupling strength between the scalar field and matter.

All weak-field deviations from GR may be described in terms of the asymptotic behaviour of α(ϕ) at spatial infinity. The theory parameters α₀ = α(ϕ₀), β₀ = β(ϕ₀) can be interpreted as asymptotic values of the function α(ϕ) and its derivative β(ϕ) = ∂α(ϕ)/∂ϕ.

In this work, without loss of generality, we use the formulation of DEF gravity with ϕ₀ = 0, two parameters {α₀,β₀} appearing in the form of A(ϕ) and A(ϕ₀) = 1. There is an equivalent formulation in terms of parameters {β₀, ϕ₀}, where A(ϕ) = exp(β₀ϕ²/2) and ϕ₀ = α₀/β₀. But the latter formulation is not well-determined for the JFBD limit β₀ = 0 (Damour 2007).

The parameterised post-Newtonian (PPN) framework allows for a wide range of metric theories of gravity to be described in their weak field approximation by ten independent PPN parameters. Damour–Esposito-Farèse gravity is a fully conservative theory and has only two PPN parameters that deviate from their GR values, which are also called Eddington parameters (Will 1993): $${\gamma }_{\text{Edd}}=1-\frac{2{\alpha }_0^2}{1+{\alpha }_0^2},$$(4) $${\beta }_{\text{Edd}}=1+\frac{{\beta }_0{\alpha }_0^2}{2{\left(1+{\alpha }_0^2\right)}^2}.$$(5)

The PPN parameter γ_Edd measures the spatial curvature induced by unit rest mass and β_Edd is a measure of the amount of non-linearity in the superposition law of gravity. Solar System experiments put limits on these parameters. The |γ_Edd − 1| ≲ 2.3 × 10⁻⁵ limit comes from the Shapiro time delay observed by the Cassini spacecraft (Bertotti et al. 2003; Will 2014). The β_Edd − 1| ≲ 8 × 10⁻⁵ limit comes from observations of the perihelion shift of Mercury (Will 2014). The Cassini experiment yields a direct limit on ${\alpha }_0^2\lesssim 1.15\times {10}^{-5}$, whereas β₀ remains unconstrained from weak-field Solar System experiments, since α₀ can be arbitrarily small.

2.2 Scalarisation of neutron stars in DEF gravity

To place firm limits on β₀ parameter, we need to explore the strong-field effects of DEF gravity. One can expect deviations from GR in the presence of highly compact massive objects, in particular for NSs, due to the scalar field sourced by the strong internal gravitational field of the star. When placing limits on β₀ it is important that for certain areas in the DEF gravity parameter space one has a significant growth of the scalar field in the NS interior, even if the deviations for weakly self-gravitating bodies are very small, or even zero (Damour & Esposito-Farèse 1992a, 1993, 1996). The origin of this spontaneous scalarisation lies in the tachyonic instability of the field equations happening when a NS reaches a sufficiently high compactness. The mechanism of this instability is similar to what happens in ferromagnets during spontaneous magnetisation (Damour & Esposito-Farèse 1996). It is a fully non-linear non-perturbative effect, thus it requires accurate techniques to be properly analysed.

For testing DEF gravity, it is crucial to introduce three gravitational form factors appearing in the expressions for PK parameters $${\alpha }_A={\frac{\partial \ln m_A}{\partial {\phi }_0}|}_{{\tilde{m}}_A},{\beta }_A={\frac{\partial {\alpha }_A}{\partial {\phi }_0}|}_{{\tilde{m}}_A},k_A=-{\frac{\partial \ln I_A}{\partial {\phi }_0}|}_{{\tilde{m}}_A},$$(6)

where all the derivatives are taken for a fixed baryonic mass (${\overline{m}}_A$). The important quantity α_A counts the effective couple strength between the ambient scalar field and a NS. The subscript A refers to the star A. α_A is a strong-field counterpart to the parameter α₀, which is its weak field approximation. If a NS mass exceeds the critical mass of m_crit then α_NS can become suddenly of order of unity α_NS ~ O(1) due to the scalarisation effect, while its weak field counterpart stays close to zero α₀ ~ 0. The value of the critical mass depends on β₀ and the chosen equation of state (EOS). The second parameter β_A reflects a change of the scalar charge when the asymptotic scalar field changes. The derivative β_A can also obtain very large values in the transition scalarisation region, where α_A changes fast with the change of a mass. The last quantity k_A describes the field dependence of the moment of inertia I_A and can become important in the eccentric systems.

2.3 Calculating grids of NSs

To test DEF gravity in precision timing experiments of binary pulsars, one must know the exact properties of NSs in that specific class of alternative gravity. In our work, we follow the procedure described in Damour & Esposito-Farèse (1996) for calculating NSs in DEF gravity. We assume the slowly rotating NS approximation with a stationary axisymmetric metric. Assuming this axial symmetry and neglecting terms of rotational velocity squared (or higher), the Einstein equations can be written as a system of eight first-order ordinary differential equations (ODEs) for a set of radial functions (see the representation by Damour & Esposito-Farèse 1996). The equations are complemented by appropriate boundary conditions placed at spatial infinity and the NS centre. Boundary conditions at spatial infinity and thus properties of NSs depend on the values of the theory parameters {α₀,β₀}.

For solving NS structures, we developed a modular program in the JULIA language (Bezanson et al. 2017). The program enables us to calculate both single NSs for selected parameters and grids of NSs for many desired parameters. The calculations are performed accurately and efficiently using the methods of parallel programming. Currently, DEF gravity is implemented as the chosen theory of gravity, but our program allows one to use any coupling function A(ϕ) to extend our tests beyond DEF gravity parametrisation.

The internal structure and NSs parameters significantly depend on the selected EOS. The strong dependence of the scalarisation phenomenon on the choice of a EOS was extensively analysed by Shao et al. (2017) and Zhao et al. (2022) in order to put constraints on this non-linear phenomenon from radio pulsars. In our work, we use a piece-wise polytropic approximation for the EOSs following the procedure in Read et al. (2009). We select 11 EOSs that allow us to cover a range from soft to stiff while still being in agreement with the maximum mass observed for a NS (for more information about selected EOSs we refer to Appendix A.2). For tests in the present paper, we select the MPA1 high-density EOS (Müther et al. 1987) in its piece-wise polytropic form (Read et al. 2009) with three polytropic pieces. MPA1 is a stiff EOS with a maximum mass of a NS equal to 2.46 M_⊙. For low densities, we use an analytic form of the SLy EOS proposed by Douchin & Haensel (2001) and approximated by four polytropic pieces (Read et al. 2009). In general stiff EOSs (e.g. MPA1) produce more conservative limits in most parts of the DEF gravity parameter space.

In our program, we utilise automatic differentiation (AutoDiff) technique (see e.g. Margossian 2018) for calculating derivatives of NSs quantities. AutoDiff is a powerful alternative to numerical differentiation deprived of numerical errors (for detailed explanation see Appendix A.1). AutoDiff algorithm utilises exact expressions of derivatives for elementary functions and a chain rule to calculate complex derivatives with the working machine-precision. Thus our program with applied AutoDiff technique calculates gravitational form factors {α_A,β_A, k_A}, presented as derivatives in Eq. (6), with very high numerical precision.

Using the developed program, we calculate gravitational form factors and masses for an extensive grid of DEF parameters and central pressures {α₀, β₀, p_c} assuming a specific EOS. Our grids cover ranges of α₀ ∈ [−10⁻¹, −10⁻⁵], β₀ ∈ [−6, +10], p_c ∈ [10³⁴, 10³⁶] dyn cm⁻² and are sampled linearly for log α₀, log p_c and β₀. The size of each grid is {101, 351, 121} points respectively. More technical details about the procedure of calculating NS properties, AutoDiff technique, and the structure of pre-calculated grids can be found in Appendix A. We stick to pre-calculating grids because it is numerically costly to calculate NS structures for a particular mass on the fly. Further one can perform interpolation over these pre-calculated grids or methods beyond interpolation, e.g reduced-order models in Guo et al. (2021a). The achieved accuracy of grid values and further interpolation over the grid is sufficient for our purposes.

The accuracy of the calculations can be checked by comparison of the scalar charge (α_A) calculated as the derivative from Eq. (6) α_{A, deriv} and its direct value from the asymptotic behaviour of the external scalar field obtained after integration α_{A, asympt} (Damour & Esposito-Farèse 1993). The relative accuracy expressed by δ_rel = |1 − α_{A, deriv}/α_{A, asympt}| lies in the range δ_re1 ~ 10⁻⁷ − 10⁻¹³ for our calculations. The calculated scalar charges are significantly more accurate compared to previous works, for instance, Anderson & Yunes (2019) with typical accuracy of δ_rel ~ 10⁻¹ −10⁻⁵. However, it is essential to mention that the accuracy of several percent is already sufficient for testing DEF gravity with binary pulsar timing.

Besides, we also find that the NS structure appears to be very sensitive to internal thermodynamic consistency. The piece-wise polytropic EOSs fulfil the thermodynamic consistency by their definition. However, often one uses tabulated EOSs instead, which must be interpolated. Unfortunately, a commonly used log-log interpolation of tabulated EOSs can lead to significant numerical errors in the scalar derivatives of Eq. (6) because of its thermodynamic inconsistency³. For this reason, we use piece-wise polytropic approximations instead of tabulated EOSs, which is absolutely sufficient for conducting our tests and for exploring the EOS dependence of pulsar constraints on DEF gravity.

3 A new timing model for scalar-tensor gravity

To place constraints on a gravity theory from the observations of a binary pulsar, one must know its predictions for the binary system’s motion and the propagation of the electromagnetic signal in the curved spacetime of the binary system. The timing formula is the tool to capture the relativistic effects predicted by the theory from a sequence of pulse arrival times on Earth. The timing formula relates the observed (topocentric) time of arrival (TOA) of the pulse and its time of emission.

3.1 Binary pulsar timing models

The first timing model ‘BT’ was introduced by Blandford & Teukolsky (1976) in order to describe the timing of the first binary pulsar, PSR B1913+16, which had been discovered earlier by Hulse & Taylor (1975). The BT model assumes that the pulsar and its companion follow a Keplerian motion with additional secular changes in the Keplerian parameters of the orbit. utilised Keplerian parameters are the orbital period (P_b), the epoch of periastron passage (T₀), the orbital eccentricity (e), the longitude of periastron (ω), and the projected semi-major axis of the pulsar’s orbit (x). The model accounts for a combination of a special-relativistic time dilation and a gravitational redshift; this periodic effect is described by an extra parameter called the Einstein delay (γ). The BT model accounts for (linear-in-time) secular changes in P_b, e, ω and x, introducing PK parameters: rates of change of the orbital period (Ṗ_b), eccentricity (ė), longitude of periastron ($\dot{\omega }$) and projected semi-major axis (ẋ). These secular changes can be caused by both relativistic and astrometric effects (Damour & Taylor 1992; Lorimer & Kramer 2004).

Later, Damour & Deruelle (1986) proved that all of the independent O(υ²/c²) timing effects could be described in a simple mathematical way for a wide range of alternative gravity theories. They developed a phenomenological (i.e. theory independent) timing model based on the full first post-Newtonian description of the two-body problem, which we refer to as the ‘DD’ model that uses a quasi-Keplerian solution to the equations of motion of a two-body problem. This model allows working within the PPK approach. Damour & Taylor (1992) then showed that the DD model could be used to constrain a wide range of conservative theories of gravity obeying the modified Einstein–Infeld–Hoffmann (mEIH) framework (see also Will 1993).

For the binary system part, the timing formula of the DD model reads as: $$t_{\text{b}}-t_0=D^{-1}\left[T+{\text{Δ}}_R(T)+{\text{Δ}}_E(T)+{\text{Δ}}_S(T)+{\text{Δ}}_A(T)\right],$$(7)

where t_b is the Solar System barycentric (infinite-frequency) arrival time, t₀ is a chosen reference epoch, T is the pulsar’s proper time, and D is a Doppler factor accounting for the relative radial motion of the centre of mass of the binary system and the Solar System barycentre. The quantities ∆_i in Eq. (7) are different time delays introducing corrections due to internal binary effects.

Splitting the timing formula into a set of different contributions is to some extent a coordinate-dependent concept, however, within the approximations used it is convenient to work with these individual expressions for timing delays. The term ∆_R (T) is called ‘Roemer delay’ and counts for the classical light travel time through the binary. ∆_E (T) is the ‘Einstein delay’ and relates the proper emission time to the coordinate time of emission. ∆_S (T) is the ‘Shapiro delay’ arising from the effect of the gravitational potential of the companion on the propagation of the pulsar signals. The ‘aberration delay’ ∆_A (T) places corrections due to the periodic changes in the direction of pulse emission (as seen in the frame of the rotating pulsar) while the pulsar follows its binary motion.

All the expressions depend on three sets of parameters. Keplerian parameters are $$\left\{p^{\text{K}}\right\}=\left\{P_{\text{b}},T_0,e_0,{\omega }_0,x_0\right\}$$(8)

and remain the same as for the BT model (the subscript 0 means a value at a given epoch). Separately measurable PK parameters are $$\left\{p^{\text{PK}}\right\}=\left\{k,\gamma ,{\dot{P}}_{\text{b}},r,s,{\delta }_{\theta },\dot{e},\dot{x}\right\}.$$(9)

The DD model introduces a periastron-shift parameter k = $\dot{\omega }{P}_{\text{b}}/(2\pi )$, ‘Shapiro shape’ (s = sin i, where i is the inclination angle) and ‘Shapiro range’ (r) parameter for the signal propagation, and a relativistic deformation of the orbit (δ_θ). Not separately measurable PK parameters are $$\left\{q^{\text{PK}}\right\}=\left\{{\delta }_r,A,B,D\right\}$$(10)

including a second parameter for the relativistic deformation (δ_r), two aberration parameters (A and B) and a Doppler factor (D).

The mentioned delays in the framework of the DD model are presented by expressions $$\begin{array}{l}{\text{Δ}}_R(T)=x\sin \omega \left[\cos u-e\left(1+{\delta }_r\right)\right]\\ +x{\left[1-{\text{e}}^2{\left(1+{\delta }_{\theta }\right)}^2\right]}^{1/2}\cos \omega \sin u\text{, }\end{array}$$(11) $${\text{Δ}}_E(T)=\gamma \sin u,$$(12) $$\begin{array}{l}{\text{Δ}}_S(T)=-2r\ln \{1-\text{e}\cos u-s[\sin \omega (\cos u-e)\\ +{\left(1-{\text{e}}^2\right)}^{1/2}\cos \omega \sin u]\},\end{array}$$(13) $$\begin{array}{l}{\text{Δ}}_A(T)=A\left\{\sin \left[\omega +A_e(u)\right]+e\sin \omega \right\}\\ +B\left\{\cos \left[\omega +A_e(u)\right]+e\cos \omega \right\},\end{array}$$(14)

where the secular changes incorporated in the projected semi-major axis x and eccentricity e are given by $$x=x_0+\dot{x}\left(T-T_0\right)$$(15) $$e=e_0+\dot{e}\left(T-T_0\right)\text{. }$$(16)

The true anomaly (A_e) and ω depend on the eccentric anomaly (u) via following the relations $$A_e(u)=2\mathrm{\hspace{0.17em}arctan\hspace{0.17em}}\left[{\left(\frac{1+e}{1-e}\right)}^{1/2}\tan \frac{u}{2}\right],$$(17) $$\omega ={\omega }_0+k{A}_e(u),$$(18)

where u is itself a function of T and defined by solving Kepler’s equation $$u-e\sin u=2\pi \left[\left(\frac{T-T_0}{P_{\text{b}}}\right)-\frac{1}{2}{\dot{P}}_{\text{b}}{\left(\frac{T-T_0}{P_{\text{b}}}\right)}^2\right].$$(19)

The parameters in Eq. (10) cannot be measured separately from the parameters in Eqs. (9) and (8) as shown by Damour & Deruelle (1986). It means that all not separately measurable parameters can be fully absorbed into the change of other parameters. More details about measurable parameters and the connection between observed parameters and the intrinsic parameters can be found in Damour & Taylor (1992). The redefinition of orbital parameters can absorb the Doppler factor (D). Such a redefinition does not affect the tests because it is only a rescaling of physical units, thus D is set to one⁴. The aberration parameters A and B can be absorbed as well by redefinition of T₀, x, e, δ_r and δ_θ.

3.2 A timing model in DEF gravity

The DD model can be applied to a wide range of fully conservative theories of gravity. Damour–Esposito–Farèse gravity is a fully conservative theory (Will 1993), thus we can investigate it in the framework of the DD model.

In different theories of gravity, the functional dependence of the PK parameters, i.e. $$p_i^{\text{PK}}=f_i^{\text{theory }}\left(m_{\text{p}},m_{\text{c}};P_{\text{b}},e_0,x_0;\text{EOS}\right),$$(20)

can differ substantially because of strong-field effects involving a highly compact NS and its companion. The PK parameters are functions of two masses, thus one can perform gravity tests if at least three PK parameters are measurable. In general, if there are n measured PK parameters, one can do n − 2 independent tests for a given gravity theory. A common procedure of making tests using expressions for PK parameters is discussed in detail in Sect. 4, devoted to applications of a new model and comparison between different techniques.

The core part of the new timing model is the predictions for PK parameters (20) in DEF gravity. We use expressions for PK parameters in DEF gravity provided in the literature (Damour & Esposito-Farèse 1992b, 1993, 1996). However, for parameters δ_θ, δ_r, A and B there are only more general expressions. A phenomenological theory-independent description of PK parameters from which we derive their expressions in DEF gravity is presented in Damour & Taylor (1992). The expressions of PK parameters and details are given in Appendix B.

The aberration parameters A and B can also be calculated in DEF gravity, provided we know the system’s geometry. The applications we use in the paper assume a special situation when the pulsar’s spin axis is aligned to the orbital angular momentum of the system. The alignment is a reasonable assumption for a recycled pulsar with a WD companion due to the preceding accretion process during the system’s evolution. In this special case, A is calculated according to Eq. (B.7) assuming angles η = −π/2 and λ = i (Damour & Taylor 1992). The second aberration parameter B = 0 in this special situation. In general, an alignment may not be the case for a particular binary pulsar (e.g. double NS systems), but the new model can also handle this misalignment. Once the system’s geometry is assumed, real A and B values can be calculated for a given epoch without the necessity to redefine other parameters⁵.

The PK parameters from Eq. (20) depend on the orbital parameters of the system and the physical properties of the pulsar. However, we also have to know properties of the companion. The PK expressions (B.1–B.13) depend on gravitational form factors of the companion {α_c, β_c}. The quantity k_c does not appear in the expressions because they take into account only leading order terms, and is therefore of no interest. These quantities, in turn, depend on the type of companion. If the companion is a NS, they are calculated in the same way as for the primary pulsar using Eq. (6). If the companion is a BH, then the gravitational form factors all vanish: $${\alpha }_{\text{BH}}=0,{\beta }_{\text{BH}}=0,k_{\text{BH}}=0.$$(21)

This is a consequence of the ‘no-hair’ theorem; the BH is not scalarised in DEF gravity (Hawking 1972; Damour & Esposito-Farèse 1992a). Section 6 is dedicated to a more thorough discussion on PSR-BH systems and the results which we can obtain from timing of a pulsar in a relativistic orbit with a BH. For a WD companion, the gravitational form factors are approximated by their weak field expressions $${\alpha }_{\text{WD}}\simeq {\alpha }_0,{\beta }_{\text{WD}}\simeq {\beta }_0,k_{\text{WD}}\simeq 0.$$(22)

Such an approximation is sufficient for our purposes because WDs are not compact enough to show the strong field effects present in NSs and BHs. The next order approximation for weakly self-gravitating objects is α_A ≃ α₀(1 − 2s_A), where s_A ≃ G_*m_A/(Rc²) is the sensitivity (R is the object’s radius) and has a typical value of s_WD ~ 10⁻⁵−10⁻³ for a WD (Damour & Esposito-Farèse 1992a, 1993; Will 2018b). Thus the usage of weak field counterparts {α₀, β₀} for a WD companion instead of precisely calculated values is justified.

To summarise, the timing model in DEF gravity is defined by two theory parameters {α₀, β₀}, the chosen EOS for a NS, and the type of the companion among {WD, NS, BH}. It concludes the theoretical description of the DDSTG model.

3.3 DDSTG implementation into TEMPO

Once we obtain the theoretical part of the new DDSTG timing model, we need to apply it to the timing data. We implemented the DDSTG model into one of the commonly used pulsar timing software – the Tempo⁶ program (Nice et al. 2015). The local implementation of the DDSTG model in TEMPO and precalculated grids of masses and gravitational form factors are supplied with the paper⁷. The authors intend to make the DDSTG model a part of the official TEMPO distribution to ensure forward compatibility.

There is a standard procedure of analysing radio pulsar timing data. The preprocessed TOAs are obtained after radio observations of the desired pulsar system. During the procedure, TEMPO reads TOAs, parameters of the binary model, and some coded instructions from supplied files. Then TEMPO fits the selected timing model accounting for the transformation to the Solar System barycentre, pulsar rotation, and its spin down for a chosen binary model.

Specifically for the DDSTG model, a user selects theory parameters {α₀, β₀}, a EOS, and a type of a companion. During the initialisation TEMPO reads pre-calculated three-dimensional (3D) grids of gravitational form factors and NS masses {α_A, β_A, k_A, m_A} (see Sect. 2.3) for the chosen EOS. Each 3D grid depends on {α₀, β₀, p_c} parameters, where p_c is the central pressure. Then gravitational form factors and masses are interpolated for the selected {α₀, β₀} values and saved into a smaller 1D grids (depending only on p_c) in the program memory. The final 1D grids remain unchanged and are used to interpolate gravitational form factors for a particular theory in the mass range during the fitting procedure of TEMPO.

TEMPO with the selected DDSTG model fits two masses: the total mass (m_tot) of the system and the companion’s mass (m_c). Every time when these masses change, the model recalculates the gravitational form factors for the pulsar and the companion. Once the gravitational form factors (Eq. (6)) are known, the model calculates all PK parameters (Eqs. (9) and (10)) using equations given in Appendix B. These calculated PK parameters are used afterwards in the timing formula (7). The details about the DDSTG implementation and the description of model parameters can be found in Appendix C.

In the timing formula (7), the Solar System barycentric arrival time (t_b) is known. However, to obtain the proper pulsar time (T), one has to perform the inversion of the timing formula, i.e. get T as a function of $t_{\text{b}}:T=t_{\text{b}}-\overline{\text{Δ}}\left(t_{\text{b}}\right)$. The original DD model utilises an approximate analytic inversion, which some-times is not accurate enough, for example, for the double pulsar PSR 0737-3039A (Kramer et al. 2021). In contrast, our DDSTG model performs accurate numerical inversion. Equation (7) is solved iteratively for T for each TOA. The discrepancy due to an inaccurate inversion may influence the test for precise timing in the future or PSR-BH systems discussed further below. Nowadays, numerical inversion is considered a new standard and implemented in the DDS model of TEMPO (Kramer et al. 2021) and generally used in PINT (Luo et al. 2021).

The DDSTG model produces the best fit to the data for DEF gravity with a selected set of parameters and a given EOS. The output of TEMPO stays the same as for other binary models, for instance, the calculated χ² value and root mean square error of fit. The χ² statistics may be applied for further tests to compare the results for different theory parameters.

3.4 Advantages of DDSTG

The tests with pulsar data allow for constraints to be placed on the gravity theory parameters. The most simple and common way to achieve this is to compare the measured PK parameters with the theoretical predictions from the theory. PK parameters are obtained by fitting a phenomenological model to the observational data, for instance, the DD model. A phenomenological model estimates PK parameters and their uncertainties.

Within a particular gravity theory, each PK parameter depends on the masses of the pulsar and the companion and thus corresponds to a curve in a mass-mass space. Together with the measurement error, each PK parameter produces a strip in the {m_p, m_c} space. If there are two measured parameters, one may find the intersection area and obtain the estimated mass values. If the number of measured parameters is more than two, one can perform a test of that particular gravity theory. The test is passed if all three curves intersect in the range of errors at one point. Generally, n measured PK parameters result in n − 2 independent tests.

If a gravity theory has arbitrary parameters, tests can be done for any fixed values of theory parameters. For DEF gravity, this procedure results in an ‘allowed’ area in {α₀, β₀} space that pass gravity tests within the measurement uncertainties. This area is bounded within the selected confidence limit by a curve, which is usually plotted. To date, the regions allowed by different tests have always included GR.

If the companion is optically bright, one can get information about the system from optical observations. For example, for a binary pulsar system with a bright WD it is possible to obtain the mass ratio and the companion mass using high-resolution optical spectroscopy of the companion (e.g. Antoniadis et al. 2012, 2013). In such cases, we do not have to measure three independent PK parameters based on the radio data, but we can combine several multi-wavelength constraints.

Another issue are possible correlations between PK and other parameters. Observed correlations can come from the theoretical correlation of parameters within the binary model (e.g. T₀ ↔ ω and $P_{\text{b}}\leftrightarrow \dot{\omega }$ in a low eccentric case) and from a non-uniform data sampling (e.g. a non-linear correlation in a parametrisation of Shapiro delay r ↔ s). Often in the past, measured PK parameters were published without any information about observed correlations. This additional information from the timing data is lost in such cases. These days it became a common practice to publish observed correlations and even provide them explicitly. Anderson & Yunes (2019) and Anderson et al. (2019) accounted for possible observed and theoretical correlations between PK parameters within DEF and Mendes-Ortiz gravity (MO, Mendes & Ortiz 2016) via computationally highly demanding Markov chain Monte Carlo (MCMC) simulations utilising the correlation matrix from the TEMPO output. However, this elaborate approach would not work if the relativistic effect is not well-measured or the dependence between the parameters is highly non-linear. If the relation is non-linear, the correlation matrix gives information only about the linear contribution. To fully account for a non-linear relation, one has to obtain a full probability density function in the parameter space.

Contrary to the PK method, the DDSTG model accounts for all even non-linear correlations naturally and breaks them by a direct fit of the masses. This was already one of the main advantages of the DDGR model; and is one of the reasons why the DDSTG model uses the same superior approach while extending it to STG theories. The model accounts for all effects internally and extracts all the information from the timing data. The information is not lost even if the correlating parameter cannot be measured but influences other parameters. Such possible correlations are relevant for PSR J1141-6545 (Venkatraman Krishnan 2019; Venkatraman Krishnan et al. 2020), where there is a correlation between the time dilation parameter (γ) of the Einstein delay and the rate of change of the projected semi-major axis due to the spin-orbit coupling caused by the fast rotating WD companion (ẋ^SO). Furthermore, the Shapiro delay cannot be measured independently but is still essential (see e.g. Bhat et al. 2008). We expect the DDSTG model to be of particular advantage for PSR J1141–6545 and there will be a dedicated paper (Venkatraman Krishnan et al., in prep.) about applying the new model to this system.

Compared to conventional methods, the constraints obtained by the DDSTG model may, in general, become either tighter or weaker depending on the particular case. However, the resulting restrictions are more reliable by the construction. This weakening of the limits can happen because of unaccounted correlations.

4 Application to PSR J2222-0137

As a demonstration of the DDSTG model, we now apply it to published timing data of a binary pulsar. We select PSR J2222–0137 because of its unusual characteristics, which we now list in detail.

4.1 About the system

We are particularly interested in recycled pulsars because of their, in general, better timing precision essential for precise gravity tests. The radio pulsar PSR J2222–0137 was discovered in the Green Bank Telescope (GBT) 350 MHz drift scan pulsar survey (Boyles et al. 2013). It is a mildly recycled pulsar which has a spin period (P) of 32.8 ms. The pulsar is in a binary system with P_b of 2.44576 days and x of 10.848 light-seconds. We also expect the pulsar’s spin axis to be aligned to the orbital angular momentum of the system. The alignment happens as a consequence of the pulsar recycling process, when the pulsar accretes matter from the companion.

PSR J2222–0137 is already known as a unique laboratory for testing gravity theories because of its special characteristics; for more detailed information, we refer the reader to Cognard et al. (2017) and Guo et al. (2021b). It is one of the closest pulsars known, with the most precise distance measured with Very Large Baseline Interferometry (VLBI, see Deller et al. 2013) and excellent timing precision. The system has a highly significant detection of the Shapiro delay as well as the measured rate of advance of periastron $(\dot{\omega })$, which yield precise mass measurements and ~1% test of the GR predictions for the Shapiro delay (Guo et al. 2021b).

The most important characteristic in the scope of this paper is the precise measurement of Ṗ_b. Given the precise masses, this can be compared with a precise prediction for the orbital decay due to GW damping, furthermore the kinematic contributions to Ṗ_b (see Sect. 4.4) can be estimated precisely because of the good distance measurement. The Ṗ_b measurement constrains dipolar GW emission in this system, which would arise in DEF gravity (Damour & Esposito-Farèse 1992a) because of the large difference in the compactness of the components (NS and WD). Finally, the mass of the pulsar (m_p = 1.831(10) M_⊙) (Guo et al. 2021b) lies in the range where spontaneous scalari-sation happens, yielding strong limits on this highly non-linear phenomenon (Shao et al. 2017; Zhao et al. 2022).

4.2 Observational data

The PSR J2222–0137 timing data used in this work are the same as used by Guo et al. (2021b). This includes observations of this pulsar going back to its original follow-up, which started on 2009 June 23 using the Green Bank Telescope (GBT); however these ended on 2011 December 26. Regular observations with the three largest European radio telescopes (3ERT) started the following years: the Nançay Radio Telescope (NC) on 2012 October 4, The Lovell Telescope (LT) at Jodrell Bank on 2012 November 20 and the dense orbital campaigns with the Effelsberg 100-m radio telescope (EFF) on 2015 October 26. These observations continue to the present day; however Guo et al. (2021b) only analysed the data obtained until 2021 May 2; their processing and how the resulting ToAs were analysed is described in detail by Cognard et al. (2017) and Guo et al. (2021b).

Observations with the MeerKAT radio telescope array and the Five hundred meter Aperture Spherical Telescope (FAST) started in 2019 September 24 and 2020 October 5, but these timing data were not used by Guo et al. (2021b). However, the latter authors did use FAST data for a detailed study of the emission properties of PSR J2222–0137. In this work, we also extend the existing data by simulated data assuming the same timing properties of all of these telescopes to simulate how the timing parameters for this system will improve in the near future (see Sect. 5).

4.3 Mass-mass diagrams in DEF gravity

The illustrative and straightforward way to explore if a given gravity theory agrees with binary pulsar observational data is to calculate mass-mass diagrams. As mentioned in Sect. 3.4, firstly, one measures values of PK parameters with their uncertainties by fitting the phenomenological DD binary model to binary pulsar data. Then the observed PK values can be compared with their theoretical predictions of a specific gravity theory. The PK parameters from Eq. (20) depend on the mass of the pulsar (m_p) and the mass of the companion (m_c) (see Appendix B). For PSR J2222–0137 we use the measured PK parameters from Guo et al. (2021b).

To illustrate how the tests are performed, we select two points in DEF gravity parameter space near the scalarisation region with large negative β₀ and calculate mass-mass diagrams for them. The test is sensitive to the predicted dipolar contribution to Ṗ_b, which rises dramatically in the region of spontaneous scalarisation. The first point with α₀ = −10⁻⁴, β₀ = −4.3 does not predict a strong enough dipolar contribution and passes the test within a 1σ limit. The corresponding mass-mass diagram is presented in the left panel of Fig. 1. In contrast, the second point with α₀ = −10⁻⁴ and a bit smaller β₀ = −4.35 is already excluded because of a rather strong scalarisation of the pulsar leading to significant dipolar GW damping (see the right panel of Fig. 1). The Ṗ_b curve nicely shows that a significant enhancement in the scalarisation happens for a particular interval in the mass range defined by the EOS.

4.4 Various contributions to Ṗ_b

At this point we need to discuss the different contributions to the observed ${\dot{P}}_{\text{b}}^{\text{obs}}$ and their influence on our results. The observed orbital decay of the system consists of many terms $${\dot{P}}_{\text{b}}^{\text{obs}}={\dot{P}}_{\text{b}}^{\text{GW}}+{\dot{P}}_{\text{b}}^{\text{Gal}}+{\dot{P}}_{\text{b}}^{\text{Shk}}+{\dot{P}}_{\text{b}}^{\text{Tid}}+{\dot{P}}_{\text{b}}^{\dot{M}}+{\dot{P}}_{\text{b}}^G$$(23)

where the first term ${\dot{P}}_{\text{b}}^{\text{GW}}$ is due to GW damping and can include dipolar and monopolar terms in DEF gravity, besides the general quadrupolar prediction of GR. The full expressions of different GW contributions are presented in Appendix (B.9). The next two terms are of kinematic origin: the Shklovskii effect (${\dot{P}}_{\text{b}}^{\text{Shk}}$) and the Galactic contribution (${\dot{P}}_{\text{b}}^{\text{Gall}}$). They are the result of a time-varying Doppler factor D due to an (apparent) radial acceleration between the pulsar binary and the Solar System (Damour & Taylor 1991). The last three terms come from tidal effects, mass loss in the system, and a possible temporal variation of the gravitational constant G.

We are mainly interested in measuring the GW emission term (${\dot{P}}_{\text{b}}^{\text{GW}}$) and comparing it with the prediction of DEF gravity. The most significant additional effects in PSR J2222–0137 come from kinematic contributions ${\dot{P}}_{\text{b}}^{\text{Gal}}$ and ${\dot{P}}_{\text{b}}^{\text{Shk}}$. These two effects arise beyond the binary system and can be combined in the overall external contribution $${\dot{P}}_{\text{b}}^{\text{ext }}={\dot{P}}_{\text{b}}^{\text{Gal}}+{\dot{P}}_{\text{b}}^{\text{Shk}}.$$(24)

The last three terms in Eq. (23) are parts of the internal contribution $${\dot{P}}_{\text{b}}^{\text{int}}={\dot{P}}_{\text{b}}^{\text{GW}}+{\dot{P}}_{\text{b}}^{\text{Tid}}+{\dot{P}}_{\text{b}}^{\text{M}}+{\dot{P}}_{\text{b}}^{\text{G}}$$(25)

and can be neglected for PSR J2222–0137. Thus the internal contribution has only one significant term left, the GW term. For the discussion about possible account of ${\dot{P}}_{\text{b}}^G$ we refer the reader to Sect. 8. Within the DDSTG model, the parameter ‘XPBDOT’ in TEMPO is used for treating both external and internal effects ${\dot{P}}_{\text{b}}^{\text{ext}}+{\dot{P}}_{\text{b}}^{\text{int}}-{\dot{P}}_{\text{b}}^{\text{GW}}$ simultaneously subtracting the GW emission term. In our case XPBDOT accounts for only external Shklovskii and Galactic effects.

There is an extensive analysis of the external terms for PSR J2222–0137 by Guo et al. (2021b) and their determined values are: ${\dot{P}}_{\text{b}}^{\text{Gal}}=-0.0142(13)\times {10}^{-12} \text{s} {\text{s}}^{-1}$ and ${\dot{P}}_{\text{b}}^{\text{Shk}}=0.2794(12)\times {10}^{-12} \text{s} {\text{s}}^{-1}$. These values correspond to the external variation of the observed orbital period of $${\dot{P}}_{\text{b}}^{\text{ext }}=0.2652(18)\times {10}^{-12} \text{s} {\text{s}}^{-1},$$(26)

which is consistent within 2σ with the total observed value $${\dot{P}}_{\text{b}}^{\text{obs}}=0.2509(76)\times {10}^{-12} \text{s} {\text{s}}^{-1},$$(27)

leaving the internal contribution ${\dot{P}}_{\text{b}}^{\text{int }}={\dot{P}}_{\text{b}}^{\text{obs }}-{\dot{P}}_{\text{b}}^{\text{ext }}=-0.0143(78)\times {10}^{-12} \text{s} {\text{s}}^{-1}$ consistent with the GR prediction for quadrupolar GW emission ${\dot{P}}_{\text{b}}^{\text{GR}}=-0.00809(5)\times {10}^{-12} \text{s} {\text{s}}^{-1}$ (Guo et al. 2021b).

According to Damour & Taylor (1991), Nice & Taylor (1995), and Lazaridis et al. (2009), the Galactic differential acceleration may be analytically approximated with the expression $$\frac{{\dot{P}}_{\text{b}}^{\text{G}al}}{P_{\text{b}}}=-\frac{K_z\left|\sin b\right|}{c}-\frac{{\text{Θ}}_0^2}{c{R}_0}\left(\cos l+\frac{\beta }{{\beta }^2+{\sin }^{2}l}\right)\cos b,$$(28)

where l is Galactic longitude, b the Galactic latitude and β = (d/R₀) cos b − cos l. The quantity K_ɀ is the vertical component of the Galactic acceleration, which for Galactic heights ɀ ≡ |d sin b| ≤ 1.5 kpc can be approximated with sufficient accuracy by $$K_z\left({10}^{-9} \text{cm} {\text{s}}^{-2}\right)\simeq 2.27z_{\text{kpc}}+3.68\left(1-{\text{e}}^{-4.31z_{\text{kpc}}}\right),$$(29)

where ɀ_kpc = ɀ(kpc) (Holmberg & Flynn 2004; Lazaridis et al. 2009). From GRAVITY Collaboration (2021) we can take a value for the Sun’s Galactocentric distance R₀ = 8275 ± 9 ± 33 pc. The Galactic circular velocity at the location of the Sun (Θ₀) is taken to be 240.5(41) km s⁻¹ (see Guo et al. 2021b, and references therein).

The Shklovskii contribution (Shklovskii 1970) can be calculated using $${\dot{P}}_{\text{b}}^{\text{Shk}}=\frac{\left({\mu }_{\alpha }^2+{\mu }_{\delta }^2\right)d}{c}{P}_{\text{b}}\text{, }$$(30)

where µ_α and µ_α are the proper motion in right ascension (RA) and declination (DEC), respectively, and d is the distance to the pulsar. The astrometric values and uncertainties for PSR J2222–0137 are taken from Guo et al. (2021b). For PSR J2222–0137, the uncertainty in ${\dot{P}}_{\text{b}}^{\text{ext}}$ is small compared to the precision of the Ṗ_b measurement. The uncertainty of ${\dot{P}}_{\text{b}}^{\text{ext}}$ in our case can be neglected and therefore we can provide the fixed XPBDOT parameter in TEMPO.

Fig. 1 Mass-mass diagrams for the PSR J2222–0137 system at the zone of high non-linearity where a fast transition to a strongly scalarised NS happens. PK parameters are obtained by fitting the DD model to the timing data from Guo et al. (2021b). Calculations are preformed within DEF gravity for the MPA1 EOS allowing NSs with maximum mass in GR of MGR ≃ 2.46 M⊙. Left panel corresponds to the point in the DEF gravity parameter space with α0 = −10−4, β0 = −4.3, which is not excluded by the test. Right panel shows the same test but for α0 = −10−4, β0 = −4.35. The prediction of a strong dipolar contribution ${\dot{P}}_{\text{b}}^D$ fails the test for the more negative β0 = −4.35. The axes correspond to the masses of the pulsar mp and the companion c. The shadowed area is the allowed region at 68% confidence level (CL) limit for a corresponding PK parameter. The solid green line corresponds to the observed value of ${\dot{P}}_{\text{b}}^{GW}$ due to the GW emission.

4.5 Results of applying DDSTG

Our main goal is to obtain limits on DEF gravity parameters {α₀, β₀} by applying the DDSTG model. TEMPO allows one to calculate a χ² value for every particularly selected pair of values (α₀, β₀). Therefore we obtain χ² values for a grid of parameters and compare them to each other. For straightforward comparison we subtract the minimum ${\chi }_{\min }^2$ over the {α₀, β₀} grid from calculated χ². The location of the minimum ${\chi }_{\min }^2$ is statistically in agreement with GR (α₀ = 0, β₀ = 0). The shifted quantity $\text{Δ}{\chi }^2={\chi }^2-{\chi }_{\min }^2$ has a χ² (d.o.f. = 2) distribution with two degrees of freedom. Finally, we construct contours of a fixed ∆χ² value corresponding to desired confidence levels. In this paper, we adhere to 90% confidence level limits with ∆χ² ≃ 4.6.

During each run, TEMPO fits spin parameters, astrometric parameters, Keplerian parameters, two masses, and other parameters of interest, for example, ẋ. A value of the external contribution to the rate of the orbital period change $\left({\dot{P}}_{\text{b}}^{\text{ext}}\right)$ is fixed and selected, so that it fully accounts for Shklovskii effect and the Galactic Ṗ_b contributions (Shklovskii 1970; Brumberg et al. 1975; Damour & Taylor 1992). The resulting ∆χ2 map in {α₀, β₀} space for PSR J2222–0137 is presented in Fig. 2. The limit on DEF gravity parameters is placed by a contour of 90% confidence level limit of χ²(d.o.f. = 2) statistics.

Generally, the uncertainty of ${\dot{P}}_{\text{b}}^{\text{ext}}$ can be important for other systems where it is not well-constrained. In this case, the DDSTG model allows a complete description with account for all the uncertainty due to external effects ${\dot{P}}_{\text{b}}^{\text{ext}}$. First we have to assume a prior distribution for the ${\dot{P}}_{\text{b}}^{\text{ext}}$ value. Then we have to calculate the χ² value for 3D map with $\left\{{\alpha }_0,{\beta }_0,{\dot{P}}_{\text{b}}^{\text{ext}}\right\}$ as parameters of the three axes. Finally, we marginalise the probability density for the ${\dot{P}}_{\text{b}}^{\text{ext}}$ axis using Bayes’ theorem and obtain the corrected 2D map of ${\chi }_{\text{corr }}^2$. The details of the procedure are presented in Appendix E.

Calculating an extensive grid in a 2D parameter space with a lot of points (e.g. −500) for each axis to obtain a finely resolved contour of the desired ∆χ2 value is too computationally demanding. In this work, we use an adaptive mesh refinement technique to trace the location of the desired contour. We start from the sparse grid with 9 × 9 points in the {α₀, β₀} space. Then, the algorithm resolves all the cells that can have the contour line inside them. The algorithm repeats the refinement of important cells until we obtain the desired curve fineness.

For a N × N grid, the naive approach utilises O(N²) iterations, whereas adaptive mesh refinement has only O(N log₂ N) complexity. For the present work we use the refinement level of six, resulting in the final grid with 513 × 513 points. We need such a large resolution in the parameter space because of the ‘horn’ feature at large α₀ values near β₀ ≃ −2. A simple analysis shows that adaptive refinement requires one to calculate 56 times fewer points to resolve a single contour for a grid with 513 × 513 points. In Fig. 3, we present an adaptive refinement map corresponding to the search of a contour from Fig. 2. A high level of refinement is performed only along the contour of 90% CL limit. In Appendix D, one can also find the mass-mass diagram showing what happens in the region of the horn in terms of PK parameters.

Fig. 2 ∆χ2 map in DEF gravity parameter space for PSR J2222–0137. The test is performed by applying the DDSTG model to the timing data of Guo et al. (2021b) and assuming the rather stiff MPA1 EOS. The red line corresponds to 90% CL limit (∆χ2 ≃ 4.6), the area above the grey line is restricted by Cassini mission. GR with α0 = 0, β0 = 0 lies beyond the plotted domain in the blue region at the bottom at the infinity.

Fig. 3 Map of the refinement level in the DEF gravity parameter space. The test is the same as in Fig. 2. Areas in the parameter space with a higher refinement level have an exponentially growing resolution. The adaptive refinement procedure resolves the contour of 90% CL limit and saves computational time dramatically.

4.6 Comparison between DDSTG and the PK method

In this section, we compare the constraining power of the newly developed approach of this paper with that of existing procedures. For this reason we perform the test with the traditional PK method based on the PK parameters from the timing data. The corresponding PK parameters are measured by means of fitting the DD model (Damour & Deruelle 1986) to the same timing data. Then for each specific choice of the theory parameters {α₀, β₀} we fit two masses m_p and m_c to minimise the χ² value. The χ² is calculated by comparison of the observed PK parameters $p_i^{\text{obs}}$ and predicted values $p_i^{\text{theory}}$ from the theory $${\chi }^2\left(p_i^{\text{obs }};m_{\text{p}},m_{\text{c}}\right)={\displaystyle \sum _i{\left({\sigma }_{p_i}^{\text{obs}}\right)}^{-2}}{\left(p_i^{\text{theory }}\left(m_{\text{p}},m_{\text{c}}\right)-p_i^{\text{obs}}\right)}^2,$$(31)

where ${\sigma }_{p_i}^{\text{obs}}$ is the standard deviation of the i-th measured PK parameter. Finally, we calculate the grid of χ² values over the desired {α₀, β₀} space with the same settings as for the DDSTG approach and shift χ² values by its minimum ${\chi }_{\min }^2$. The detailed explanation of this common PK method can be found in Damour & Esposito-Farèse (1998). For both methods the minimum ${\chi }_{\min }^2$ is in statistical agreement with the GR value ${\chi }_{\text{GR}}^2$.

Figure 4 shows the comparison between the two methods. Each point on the plot presents a particular pair of {α₀, β₀} parameters. The analysis shows that, for PSR J2222–0137, the DDSTG model produces higher or equal values of ∆χ2 compared to the traditional PK method. Despite the absence of any strong correlations between the PK parameters in PSR J2222–0137, the new approach produces slightly more restrictive results for the area with negative β₀. For certain areas in {α₀, β₀} space the difference in the shifted $\text{Δ}{\chi }^2={\chi }^2-{\chi }_{\text{min}}^2$ values becomes statistically significant when we calculate contours of 90% confidence level limit (∆χ2 ≃ 4.6).

Fig. 4 Comparison in the constraining power between the DDSTG model and the method based on measured PK parameters with the DD model. Both methods calculate χ2 values for a grid in {α0, β0} space, each point corresponds to a unique theory. β0 values are presented by the colour, while α0 values cover the [−10−1, −10−4] range. Red vertical line corresponds to the 90% CL limit, all points to the left are allowed by the test.

5 Predictions for PSR J2222–0137

Our next step is to estimate what enhancement in limits we can expect with future observations of PSR J2222–0137. For this purpose we simulate fake TOAs for a set of radio observatories, assuming realistic timing precision estimated from real timing data.

5.1 Simulated datasets for FAST, MeerKAT, and 3ERT

The PSR J2222–0137 timing data used up to this point are described in Sect. 4.2. We now describe our simulations, which show what we might be able to achieve in the near and foreseeable future with timing from this system.

We simulate TOAs spanning 10 yr from 2021 to 2030 based on the current precision of the TOAs obtained with the 3ERT telescopes, as well as the TOA precisions from ongoing observations from MeerKAT and FAST. These simulations are conservative since they assume that there will be no improvement in the existing capabilities at these telescopes.

For TOAs from FAST, the radiometer noise reduces significantly thanks to its large collecting area, while the jitter noise becomes the primary limitation of timing precision. We find, however, that increasing the integration time to 15 min can largely reduce the jitter noise and eliminate its contribution to the timing precision (as it scales with the number of averaged pulses (N_p) as ${\sigma }_{\text{J}}\propto 1/\sqrt{N_p}$ (Lorimer & Kramer 2004). Therefore, the median TOA uncertainty from 15-min TOAs are adopted in the simulation.

Table 1 lists the telescopes assumed in our simulation, with the information on their effective diameter (Ø_eff), observing bandwidth, and TOA uncertainties at L-band. All TOA uncertainties are scaled to an integration time of 15 min over the full bandwidth. For each telescope, we assume one full orbit observation (~60 h) per year, and split the observations into a monthly cadence, i.e. 5 h per month⁸, to allow a good estimation of timing parallax (which requires a good coverage of Earth’s orbit). This is important for the estimation of uncertainties in the external Ṗ_b contributions shown in the next section.

The simulations are performed using the program developed in Hu et al. (2020). First, we simulate TOAs based on the above assumptions, and add the TOAs from Effelsberg, Nançay, and Lovell telescopes together to be compared with MeerKAT and FAST. For 3ERT, the simulated TOAs are combined with the existing TOAs in Guo et al. (2021b). We then adjust the TOAs to perfectly match the timing parameters measured in Guo et al. (2021b), and add a Gaussian white noise to each TOA based on its σ_TOA. Finally, we fit for timing parameters and obtain their uncertainties, among which Ṗ_b and timing parallax are of most important here. The whole process is done with the TEMPO DDSTG model.

5.2 Contributions to the uncertainty of Ṗ_b

The predicted uncertainties of Ṗ_b are shown in Fig. 5, where the pink, orange, and blue lines show the improvement in time with the simulated data from 3ERT, MeerKAT, and FAST, respectively. As discussed in Sect. 4.4, we also need to account for uncertainties from external effects, ${\dot{P}}_{\text{b}}^{\text{ext}}$. They are also expected to become more precise in the future because of anticipated improvements in the model for the Galactic gravitational potential, distance, and proper motion from future observations. Correspondingly improved values for ${\dot{P}}_{\text{b}}^{\text{Gal}}$ and ${\dot{P}}_{\text{b}}^{\text{Shk}}$ may then be estimated using Eqs. (28) and (30). For the Shklovskii effect, we find that it is mostly limited by the uncertainty in the distance, which comes from VLBI parallax or timing parallax measurement, whichever is better. With the simulated FAST data, the measurements of timing parallax and proper motion improve quickly with time. In particular, the uncertainty of timing parallax will soon surpass the VLBI parallax and hence improve the distance measurement. The corresponding uncertainty in ${\dot{P}}_{\text{b}}^{\text{Shk}}$ will decrease with time as shown by the brown line in Fig. 5. This line is below the predicted uncertainties from observed Ṗ_b most of the time, except at the very end of the simulation when compared to $\text{Δ}{\dot{P}}_{\text{b}}^{\text{obs }}$ with FAST (blue line). In addition, future VLBI parallax measurements will likely be improved so that Shklovskii effect will not be a limiting factor for Ṗ_b.

As for the contribution from the Galactic acceleration, a typical uncertainty in its vertical component (ΔK_z) is about 10%, which contributes the most in $\text{Δ}{\dot{P}}_b^{Gal}$, shown as the green dashed line in Fig. 5. With FAST, Ṗ_b will then be limited by ΔK_ɀ from 2024 onwards, if there is no improvement for this quantity. In fact, we find that ΔK_ɀ needs to be improved to ≲3% (see the green dash-dotted line) to not limit the precision of Ṗ_b before 2030. The uncertainties in the Galactic potential do not limit the precision in this case. For the scope of further analysis, we assume that with future observations on pulsar timing, VLBI parallax, and Galactic acceleration, the precision of Ṗ_b will not be limited by Shklovskii and Galactic effects. For Shklovskii contribution it is a reasonable assumption at least up to 2030. We also expect our knowledge about the Galactic potential to improve in the future especially in the proximity of the Solar System, thus we assume the improvement of ΔK_ɀ for the selected very close pulsar to be 3%.

Fig. 5 Comparison of different contributions in the uncertainty of Ṗb for simulated data from 2021 to 2030. The pink, orange and blue lines show the uncertainty of Ṗb using the simulated data from 3ERT (EFF, NC, LT), MeerKAT (MK), and FAST, respectively. The brown line indicates the uncertainty of ${\dot{P}}_{\text{b}}^{\text{Shk}}$ when using the timing parallax (πx) measured from simulated FAST data. The green lines indicate the uncertainties of ${\dot{P}}_{\text{b}}^{\text{Gal}}$ when assuming 10% (dashed) and 3% (dash-dotted) uncertainties in the vertical component of Galactic acceleration Kɀ.

5.3 Potential future constraints on DEF gravity from PSR J2222—0137

We apply the same techniques as described in Sect. 4 including the adaptive mesh refinement for contours. The map of Δχ² for a combined simulated data is presented in Fig. 6. By 2030, we can expect a significant improvement in the limits for large positive β₀. This region is susceptible to dipolar gravitational emission. The significant improvement in Ṗ_b measurement pushes the limit below the Cassini limit.

In Fig. 6, we also present the comparison of constraining power for different observatories. The tightest constraints may be obtained for the combination of all three observatories. However, the precision of FAST is high enough to be significantly constraining on its own.

Fig. 6 Δχ2 map for PSR J2222–0137 from simulated data (2021–2030) combining different observatories and using the MPA1 EOS. Assumed one orbit per year for FAST, MeerKAT (MK), and 3ERT (EFF, NC, LT) based on their observation precision, appending to the existing dataset (Guo et al. 2021b). Each contour corresponds to 90% CL limits on DEF gravity parameters from different datasets, grey line depicts the limit from the Cassini mission.

6 Simulations for PSR-BH systems

As a further application of the new timing model, we investigate what limits on DEF gravity we can obtain from a PSR-BH system. As mentioned above, our DDSTG implementation already allows for a BH to be specified as a companion, even though such a system has so far not been found. We apply the same technique used in Sects. 4 and 5 to synthetic TOAs for three different hypothetical PSR-BH systems.

6.1 Simulated FAST datasets

To investigate possible limits on DEF gravity from PSR-BH systems, we select presumed realistic parameters for simulating fake TOAs. These systems are assumed to comprise a 1.4 M_⊙ pulsar and a 10.0 M_⊙ BH in a highly eccentric orbit (e = 0.6). As PSR-BH systems are more likely to reside in GCs (see the discussion in Sect. 6.4 and references therein), we assume that our hypothetical systems are located in the GC M15 and take the position (RA and Dec) and dispersion measure (DM) of PSR B2127+11C (M15C) for our simulations. To investigate how the limits on DEF gravity depend on the orbital period, we consider three cases with orbital periods of 3 days, 1 days, and 8 h, respectively. The selected parameters for PSR-BH systems are presented in Table 2.

We simulate datasets with two different TOA uncertainties (for 15-min integration time): moderate 10 µs (which is typical for pulsars in GCs), and precise 1 µs. Two TOA uncertainties help us to explore the dependence of the DEF gravity test on the TOA precision. We effectively assume a recycled millisecond pulsar (MSP), because they in general produce a better timing precision. We assume 6 h of timing observation on these systems every two weeks, each session starts at a random orbital phase. All simulations cover a time span of 5 yr and have the same number of TOAs (n_TOA = 3144).

6.2 Tests of DEF gravity in case of a BH companion

Binary systems consisting of a pulsar and a white dwarf (PSR-WD) are particularly interesting for constraining STG theories due to their high asymmetry in compactness (α_WD ≃ α₀). The theory predicts a higher rate of orbital energy loss due to dipolar radiation from such asymmetric systems (Will 1993; Damour & Esposito-Farèse 1996). However, generally, PSR-BH systems are expected to be even more asymmetric up to a certain positive β₀, as was pointed out by Damour & Esposito-Farèse (1998). As a result of the no-scalar-hair theorem to BHs in DEF gravity one has (Hawking 1972; Damour & Esposito-Farèse 1992a) $${\alpha }_{\text{BH}}=0,{\beta }_{\text{BH}}=0.$$(32)

This approximation is only valid for stationary BHs with an asymptotically flat spacetime and an asymptotically constant scalar field (Berti et al. 2013), which is not true anymore in the presence of a compact scalarised companion. However, a justification that it is still an excellent approximation can be found in Liu et al. (2014). The pulsar orbiting in an eccentric orbit will eventually induce a time-dependent scalar field at the location of the BH. This scalar field results in an induced effective scalar charge $\left({\alpha }_{\text{BH}}^{\text{induced }}\right)$ which, however, is totally negligible (≲8 × 10⁻¹⁴ α_p; cf. Eq. (50) in Liu et al. 2014).

The absence of scalar charges for the BH results in simplified relations of the PK parameters. The main consequence of a BH presence is that all PK parameters except Ṗ_b are identical to their GR expressions with an appropriate rescaling of the masses. In the GW damping sector, all the multipoles, for example, monopolar, dipolar and quadrupolar modes, are still present. Thus PSR-BH systems are susceptible but only to the test of GW damping via Ṗ_b (Damour & Esposito-Farèse 1992a; Mirshekari & Will 2013; Liu et al. 2014).

Fig. 7 Δχ2 map in the DEF gravity parameter space from the simulated 5-yr timing data for three different PSR-BH systems. The left panel corresponds to 10 µs TOA uncertainty, the right panel to 1 µs. Solid lines show the 90% CL limits for PSR-BH systems with different orbital periods: 3 days, 1 day and 8 h. The grey dashed line depicts the limit from the Cassini mission and green dash-dotted line is the current limit from PSR J2222–0137 discussed in Sect. 4 (see Fig. 2).

6.3 Results

The limits obtained from the DDSTG model for the simulated PSR-BH systems are shown in Fig. 7. For the tests presented here, we do not include any potential uncertainty in the observed ${\dot{P}}_{\text{b}}^{\text{obs}}$ value due to the external effects. The restrictive power of the test increases for more relativistic systems with shorter orbital periods. To place new limits on DEF gravity with a moderate TOA precision of 10 µs, the orbital period should be a fraction of a day.

Another way to improve the restrictive power is to increase the precision of the TOAs (compare left and right panels of Fig. 7). The obtained limits strongly depend on the accuracy of Ṗ_b measurement as most of the deviations from GR come from the predicted dipolar GW emission. Table 3 shows the predicted GR values of Ṗ_b and compares them to the corresponding ∆Ṗ_b uncertainties. The lower uncertainty and the higher the predicted value the better the test.

We expect a dramatic improvement of limits on DEF gravity from a relativistic PSR-BH system. Due to higher asymmetry, the test is extremely sensitive to the precision of Ṗ_b measurement. To answer the possibility of obtaining accurate enough TOAs, we have to argue where we can find such a system.

6.4 Origin of PSR-BH systems

A binary pulsar system with a stellar-mass BH may originate from several different evolutionary scenarios. The first way is a standard evolution of a massive binary system (Voss & Tauris 2003). It results in a PSR-BH system with a young slow spin-period (~0.1–1 s) pulsar and a wide, eccentric orbit. The second path is a reversal mechanism (Sipior et al. 2004) taking place under a specific set of circumstances. The pulsar is formed first and later spun-up by accretion during the red-giant phase of the companion (Pfahl et al. 2005), later becomes a BH. This is similar to the origin of PSR J1141–6545, a system where a massive WD star formed before a more massive NS (Tauris & Sennels 2000). The reversal mechanism may give a more desirable result of a recycled pulsar in a system with a BH, because recycled pulsars have generally more precise timing and a more stable rotation than the slow ‘normal’ pulsars (Verbiest et al. 2009).

Moreover, the third possible way to form a PSR-BH system is a multiple body encounter, happening in regions of high stellar density, for example, GCs and the Galactic Centre region (Verbiest et al. 2009; Clausen et al. 2014). Such encounters are the reason why there are ~10³ times more low-mass X-ray binaries (LMXBs) per unit stellar mass in GCs than in the Galactic disk (Clark 1975), which results in a similarly enhanced proportion of MSPs. In these encounters, an old, inactive NS approaches a main sequence (MS) binary so closely that a chaotic interaction ensues. The most likely result of such an interaction is that the two most massive objects (in this case the old, recycled NS and the more massive MS star) will form a more compact binary system, with the lighter MS star component being ejected at high velocity (see review by Phinney 1992). That MS will evolve, fill its Roche lobe, and start transferring matter onto the NS, which is spun up in the process, this is the aforementioned LMXB stage.

If left undisturbed, many of these then evolve into binary MSP systems. A consequence of this is that GCs not only have many MSPs, but that the number of MSPs in each cluster appears to be roughly proportional to its rate of stellar encounters (Γ) (Verbunt & Hut 1987).

However, in some GCs – especially those with collapsed cores – stellar densities are so high that only few binaries evolve without being disturbed; these GCs have a large interaction rate per binary (γ_GC) (Verbunt & Freire 2014). This means that, even after being recycled, an MSP has a high probability of undergoing further (‘secondary’) exchanges.

In such exchanges, an incoming massive star interacts chaotically with the components of either a LMXB or a MSP binary. Again, the least massive object is the most likely to be ejected, in this case that will be the low-mass star that recycled (or was still recycling) the NS. A new binary will form, consisting of the NS and the intruding massive star. If the latter is degenerate, then the system will not undergo accretion and the orbital circulari-sation that comes with it, but will keep the orbital eccentricity it acquired after formation.

Several such obvious products of secondary exchange encounters have been found in GCs, invariably with a large γ_GC: PSR B2127+11C, in the core-collapsed M15 GC (Prince et al. 1991; Jacoby et al. 2006), PSR J0514–4002A, in NGC 1851 (Freire et al. 2004; Ridolfi et al. 2019; recently two additional such systems, D and E, were found in the same GC, see Ridolfi et al. 2022), PSR J1807-2500B, in the core-collapsed NGC 6544 (Lynch et al. 2012), PSR J835-3259A, in the core-collapsed NGC 6652 (DeCesar et al. 2015) and PSR J1823–3021G, in the core-collapsed NGC 6624 (Ridolfi et al. 2021). These discoveries suggest (but by no means assure us) of the possibility that similar secondary exchange encounters might produce MSP binaries with even more massive companions, such as BHs. As for the massive MSP binaries above, such MSP-BH systems would be preferentially produced in the GCs with large γ_GC; this is one of the reasons why they are targeted by the MeerKAT TRAPUM survey (Ridolfi et al. 2021).

These secondary exchange products are that their orbital periods vary between 8 h in the case of B2127+11C and 18.8 d in the case of PSR J0514–4002A, their orbital eccentricities are between 0.38 and 0.90. Thus the simulated PSR-BH systems listed in Table 2 have realistic orbital periods and eccentricities.

6.5 Influence of a globular cluster origin

For the pulsars observed in GCs, the derivatives of the spin period, P, and in the case of binary pulsars, the derivative of the orbital period, Ṗ_b are contaminated (and in most cases dominated) by the line-of-sight component of the acceleration of the pulsar (or binary) in the gravitational potential of the GC (a_GC); this enters Eq. (23) as an additional term that is similar to ${\dot{P}}_{\text{b}}^{\text{Gal}}$: $${\dot{P}}_{\text{b}}^{\text{GC}}=\frac{a_{\text{GC}}}{c}{P}_{\text{b}}.$$(33)

For pulsars in GCs, the only radiative test done to date was with PSR B2127+11C (Jacoby et al. 2006). The orbital decay observed in this system is −3.95 ± 0.13 × 10⁻¹² s s⁻¹, which is within ~3% of the predicted value. That test, however, is limited by the fact that the maximum value for |a_GC/c| ~ 6 × 10⁻¹⁸ s⁻¹ (Phinney 1993), i.e., in this case, the maximum value for $\left|{\dot{P}}_{\text{b}}^{\text{GC}}\right|~0.17\times {10}^{-12} \text{s} {\text{s}}^{-1}$. This is of the same order as their measurement precision, which means that this radiative test cannot be improved, unless one could measure the acceleration of the pulsar independently. The situation would be much worse if PSR B2127+11C were closer to the cluster centre, where accelerations are much larger. For instance, PSR B2127+11A, a couple of arcseconds from the centre, has P = 0.1106 s and Ṗ=−2.107×10⁻⁷; this implies that |a_GC/c| > 1.9 × 10⁻¹⁶ s⁻¹ and thus $\left|{\dot{P}}_{\text{b}}^{\text{GC}}\right|>5.5\times {10}^{-12} \text{s} {\text{s}}^{-1}$, a value larger than the orbital decay predicted by GR.

For the three PSR-BH systems in Table 3, the situation would be similar. If they were at the locations of PSR B2127+11C, their values of P_b would imply that the maximum values for ${\dot{P}}_{\text{b}}^{\text{GC}}$ would be 9, 3 and 1 times larger respectively, i.e., 1.53, 0.51 and 0.17 × 10⁻¹² s s⁻¹. This corresponds to, respectively, ~6.9, ~0.37 and 0.02 times the values of ${\dot{P}}_{\text{b}}^{\text{GR}}$ in Table 3. That means that, for the latter system, a ~50-σ test of the radiative properties of a PSR-BH system would be possible. These numbers illustrate the immense advantage of an increasingly shorter P_b, either for GW tests in GCs or in the Galaxy: on one hand, the ${\dot{P}}_{\text{b}}^{\text{GR}}$ increases with $P_{\text{b}}^{-5/3}$, the ‘polluting’ part ${\dot{P}}_{\text{b}}^{\text{GC}}$ decreases as P_b. Thus, the significance of a particular GW test grows, everything else being identical, with $P_{\text{b}}^{-8/3}$.

However, if the 8-h PSR-BH system is placed at the location of PSR B2127+11A, the test would lose its significance almost entirely. Fortunately, even in this case, we can get a firm upper limit for $\left|{\dot{P}}_{\text{b}}^{\text{GR}}\right|$; the reason is that we also measure the spin period derivative. That is also affected by the acceleration in the cluster and by other terms, such as the Shklovskii effect. Adding the equations for Ṗ and Ṗ_b, we obtain: $${\dot{P}}_{\text{b}}^{\text{GR}}={\dot{P}}_{\text{b}}^{\text{obs}}-P_{\text{b}}{\left(\frac{\dot{P}}{P}\right)}^{\text{obs}}+\frac{1}{2{\tau }_c}{P}_{\text{b}}.$$(34)

We note that all the terms on the right can be measured precisely, except for the characteristic age of the pulsar, τ_c, which cannot be measured independently. Therefore, the last term then quantifies the uncertainty of the test, which is, again, proportional to P_b, which implies that the significance of the test is again proportional to $P_{\text{b}}^{-8/3}$.

However, that last term is necessarily positive. If we assume the pulsar is extremely old, then this term will be very small and we obtain, from the other two terms, a hard lower limit for ${\dot{P}}_{\text{b}}^{\text{GR}}$. Since the latter is negative, this represents a hard upper limit of its magnitude. This means that such a test could still, in principle, falsify GR. An upper limit on the last term can be obtained from the lowest likely value for τ_c; for most MSPs τ_c < 1 Gyr. This means that, for the 8-h PSR-BH system, this unknown term would be at most 0.46 × 10⁻¹² s s⁻¹, implying a ~20-σ test of GR. Thus, this test will be the more significant the larger τ_c is compared to the orbital decay timescale for the binary.

Despite such possible mitigation, it is clear that the location in a GC always degrades the quality of radiative tests. For the 8-h PSR-BH system discussed above, significances of 50 and 20 σ in the measurement of Ṗ_b represent a significant degradation relative to the tests listed in Table 3, where for a 8-h PSR-BH system timed with 10 µs the significance of the Ṗ_b measurement is larger than 1600.

7 DDSTG and other gravity theories

The approach developed in this work is not restricted to DEF gravity. It can be extended straightforwardly to investigate a larger set of alternative gravity theories. Broadly speaking, every theory that maps to the phenomenological DD model (Damour & Deruelle 1986; Damour & Taylor 1992) can be put into the DDSTG framework with a new implementation of appropriate PK formulae. The DD model is a quasi-Keplerian solution in the first post-Newtonian approximation to the dynamics of a two-body system within the modified Einstein-Infeld-Hoffmann (mEIH) framework (Damour & Taylor 1992; Will 1993, 2018b). The mEIH formalism covers a large set of fully conservative gravity theories without the ‘Whitehead term’ in the post-Newtonian limit. Later it was extended by Will (2018a) from fully conservative to semi-conservative theories of gravity. However, this extension requires additional terms in the mEIH Lagrangian, which are not accounted for in the DD model.

The DDSTG model can be directly applied to a certain class of gravity theories without a modification of the PK parameter equations. For example, DEF gravity is a specific case of a mass-less mono-scalar tensor gravity theory with a particular expression of the conformal coupling A(φ). The DDSTG model covers STG theories with any conformal coupling function depending on two arbitrary parameters {α₀, β₀} and a single massless scalar field. To work, the model only requires pre-calculated gravitational form factors for the new coupling function. In case of a more complex coupling function depending on a higher number of parameters the TEMPO code needs to be adapted.

Recently, Mendes and Ortiz (MO, Mendes & Ortiz 2016) introduced an extension of DEF gravity. Mo gravity is an example of a theory that can be easily incorporated into the DDSTG model. Its difference from DEF gravity is in the form of the conformal coupling $$A(\phi )={\left[\text{\hspace{0.17em}cosh\hspace{0.17em}}\left(\sqrt{3}{\beta }_0\phi \right)\right]}^{1/\left(3{\beta }_0\right)},$$(35)

where β₀ is a free parameter. The second parameter α₀ is hidden in MO gravity to the scalar field at infinity. MO theory received attention in recent years and was originally introduced as an analytical approximation to a more fundamental theory, where the action includes quadratic terms of the scalar field coupled to curvature. The developed framework can therefore be extended straightforwardly to MO gravity without the change of the TEMPO implementation.

8 DDSTG and a time-varying gravitational constant

Another interesting extension of the DDSTG model – planned for a future release – is the inclusion of the effects of a temporal variation of the gravitational constant (G). A time-evolving asymptotic scalar field (φ₀) of a gravitating system, generally, leads to a temporal variation of the local gravitational constant. In Damour–Esposito–Farèse gravity, such a change of the Newtonian gravitational constant as measured in the Solar System reads $$\frac{{\dot{G}}_{\text{Cav }}}{G_{\text{Cav }}}=2\left(1+\frac{{\beta }_0}{1+{\alpha }_0^2}\right){\alpha }_0{\dot{\phi }}_0$$(36)

(see e.g. Uzan 2011). Generally, one expects a temporal change in φ₀ to arise from the expansion of the universe and ${\dot{\phi }}_0$ to result from a cosmological model based on DEF gravity (see e.g. Damour & Nordtvedt 1993). However, as part of a more agnostic approach, ${\dot{\phi }}_0$ can be treated as an additional, independent timing parameter.

The main impact on the orbital motion of a binary system that arises from a time-varying gravitational constant is a secular change in the orbital period (Damour et al. 1988). For two weakly self-gravitating masses, one simply has ${\dot{P}}_{\text{b}}^{\dot{G}}/P_{\text{b}}=-2\dot{G}/G$. However, for binary pulsar systems, this simple expression needs to be extended by body-dependent contributions (Nordtvedt 1990, 1993). One then has $$\frac{{\dot{P}}_{\text{b}}^G}{P_{\text{b}}}=-2\left(\frac{{\dot{G}}_{\text{Cav}}}{G_{\text{Cav}}}\right){ℱ}_{AB},$$(37)

where the factor 𝓕_AB accounts for all the corrections related to the strong gravitational fields of the pulsar and its companion, if the latter is also a NS. More specifically, following Nordtvedt (1993), 𝓕_AB accounts for a change in the body-dependent part in the effective gravitational constant G_AB as well as a change in the masses resulting directly from ${\dot{\phi }}_0$⁹.

It has been demonstrated by Wex (2014) that, depending on a parameter space, a pulsar mass, and a EOS, strong-field effects can considerably enhance the effect of a time-varying gravitational constant, i.e. 𝓕_AB ≫ 1. Consequently, accounting for ${\dot{P}}_{\text{b}}^{\dot{G}}$ in binary pulsar tests not only provides an independent test for a varying gravitational constant but also probes strong-field aspects related to a time-varying gravitational constant.

9 Discussion and conclusions

In this work, we developed a new, improved approach for testing STG theories. We examined a specific class of STG theories known as DEF gravity. This approach is based on a new timing model, called the DDSTG model, which is an extension of the DD model to work within DEF gravity. Analysis of pulsar timing data with this model overcomes some of the problems of conventional methods, which we have discussed in detail. The DDSTG timing model uses theoretical predictions for PK parameters in DEF gravity and therefore uses a minimal set of binary parameters. For that reason, it accounts for all the possible correlations between these parameters. All the information from the observational data is used to provide the most reliable tests of an alternative theory directly, without intermediate steps with phenomenological parameters of the DD model.

As a demonstration of the DDSTG model, we applied it to the most recently published timing data of the binary pulsar, PSR J2222–0137, described and used by Guo et al. (2021b). This system is of great importance for testing alternative gravity theories because it is very close to us (resulting in the most precise VLBI distance for any pulsar), it has precise timing, and it shows a set of well-measured relativistic effects. The system has a high asymmetry in the compactness between the components: it comprises a massive NS (m_NS ~ 1.82 M_⊙) and a massive WD (m_WD ~ 1.31 M_⊙). This high asymmetry results, for some areas in the DEF gravity parameter space, in the prediction of a very strong dipolar GW contribution to the rate of orbital decay $\left({\dot{P}}_{\text{b}}^D\right)$. The non-detection of dipolar GWs in this system is used to constrain the DEF gravity parameter space. Moreover, the mass of the pulsar lies in the ‘scalarisation gap’ (m_NS ≳ 1.5 M_⊙); this means that strict limits can be placed on the occurrence of spontaneous scalarisation, a highly non-linear phenomenon (Zhao et al. 2022).

The results from the new method confirmed and improved the existing limits on DEF gravity parameters from this system. The DDSTG model appeared to be more constraining in the area near spontaneous scalarisation (β₀ ≲ −4.0) when compared to the commonly used PK method. This suggests that the combination of the DDSTG model with a EOS agnostic approach can improve limits placed on the spontaneous scalarisation.

Moreover, we applied the DDSTG timing model to the simulated TOA for PSR J2222–0137 covering the period of 2021–2030 to see what improvement we can expect from that system in the future. The mock timing datasets were simulated for several large observatories: FAST, 3ERT (EFF, NC, LT), and MeerKAT with TOA uncertainties based on the real timing data. We have discussed the future importance of kinematic contributions (Shklovskii and Galactic) to Ṗ_b, and consequently to the precision of these tests. The main limiting factor to Ṗ_b comes from the uncertainty in the Galactic contribution $\left({\dot{P}}_{\text{b}}^{\text{Gal}}\right)$. In particular, the limit comes from the current uncertainty in the vertical component of the Galactic acceleration (ΔK_ɀ). Our analysis predicts that this limitation will disappear if the current uncertainty of ΔK_ɀ ~ 10% can be improved to ΔK_ɀ ≲ 3%, which is likely to be the case in the future with improvement of models for the gravitational potential of our Galaxy. Future observations are expected to significantly improve the limits on DEF gravity, especially with the use of FAST data.

One of the most promising systems for testing gravity, which we hope to have in the near future, are binary PSR-BH systems. We simulated artificial timing datasets for three eccentric PSR-BH systems with reasonable orbital parameters for three different orbital periods. The results of applying the DDSTG model to the simulated PSR-BH data strongly depended on the precision of the Ṗ_b measurement. We have briefly discussed possible evolution scenarios leading to the formation of the PSR-BH system, such as GC origin and a reversal mechanism. Depending on the place of origin, there might be issues in obtaining a precise intrinsic Ṗ_b, that is to say in accounting for the contamination of Ṗ_b from a kinematic contribution due to the acceleration of the system in the gravitational field of the GC. Depending on the timing precision and orbital properties, PSR-BH systems can place stringent limits on DEF gravity.

In the future, the DDSTG model can be applied to a range of different binary pulsar systems to improve the limits on DEF gravity. We expect especially interesting results from PSR J1141–6545 where the DDSTG nodel is expected to be superior to standard approaches (Venkatraman Krishnan 2019). PSR J1141–6545 is an asymmetric PSR-WD system in a 4.7 h orbit with significant spin-orbit coupling due to the fast rotating WD. Due to the spin-orbit coupling, the system shows a change in the projected semi-major axis ẋ^SO, which has a strong correlation with the time dilation parameter γ. The latter parameter is caused by the precession which cannot be calculated because we do not know the exact spin properties of the WD. Moreover, PSR J1141–6545 shows a weak Shapiro delay in the timing data which the DDSTG model can fully exploit, resulting in a significant improvement to the test. The high asymmetry in the compactness between the components (m_p ~ 1.26 M_⊙ NS and m_c ~ 1.02M_⊙ WD) makes this system a perfect tool for radiative tests of gravity. We expect the DDSTG model to be particularly advantageous because it accounts for both possible correlations and weak relativistic effects (Venkatraman Krishnan et al., in prep.).

Another perspective system to perform DEF gravity tests with the DDSTG model is the double pulsar, PSR 0737-3039A which shows the largest number of PK parameters in the timing data (Kramer et al. 2021). The system consists of two radio pulsars with masses of m_p ≃ 1.34 M_⊙ and m_c ≃ 1.25 M_⊙. Properties of NSs in DEF gravity (gravitational form factors) strongly depend on the choice of a EOS, which in turn affects the test. Thus to put reliable limits on DEF gravity from PSR 0737-3039A independent of a choice of a EOS, one must perform a EOS-agnostic analysis. An EOS-agnostic test means that it is performed with a set of EOSs which are diverse in their properties (see Voisin et al. 2020, who used such a EOS-agnostic approach to constrain DEF gravity with a pulsar in a stellar triple system). With this paper, we have included a set of 11 EOSs varying from soft to stiff (see Appendix A.2) which can be used for such agnostic tests in the future.

Acknowledgements

We thank Lijing Shao for carefully reading the manuscript and giving important comments. We further thank Thomas Tauris and Selma de Mink for useful discussions and suggestions. AB and HH are members of the International Max Planck Research School for Astronomy and Astrophysics at the Universities of Bonn and Cologne. HH acknowledges the support by the Max-Planck Society as part of the ‘LEGACY’ collaboration with the Chinese Academy of Sciences on low-frequency gravitational wave astronomy. This study is partly based on observations with the 100-m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). LG, IC, and GT acknowledge financial support from the ‘Programme National Gravitation, Références, Astronomie, Métrologie (PNGRAM)’ of CNRS/INSU, France. J.W. McKee gratefully acknowledges support by the Natural Sciences and Engineering Research Council of Canada (NSERC), [funding reference #CITA 490888-16]. This publication made use of open source libraries including DIFFERENTIALEQUATIONS.JL (https://diffeq.sciml.ai/) (Rackauckas & Nie 2017), OPTIM.JL (Mogensen & Riseth 2018), FORWARDDIFF.JL (Revels et al. 2016), MATPLOTLIB (Hunter 2007), along with a pulsar analysis package TEMPO (http://tempo.sourceforge.net/) (Nice et al. 2015).

Appendix A Details on calculating grids of NSs

We select slowly rotating axisymmetric approximation of NSs following Damour & Esposito-Farèse (1996). The internal structure of such a NS in DEF gravity is described by a set of 8 ordinary differential equations depending on the radial variable (r). To calculate a single NS the program solves the system of ODEs by applying appropriate boundary conditions at the centre of a NS and spatial infinity. We use a shooting technique to match initial internal parameters with the external structure of the spacetime. The differential equations are solved by means of the Julia library DIFFERENTIALEQUATIONS.JL¹⁰ (Rackauckas & Nie 2017).

A particular solution is determined by the choice of two theory parameters {α₀,β₀}, a EOS and a central pressure (p_c). The central value of the scalar field (φ_c) is determined by the shooting procedure to correspond to the scalar field value at spatial infinity φ₀ = φ(r = ∞), which in our case is, without loss of generality, set to zero. Thus φ_c is not an arbitrary parameter. The central pressure uniquely corresponds to the gravitational mass m_A for most masses of interest and fixed theory parameters {α₀,β₀}. In general, this is not the case – there may be several stable solutions for a fixed mass in the area of strong scalarisation. However, multiple solutions happen for large NS masses and very negative β₀ < −4.5 which is already ruled out and thus of no interest.

In the recent papers devoted to the calculation of NSs in STG theories (Anderson & Yunes 2019), the gravitational form factors from Eq. (6) (they are also often called ‘scalar charges’) are calculated using numerical differentiation. Numerical differentiation is performed on a grid of dependent variables. The precision of this approach is determined by two factors: a) the error in the calculation of the desired quantity and b) the error due to the numerical differentiation formula. The first error depends on the precision of the numerical integration of the structure equations. The second error strongly depends on the fineness of the grid on which the derivatives are calculated. If the grid is too sparse the formula such as the central difference formula ƒ′(x) ≃ [ƒ(x + δx) − ƒ(x − δx)] / [2δx] is not accurate enough. On the other hand, if the grid is too fine, the error arises because of the subtraction of two close numbers in the computer memory. The total relative errors for calculating gravitational form factors usually lie in the range of δ_re1 ~ 10⁰–10⁻⁵. In our work, we propose a more precise way to calculate all the quantities.

A.1 Automatic differentiation

In our program, we utilise automatic aifferentiation (AutoDiff) technique. It is a special technique to calculate derivatives, when the algorithm knows the exact expressions for derivatives of all elementary functions and uses the chain rule to unfold complex derivatives. A review on AutoDiff can be found, for instance, in Margossian (2018). We use an implementation of AutoDiff technique in Julia from ForwardDiff.jl library, which is a part of a broad JuliaDiff framework (Revels et al. 2016). AutoD-iff calculates complex derivatives of quantities with respect to their arguments with the machine-precision (for Double Float numbers the corresponding precision is δ ~ 2 × 10⁻¹⁶ on 64-bit systems). It is done by internal implementation of sophisticated arithmetic on a special type of numbers (dual numbers) in the programming language. A derivative is calculated exactly in the desired point and does not require points in the vicinity. The precision of the method is constant and thus does not depend anymore on the fineness of the grid.

Automatic differentiation enables us to simultaneously calculate both the function F = F (x, y) and its derivatives with respect to arguments $\frac{\partial F}{\partial x},\frac{\partial F}{\partial y}$.It takes twice as long as calculating a mere function itself and both are done with the machine-precision . To calculate a complex derivative, then the arguments are also functions B = B(x, y), C = C(x, y) we apply the chain rule: $${\frac{\partial A(B,C)}{\partial B}|}_C=\left(\frac{\partial A}{\partial x}\frac{\partial C}{\partial y}-\frac{\partial A}{\partial y}\frac{\partial C}{\partial x}\right)/\left(\frac{\partial B}{\partial x}\frac{\partial C}{\partial y}-\frac{\partial B}{\partial y}\frac{\partial C}{\partial x}\right).$$(A.1)

In our program, we calculate the gravitational form factors from Eq. (6), which are complex derivatives of quantities ( e.g. a mass (m_A) and a moment of inertia (I_A) of a NS) taken with respect to external parameters (the value of the scalar field at spatial infinity (φ₀). The properties of NSs are obtained through numerical integration of the structure equations and depend on the central scalar field (φ_c) and central pressure (p_c). For example for α_A using Eq. (A.1) we simply obtain: $$\begin{array}{l}{\alpha }_A={\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}{m}_A}{\partial {\phi }_0}|}_{{\overline{m}}_A}=\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}{m}_A}{\partial {\phi }_{\text{c}}}\frac{\partial {\overline{m}}_A}{\partial p_{\text{c}}}-\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}{m}_A}{\partial p_{\text{c}}}\frac{\partial {\overline{m}}_A}{\partial {\phi }_{\text{c}}}\right)/\\ \left(\frac{\partial {\phi }_0}{\partial {\phi }_{\text{c}}}\frac{\partial {\overline{m}}_A}{\partial p_{\text{c}}}-\frac{\partial {\phi }_0}{\partial p_{\text{c}}}\frac{\partial {\overline{m}}_A}{\partial {\phi }_{\text{c}}}\right),\end{array}$$(A.2)

where every simple derivative is calculated with the machine-precision due to AutoDiff.

As a result, the total precision of calculated gravitational form factors is determined only by the accuracy of the differential equation solver and may be easily enhanced. We use the Julia library DIFFERENTIALEQUATIONS.JL which allows us to use AutoDiff on the results of the numerical integration. The numerical integration can take advantage of AutoDiff itself because the Jacobian of the system can be calculated more accurately and used to integrate the ODE system. We obtain a higher precision if we select a finer adaptive step in the radial variable r. For our grids we set the integrator to the relative error of δ_Intrel ~ 10⁻¹¹. The obtained relative precision for gravitational form factors typically lies in δ_rel ~ 10⁻⁷ − 10⁻¹³ range, which is more precise than previous calculations.

A.2 Pre-calculated grids of neutron stars

To place limits on the DEF gravity parameters {α₀,β₀} we need to know the gravitational form factors {α_A,β_A,k_A} in Eq. (6) as functions of the mass m_A. To approximate these relations we have to calculate grids of NSs. We calculate grids for several EOSs, including MPA1 used in this paper. Each grid contains 4 calculated values: masses mA and gravitational form factors {α_A, β_A,k_A} for a range of values of {α₀,β₀, p_c} with number of points in each axis respectively {101, 351, 121}. The grid properties are shown in the Table A.1. α₀ and p_c are selected to be logarithmically distributed. Whereas, the parameter β₀ is piece-wise linearly distributed. The grid for p_c is more dense near the common masses of NSs 1.4–2.0 M_⊙. We also include more points in the region of spontaneous scalarisation with β₀ ∈ [−5.0, −4.0]. The resulting calculated gravitational form factors are saved in files which then can be used in DEF gravity tests, the DDSTG model in particular.

With the developed DDSTG model we supply grids for 11 different EOSs ranging from soft to stiff. Each EOS satisfies the requirement on a maximal mass in GR: $M_{\text{NS},\max }^{\text{GR}}>2M_{\odot }$ (Antoniadis et al. 2013; Fonseca et al. 2021). Table B.1 presents the selected EOSs with their maximum masses in GR. In our work we use the piece-wise polytropic approximation described in Read et al. (2009).

Appendix B PK parameters in DEF gravity

In this Appendix we present the expressions for the PK parameters in DEF gravity as functions of the masses of pulsar (m_p) and companion (m_c), and the three Keplerian timing parameters P_b, e, and x. The formulae for the PK parameters {γ, k, r, s} can be found in Damour & Esposito-Farèse (1992a) and Damour & Esposito-Farèse (1996). For the PK parameters {δ_r,δ_θ, A, B} there are only more general theory-independent presentations in Damour & Taylor (1992), from which the expressions in DEF gravity, however, can be derived straightforwardly. The total set of PK parameters in DEF gravity reads: $$\gamma =\frac{e}{n}\frac{X_{\text{c}}{\beta }_O^2}{1+{\alpha }_{\text{p}}{\alpha }_{\text{c}}}\left[X_{\text{c}}\left(1+{\alpha }_{\text{p}}{\alpha }_{\text{c}}\right)+1+k_{\text{p}}{\alpha }_{\text{c}}\right],$$(B.1) $$k=\frac{3{\beta }_O^2}{1-{\text{e}}^2}\left[\frac{1-\frac{1}{3}{\alpha }_{\text{p}}{\alpha }_{\text{c}}}{1+{\alpha }_{\text{p}}{\alpha }_{\text{c}}}-\frac{X_{\text{p}}{\beta }_{\text{c}}{\alpha }_{\text{p}}^2+X_{\text{c}}{\beta }_{\text{p}}{\alpha }_{\text{c}}^2}{6{\left(1+{\alpha }_{\text{p}}{\alpha }_{\text{c}}\right)}^2}\right],$$(B.2) $$r=\frac{G_{*}{m}_{\text{c}}}{c^3},$$(B.3) $$s=\frac{nx}{X_{\text{c}}{\beta }_O}\left[1+\frac{1}{3}\left(\frac{9-{\alpha }_{\text{p}}{\alpha }_{\text{c}}}{1+{\alpha }_{\text{p}}{\alpha }_{\text{c}}}-X_{\text{p}}{X}_{\text{c}}\right){\beta }_O^2\right],$$(B.4) $$\begin{array}{l}{\delta }_r=\frac{{\beta }_O^2}{m_{\text{tot}}^2\left(1+{\alpha }_{\text{p}}{\alpha }_{\text{c}}\right)}[\left(3-{\alpha }_{\text{p}}{\alpha }_{\text{c}}\right)m_{\text{p}}^2\\ +\left(6-{\alpha }_{\text{p}}{\alpha }_{\text{c}}-k_{\text{p}}{\alpha }_{\text{c}}\right)m_{\text{p}}{m}_{\text{c}}+\left(2-{\alpha }_{\text{p}}{\alpha }_{\text{c}}-k_{\text{p}}{\alpha }_{\text{c}}\right)m_{\text{c}}^2],\end{array}$$(B.5) $$\begin{array}{l}{\delta }_{\theta }=\frac{{\beta }_O^2}{m_{\text{tot}}^2\left(1+{\alpha }_{\text{p}}{\alpha }_{\text{c}}\right)}[\left(\frac{7}{2}-\frac{1}{2}{\alpha }_{\text{p}}{\alpha }_{\text{c}}\right)m_{\text{p}}^2\\ +\left(6-{\alpha }_{\text{p}}{\alpha }_{\text{c}}-k_{\text{p}}{\alpha }_{\text{c}}\right)m_{\text{p}}{m}_{\text{c}}+\left(2-{\alpha }_{\text{p}}{\alpha }_{\text{c}}-k_{\text{p}}{\alpha }_{\text{c}}\right)m_{\text{c}}^2],\end{array}$$(B.6) $$A=-\frac{{\beta }_O{X}_{\text{c}}}{2\pi v_{\text{p}}{\left(1-{\text{e}}^2\right)}^{1/2}}\frac{\sin \eta }{\sin \lambda },$$(B.7) $$B=-\frac{{\beta }_O{X}_{\text{c}}}{2\pi v_{\text{p}}{\left(1-{\text{e}}^2\right)}^{1/2}}\frac{\cos i\cos \eta }{\sin \lambda },$$(B.8) $${\dot{P}}_{\text{b}}={\dot{P}}_{\text{b}}^{\phi ,\text{ mon }}+{\dot{P}}_{\text{b}}^{\phi ,\text{ dip }}+{\dot{P}}_{\text{b}}^{\phi ,\text{ quad }}+{\dot{P}}_{\text{b}}^{g,\text{ quad }}$$(B.9)

where we used variables β_o = (G_pcMN/c³)^1/3, G_pc = G_∗(1 + α_p α_c), m_tot = m_p + m_c, X_p = m_p/m_tot, and X_c = m_c/m_tot = 1 − X_p. Moreover, n = 2π/P_b is the orbital circular frequency and ν_p = 1/P is the pulsar’s rotational frequency. The angle i is the inclination of the orbital plane with respect to the plane of sky, while angles λ and η are the polar angles of the spin axis. In the special case of alignment between the orbital and spin axes λ = i and η = −π/2 forcing B = 0.

The expression for Ṗ_b is composed from four different contributions (Damour & Esposito-Farèse 1992a) $$\begin{array}{l}{\dot{P}}_{\text{b}}^{\phi ,\text{ mon }}=-\frac{3\pi X_{\text{p}}{X}_{\text{c}}}{1+{\alpha }_{\text{p}}{\alpha }_{\text{c}}}{\beta }_O^5\frac{{\text{e}}^2\left(1+{\text{e}}^2/4\right)}{{\left(1-{\text{e}}^2\right)}^{7/2}}\\ \times {\left[{\alpha }_{\text{p}}+{\alpha }_{\text{c}}+\frac{2}{3}\left({\alpha }_{\text{p}}{X}_{\text{c}}+{\alpha }_{\text{c}}{X}_{\text{p}}\right)+\frac{{\beta }_{\text{p}}{\alpha }_{\text{c}}+{\beta }_{\text{c}}{\alpha }_{\text{p}}}{1+{\alpha }_{\text{p}}{\alpha }_{\text{c}}}\right]}^2,\end{array}$$(B.10) $$\begin{array}{l}{\dot{P}}_{\text{b}}^{\phi ,\text{ dip }}=-\frac{2\pi X_{\text{p}}{X}_{\text{c}}}{1+{\alpha }_{\text{p}}{\alpha }_{\text{c}}}{\beta }_O^3\frac{1+{\text{e}}^2/2}{{\left(1-{\text{e}}^2\right)}^{5/2}}{\left({\alpha }_{\text{p}}-{\alpha }_{\text{c}}\right)}^2\\ -\frac{4\pi X_{\text{p}}{X}_{\text{c}}}{1+{\alpha }_{\text{p}}{\alpha }_{\text{c}}}{\beta }_O^5\frac{1}{{\left(1-{\text{e}}^2\right)}^{7/2}}[\frac{8}{5}\left(1+\frac{31{\text{e}}^2}{8}+\frac{19{\text{e}}^4}{32}\right)\left({\alpha }_{\text{p}}-{\alpha }_{\text{c}}\right)\\ \times \left({\alpha }_{\text{p}}{X}_{\text{p}}+{\alpha }_{\text{c}}{X}_{\text{c}}\right)\left(X_{\text{p}}-X_{\text{c}}\right)+\left(1+3{\text{e}}^2+\frac{3{\text{e}}^4}{8}\right)\\ \times \frac{\left({\alpha }_{\text{p}}-{\alpha }_{\text{c}}\right)\left({\beta }_{\text{c}}{\alpha }_{\text{p}}{X}_{\text{p}}-{\beta }_{\text{p}}{\alpha }_{\text{c}}{X}_{\text{c}}\right)}{1+{\alpha }_{\text{p}}{\alpha }_{\text{c}}}],\end{array}$$(B.11) $$\begin{array}{l}{\dot{P}}_{\text{b}}^{\phi ,\text{ quad }}=-\frac{32\pi X_{\text{p}}{X}_{\text{c}}}{5\left(1+{\alpha }_{\text{p}}{\alpha }_{\text{c}}\right)}{\beta }_O^5\frac{1+73{\text{e}}^2/24+37{\text{e}}^4/96}{{\left(1-{\text{e}}^2\right)}^{7/2}}\\ \times {\left({\alpha }_{\text{p}}{X}_{\text{c}}+{\alpha }_{\text{c}}{X}_{\text{p}}\right)}^2,\end{array}$$(B.12) $${\dot{P}}_{\text{b}}^{g,\text{ quad }}=-\frac{192\pi X_{\text{p}}{X}_{\text{c}}}{5\left(1+{\alpha }_{\text{p}}{\alpha }_{\text{c}}\right)}{\beta }_O^5\frac{1+73{\text{e}}^2/24+37{\text{e}}^4/96}{{\left(1-{\text{e}}^2\right)}^{7/2}},$$(B.13)

where scalar field φ gives rise to a monopolar ${\dot{P}}_{\text{b}}^{\phi ,\text{ mon }}$ dipolar ${\dot{P}}_{\text{b}}^{\phi ,\text{ dip }}$ and a quadrupolar ${\dot{P}}_{\text{b}}^{\phi ,\text{ quad }}$ contribution. Another quadrupo-lar term ${\dot{P}}_{\text{b}}^{ɡ,\text{ quad }}$ is associated with metric ${ɡ}_{\mu v}^{\ast }$. Following Damour & Esposito-Farèse (1992a), we have ignored a term of order ${\left({\alpha }_{\text{p}}-{\alpha }_{\text{c}}\right)}^2{\beta }_O^5$ in Eq.(B.11).

Appendix C DDSTG implementation into TEMPO

As part of this work, we implement a new DDSTG timing model into an independent version of the standard timing software TEMPO⁷. The new model is based on the DD (Damour & Deruelle 1986) and DDGR (Taylor & Weisberg 1989) models but modified to work within STG theories. In addition to the parameters present in the DDGR model, the DDSTG model utilises four additional parameters, added to the corresponding parameter input file. New parameters and their possible values are presented in Table C.1.

The first two parameters ‘ALPHA0’ and ‘BETA0’ correspond to the values of two arbitrary DEF gravity parameters {α₀, β₀}. They can be selected to be zero to recover GR. Otherwise, the applicable range of values is determined by the range of the supplementary data files, containing grids of the calculated gravitational form factors and masses. The grids are calculated as part of this work and supplied with the model. The supplied grids cover α₀ ∈ [−10⁻¹, −10⁻⁵] and β₀ ∈ [−6, +10] in the DEF gravity parameter space.

The third parameter ‘EOS’ selects the equation of state to be used to perform the test. In the present work, we generally select the rather stiff EOS MPA1. TEMPO reads the supplied pre-calculated grids of the gravitational form factors and masses corresponding to the selected EOS (see also Table B.1).

The last parameter ‘COMP_TYPE’ is used to select the type of the companion to the pulsar, among a white dwarf (‘WD’), a neutron star (‘NS’) and a black hole (‘BH’). These three compact objects have different properties in terms of gravitational form factors. NSs in DEF gravity can be moderately |α_NS| ≲ |α₀| or strongly |α_NS| ~ O(1) scalarised depending on the mass m_A and the chosen DEF gravity parameters. White dwarfs are not enough compact to significantly deviate from the weak-field approximation |α_WD| ≲ |α₀|. Thus for WDs we apply α_WD = α₀, β_WD = β₀.¹¹ In contrast, BHs in DEF gravity fulfil a no-hair theorem and they are completely descalarised α_BH = 0 and β_BH = 0. The choice of the companion type plays a significant role for the calculated values of PK parameters.

During initialisation, TEMPO reads the parameter file and the file with the TOAs. When the {α₀,β₀} values and a EOS are provided, it reads the supplementary grids for the selected EOS. The supplied grids are 3D and consist of the masses of NSs m_A and their gravitational form factors {α_A, β_A, k_A}. Two axes are {α₀,β₀} and the third axis is a fixed range of central pressure (p_c) of a NS.

Bilinear interpolation is performed in the {α₀,β₀} parameter space for the selected parameters. The result of the interpolation are four 1D arrays for masses and gravitational form factors depending on a range of central pressures. Interpolated 1D arrays are saved in the computer memory and used further without being changed. From this point gravitational form factors are ready to be interpolated in the mass range during the fitting procedure.

TEMPO with the selected DDSTG model fits two masses in addition to spin, astrometric and Keplerian parameters. The masses are iteratively changed during the fitting procedure. Only when the masses change TEMPO interpolates the gravitational form factors from the pre-calculated 1D grids. Gravitational form factors are calculated for both the pulsar and the companion and are kept constant for all TOAs within a single iteration.

TEMPO tries to find the companion mass (m_c) and the total mass (m_tot = m_p + m_c) corresponding to the best fit. During the fitting, it utilises the derivatives of the timing formula (7) with respect to these two masses. We calculate derivatives separately for each delay in Equations (11–14). We used the approximated expressions for the derivatives of the delays with respect to PK parameters provided by Damour & Deruelle (1986). After applying the chain rule for our derivatives we come to expressions depending on the derivatives of PK parameters with respect to the two masses. The model calculates the derivatives of the PK parameters of DEF gravity (Eqs. B.1 – B.9) with respect to the two masses of the timing model, $${\left(\frac{\partial p_i^{\text{PK}}\left(m_{\text{c}},m_{\text{tot}}\right)}{\partial m_{\text{c}}}\right)|}_{m_{\text{tot}}},{\left(\frac{\partial p_i^{\text{PK}}\left(m_{\text{c}},m_{\text{tot}}\right)}{\partial m_{\text{tot}}}\right)|}_{m_{\text{c}}},$$(C.1)

which are certainly different to those in GR. Moreover, for proper convergence of the model in highly non-linear areas the additional corrections to the derivatives related to variation of gravitational form factors with respect to masses (i.e. ∂α_A/∂m_A, ∂β_A/∂m_A) are required.

The formulae in Eq. (C.1) are large but can be obtained by straightforward differentiation of the PK formulae in Appendix B. The expressions are incorporated in the model and depend on gravitational form factors, masses and Keplerian parameters which are known on each iteration. Once the derivatives of PK parameters are calculated, the model calculates the derivatives of the timing delays with respect to the two masses and proceeds the iteration. The outcome of the DDSTG model in TEMPO is the same as for the DDGR model. The most important result is the χ² value, which is used further to perform tests of gravity.

The DDSTG timing model converges well if the eccentricity of the orbit is not too small and TEMPO is supplied with a reasonable initial guess of the orbital parameters. The model can run into difficulties in the convergence for solutions if those parameters are far away from GR or if the orbit is almost circular. The problem can happen if TEMPO changes masses significantly during the iteration and reaches a point where the Shapiro shape parameter calculated with new masses using Eq. (B.4) becomes greater than one, i.e. s > 1. In this case the expression of the Shapiro delay from Eq. (13) does not have a physical meaning and TEMPO breaks with an error for TOAs near conjunction, as the argument of the logarithm becomes negative. The convergence can be enforced by applying the special ‘GAIN’ parameter, which however comes at the cost of increased computational time. A GAIN parameter of less than one forces TEMPO to use smaller steps while iterating. Another way which ensures the convergence is to supply TEMPO with a good initial guess for the masses in DEF gravity. The rough estimate of the masses in DEF gravity can be obtained from the traditional PK method. The combination of these two methods helps to converge systems with a small eccentricity and Shapiro shape parameter near one (s ≃ 1) even in the strongly non-linear part of the {α₀,β₀} plane.

Appendix D Mass-mass diagram near the horn

Many binary pulsar tests are not sensitive at high |α₀| values and β₀~ − 2. In this region near the ‘horn’, α_NS ≈ α₀ and the overall ${\dot{P}}_{\text{b}}$ is close to the GR value, as the dipolar contribution is greatly suppressed. As a consequence, the test is passed. In Fig. D.1 we show a mass-mass diagram in DEF gravity with α₀ = −0.1, β₀ = −1.9. This point in the DEF gravity parameter space passes the test with timing data but is excluded by Cassini experiment (|α₀| ≲ 3.4 × 10⁻³). The picture is dramatically different from what we see in the zone of scalarisation in Fig. 1.

Fig. D.1 Mass-mass diagram in DEF gravity corresponding to a point near the ‘horn’ with α0 = −10−1, β0 = −1.9. The performed test is the same as for Fig. 1. The shadowed area is the allowed region at 68% CL limit for a corresponding PK parameter. The solid green line corresponds to the observed value of ${\dot{P}}_{\text{b}}^{GW}$.

Appendix E Proper account of ${\dot{P}}_b^{ext}$ uncertainty

If the observed ${\dot{P}}_{\text{b}}$ is measured more precisely than the external contribution ${\dot{P}}_{\text{b}}^{\text{ext}}$ we have to take the uncertainty of the external contribution into account. In most situations ${\dot{P}}_{\text{b}}^{\text{ext}}$ consists of the Shklovskii and the Galactic contribution. The calculation of the Shklovskii contribution is limited by the uncertainty in the distance and the one in the proper motion. The uncertainty in the estimation of the Galactic contribution is determined by the uncertainty in the distance and our imperfect knowledge of the Galactic gravitational potential. The latter is partly systematic in nature and therefore somewhat more difficult to quantify. If the uncertainty $\text{Δ}{\dot{P}}_{\text{b}}^{\text{ext}}$ is the limiting factor we no longer can have ${\dot{P}}_{\text{b}}^{\text{ext}}$ fixed on one value for the whole experiment (i.e. provide a fixed value for XPBDOT parameter in TEMPO) because its uncertainty can affect the derived limits.

To properly account for the uncertainty in ${\dot{P}}_{\text{b}}^{\text{ext}}$ we follow the method described in Splaver et al. (2002) based on Bayes inference. On the first stage TEMPO calculates χ² grid for three independent parameters $\left\{{\alpha }_0,{\beta }_0,{\dot{P}}_{\text{b}}^{\text{ext}}\right\}$. Then the shifted value $\text{Δ}{\chi }^2={\chi }^2-{\chi }_{\min }^2$ is described by χ² distribution with 3 degrees of freedom. The number of degrees of freedom is important when we calculate confidence level limits. The ∆χ² maps as usual to a Bayesian likelihood function, $$p\left(\left\{t_j\right\}\mid {\alpha }_0,{\beta }_0,{\dot{P}}_{\text{b}}^{\text{ext }}\right)=\frac{1}{2}{\text{e}}^{-\text{Δ}{\chi }^2/2},$$(E.1)

where {t_j} refers to the used timing data. We treat the obtained probability density $p\left(\left\{t_j\right\}\mid {\alpha }_0,{\beta }_0,{\dot{P}}_{\text{b}}^{\text{ext}}\right)$ as the likelihood for Bayes’ theorem $$p\left({\alpha }_0,{\beta }_0,{\dot{P}}_{\text{b}}^{\text{ext}}\mid \left\{t_j\right\}\right)=\frac{p\left(\left\{t_j\right\}\mid {\alpha }_0,{\beta }_0,{\dot{P}}_{\text{b}}^{\text{ext}}\right)}{p\left(\left\{t_j\right\}\right)}p\left({\alpha }_0,{\beta }_0,{\dot{P}}_{\text{b}}^{\text{ext}}\right),$$(E.2)

where $p\left({\alpha }_0,{\beta }_0,{\dot{P}}_{\text{b}}^{\text{ext}}\mid \left\{t_j\right\}\right)$ is the desired joint posterior probability function and $p\left({\alpha }_0,{\beta }_0,{\dot{P}}_{\text{b}}^{\text{ext}}\right)$ is the prior. The Bayesian evidence p({t_j}) is obtained from normalisation of posterior probability over the whole grid {α₀,β₀}.

In the next stage we select the prior $p\left({\alpha }_0,{\beta }_0,{\dot{P}}_{\text{b}}^{\text{ext}}\right)$. We assume no prior information about a0 and β₀, so their distributions are taken uniformly for β₀ ∈ (−∞, +∞) and α₀ ∈ (−∞, 0]. At this step, we utilise the information about the uncertainties in ${\dot{P}}_{\text{b}}^{\text{ext}}$. For example, we can select the prior $p\left({\dot{P}}_{\text{b}}^{\text{ext}}\right)$ to be normally distributed with the mean measured value ${\dot{P}}_{\text{b}}^{\text{ext}}$ and its standard deviation $\text{Δ}{\dot{P}}_{\text{b}}^{\text{ext}}$. Or we can use more elaborate probability distribution accounting for our uncertainty in the measurements of distance, proper motion, and Galactic gravitational potential, $$p\left({\alpha }_0,{\beta }_0,{\dot{P}}_{\text{b}}^{\text{ext }}\right)=p\left({\alpha }_0\right)\times p\left({\beta }_0\right)\times p\left({\dot{P}}_{\text{b}}^{\text{ext}}\right).$$(E.3)

The last step is to marginalise over the ${\dot{P}}_{\text{b}}^{\text{ext}}$ variable using the correct joint posterior probability with desired priors $$p\left({\alpha }_0,{\beta }_0\mid \left\{t_j\right\}\right)={\displaystyle \int d}{\dot{P}}_{\text{b}}^{\text{ext }}p\left({\alpha }_0,{\beta }_0,{\dot{P}}_{\text{b}}^{\text{ext }}\mid \left\{t_j\right\}\right),$$(E.4)

where the integral is calculated on the grid of ${\dot{P}}_{\text{b}}^{\text{ext}}$ parameter. The marginalised probability than can be used to derive limits, for instance going back to χ² representation with 2 degrees of freedom.

References