Binary Pulsar constraints on massless scalar-tensor theories using Bayesian statistics

TL;DR

This paper uses Bayesian statistics and pulsar timing data to tightly constrain massless scalar-tensor theories of gravity, highlighting the importance of orbital decay measurements.

Contribution

It compares Bayesian and traditional methods for constraining scalar-tensor theories using binary pulsar data, demonstrating Bayesian analysis's ability to explore parameter posteriors.

Findings

Both methods yield similar constraints on theories.

Orbital period decay measurements dominate the constraints.

Bayesian approach provides posterior distributions of parameters.

Abstract

Binary pulsars provide some of the tightest current constraints on modified theories of gravity and these constraints will only get tighter as radio astronomers continue timing these systems. These binary pulsars are particularly good at constraining scalar-tensor theories in which gravity is mediated by a scalar field in addition to the metric tensor. Scalar-tensor theories can predict large deviations from General Relativity due to the fact that they allow for violation of the strong-equivalence principle through a phenomenon known as scalarization. This effect appears directly in the timing model for binary pulsars, and as such, it can be tightly constrained through precise timing. In this paper, we investigate these constraints for two scalar-tensor theories and a large set of realistic equations of state. We calculate the constraints that can be placed by saturating the current $1\sigma$ bounds on single post-Keplerian parameters, as well as employing Bayesian methods through Markov-Chain-Monte-Carlo simulations to explore the constraints that can be achieved when one considers all measured parameters simultaneously. Our results demonstrate that both methods are able to place similar constraints and that they are both indeed dominated by the measurements of the orbital period decay. The Bayesian approach, however, allows one to simultaneously explore the posterior distributions of not only the theory parameters but of the masses as well.

Tables

Table 1. Table 1: The Keplerian and post-Keplerian parameters measured for the systems we consider in this paper. Values with and asterisks are derived quantities assuming GR as the underlying theory. The values in parentheses represent the 1 ​ σ 1 𝜎 1\sigma errors associated with each quantity. Note that the Shaprio parameter r 𝑟 r is measured in units of T ⊙ = G ​ M ⊙ / c 3 = 4.925490947 ​ μ subscript 𝑇 direct-product 𝐺 subscript 𝑀 direct-product superscript 𝑐 3 4.925490947 𝜇 T_{\odot}=GM_{\odot}/c^{3}=4.925490947\mu s for J2222-0137.

	J1738+0333	J1012+5307	J2222-0137	J1909-3744
$P_b$ (d)	0.3547907398724(13)	0.60467271355(3)	2.44576469(13)	1.533449474406(13)
$x$ (s)	0.343429130(17)	0.5818172(2)	10.8480239(6)	1.89799118
$e$	$3.5\left(1.1\right)\times {10}^{-7}$	$1.2\left(3\right)\times {10}^{-6}$	0.000380967(30)	$1.14\left(10\right)\times {10}^{-7}$
$m_A$ ($M_{\odot }$)	1.47(7)*	1.64(22)*	1.76(6)*	1.47(3)*
$m_B$ ($M_{\odot }$)	0.181(8)	0.16(2)	1.293(25)*	0.208(2)
$q\equiv m_A/m_B$	8.1(2)	10.5(5)	—	—
${\dot{P}}_b^{\text{int}}$ (fs s-1)	-25.9(3.2)	-15(15)	-60(90)	—
$\dot{\omega }$ ($\mathrm{deg}$ yr-1)	—	—	0.1033(29)	—
$s\equiv \sin \iota$	—	—	0.99559	0.99771
$r$ ($T_{\odot }$)	—	—	1.293(25)	—

Equations

$$S=\int \frac{d^{4}x}{c}\frac{\sqrt{-g}}{4\kappa }\left[R-2g^{\mu \nu }{\partial }_{\mu }\phi {\partial }_{\nu }\phi \right]+S_{\text{mbox}mat}\left[\chi ,A^2(\phi )g_{\mu \nu }\right],$$

$$R_{\mu \nu }$$

$$\square \phi$$

$$T_{\mu \nu }^{\text{mbox}mat}\equiv \frac{2c}{\sqrt{-g}}\left(\frac{\delta S_m}{\delta g^{\mu \nu }}\right),$$

$$\alpha (\phi )=\frac{\partial \ln A(\phi )}{\partial \phi },$$

$$\beta (\phi )=\frac{\partial \alpha (\phi )}{\partial \phi },$$

$$A(\phi )$$

$$\alpha (\phi )$$

$$\beta (\phi )$$

$$A(\phi )$$

$$\alpha (\phi )$$

$$\beta (\phi )$$

$${\alpha }_{\infty }^{\text{mbox}DEF}$$

$${\beta }_{\infty }^{\text{mbox}DEF}$$

$${\alpha }_{\infty }^{\text{mbox}MO}$$

$${\beta }_{\infty }^{\text{mbox}MO}$$

$$G_N=G\left[A_{\infty }^2\left(1+{\alpha }_{\infty }^2\right)\right],$$

$$\bar{\gamma }=\left|1-{\gamma }_{\text{mbox}PPN}\right|=\frac{2{\alpha }_{\infty }^2}{1+{\alpha }_{\infty }^2},$$

$$\bar{\beta }=\left|1-{\beta }_{\text{mbox}PPN}\right|=\frac{|{\beta }_{\infty }|{\alpha }_{\infty }^2}{2{\left(1+{\alpha }_{\infty }^2\right)}^2},$$

$${\alpha }_{\infty }^2\lesssim \frac{\bar{\gamma }}{2-\bar{\gamma }},$$

$${\alpha }_{\infty }^2({\beta }_0)\lesssim \frac{|{\beta }_{\infty }|}{\bar{\beta }}-1-\frac{1}{4\bar{\beta }}\sqrt{|{\beta }_{\infty }|(|{\beta }_{\infty }|-8\bar{\beta }}),$$

$${\alpha }_A={\frac{\partial \ln m_A}{\partial {\phi }_{\infty }}|}_{{\bar{m}}_A},$$

$${\beta }_A={\frac{\partial {\alpha }_A}{\partial {\phi }_{\infty }}|}_{{\bar{m}}_A},$$

$$k_A={\frac{\partial \ln I_A}{\partial {\phi }_{\infty }}|}_{{\bar{m}}_A},$$

$$D{\tau }_a=T_e+{\Delta }_R+{\Delta }_E+{\Delta }_S+\mathcal{O}(c^{-4}),$$

$${\Delta }_E=\gamma \sin u,$$

$$u-e\sin u=n\left[(T_e-T_0)-\frac{{\dot{P}}_b}{2p_b}(T_e-T_0{)}^2\right].$$

$$x=x_0+\dot{x}(T_e-T_0),$$

$$e=e_0+\dot{e}(T_e-T_0),$$

$$\omega ={\omega }_0+\dot{\omega }(T_e-T_0),$$

$$\gamma$$

$$\dot{\omega }$$

$$r$$

$$s$$

$${\dot{P}}_b={\dot{P}}_b^{\text{mbox}Acc}+{\dot{P}}_b^{\text{mbox}Shk}+{\dot{P}}_b^{\dot{M}}+{\dot{P}}_b^T+{\dot{P}}_b^{\dot{G}}+{\dot{P}}_b^{\text{mbox}int}.$$

$${\dot{P}}_b^{\text{mbox}int}={\dot{P}}_b^{\phi ,\text{mbox}mon}+{\dot{P}}_b^{\phi ,\text{mbox}dip}+{\dot{P}}_b^{\phi ,\text{mbox}quad}+{\dot{P}}_b^{g,\text{mbox}quad}.$$

Full text

I Introduction

General Relativity (GR) has been the most successful theory of gravity over the last century, passing all current tests with great success. Our ability to test gravitational theories is becoming even stronger Yunes and Siemens (2013); Yunes et al. (2016) with the recent observations of gravitational waves (GW) Abbott et al. (2016, 2017a, 2017b) and the continued monitoring of binary pulsars (PSRs) Freire et al. (2012); Kramer and Wex (2009); Kramer (2016). As more neutron star (NS) systems are discovered through gravitational wave observations we will be able to pin down the nuclear equation of state (EOS) and other NS properties Abbott et al. (2018a, b). Such observations can help break the degeneracy that arise between EOS effects and those introduced by modified theories of gravity.

Even with the success of GR, there are many theoretical reasons to consider modified theories. Some of the most well-motivated modified theories are scalar-tensor theories (STTs) of gravity, in which gravity is mediated by the metric tensor and a dynamical scalar field. These theories are defined by a choice of conformal coupling function that dictates the manner in which the scalar field couples to matter, and ultimately the level at which the strong-equivalence principle (SEP) is violated. Originally proposed by Jordan Jordan (1955, 1959), Fierz Fierz (1956), Brans Brans and Dicke (1961), and Dicke (JFBD) as some of the most natural alternatives to GR, STTs were later extended by Damour and Esposito-Farése (DEF) to include higher order terms in the conformal coupling, as well as multiple scalar fields Damour and Esposito-Farese (1992). Recently, a slight variation on this theory, proposed by Mendes and Ortiz (MO) Mendes and Ortiz (2016), has gained some attention as it arises from more fundamental considerations.

While Solar System observations have the ability to tightly constrain STTs, these weak field constraints do not always translate to tight restrictions in the strong field regime. In particular, STTs give rise to a phenomenon in NSs known as scalarization in which the scalar field can become excited far above its background (weak field) value. Neutron stars in STTs can acquire a so-called scalar charge which quantifies the $1/r$ behavior of the scalar field in a far field expansion from the NS. Such modifications to NS spacetimes directly affect observables, particularly those that can be probed with binary PSR observations. The scalar charges that NSs can develop appear directly in the timing model that is used to predict when we should observe pulses from binary PSRs. Astronomers are able to time PSRs so precisely that even slight deviations from the GR predictions are highly constrained, thus providing some of the best constraints on STTs.

Typically, constraints on STTs are placed through the observed rate of decay of the orbital period of the binary. STTs predict dipolar gravitational radiation, which enters at lower post Newtonian (PN) order relative to the typical quadrupole radiation predicted by GR, and thus, speeds up the orbital decay rate. Since observations seem to suggest the absence of this extra dipole effect, STTs can be constrained from the observation of the orbital decay rate with precision that can exceed that of Solar System observations. However, STTs predict that all post-Keplerian (PK) parameters Damour and Taylor (1992) are modified from the GR prediction. Thus, observations of post-Keplerian parameters other than the orbital period decay can in principle be combined to allow for even tighter constraints.

In this paper we investigate the type of constraints that can be placed on STTs, particularly the theories proposed by DEF and MO, by using multiple post-Keplerian parameters, as well as multiple PSR systems. We use a combination of PSR-white-dwarf systems and PSR-NS systems with a wide range of PSR masses to explore the effects of scalarization and place constraints in a self-consistent manner. We also explore the effects of the NS EOS and how it affects the strength of such constraints. Such studies are carried out with Bayesian methods through Markov-Chain Monte-Carlo simulations that explore the full parameter space and determine which values of the STT parameters are most consistent with observations. This approach allows us to accurately take into account any correlations between observed post-Keplerian parameters and ensures that the constraints we place are self-consistent.

The study and results summarized above are a comprehensive extension of other work that also employed a Markov-Chain Monte-Carlo methodology for binary pulsar tests of GR. In 2017, Shao et al. (2017) carried out a Markov-Chain Monte-Carlo exploration of the posterior probability distribution of the two main STT parameters that characterize DEF theory. This analysis used 5 PSRs that currently have good measurements of the component masses (either through optical counterparts or through measurements of the PPK parameters associated with the Shapiro time delay) and of one additional PPK parameter (the orbital period decay). A Bayesian analysis was then employed to obtain posteriors for the STT parameters for 11 different EOSs, each marginalized over the masses.

Our work differs from this analysis in various ways. First, we do not marginalize only over the component masses, but also over the equations of state (treating the latter as a discrete parameter). Second, we analyze binaries with measurements of any PPK parameter and not just the orbital period decay, which is possible due to a prior deep exploration of scalar charges Anderson . Third, we greatly extend the exploration of the STT parameter space, not only considering the region where spontaneous scalarization activates. Fourth, we do not only focus on the DEF theory, which is already stringently constrained by cosmological and Solar System observations Sampson et al. (2014); Anderson et al. (2016); Anderson and Yunes (2017), but we also study the MO theory.

Using the data for scalar charges provided in Anderson , we place constraints on DEF and MO theory using all measured post-Keplerian parameters available for a set of 7 PSRs with 11 different EOSs. We show that the EOS has relatively little impact on these types of parameters and does not change the relative strength of any constraints placed from different PSRs. We find that other post-Keplerian parameters beyond the orbital period decay, like the rate of periastron advance, are, in some cases, able to place tighter constraints on STTs for PSR-NS systems. Finally, when we take a Bayesian approach to place constraints, we find similar results to what one finds when using the approach taken in the past Esposito-Farese (2011); Freire et al. (2012), except that a Bayesian approach allows us to (i) easily stack all constraints from all PSR observations and (ii) easily marginalize over masses and EOSs.

We start with a description of STTs and the timing model in Sec. II in order to lay down the foundation for the rest of the paper. Here we discuss the current constraints on STTs from Solar System observations, as well as the details of the various post-Keplerian parameters appearing in the timing model. In Sec. III we investigate the various constrains that can be placed on STTs, first through the standard method that is employed in the literature, and then through MCMC methods in §IV. We then conclude in Sec. V with a discussion of our results and future work. Throughout this paper, we follow the conventions of Misner et al. (1973).

II Timing binary PSRs in STTs

In this section we will introduce the basics of STTs and the post-Keplerian parameters that appear in the timing model. We will discuss the particular theories we consider in this paper and the current constraints that Solar System observations place on these theories. We give a brief discussion of the scalar charges and how they are calculated, and discuss their importance for timing binary PSRs. For completeness we provide a summary of the timing model and present the various post-Keplerian parameters in the context of STTs.

II.1 Scalar-tensor theories

In this paper we will focus explicitly on massless STTs in which there exists a single scalar field non-minimally coupled to the metric tensor $g_{\mu\nu}$. These theories can be derived from an action in the so called Einstein frame given by Damour and Esposito-Farese (1992, 1993, 1996)

[TABLE]

where $g$ and $R$ are the determinant and Ricci scalar associated with the metric $g_{\mu\nu}$, $\kappa=4\pi G/c^{4}$, $\chi$ are any matter fields, and $A(\varphi)$ is a conformal factor that determines how the scalar field $\varphi$ couples to matter.

The field equations resulting from variation of action above with respect to $g_{\mu\nu}$ and $\varphi$ are given by

[TABLE]

where the stress-energy tensor is defined by

[TABLE]

and $T^{\mbox{\tiny mat}}\equiv g_{\mu\nu}T_{\mu\nu}^{\mbox{\tiny mat}}$ is its trace.

In addition to the $A(\varphi)$ function that appears in the action, it is also convenient to introduce a couple of additional quantities. Let us then define

[TABLE]

which we designate the conformal coupling, as it appears directly in the equations of motion for the scalar field. Another quantity that appears directly in the parameterized-post-Newtonian (PPN) and parameterized-post-Keplerian (PPK) formalisms is the derivative of this quantity, namely,

[TABLE]

and quantifies the non-linear behavior of $A(\varphi)$.

A particular STTs is defined by one’s choice of the conformal factor $A(\varphi)$, or likewise $\alpha(\varphi)$, and will ultimately determine how the scalar field reacts to the presence of matter, or likewise the gravitational potential generated by such matter. In this paper, we consider two models that exhibit similar behavior when $\varphi\ll 1$ but allows for much different behavior elsewhere. The first model we consider is DEF theory Damour and Esposito-Farese (1993, 1992), which is defined by

[TABLE]

where $\beta_{0}$ enters directly as a free parameter. The other model we consider is MO theory Mendes and Ortiz (2016), which is defined by

[TABLE]

where, again, $\beta_{0}$ enters directly as a free parameter.

Both theories are subject to a boundary condition at spatial infinity such that the scalar field has a cosmologically determined background value $\varphi_{\infty}$. This quantity then becomes another free parameter of these theories and must be chosen to satisfy weak-field constraints, like those from Solar System observations. However, we choose the parameterization $\varphi_{\infty}=\alpha_{0}/\beta_{0}$ and let $\alpha_{0}$ become our second free parameter appearing in both theories. A more detailed description of this parameterization of STTs can be found in Refs. Anderson ; Damour (2007). It is worth pointing out, however, the this parameterization enforces certain relations when one inserts $\varphi_{\infty}$ into Eqs. (8)-(9), namely

[TABLE]

in DEF theory, and when using Eqs. (11)-(12)

[TABLE]

in MO theory.

The quantities in Eqs. (13)-(16) appear directly in the weak field predictions made by these theories. The choice of parameters $(\alpha_{0},\,\beta_{0})$ will determine the local value of the gravitational constant via the relation

[TABLE]

as well as the PPN parameters $\gamma_{\mbox{\tiny PPN}}$ and $\beta_{\mbox{\tiny PPN}}$ Will and M (1993); Will (2014). The former is given by

[TABLE]

while the latter is given by

[TABLE]

The $\gamma_{{\mbox{\tiny PPN}}}$ parameter is a measure of the spatial curvature induced by a unit rest mass and has been measured from the Shapiro time delay observed by the Cassini spacecraft Bertotti et al. (2003); Will (2014), being constrained to $\left|1-\gamma_{\mbox{\tiny PPN}}\right|\lesssim 2.3\times 10^{-5}$. The $\beta_{{\mbox{\tiny PPN}}}$ parameter is a measure of the amount of non-linearity in the superposition law for gravity, and it is constrained from observations of the perihelion shift of Mercury Will (2014) to be $\left|1-\beta_{\mbox{\tiny PPN}}\right|\lesssim 8\times 10^{-5}$. From these relations we infer constraints on $\alpha_{\infty}$ and $\beta_{\infty}$, namely

[TABLE]

from Eq. (18), and

[TABLE]

from Eq. (19), which only places tighter constraints for $|\beta_{\infty}|\gtrsim 14$.

II.2 Scalarization of neutron stars

While STTs can be constrained rather tightly from Solar System observations they are still able to produce significant deviations from GR near strongly self-gravitating matter like NSs. These strong field effects are a result of a well-studied phenomenon known as scalarization Damour and Esposito-Farese (1993, 1992, 1996) . When the coupling parameter $\beta_{0}$ is sufficiently negative ($\lesssim-4.3$) the scalar field can grow rapidly inside the NS even when the cosmologically determined background value of the field approaches zero. This is precisely how STTs satisfy Solar System constraints but produce observable deviations from GR in the strong field regime.

The effects of scalarization can be quantified by quantities known as scalar charges and they enter directly into the timing model presented in the next section. There are three scalar charges that appear, the first of which is defined via

[TABLE]

for the $A$th NS of a system and where the derivative must be taken with the baryonic mass held constant. This scalar charge measures the “sensitivity” of the NS’s mass to variations in the background scalar field, and thus, it represents an effective coupling between the NS and the scalar field. This quantity is the strong field counterpart to the weak field parameter $\alpha_{\infty}$. The second relevant scalar charge is the strong field equivalent to $\beta_{\infty}$, and thus, it is a derivative of Eq. (22), i.e.

[TABLE]

where again the baryonic mass must be held constant. This scalar charge encodes non-linear interactions between the binary component. The final scalar charge of interest is linked to the NS’s moment of inertia $I_{A}$, and it is defined as

[TABLE]

with the baryonic mass once more held fixed. Similar to the relations between $\alpha_{A}$ and the NS’s mass, $k_{A}$ quantifies the sensitivity of the NS’s moment of inertia to the scalar field. This quantity becomes most relevant when a NS is in a binary system with another NS. As the NSs orbit one another, they will move through each other’s gravitational potential, effectively altering the local value of the scalar field. The NS’s moment of inertia will then fluctuate throughout the orbit and produce an additional time delay that is measurable in principle.

Calculating these scalar charges is not very difficult but can be numerically expensive over large regions of the $(\alpha_{0},\,\beta_{0})$ parameter space. For the Bayesian methods used in this paper, it would be nearly impossible to compute these quantities as often as they would be needed for millions of likelihood calculations. However, these quantities have recently been calculated over a large region of parameter space and for a wide range of EOS Anderson , consistent with the PSR mass in J0348+0432 and the recent constraints placed from LIGO and VIRGO Abbott et al. (2018a). In this paper, we make use of this data in order to avoid solving for the scalar charges on the fly during our MCMC simulations.

II.3 Timing model in STTs

In order to use binary PSRs to place constraints on a gravitational theory, one must work out how that theory predicts modifications to the motion of a binary system. Therefore, one must develop a timing formula that captures the relativistic effects of the theory by relating the observed time of arrival (TOA) of the pulse to the time the pulse was emitted. Blandford and Teukolsky (BT) Blandford and Teukolsky (1976) originally developed such a formula to explain the, then recent, discovery by Hulse and Taylor of the the first binary PSR B1913+16 Hulse and Taylor (1975). Later, Damour and Deruelle (DD) developed a full 1PN description of the two-body problem in a way that allowed for a parameterized description of the timing model, encapsulating all relativistic effects in GR. This model was then extended by Damour and Taylor Damour and Taylor (1992) in a phenomenological manner such that the DD model could be used generically to constrain any conservative theory of gravity. The parameters of the Damour-Taylor model are the two masses ($m_{A}$ and $m_{B}$) of the binary, the standard Keplerian orbital parameters $\{P_{b},n=2\pi/P_{b},\,T_{0},\,\omega_{0},\,e_{0},\,x_{0}\}$ (orbital period, orbital frequency, time of periastron passage, location of periastron, eccentricity, and projected semi-major axis), and a set of post-Keplerian parameters. In this section we will review the timing formula of DD and present the relevant post-Keplerian parameters of Damour and Taylor in the context of scalar-tensor theories.

The timing model is traditionally (BT and DD) written as

[TABLE]

where $\tau_{a}$ is the infinite frequency barycenter arrival time, $T_{e}$ is the proper time of emission, and $D$ is a Doppler factor accounting for center of mass motion between the binary and the Solar System barycenter. Using $\tau_{a}$ here means that one must have accurately accounted for the time delays associated with dispersion from the interstellar medium and corrections that transform this barycenter arrival time to the time kept at the observatory on Earth’s surface Edwards et al. (2006); Lorimer and Kramer (2005). The Doppler factor $D$ is an inconsequential constant that can be transformed out of the timing formula through a redefinition of units, and restored later if needed Damour and Taylor (1992). The last three terms in Eq. (25) account for all other time delays occurring between the binary and the Solar System barycenter. The Römer delay $\Delta_{R}$ is simply the classical light travel across the binary and depends on the PSR’s position in its orbit. The magnitude of this time delay is determined by the projected semi-major axis along the line of sight.

The Einstein delay $\Delta_{E}$ arises from relating coordinate time of emission to the proper time of emission, and thus, it depends on the metric component $g_{00}$. This time delay takes the form

[TABLE]

where $\gamma$ is a post-Keplerian parameter and $u$ is an eccentric anomaly like variable that is the function of $T_{e}$ found by solving

[TABLE]

with $\dot{P}_{b}$ being another post-Keplerian parameter describing the change in orbital period, which we will discuss later. The coefficient $\gamma$ is determined when one accounts for the gravitational redshift of the companion and the second-order Doppler shift from the PSR’s motion in the line element, using Kepler’s equations to integrate.

There is also an additional effect described by the Einstein delay that is related to violation of the SEP. Theories that violate the SEP will allow the PSR’s moment of inertia to change throughout the orbit as it moves through the gravitational potential of the companion. This means that the PSR’s rotational frequency will vary throughout its orbit in order to conserve angular momentum. Such effects add another time delay to the timing formula and are directly proportional to $\sin u$, and therefore, they are included in the definition of $\gamma$ Will and M (1993); Damour and Taylor (1992).

The final time delay appearing in Eq. (25) is the Shapiro delay $\Delta_{S}$ resulting from the fact that light must travel in the curved background of the binary. This term can be thought of as the first relativistic correction to $\Delta_{R}$ and is most easily measured when one observes the binary edge on. The magnitude of the Shapiro delay is denoted by a post-Keplerian parameter $r$ and is labeled the “range” of the Shapiro delay. Another post-Keplerian parameter, $s\equiv\sin\iota$, representing the “shape” of the Shapiro delay, and enters here to characterize how this time delay is effected by the orbit’s inclination relative to the line of sight.

While Eq. (25) captures the delays associated with the pulse traveling in the curved spacetime between the PSR and the Solar System barycenter, it does not account for any secular variations in the Keplerian parameters of the binary. Thus, in principle, one must also account for three other post-Keplerian parameters, i.e. $\dot{x}$, $\dot{e}$, and $\dot{\omega}$, via

[TABLE]

and insert these back into the timing formula of Eq. (25). In practice, however, the rate of periastron advance $\dot{\omega}$ is usually the only one of these parameters that can be measured consistently and its theoretical value can be determined by the methods in Damour and Deruelle (1985); Damour and Taylor (1992).

Thus far, we have worked in a theory independent framework and we have parameterized the various time delays with a set of measurable post-Keplerian parameters $\{\dot{\omega},\,\gamma,\,r,\,s\}$. In GR, these four parameters take on simple forms that are a function of the two masses and the Keplerian parameters. In other theories of gravity, however, one must determine their functional form using the formalism in Damour and Taylor (1992). Thus, in the context of STTs Damour and Taylor (1992); Damour and Esposito-Farese (1992); Damour (2007); Horbatsch and Burgess (2012), the post-Keplerian parameters that are typically measured are given by

[TABLE]

in which we have made use of the notation in Damour and Taylor (1992), where $M$ is the total mass and $X_{A,B}=m_{A,B}/M$. One will notice that the scalar charge $(\alpha_{A},\,\beta_{A},\,k_{A})$ described earlier appear directly in these definitions.

The effects we have incorporated thus far do not capture radiative effects such as gravitational wave emission. The orbital period derivative, $\dot{P}_{b}$, accounts for these radiative effects, along with others, and in general has the contributions Damour and Taylor (1992)

[TABLE]

The quantity $\dot{P}_{b}^{\mbox{\tiny Acc}}$ is due to the relative acceleration between the binary and the Solar System barycenter along the line of sight. The quantity $\dot{P}_{b}^{\mbox{\tiny Shk}}$ is the so-called Shklovskii effect and is due to centrifugal acceleration between the two binary components. The quantity $\dot{P}_{b}^{\dot{M}}$ accounts for mass loss of the system, and finally $\dot{P}_{b}^{T}$ accounts for any tidal effects. The term $\dot{P}_{b}^{\dot{G}}$ accounts for a possibly varying gravitational constant which generally happens on cosmological time scales, if for example the scalar field evolves in a cosmological potential and influences the gravitational constant.

The intrinsic orbital period derivative, $\dot{P}_{b}^{\mbox{\tiny int}}$, is of most interest for constraining STTs and can be found directly from the orbital energy that is lost due to gravitational radiation. In the context of STTs, $\dot{P}_{b}^{\mbox{\tiny int}}$ has the following contributions

[TABLE]

Each of these pieces can be calculated from a multipole expansion of the energy fluxes for the scalar field and the metric, which yield Damour and Esposito-Farese (1992)

[TABLE]

In GR only the quadrupolar term survives, but in STTs, both monopole and dipole radiation exist as well, the latter entering at -1 PN order relative to the GR term. In practice, the dipolar contribution from the scalar field will dominate the energy loss, followed by the quadrupolar terms. Thus, we will neglect the monopole contribution and the higher PN order contributions to the dipolar radiation in our calculations. The set of parameters $\{\dot{P}_{b},\gamma,\dot{\omega},r,s\}$ will be referred to as the post-Keplerian parameters from now on and will be the main set of parameters we use to constrain STTs.

II.4 Systems Considered

Typically, PSR-WD systems provide the tightest constraints on STTs due to the predicted dipolar contribution to $\dot{P}_{b}$, which has the largest effect when the binary components have very different scalar charges (as is the case with PSR-WD systems). However, one of the goals of this paper is to investigate how other post-Keplerian parameters affect the $(\alpha_{0},\,\beta_{0})$ parameter space. Therefore, we are also particularly interested in PSR-NS systems that have multiple post-Keplerian parameters measured to high precision. While the dipolar radiation is generally suppressed in such systems, $\dot{\omega}$ and $\gamma$ can become very large (and therefore inconsistent with observations) due to the presence of the scalar charges $\beta_{A,B}$ and $\kappa_{A}$, which can take values $\sim 10^{2}$ and even $\sim 10^{3}$ in some cases Anderson .

The set of PSRs in Shao et al. (2017) satisfy our first criteria for PSR-WD systems with mass measurements. These PSRs include J1738+0333 Antoniadis et al. (2012); Freire et al. (2012), J2222-0137 Kaplan et al. (2014); Cognard et al. (2017), J1012+5307 Lazaridis et al. (2009), and J1909-3744 Reardon et al. (2016)111We exclude J0348+0432 due to the large mass ($\sim 2M_{\odot}$) of the PSR in this system. The scalar charges we are using, a subset of what is in Anderson , only go up to 2.1$M_{\odot}$ ($2M_{\odot}$ for some EOS), and thus, the priors we use limit our exploration of the posterior in the neighbourhood of $2M_{\odot}$. This can be easily relaxed in the future if desired.. In one way or another, there have been measurements of $\dot{P}_{b}$, the companion mass $m_{B}$, and the mass ratio $q\equiv m_{A}/m_{B}$, thus allowing us to pin down the masses of the system. In the case of J1738+0333 and J1012+5307, there have been independent mass measurements from optical observations of the white dwarf companions, and thus, there are no correlations between any of the parameters. For J2222-0137 and J1909-3744, there exist measurements of the Shapiro parameters $r$ and $s$, which provide the extra information needed to pin down the masses. These systems also contain a diverse set of PSR masses that allow us to probe the $(\alpha_{0},\,\beta_{0})$ parameter space of STTs over a wide range of masses in which strong field effects like scalarization behave differently.

In terms of PSR-NS systems, we study J0737-3039A Kramer et al. (2006); Kramer and Wex (2009) and B1913+16 Weisberg et al. (2010), as these PSRs have been precisely timed and there exists measurements of multiple post-Keplerian parameters. Table 1 summarizes the Keplerian and post-Keplerian parameters associated with each of these PSRs. For B1913+16 there are measurements of all 5 post-Keplerian parameters introduced in Sec. II.3. J0737-3039A is part of a double PSR system, and there are measurements of all 5 post-Keplerian as well. But here there is also an independent measurement of the mass ratio from the orbital velocities of the two PSRs, providing a total of 6 parameters that can be used to place constraints.

III Constraints on STTs from binary PSRs: $1\sigma$ Constraints

In this section we will focus on placing constraints on the STT parameters $\alpha_{0}$ and $\beta_{0}$ by saturating the bounds on individual post-Keplerian parameters. As we have mentioned, some of the tightest constraints on STTs come from the lack of observed dipolar radiation in binary PSRs. In principle, if one has measured the orbital decay rate to be $\dot{P}_{b}{}_{,{\mbox{\tiny obs}}}$ up to some uncertainty $\delta\dot{P}_{b}$ and found that it is consistent with the GR prediction, then Eq. (36) can be used to place constraints on the scalar charges $\alpha_{A,B}$, and therefore, the theory parameters $(\alpha_{0},\,\beta_{0})$. The connection between $\alpha_{A,B}$ and the theory parameters, however, is dependent on the equation of state used to calculate the scalar charges. These constraints, therefore, must assume a particular EOS, and then be repeated for all possible choices. We investigate how these constraints behave when using a large number of EOSs, and how they can be improved when using measurements of other post-Keplerian parameters, considering both DEF and MO theories.

III.1 Constraints from $\dot{P}_{b}^{\mbox{\tiny int}}$

To place constraints using measurements of $\dot{P}_{b}$, we evaluate Eq. (36) for the entire region of parameter space investigated in Anderson , i.e. $-5.5\lesssim\log_{10}(\alpha_{0})\lesssim 0$ and $-5\lesssim\beta_{0}\lesssim 5$. While the masses are measured to a finite precision, we evaluate the dipole term at the best fit value of the masses for the moment, and returning to this point later. Then, by determining if these predicted values of $\dot{P}_{b}^{{\mbox{\tiny int}}}$ lie with the range $\dot{P}_{b}{}_{,{\mbox{\tiny obs}}}\pm\delta\dot{P}_{b}$ we can determine if that point in parameter space is consistent or inconsistent with observations.

Figure 1 shows such constraints for multiple PSR systems that have accurate measurements of $\dot{P}_{b}$ for both DEF theory (left panel) and MO theory (right panel) with the AP3 EOS Akmal et al. (1998). Let us first focus on the constraints on DEF theory. These constraints confirm the results first presented in Freire et al. (2012). The only difference between those results and the ones found here arises because the scalar charges have been here calculated much more finely in $(\alpha_{0},\beta_{0})$ space, allowing us to resolve the structure of the “horn” constraints more accurately. These horns arise because there are certain values of $\alpha_{0}$ and $\beta_{0}$ for which the dipole term is significantly suppressed, thus preventing any constraint.

The constraints presented on the right panel of Fig. 1 on MO theory are new. Observe that the strength of the constraints in the two theories is roughly the same. This is mostly because these theories are nearly identical to each other in the limit $\varphi\ll 1$, and therefore, they only differ substantially from each other when $\beta_{0}$ is very negative and/or $\alpha_{0}\gtrsim 0.01$. Indeed, we see that the main difference between the constraints on the different theories is the “horns” that appear in Fig. 1 for large values of $\alpha_{0}$ and $\beta_{0}<0$.

The EOS affects the magnitude of the scalar charges, and thus, the mass at which spontaneous scalarization occurs, and the constraints that can be obtained from observations of $\dot{P}_{b}$. Let us then repeat the analysis presented above, but this time using 11 different EOSs, namely the ones previously studied in Refs. Shao et al. (2017); Anderson : AP3-4 Akmal et al. (1998), ENG Engvik et al. (1996), H4 Lackey et al. (2006), MPA1 Muther et al. (1987), MS0 Mueller and Serot (1996), MS2 Mueller and Serot (1996), PAL1 Prakash et al. (1988), SLy4 Douchin and Haensel (2001), and WFF1-2 Wiringa et al. (1988). Figure 2 illustrates the effect that the EOS has on the constraints placed by J1738+0333. As one can see, the EOS simply shifts the location of the horns left and right. This is because different EOS predict different masses at which scalarization kicks in, modifying the values of $(\alpha_{0},\beta_{0})$ at which the dipole term is suppressed. Aside from this modification, the variation of the EOS does not affect the overall strength of the constraints much, especially in regions where spontaneous scalarization does not occur, confirming previous results from Freire et al. (2012); Shao et al. (2017). Although one ought to marginalize over our ignorance of the EOS, as we will discuss in the next section, the above analysis suggests that the constraints placed on DEF theory in the past Freire et al. (2012) are robust to this ignorance.

III.2 Constraints from other Post-Keplerian parameters

For most systems, measurements of $\dot{P}_{b}$ place the tightest constraints on STTs since the dipolar contribution to Eq. (36) enters at -1 PN order relative to the quadrupole terms, and is therefore dominant. However, for double NS systems, the other post-Keplerian parameters can be significantly affected by the scalar charges, even though the dipolar contribution to Eq. (36) can be negligible. Thus, using the same methods used above for $\dot{P}_{b}^{\mbox{\tiny int}}$, we here investigate the constraints that can be placed from other post-Keplerian parameters in these double NS systems. More precisely, for any one post-Keplerian parameter $d^{i}_{{\mbox{\tiny obs}}}$ measured with $1\sigma$ accuracy $\delta d^{i}$, we calculate the values of $(\alpha_{0},\beta_{0})$ for which the predicted value of $d^{i}_{{\mbox{\tiny th}}}$ lies inside $d^{i}_{{\mbox{\tiny obs}}}\pm\delta d^{i}$.

Figure 3 shows the relative strength of the constraints on $(\alpha_{0},\beta_{0})$ from B1913+16 and J0737-3039 using the post-Keplerian parameters $\dot{P}_{b}$, $\dot{\omega},$ and $\gamma$. In both cases, the observations of $\dot{\omega}$ place significantly tighter constraints on STTs than measurements of $\dot{P}_{b}$, especially for $\beta_{0}<0$. This may come as a surprise but it can be understood by the functional form of $\dot{\omega}$. The equation for $\dot{\omega}$ contains the higher-order scalar charge $\beta_{A}$, which can take on significantly large values when $\beta_{0}<0$. Typically, however, these effects are suppressed by a factor of $\alpha_{0}^{2}$ for PSR-WD systems, which is enough to make these effects completely negligible Anderson . However, when the companion is another NS, $\alpha_{B}^{2}$ can be of order unity, making these contributions large and lead to confrontation with observations as we are seeing in Fig. 3.

These constrains, particularly the ones places by $\dot{\omega}$, are very sensitive to the masses used in Eq. (32), especially when the error is small. Consider B1913+16 in GR, for instance, in which $\dot{\omega}$ has been measured to very high precision. The current measurements of $\dot{\omega}$ are so precise that the total mass of the binary in GR can be determined to one part in $10^{6}$. However, the individual masses cannot be constrained this well from the other parameters, which means that one is not well-justified in using the best-fit masses into the equations to place the above constraints. Moreover, these plots should not be taken as a hard upper limit on $(\alpha_{0},\beta_{0})$, but rather more as a guideline of what constraints could be placed.

There are two ways to address these types of constraints in a more consistent manner and that are not prone to issues described above. One such way is to perform an analysis similar to the one in Ref. Esposito-Farese (2011) in which one calculates the theoretical values of the post-Keplerian parameters and their errors for every combination of $\alpha_{0}$ and $\beta_{0}$. Then, that point in parameter space is only excluded if there does not exist a pair of masses $(m_{A},\,m_{B})$ that lies in the intersection of these curves. Another method, however, involves using an MCMC to investigate these constraints and it not limited to only exploring the $1\sigma$ constraints of this section. Thus, in the next section, when we use Bayesian methods to investigate these constraints in a more informative and self-consistent manner.

IV Bayesian inference

Thus far we have considered the tightest possible $1\sigma$ constraints that can be placed on STTs using individual PSRs. However, we have not investigated how these constraints are affected by possible covariances between the post-Keplerian parameters, or by marginalizing over our ignorance on the EOS. Thus, in this section we will employ Bayesian methods through MCMC simulations to address these short comings.

IV.1 Basics of Bayesian Inference

In Bayesian statistics, given a data set $d_{n}$ and some hypothesis $H$ (playing the role of the theory in this case), the posterior distribution on the parameters of the hypothesis are determined by Bayes’ theorem

[TABLE]

where $P(\vec{\lambda}|H)$ is the prior probability density of the parameters $P(d_{n}|\vec{\lambda},H)$ is the likelihood of the data given the model $H$ and parameters $\vec{\lambda}$, and $P(d_{n}|H,I)$ is the model evidence that plays the role of a normalization factor here. Therefore, if one has the priors and the likelihood, then in principle the posterior distribution of $\vec{\lambda}$ can be calculated. In order to explore these posterior distributions efficiently, we use MCMC simulations for each data set $d_{n}$, corresponding to the independent observations of each PSR.

For our purposes, $\vec{\lambda}=\{\xi_{i},m_{A},\,m_{B},\,\log_{10}(\alpha_{0}),\,\beta_{0}\}$, and we enforce uniform priors for all parameters. The discrete parameter $\xi_{i}$ represents a given EOS and can take integer values from 0 to 10, with 0 corresponding to AP3 and 10 to WFF2 (ordered alphabetically). For the masses, the priors ranges are $0.01M_{\odot}<m<1.44M_{\odot}$ if the mass corresponds to a WD and $1M_{\odot}<m<2M_{\odot}$ if the mass corresponds to a NS222We have restricted our investigation to this range of NS masses because it becomes numerically complicated to deal with the scalar charges above this range. For our study, however, this does not limit our investigation because we do not include PSR systems with large NS masses and indeed these priors do not affect our final results.. The STTs parameters have the prior ranges $-5.5<\log_{10}\alpha_{0}<-1$ and $-5<\beta_{0}<5$. The priors we have chosen are primarily constrained by the information we have available about the scalar charges today. Despite being able to, we do not allow $\log_{10}\alpha_{0}$ to reach 0, nor do we allow the NS mass to exceed $2M_{\odot}$. We restrict $\log_{10}\alpha_{0}$ simply because we find that our MCMC chains never explore these regions (as one might expect given the results of Sec. III.1. We restrict the masses of the NSs to be less than $2M_{\odot}$ because not all of the EOS we consider can achieve masses much larger than this.

For our logarithmic likelihood we use a multivariate Gaussian distribution given by

[TABLE]

where $d_{n}^{i}$ are the post-Keplerian parameters of the $n$th data set (or $n$th PSR), $d_{{\mbox{\tiny th}}}^{i}(\vec{\lambda})$ are the theoretical prediction for the same post-Keplerian parameters given the parameters $\vec{\lambda}$, and $\mathbb{C}^{-1}_{ij}$ is the inverse of the correlation matrix associated with the data set $d_{n}$. This likelihood has a maximum when the theoretical post-Keplerian parameters are precisely equal to their observed values and adequately handles the covariances that exist between the observed post-Keplerian parameters. The likelihood contains the theoretical predictions of the post-Keplerian parameters $d^{i}_{{\mbox{\tiny th}}}$, given in Eqs. (31)-(34) and Eqs. (36)-(40), which in turn all depend explicitly on the scalar charges $(\alpha_{A,B},\,\beta_{A,B},\,k_{A})$. For any WD companions, the charges reduce to $\alpha_{B}=\alpha_{\infty}$ and $\beta_{A}=\alpha_{A}/\varphi_{\infty}\,=\,\beta_{\infty}$ since WD are weakly self gravitating. For NSs on the other hand, the scalar charges are functions of the NS mass, $\alpha_{0}$, and $\beta_{0}$ and must be solved numerically. To avoid this last step, we linearly interpolate the data from Anderson for the scalar charges, allowing the quick computation of the likelihood over the entire prior range of the parameters.

We start our MCMC chains near the binary masses predicted by GR, with EOS AP3, $\log_{10}\alpha_{0}=-2$, $\beta_{0}=-2$. While we start the chains near the expected peaks of the posterior distributions, the chain are allowed to explore the entire prior range, and we find that, after a burn-in phase, the chains always find these peaks regardless of where we start them. The parameter space is then explored through a proposal distribution that is equivalent to a random draw from a Gaussian distribution, with a variance that is different for each of the parameters. The variances are chosen to ensure decent acceptance ratios and autocorrelation lengths, and in practice they are different for each PSR..

Proposed jumps are accepted or rejected based on the Metropolis-Hastings algorithm. The Metropolis ratio for symmetric proposal distributions is just the ratio of the posterior distributions, namely

[TABLE]

where we used Bayes theorem from Eq. (41) to rewrite this ratio in terms of the priors and likelihoods. The acceptance probability is then defined via $a(\vec{\lambda}_{old},\vec{\lambda}_{new})=\text{min}(1,r)$, which determines whether or not the proposed set of parameters is accepted by the chain or not. This is accomplished by drawing a random number $u$ between 0 and 1 and comparing it to the acceptance probability, i.e. if $a\geq u$ the jump is made, and otherwise it is not. If the proposed jump leads to a larger likelihood, the jump is always accepted, but even if the new likelihood is smaller, there is still a chance that the proposal will be accepted. We repeat this process millions of time for each PSR to ensure convergence and good exploration of the posterior distributions.

The result of an MCMC is the joint posterior distribution on all parameters $\vec{\lambda}$. To obtain information about a single parameter, say $\beta_{0}$, we marginalize, or integrate, over all other parameters $\vec{\theta}$, for instance

[TABLE]

In practice, this is equivalent to creating a histogram of the likelihood evaluations over the parameter that is being marginalized.

IV.2 Bayesian Results

Figure 4 shows the marginalized posteriors on the parameters $(m_{A},\,m_{B},\,\log_{10}\alpha_{0},\,\beta_{0})$ for J1738+0333. The marginalized posteriors on $\xi_{i}$ for the EOS is uniform, meaning that the MCMC showed no preference toward any particular EOS and thus we do not show it because it is uninformative. We recover mass distributions that are consistent with observations and the predictions of GR. The posteriors on $\alpha_{0}$ and $\beta_{0}$ behave as one would expect. Lower values of $\alpha_{0}$ (towards GR) are preferred and the posteriors on $\beta_{0}$ have a nearly identical form as the constraints placed by $\dot{P}_{b}$ in the previous section. The latter makes sense since we are not enforcing Solar System priors on the parameters, and thus, the MCMC spends more time exploring regions that are less constrained, i.e. near $\beta_{0}=-2$.

The joint posterior on $\alpha_{0}$ and $\beta_{0}$ are presented in Fig. 5, where we also include an overlay of the constraints placed from the previous section. This figure is really a heat plot of the posterior distribution, with brighter (darker) colors representing low (high) posterior probabilities. Therefore, the dark regions are unconstrained by our Bayesian analysis, while the bright regions are excluded. We see great consistency between the Bayesian constraints and the constraints from the previous subsection.

Thus far we have only discussed results involving single PSRs with associated data $d_{n}$, so let us now combine the constraints from multiple observations in two different ways. The first involves constructing a new log likelihood that contains information from all of the PSRs we consider. This is equivalent to simply adding the log likelihoods (multiplying the actual likelihoods) associated with each individual PSR and running a new MCMC with a new parameter vector given by $\vec{\lambda}_{\text{total}}=\{\sum_{n}m_{A,n},\sum_{n}\,m_{B,n},\,\xi_{i}\,,\alpha_{0},\,\beta_{0}\}$ where $n$ stands for the $n$th binary system. One would then marginalize the likelihood over all parameters that are not $(\alpha_{0},\beta_{0})$ to find the new joint posterior. This is quite possible but it requires us to effectively recompute the posteriors we have already found.

Another way to combine information is to carefully “stack” the joint posteriors on the parameters $\vec{\sigma}\,=\,(\xi_{i},\,\alpha_{0},\,\beta_{0})$ that we have already calculated333It is impossible to stack posteriors on single parameters simply because the operations of integration and multiplication do not commute. Thus, we are only able to stack marginalized posteriors if we have marginalized over parameters that are not shared in different observations, like the masses in this case.. In this case the total posterior distribution is just the product of the individual posteriors. To do this we start with Bayes theorem in Eq. (41) with total data $D$ from $N$ systems, and we note that the priors are identical with the likelihood following the relation

[TABLE]

Putting this information back into Bayes theorem for the total posterior we find the relation

[TABLE]

Therefore, we are able to take the individual joint posteriors on $\vec{\sigma}$ that we have already calculated and multiply them in this manner to find the total posterior we would have found had we done a new MCMC simulation with all PSRs included.

Figure 6 shows the total joint posterior on $(\alpha_{0},\,\beta_{0})$, marginalized over the EOS, after combining the posteriors from the PSRs in Table 1 via Eq. (46). The overall effect of this “stacking” of the posteriors does not produce a significant effect but it does visibly reduce the probability of the larger values of $\alpha_{0}$ near $\beta_{0}\sim-2$. This results makes sense when compared to results found in Sec. III.1. Each PSR tends to place different constraints in this region of parameter space, and the Bayesian methods we have employed here allow us to take these multiple constraints into consideration simultaneously.

Rather than assuming flat priors for the STTs parameters we can enforce constraints that have been placed from other tests, like Solar System tests, through different priors. The Solar System constraint on the post-Keplerian parameters appear in Eqs. (18)-(21), and numerically, these constraints turn into $\alpha_{\infty}\leq 0.003391$. As we have pointed out in Sec. II, this constraint on $\alpha_{\infty}$ technically translates to different constraints on $\alpha_{0}$ in each theory. However, for practical purposes the constraints on $\alpha_{0}$ are identical in both DEF and MO theory because of the small magnitude of $\alpha_{\infty}$.

The effect of more stringent priors is to tighten the posteriors on all parameters. These priors effectively remove the posterior distributions on $\alpha_{0}$ larger than $0.003391$ that we have presented thus far and therefore reduce the parameter space that the MCMC explores. This results in tighter posteriors for the masses of the systems, as well as more stringent constraints on $(\alpha_{0},\beta_{0})$ as shown in Fig. 7.

V Conclusion

We have investigated the constraints that binary PSRs can place on the parameters $\alpha_{0}$ and $\beta_{0}$ that appear in scalar-tensor theories, particular in DEF theory and in MO theory. We show for the first time in MO theory the type of constraints that one can place on STTs using observed orbital period decay $\dot{P}_{b}$. Then, using other post-Keplerian parameter we demonstrate that for certain systems, like double NS systems, even tighter constraints can be place on the $(\alpha_{0},\,\beta_{0})$ parameter space of STTs. We find consistency with these results when we use Bayesian methods to investigate the posterior distributions of the parameters and find that the dominant effects the dictate the overall constraints is indeed $\dot{P}_{b}$.

This work has been one of the first attempts to apply an MCMC to explore the constraints on STT parameter space from binary PSRs, and therefore, it should be considered a first step towards a more in depth Bayesian study. For example, we have only used the quoted values of the post-Keplerian parameters, along with the correlations, to construct our likelihood functions and explore the posterior distributions of $\alpha_{0}$ and $\beta_{0}$. Such a method introduces systematic errors because we are not able to incorporate covariances between the normal Keplerian parameters, even though these parameters are measured to extremely high precision. Our analysis also potentially suffers from having a limited number of data points, i.e. the post-Keplerian parameters here.

A more complete study, and one we hope to perform in the future, would involve using more sophisticated MCMC techniques like parallel tempering, to allow us to explore potentially multi-modal posteriors, or like those implemented in the emcee_hammer package available for Python. Moreover, one could also re-analyze the PSR data to incorporate the TOAs directly into the likelihood and allow an MCMC to perform the fits for the Keplarian and the post-Keplerian parameters simultaneously, essentially performing the computations that timing packages like TEMPO and TEMPO2 handle. This would allow one to easily fit the TOAs given any gravitational theory, like GR and STTs, that have a well developed timing formula. The benefit of handling the TOA directly is that it removes any systematics and inconsistencies that arise in our study, i.e. trying to analyze data that has already been processed.

Another possible avenue for future work is to consider more complicated binary pulsar systems in a Bayesian approach to test GR. A particular interesting system is the triple system recently observed . This observation allows for a measurement of a parameter that constraints a certain type of SEP violation. This parameter depends on the scalar charges in STTs, as well as on the component masses of inner binary. One could thus repeat the analysis presented here for this system, marginalizing over the component masses and over the EOSs to obtain new stringent constraints on DEF and MO theory.

Acknowledgements.

We would like to thank Norbert Wex for useful discussions on constraining STTs and Neil Cornish and Travis Robson for useful discussions on Bayesian statistics. DA would like to thank Paulo Freire and the Max-Planck-Institut für Radioastronomie for their hospitality during part of this work. NY and DA acknowledge support from NSF grant PHY-1759615 and NASA grants NNX16AB98G and 80NSSC17M0041.