New binary pulsar constraints on Einstein-æther theory after GW170817

Lorentz symmetry has been the foundation of the magnificent edifice of theoretical physics for more than a century, playing a central role in special and general relativity (GR), as well as in the quantum theory of fields. Because of its special status, Lorentz invariance has been tested to exquisite precision in the matter sector via particle physics experiments [45, 49, 50, 54]. More recently, this experimental program has been extended to the matter-gravity [48], dark matter [13, 16], and pure-gravity sectors [42, 53], where bounds on Lorentz violations (LVs) have been historically looser (because of the intrinsic weakness of the gravitational interaction).

Compelling theoretical reasons to seriously consider the possibility of LVs in the purely gravitational sector were provided by the realization that they could generate a better behavior in the ultraviolet (UV) limit. In particular, Hořava [39] showed that by allowing for a non-isotropic scaling between space and time, one can construct a theory that is power-counting renormalizable in the UV. Renormalizability beyond power counting (i.e. perturbative renormalizability) in special ('projectable') versions of Hořava gravity has also been proven [10].

The low-energy limit of Hořava gravity reduces to 'khronometric theory' [15, 43], which consists of GR plus an additional hypersurface-orthogonal and timelike vector field, often referred to as the 'æther'. Because this vector field is hypersurface orthogonal, it selects a preferred spacetime foliation, which makes LVs manifest. A more general boost-violating low-energy gravitational theory, however, can be obtained by relaxing the assumption that the æther be hypersurface-orthogonal, in which case it selects a preferred time threading of the spacetime rather than a preferred foliation. The resulting theory is known as Einstein-æther theory [46].

Despite allowing for an improved UV behavior, LVs in gravity face long-standing experimental challenges, particularly when it comes to their percolation into the matter sector, where particle physics experiments are in excellent agreement with Lorentz symmetry. While some degree of percolation is inevitable, because of the coupling between matter and gravity, mechanisms suppressing it have been put forward, including suppression by a large energy scale [58], or the effective emergence of Lorentz symmetry at low energies as a result of renormalization group flows [11, 12, 23] or accidental symmetries [38].

At the same time, purely gravitational bounds on LVs are becoming increasingly compelling. The parameters ('coupling constants') of both Einstein-æther and khronometric theory have been historically constrained by theoretical considerations (absence of ghosts and gradient instabilities [18, 37, 41], well-posedness of the Cauchy problem [60]), by the absence of vacuum Cherenkov cascades in cosmic-ray experiments [30]), by solar-system tests [18, 19, 34, 55, 69], by observations of the primordial abundances of elements from Big-Bang nucleosynthesis [22], by other cosmological tests [8], and by precision timing of binary pulsars (where LVs generically predict violations of the strong equivalence principle) [9, 32, 33, 73, 74]. More recently, the coincident detection [2, 3] of gravitational waves (GW170817) and gamma rays (GRB170817A) emitted by the coalescence of two neutron stars and the subsequent kilonova explosion has allowed extremely strong constraints on the propagation speed of gravitational waves, which must equal that of light to within ⁸ 10⁻¹⁵ [1], which in turn places even more stringent bounds on the couplings of both theories [31, 56, 59, 60].

The bounds from the coincident GW170817/GRB170817A observations force us to rethink the parameter spaces of both Einstein-æther and khronometric theory, as the only currently allowed regions appear to be ones that were previously thought to be of little interest, and which were not explored extensively. In the case of khronometric theory, references [9, 59] found that the couplings that remain viable after GW170817 and GRB170817 produce exactly no deviations away from the predictions of GR, not only in the Solar System, but also in binary systems of compact objects, be they black holes (BHs) or neutron stars (NSs), to leading post-Newtonian (PN) order. Reference [35] extended this result to the quasinormal modes of spherically symmetric BHs and to fully non-linear (spherical) gravitational collapse, where again no deviations from the GR predictions are found. It would therefore seem that the most promising avenue to further test khronometric theory may be provided by cosmological observables (e.g. Big-Bang nucleosynthesis abundances or CMB physics), where the viable couplings do produce non-vanishing deviations away from the GR phenomenology.

Like for khronometric theory, the parameter space where detailed predictions for isolated/binary pulsars were obtained in Einstein-æther theory [73, 74] does not include the region singled out by the combination of the GW170817/GRB170817A bound and existing solar-system constraints (see reference [60] for a discussion). The goal of this paper is therefore to extend the previous analysis of binary/isolated-pulsar data by some of us [73, 74] to this region of parameter space. This will require a significant modification of the formalism that references [73, 74] utilized to calculate pulsar 'sensitivities', i.e. the parameters that quantify violations of the strong-equivalence principle in these systems. Moreover, we will extend our analysis to include additional data over that considered in references [73, 74], namely the triple system PSR J0337+1715 [7]. Overall, we find that observations of the damping of the period of quasi-circular binary pulsars, and that of the triple system PSR J0337+1715, reduce the viable parameter space of Einstein-æther theory by about an order of magnitude over previous constraints.

We will also amend an error (originally pointed out in reference [70]) in the calculation of the strong-field preferred-frame parameters ${\hat{\alpha }}_{1}$ and ${\hat{\alpha }}_{2}$ for isolated pulsars, which were presented in references [73, 74]. While we have checked that this error does not impact the bounds presented in references [73, 74], we present in appendix A a detailed derivation of ${\hat{\alpha }}_{1}$ and ${\hat{\alpha }}_{2}$ for possible future applications, also correcting a few typos present in the original calculation of reference [70].

This paper is organized as follows. In section 2 we give a succinct introduction to Einstein-æther theory, including the modified field equations and the current observational bounds on the coupling constants. In section 3 we introduce the concept of stellar sensitivities as parameters regulating violations of the strong equivalence principle. Solutions describing slowly moving stars are derived in section 4, and they are used in section 5 to compute the sensitivities. Section 6 uses the sensitivities to obtain the constraints on Einstein-æther theory resulting from observations of binary and triple pulsar systems. We summarize our conclusions in section 7. Appendix A contains a calculation of the strong-field preferred-frame parameters ${\hat{\alpha }}_{1}$ and ${\hat{\alpha }}_{2}$ in Einstein-æther theory, fixing an oversight in [73], which was pointed out by [70], and correcting also a few typos present in [70] itself. We will adopt units where c = 1 and a signature + −−−, in accordance with most of the literature on Einstein-æther theory.

In order to break boost (and thus Lorentz) symmetry, Einstein-æther theory introduces a dynamical threading of the spacetime by a unit-norm, time-like vector field U . This vector field, often referred to as the æther, physically represents a preferred 'time direction' at each spacetime event. Requiring the action to also include the usual spin-2 graviton of GR, to be quadratic in the æther derivatives, and to feature no direct coupling between the matter and the æther (so as to enforce the weak equivalence principle, i.e. the universality of free fall, and the absence of matter LVs at tree level), one obtains the action [44, 46]

$$\begin{align}S& =-\frac{1}{16\pi G}\int \left[R+\frac{1}{3}{c}_{\theta }{\theta }^{2}+{c}_{\sigma }{\sigma }_{\mu \nu }{\sigma }^{\mu \nu }+{c}_{\omega }{\omega }_{\mu \nu }{\omega }^{\mu \nu }+{c}_{a}{A}_{\mu }{A}^{\mu }\right.\\ & \left.\quad +\lambda ({U}^{\mu }{U}_{\mu }-1)\right]\sqrt{-g}\enspace {\mathrm{d}}^{4}x+{S}_{\text{mat}}(\psi ,{g}_{\mu \nu }),\end{align} \tag{ 1 }$$

where R is the four-dimensional Ricci scalar, g the determinant of the metric, G the bare gravitational constant (related to the value G_N measured locally by G_N = G/(1 − c_a /2) [22, 42]), ψ collectively denotes the matter degrees of freedom, λ is a Lagrange multiplier enforcing the æther's unit norm, c_θ , c_σ , c_ω and c_a are dimensionless constants ⁹ , and we have decomposed the æther congruence into the expansion θ, the shear σ_μν , the vorticity ω_μν and the acceleration A_μ as follows:

$$\begin{equation}{A}^{\mu }={U}^{\nu }{\nabla }_{\nu }{U}^{\mu },\end{equation} \tag{ 2 }$$

$$\begin{equation}\theta ={\nabla }_{\mu }{U}^{\mu },\end{equation} \tag{ 3 }$$

$$\begin{equation}{\sigma }_{\mu \nu }={\nabla }_{\left(\right.\nu }{U}_{\mu \left.\right)}+{A}_{\left(\right.\mu }{U}_{\nu \left.\right)}-\frac{1}{3}\theta {h}_{\mu \nu },\end{equation} \tag{ 4 }$$

$$\begin{equation}{\omega }_{\mu \nu }={\nabla }_{\left[\right.\nu }{U}_{\mu \left.\right]}+{A}_{\left[\right.\mu }{U}_{\nu \left.\right]},\end{equation} \tag{ 5 }$$

with h_μν = g_μν − U_μ U_ν the projector onto the hyperspace orthogonal to U .

By varying the action with respect to the metric, the æther and the Lagrange multiplier, and by eliminating the latter from the equations, one obtains the generalized Einstein equations

$$\begin{equation}{E}_{\alpha \beta }\equiv {G}_{\alpha \beta }-{T}_{\alpha \beta }^{\unicode{198}}-8\pi G{T}_{\alpha \beta }^{\text{mat}}=0\end{equation} \tag{ 6 }$$

and the æther equations

$$\begin{equation}{\unicode{198}}_{\mu }=\left[{\nabla }_{\alpha }{J}^{\alpha \nu }-\left({c}_{a}-\frac{{c}_{\sigma }+{c}_{\omega }}{2}\right){A}_{\alpha }{\nabla }^{\nu }{U}^{\alpha }\right]{h}_{\mu \nu }=0,\end{equation} \tag{ 7 }$$

where G_αβ is the Einstein tensor, the æther stress–energy tensor is

$$\begin{align}{T}_{\alpha \beta }^{\unicode{198}}& ={\nabla }_{\mu }\left({J}_{\left(\right.\alpha }{}^{\mu }U_{\beta \left.\right)}-{J}^{\mu }{}_{\left(\right.\alpha }U_{\beta \left.\right)}-{J}_{(\alpha \beta )}{U}^{\mu }\right)+\frac{{c}_{\omega }+{c}_{\sigma }}{2}\left[({\nabla }_{\mu }{U}_{\alpha })({\nabla }^{\mu }{U}_{\beta })\right.\\ & \left.\quad -({\nabla }_{\alpha }{U}_{\mu })({\nabla }_{\beta }{U}^{\mu })\right]+{U}_{\nu }({\nabla }_{\mu }{J}^{\mu \nu }){U}_{\alpha }{U}_{\beta }-\left({c}_{a}-\frac{{c}_{\sigma }+{c}_{\omega }}{2}\right)\\ & \quad {\times}\left[{A}^{2}{U}_{\alpha }{U}_{\beta }-{A}_{\alpha }{A}_{\beta }\right]+\frac{1}{2}{M}^{\sigma \rho }{}_{\mu \nu }\nabla _{\sigma }{U}^{\mu }{\nabla }_{\rho }{U}^{\nu }{g}_{\alpha \beta },\end{align} \tag{ 8 }$$

with

$$\begin{align*}J_{\enspace \mu }^{\alpha }& \equiv {M}^{\alpha \beta }{}_{\mu \nu }\nabla _{\beta }{U}^{\nu },\\ M_{\quad \mu \nu }^{\alpha \beta }& =\left(\frac{{c}_{\sigma }+{c}_{\omega }}{2}\right){h}^{\alpha \beta }{g}_{\mu \nu }+\left(\frac{{c}_{\theta }-{c}_{\sigma }}{3}\right){\delta }_{\mu }^{\alpha }{\delta }_{\nu }^{\beta }+\left(\frac{{c}_{\sigma }-{c}_{\omega }}{2}\right){\delta }_{\nu }^{\alpha }{\delta }_{\mu }^{\beta }+{c}_{a}{U}^{\alpha }{U}^{\beta }{g}_{\mu \nu },\end{align*}$$

and the matter stress–energy tensor is defined as usual by

$$\begin{equation}{T}_{\text{mat}}^{\alpha \beta }\equiv -\frac{2}{\sqrt{-g}}\frac{\delta {S}_{\text{mat}}}{\delta {g}_{\alpha \beta }}.\end{equation} \tag{ 9 }$$

As already mentioned, a number of experimental and theoretical results constrain Einstein-æther theory and the couplings c_i . In more detail, perturbing the field equations about Minkowski space yields propagation equations for spin-0 (i.e. scalar), spin-1 (i.e. vector) and spin-2 (i.e. tensor gravitons) modes. Their propagation speeds are respectively given by [41]

$$\begin{equation}{c}_{T}^{2}=\frac{1}{1-{c}_{\sigma }},\end{equation} \tag{ 10 }$$

$$\begin{equation}{c}_{V}^{2}=\frac{{c}_{\sigma }+{c}_{\omega }-{c}_{\sigma }{c}_{\omega }}{2{c}_{a}(1-{c}_{\sigma })},\end{equation} \tag{ 11 }$$

$$\begin{equation}{c}_{S}^{2}=\frac{({c}_{\theta }+2{c}_{\sigma })(1-{c}_{a}/2)}{3{c}_{a}(1-{c}_{\sigma })(1+{c}_{\theta }/2)}.\end{equation} \tag{ 12 }$$

In order to ensure stability at the classical level (i.e. no gradient instabilities) and at the quantum level (i.e. no ghosts) one needs to have ${c}_{T}^{2}{ >}0$, ${c}_{V}^{2}{ >}0$ and ${c}_{S}^{2}{ >}0$ [37, 41]. If we also require the modes to carry positive energy, we get c_a > 0 and c_ω > 0 [29]. Furthermore, significantly subluminal graviton propagation would cause ultrarelativistic matter to lose energy to gravitons via a Cherenkov-like process [30]. Since this effect is not observed e.g. in ultrahigh energy cosmic rays, one must have ${c}_{I}^{2}\;\gtrsim \;1-\mathcal{O}(1{0}^{-15})$ (with I = T, V, S). More recently, the coincident detection of a neutron-star merger in GW170817 (gravitational waves) and GRB170817A (gamma rays) had led to the bound −3 × 10⁻¹⁵ < c_T − 1 < 7 × 10⁻¹⁶ [1].

Expanding the field equations through 1PN order leads to the conclusion that the 1PN dynamics is well described (like in GR) by the parametrized PN (PPN) expansion [55, 69]. However, unlike in GR, the preferred frame parameters α₁ and α₂ appearing in the PPN expansion do not vanish, but are given by [34]

$$\begin{equation}{\alpha }_{1}=4\frac{{c}_{\omega }({c}_{a}-2{c}_{\sigma })+{c}_{a}{c}_{\sigma }}{{c}_{\omega }({c}_{\sigma }-1)-{c}_{\sigma }},\end{equation} \tag{ 13 }$$

$$\begin{equation}{\alpha }_{2}=\frac{{\alpha }_{1}}{2}+\frac{3({c}_{a}-2{c}_{\sigma })({c}_{\theta }+{c}_{a})}{(2-{c}_{a})({c}_{\theta }+2{c}_{\sigma })}.\end{equation} \tag{ 14 }$$

Solar system experiments require |α₁| ≲ 10⁻⁴ and |α₂| ≲ 10⁻⁷ [55, 69]. By saturating these bounds (i.e. requiring in particular that |α₁| ≲ 10⁻⁴ but not |α₁| ≪ 10⁻⁴) and combining them with the constraints on the propagation speeds, one finds ${c}_{\sigma }\approx \mathcal{O}(1{0}^{-15})$, ${c}_{a}\approx \mathcal{O}(1{0}^{-4})$, and ${c}_{\theta }\approx 3{c}_{a}[1+\mathcal{O}(1{0}^{-3})]$. The resulting experimentally viable parameter space, therefore, is effectively (i.e. to within a fractional width of 10⁻⁴ or better in the parameters) one-dimensional: c_σ , c_a , c_θ ≈ 0, but c_ω is essentially unconstrained [60].

Another viable region of the parameter space can be obtained by not saturating the PPN constraints [60]. In more detail, one may require |α₁| be much smaller than its upper limit, so as to automatically satisfy the bound on α₂ (since α₂ ∝ α₁ if c_σ ≈ 0, as imposed by GW170817 and GRB170817A). This leads to |c_a | ≲ 10⁻⁷ and thus to an effectively two-dimensional experimentally viable parameter space (c_θ , c_ω ), with the only additional requirement that |c_θ | ≲ 0.3 to ensure that the production of light elements during Big Bang nucleosynthesis gives predictions in agreement with observations [22].

Both of the viable regions of parameter space identified above were not considered in references [73, 74], where neutron-star sensitivities in Einstein-æther theory were first computed. This is because back when references [73, 74] were written, the strongest constraints available were the Solar System ones (since the GW170817/GRB170817A constraint was not yet available). References [73, 74] solved for c_a and c_θ in terms of c_σ , c_ω , α₁ and α₂, and fixed the latter two to their largest allowed values (respectively α₁ = 10⁻⁴ and α₂ = 10⁻⁷). Therefore, (i) this does not include the first region listed above ((c_θ , c_σ , c_a ) ≲ 10⁻⁴ with c_ω kept free), which is obtained by choosing α₁ ∝ α₂; and (ii) references [73, 74] did not sample accurately enough the sub-region c_σ ≈ 0, α₁ ≲ 10⁻⁴, α₂ ≲ 10⁻⁷, since there was no reason to do so at that time and because the techniques employed broke down in that sub-region (as we will show in the following).

Most theories extending/modifying GR involve additional degrees of freedom besides the massless tensor graviton of GR. These additional gravitational polarizations cannot directly couple with matter significantly, to avoid introducing unwanted fifth forces in particle physics experiments, and to prevent violations of the weak equivalence principle (and particularly violations of the universality of free fall for weakly gravitating objects). Nevertheless, effective couplings between the extra gravitons and matter may be mediated by the metric perturbations (i.e. by the tensor gravitons present also in GR), which are typically coupled non-minimally to the extra gravitational degrees of freedom. These effective couplings become important when the metric perturbations are 'large', which is the case for strongly gravitating systems such as those involving NSs and/or BHs.

A useful way to parametrize this effective coupling is provided by the sensitivity parameters. Because of the aforementioned effective couplings, the mass of strongly gravitating objects will be comprised not only of the contributions from matter and the metric (like in GR), but it will also generally depend on the additional gravitational fields. We can thus describe isolated objects, and members of a widely separated binary, by a point particle model (like in GR), but with a non-constant mass depending on the extra fields. Because the mass is a scalar quantity, it must depend on a scalar constructed from the æther field U , the simplest of which is the Lorentz factor γ ≡ u ⋅ U , where u is the particle's (i.e. the body's) four-velocity.

In many practical situations (including the long inspiral of a binary system of compact objects) one may assume that the relative speed between the æther and the object is small compared to the speed of light, and thus Taylor-expand the mass μ(γ) around γ = 1:

$$\begin{equation}\mu (\gamma )=~{m}\left[1+\sigma (1-\gamma )+\frac{1}{2}{\sigma }^{\prime }{(1-\gamma )}^{2}+\cdots \enspace \right],\end{equation} \tag{ 15 }$$

where $~{m}$, σ and σ' are constant parameters. In particular, the latter two are often referred to as the 'sensitivities' and their derivatives:

$$\begin{equation}\sigma \equiv -\frac{\mathrm{d}\mathrm{ln}\enspace \mu (\gamma )}{\mathrm{d}\mathrm{ln}\enspace \gamma }{\left\vert \right.}_{\gamma =1},\end{equation} \tag{ 16 }$$

$$\begin{equation}{\sigma }^{\prime }\equiv \sigma +{\sigma }^{2}+\frac{{\mathrm{d}}^{2}\mathrm{ln}\enspace \mu (\gamma )}{\mathrm{d}{(\mathrm{ln}\enspace \gamma )}^{2}}{\left\vert \right.}_{\gamma =1}.\end{equation} \tag{ 17 }$$

In order to understand the effect of the sensitivities and their derivatives on the dynamics of binary systems, one can derive the equations of motion simply by varying the point particle action

$$\begin{equation}{S}_{\text{pp}}=-\sum\limits _{A}\int {\mu }_{A}({\gamma }_{A})\mathrm{d}{\tau }_{A},\end{equation} \tag{ 18 }$$

where A is an index identifying the objects, and τ is the proper time. This yields the equation of motion

$$\begin{align}[{\mu }_{A}({\gamma }_{A})-\mu _{A}^{\prime }({\gamma }_{A}){\gamma }_{A}]{a}_{\beta }^{A}& =-\mu _{A}^{\prime }({\gamma }_{A})(-{u}_{A}^{\mu }{\nabla }_{\beta }{U}_{\mu }^{A}+{u}_{A}^{\mu }{\nabla }_{\mu }{U}_{\beta }^{A})\\ & \quad -{\mu }_{A}^{\prime \prime }({\gamma }_{A}){\dot {\gamma }}_{A}({U}_{\beta }^{A}-{\gamma }_{A}{u}_{\beta }^{A}),\end{align} \tag{ 19 }$$

where again the index A identifies the particle under consideration (when used in μ, γ and u ) or at which position the æther field U and its acceleration A are to be computed, the prime denotes a derivative with respect to the function's argument, and the overdot represents a derivative along u (i.e. with respect to the proper time).

Reinstating the dependence on the speed of light c and expanding in PN orders (i.e. for c → ∞), one obtains the 1PN equations of motion for a binary as

$$\begin{align}\frac{\mathrm{d}{\boldsymbol{v}}_{A}}{\mathrm{d}t}& =-\frac{{m}_{B}\boldsymbol{n}}{{r}^{2}}\left\{{\mathcal{G}}_{AB}-\left(3{\mathcal{G}}_{AB}{\mathcal{B}}_{AB}+{\mathcal{D}}_{ABB}\right)\frac{{m}_{B}}{r}\right.\\ & \left.\quad -\frac{1}{2}\left[2{\mathcal{G}}_{AB}^{2}+6{\mathcal{G}}_{AB}{\mathcal{B}}_{(AB)}+2{\mathcal{D}}_{B\hspace{0.5pt}AA}+{\mathcal{G}}_{AB}({\mathcal{C}}_{AB}+{\mathcal{E}}_{AB})\right]\frac{{m}_{A}}{r}\right.\\ & \left.\quad +\frac{1}{2}\left[3{\mathcal{B}}_{AB}-{\mathcal{G}}_{AB}(1+{\mathcal{A}}_{A})\right]{v}_{A}^{2}+\frac{1}{2}(3{\mathcal{B}}_{21}+{\mathcal{G}}_{AB}+{\mathcal{E}}_{AB}){v}_{B}^{2}\right.\\ & \left.\quad -\frac{1}{2}\left(6{\mathcal{B}}_{(AB)}+2{\mathcal{G}}_{AB}+{\mathcal{C}}_{AB}+{\mathcal{E}}_{AB}\right){\boldsymbol{v}}_{A}\cdot {\boldsymbol{v}}_{B}-\frac{3}{2}\left({\mathcal{G}}_{AB}+{\mathcal{E}}_{AB}\right){(\boldsymbol{n}\cdot {\boldsymbol{v}}_{B})}^{2}\right\}\\ & \quad +\frac{{m}_{B}{\boldsymbol{v}}_{A}}{{r}^{2}}\boldsymbol{n}\cdot \left\{\left[3{\mathcal{B}}_{AB}+{\mathcal{G}}_{AB}(1+{\mathcal{A}}_{A})\right]{\boldsymbol{v}}_{A}-3{\mathcal{B}}_{AB}{\boldsymbol{v}}_{B}\right\}\\ & \quad -\frac{1}{2}\frac{{m}_{B}{\boldsymbol{v}}_{B}}{{r}^{2}}\boldsymbol{n}\cdot \left[\left(6{\mathcal{B}}_{(AB)}+2{\mathcal{G}}_{AB}+{\mathcal{C}}_{AB}+{\mathcal{E}}_{AB}\right){\boldsymbol{v}}_{A}\right.\left.-\left(6{\mathcal{B}}_{(AB)}+{\mathcal{C}}_{AB}-{\mathcal{E}}_{AB}\right){\boldsymbol{v}}_{B}\right]\\ & \quad -\frac{1}{2}\frac{{m}_{B}\boldsymbol{w}}{{r}^{2}}\boldsymbol{n}\cdot \left[\left({\mathcal{C}}_{AB}-6{\mathcal{B}}_{[AB]}\right.\right.\left.\left.+{\mathcal{E}}_{AB}-2{\mathcal{G}}_{AB}{\mathcal{A}}_{A}\right){\boldsymbol{v}}_{A}-\left({\mathcal{C}}_{AB}-6{\mathcal{B}}_{[AB]}-{\mathcal{E}}_{AB}\right){\boldsymbol{v}}_{B}\right],\\ \end{align} \tag{ 20 }$$

where indices A ≠ B, r ≡ | x _A − x _B |, n ≡ ( x _A − x _B )/r, v _A ≡ d x /dt, and

$$\begin{align}\begin{aligned}{\mathcal{G}}_{AB}& =\frac{{G}_{N}}{(1+{\sigma }_{A})(1+{\sigma }_{B})},\\ {\mathcal{A}}_{A}& =-\frac{{\sigma }_{A}^{\prime }}{1+{\sigma }_{A}},\\ {\mathcal{B}}_{AB}& ={\mathcal{G}}_{AB}(1+{\sigma }_{A}),\\ {\mathcal{D}}_{ABB}& ={{\mathcal{G}}_{AB}}^{2}(1+{\sigma }_{A}),\\ {\mathcal{C}}_{AB}& ={\mathcal{G}}_{AB}\left[{\alpha }_{1}-{\alpha }_{2}-3\left({\sigma }_{A}+{\sigma }_{B}\right)-{\mathcal{Q}}_{AB}-{\mathcal{R}}_{AB}\right],\\ {\mathcal{E}}_{AB}& ={\mathcal{G}}_{AB}\left[{\alpha }_{2}+{\mathcal{Q}}_{AB}-{\mathcal{R}}_{AB}\right],\\ {\mathcal{Q}}_{AB}& =-\frac{1}{2}\left(\frac{2-{c}_{a}}{2{c}_{\sigma }-{c}_{a}}\right)({\alpha }_{1}-2{\alpha }_{2})({\sigma }_{A}+{\sigma }_{B})+3\left(\frac{2-{c}_{a}}{2{c}_{\sigma }+{c}_{\theta }}\right){\sigma }_{A}{\sigma }_{B},\\ {\mathcal{R}}_{AB}& =\frac{1}{2}\left(\frac{8+{\alpha }_{1}}{{c}_{\omega }+{c}_{\sigma }}\right)\left[-{c}_{\omega }({\sigma }_{A}+{\sigma }_{B})+(1-{c}_{\omega }){\sigma }_{A}{\sigma }_{B}\right].\end{aligned}\end{align} \tag{ 21 }$$

Note that we have defined the 'active' masses ${m}_{A}\equiv {~{m}}_{A}(1+{\sigma }_{A})$, in terms of which the Newtonian acceleration matches the GR result, albeit with a rescaled gravitational constant ${\mathcal{G}}_{AB}$. To derive equation (20) we have also used the PN-expanded solutions for the metric and æther found in [33, 73] (dropping divergent terms due to the point particle approximation, as usual in PN calculations):

$$\begin{align}{g}_{00}& =1-\frac{2{G}_{N}{~{m}}_{1}}{{r}_{1}{c}^{2}}+\frac{1}{{c}^{4}}\left[\frac{2{G}_{N}^{2}{~{m}}_{1}^{2}}{{r}_{1}^{2}}+\frac{2{G}_{N}^{2}{~{m}}_{1}{~{m}}_{2}}{{r}_{1}{r}_{2}}+\frac{2{G}_{N}^{2}{~{m}}_{1}{~{m}}_{2}}{{r}_{1}{r}_{12}}\right.\\ & \left.\quad -\frac{3{G}_{N}{~{m}}_{1}}{{r}_{1}}{v}_{1}^{2}\left(1+{\sigma }_{1}\right)\right]+1{\leftrightarrow}2+\mathcal{O}(1/{c}^{6}),\end{align} \tag{ 22 }$$

$$\begin{equation}{g}_{0i}=-\frac{1}{{c}^{3}}\left[{B}_{1}^{-}\frac{{G}_{N}{~{m}}_{1}}{{r}_{1}}{v}_{1}^{i}+{B}_{1}^{+}\frac{{G}_{N}{~{m}}_{1}}{{r}_{1}}{v}_{1}^{j}{n}_{1}^{j}{n}_{1}^{i}\right]+1{\leftrightarrow}2+\mathcal{O}(1/{c}^{4}),\end{equation} \tag{ 23 }$$

$$\begin{equation}{g}_{ij}=-\left(1+\frac{1}{{c}^{2}}\frac{2{G}_{N}{~{m}}_{1}}{{r}_{1}}\right){\delta }_{ij}+1{\leftrightarrow}2+\mathcal{O}(1/{c}^{4}),\end{equation} \tag{ 24 }$$

$$\begin{equation}{U}^{0}=1+\frac{1}{{c}^{2}}\frac{{G}_{N}{~{m}}_{1}}{{r}_{1}}+1{\leftrightarrow}2+\mathcal{O}(1/{c}^{4})\enspace ,\end{equation} \tag{ 25 }$$

$$\begin{equation}{U}^{i}=\frac{1}{{c}^{3}}\frac{{G}_{N}{~{m}}_{1}}{{r}_{1}}\left({C}_{1}^{-}{v}_{1}^{i}+{C}_{1}^{+}{v}_{1}^{j}{n}_{1}^{j}{n}_{1}^{i}\right)+1{\leftrightarrow}2+\mathcal{O}(1/{c}^{5}),\end{equation} \tag{ 26 }$$

$$\begin{align}{B}_{A}^{{\pm}}& \equiv {\pm}\frac{3}{2}-2{\pm}\frac{1}{4}({\alpha }_{1}-2{\alpha }_{2})\left(1+\frac{2-{c}_{a}}{2{c}_{\sigma }-{c}_{a}}{\sigma }_{A}\right)\\ & \quad -\frac{2{c}_{\omega }}{{c}_{\omega }+{c}_{\sigma }}{\sigma }_{A}-\frac{1}{4}{\alpha }_{1}\left(1+\frac{{c}_{\omega }}{{c}_{\omega }+{c}_{\sigma }}{\sigma }_{A}\right),\end{align} \tag{ 27 }$$

$$\begin{equation}{C}_{A}^{{\pm}}\equiv \frac{1}{4}\left(\frac{8+{\alpha }_{1}}{{c}_{\omega }+{c}_{\sigma }}\right)\hspace{-1.5pt}[{c}_{\omega }-(1-{c}_{\omega }){\sigma }_{A}]{\pm}\frac{2-{c}_{a}}{2}\hspace{-1.5pt}\left(\frac{2{\alpha }_{2}-{\alpha }_{1}}{2(2{c}_{\sigma }-{c}_{a})}\hspace{-1.5pt}+\hspace{-1.5pt}\frac{3{\sigma }_{A}}{2{c}_{\sigma }+{c}_{\theta }}\right)\hspace{-1pt},\end{equation} \tag{ 28 }$$

where r_A ≡ | x − x _A | and n _A ≡ ( x − x _A )/r_A .

Note that the æther solution (25)–(26) has space components Uⁱ vanishing at large distances from the binary, i.e. the equations of motion are valid in a preferred reference frame in which the æther is asymptotically at rest. The dependence of the dynamics on the velocity w of the binary's center of mass with respect to the preferred frame can be reinstated by performing a boost. If w ≪ c, one then obtains [70]

$$\begin{align}\frac{\mathrm{d}{\boldsymbol{v}}_{A}}{\mathrm{d}t}& =-\frac{{m}_{B}\boldsymbol{n}}{{r}^{2}}\left\{{\mathcal{G}}_{AB}-\left(3{\mathcal{G}}_{AB}{\mathcal{B}}_{AB}+{\mathcal{D}}_{ABB}\right)\frac{{m}_{B}}{r}\right.\\ & \left.\quad -\frac{1}{2}\left[2{\mathcal{G}}_{AB}^{2}+6{\mathcal{G}}_{AB}{\mathcal{B}}_{(AB)}+2{\mathcal{D}}_{B\hspace{0.5pt}AA}+{\mathcal{G}}_{AB}({\mathcal{C}}_{AB}+{\mathcal{E}}_{AB})\right]\frac{{m}_{A}}{r}\right.\\ & \left.\quad +\frac{1}{2}\left[3{\mathcal{B}}_{AB}-{\mathcal{G}}_{AB}(1+{\mathcal{A}}_{A})\right]{v}_{A}^{2}+\frac{1}{2}(3{\mathcal{B}}_{21}+{\mathcal{G}}_{AB}+{\mathcal{E}}_{AB}){v}_{B}^{2}\right.\\ & \left.\quad -\frac{1}{2}\left(6{\mathcal{B}}_{(AB)}+2{\mathcal{G}}_{AB}+{\mathcal{C}}_{AB}+{\mathcal{E}}_{AB}\right){\boldsymbol{v}}_{A}\cdot {\boldsymbol{v}}_{B}-\frac{3}{2}\left({\mathcal{G}}_{AB}+{\mathcal{E}}_{AB}\right){(\boldsymbol{n}\cdot {\boldsymbol{v}}_{B})}^{2}\right.\\ & \left.\quad +\frac{1}{2}\left({\mathcal{C}}_{AB}+{\mathcal{G}}_{AB}{\mathcal{A}}_{A}\right){w}^{2}+\frac{1}{2}\left({\mathcal{C}}_{AB}-6{\mathcal{B}}_{[AB]}+{\mathcal{E}}_{AB}\right.\right.\left.\left.+2{\mathcal{G}}_{AB}{\mathcal{A}}_{A}\right){\boldsymbol{v}}_{A}\cdot \boldsymbol{w}\right.\\ & \quad +\frac{1}{2}\left({\mathcal{C}}_{AB}+6{\mathcal{B}}_{[AB]}-{\mathcal{E}}_{AB}\right){\boldsymbol{v}}_{B}\cdot \boldsymbol{w}\\ & \left.\quad +\frac{3}{2}{\mathcal{E}}_{AB}\left[{(\boldsymbol{w}\cdot \boldsymbol{n})}^{2}+2(\boldsymbol{w}\cdot \boldsymbol{n})({\boldsymbol{v}}_{B}\cdot \boldsymbol{n})\right]\right\}\\ & \quad +\frac{{m}_{B}{\boldsymbol{v}}_{A}}{{r}^{2}}\boldsymbol{n}\cdot \left\{\left[3{\mathcal{B}}_{AB}+{\mathcal{G}}_{AB}(1+{\mathcal{A}}_{A})\right]{\boldsymbol{v}}_{A}-3{\mathcal{B}}_{AB}{\boldsymbol{v}}_{B}+{\mathcal{G}}_{AB}{\mathcal{A}}_{A}\boldsymbol{w}\right\}\\ & \quad -\frac{1}{2}\frac{{m}_{B}{\boldsymbol{v}}_{B}}{{r}^{2}}\boldsymbol{n}\cdot \left\{\left(6{\mathcal{B}}_{(AB)}+2{\mathcal{G}}_{AB}+{\mathcal{C}}_{AB}+{\mathcal{E}}_{AB}\right){\boldsymbol{v}}_{A}\right.\\ & \left.\quad -\left(6{\mathcal{B}}_{(AB)}+{\mathcal{C}}_{AB}-{\mathcal{E}}_{AB}\right){\boldsymbol{v}}_{B}+2{\mathcal{E}}_{AB}\boldsymbol{w}\right\}\\ & \quad -\frac{1}{2}\frac{{m}_{B}\boldsymbol{w}}{{r}^{2}}\boldsymbol{n}\cdot \left\{\left({\mathcal{C}}_{AB}-6{\mathcal{B}}_{[AB]}+{\mathcal{E}}_{AB}-2{\mathcal{G}}_{AB}{\mathcal{A}}_{A}\right){\boldsymbol{v}}_{A}\right.\\ & \left.\quad -\left({\mathcal{C}}_{AB}-6{\mathcal{B}}_{[AB]}-{\mathcal{E}}_{AB}\right){\boldsymbol{v}}_{B}-2\left({\mathcal{G}}_{AB}{\mathcal{A}}_{A}-{\mathcal{E}}_{AB}\right)\boldsymbol{w}\right\},\\ \end{align} \tag{ 29 }$$

with which equation (20) of course agrees for w = 0.

The sensitivities and their derivatives enter into the conservative dynamics of the binary system at Newtonian and 1PN order, as can be seen explicitly in equations (20) and (29). In appendix A we will use these equations as starting point for studying in more detail the 1PN dynamics of binaries in Einstein-æther theory. In doing so, we will also amend the calculation of the strong-field preferred-frame parameters ${\hat{\alpha }}_{1}$ and ${\hat{\alpha }}_{2}$ performed in [73], fixing an oversight pointed out by [70] and correcting also a few typos present in [70] itself ¹⁰ .

The sensitivities and their derivatives, however, also enter the dissipative dynamics. In more detail the total energy emitted in GWs (including not only tensor but also scalar and vector æther modes) by a binary in quasi-circular orbits was derived in [33, 73] via a standard multipole expansion and reads

$$\begin{align}\frac{{\dot {E}}_{b}}{{E}_{b}}& =2\left(\frac{{\mathcal{G}}_{12}G{m}_{1}{m}_{2}}{{r}^{3}}\right)\left\{\frac{32}{5}({{\Psi}}_{1}+\mathcal{S}{{\Psi}}_{2}+{\mathcal{S}}^{2}{{\Psi}}_{3}){v}_{21}^{2}\right.\\ & \quad +{\left({s}_{1}-{s}_{2}\right)}^{2}\left.\left[{\zeta }_{2}+2{\zeta }_{3}[{w}^{2}-{(\boldsymbol{w}\cdot \boldsymbol{n})}^{2}]+\frac{18}{5}{{\Psi}}_{3}\enspace {w}^{2}+\left(\frac{6}{5}{{\Psi}}_{3}+36{\zeta }_{1}\right){(\boldsymbol{w}\cdot \boldsymbol{n})}^{2}\right]\right\}\end{align} \tag{ 30 }$$

where s_A ≡ σ_A /(1 + σ_A ) is the rescaled sensitivity for the Ath body; v ₂₁ = v ₂ − v ₁ is the relative velocity of the two bodies; the total (potential and kinetic) energy of the binary is

$$\begin{equation}{E}_{b}=-\frac{{\mathcal{G}}_{12}{m}_{1}{m}_{2}}{2r};\end{equation} \tag{ 31 }$$

w is the velocity of the binary's center of mass with respect to the preferred frame; and we have introduced the definitions

$$\begin{equation}{{\Psi}}_{1}\equiv \frac{1}{{c}_{T}}+\frac{2{c}_{a}{c}_{\sigma }^{2}}{{({c}_{\sigma }+{c}_{\omega }-{c}_{\sigma }{c}_{\omega })}^{2}{c}_{V}}+\frac{3{c}_{a}{(Z-1)}^{2}}{2{c}_{S}(2-{c}_{a})},\end{equation} \tag{ 32 }$$

$$\begin{equation}{{\Psi}}_{2}\equiv \frac{2(Z-1)}{({c}_{a}-2){c}_{S}^{3}}-\frac{2{c}_{\sigma }}{({c}_{\sigma }+{c}_{\omega }-{c}_{\sigma }{c}_{\omega }){c}_{V}^{3}},\end{equation} \tag{ 33 }$$

$$\begin{equation}{{\Psi}}_{3}\equiv \frac{1}{2{c}_{V}^{5}{c}_{a}}+\frac{2}{3{c}_{a}(2-{c}_{a}){c}_{S}^{5}},\hspace{25.0pt}{\zeta }_{1}\equiv \frac{1}{9{c}_{a}{c}_{S}^{5}(2-{c}_{a})},\end{equation} \tag{ 34 }$$

$$\begin{equation}{\zeta }_{2}\equiv \frac{4}{3{c}_{S}^{3}{c}_{a}(2-{c}_{a})}+\frac{4}{3{c}_{a}{c}_{V}^{3}},\hspace{25.0pt}{\zeta }_{3}\equiv \frac{1}{6{c}_{V}^{5}{c}_{a}},\end{equation} \tag{ 35 }$$

$$\begin{equation}Z\equiv \frac{({\alpha }_{1}-2{\alpha }_{2})(1-{c}_{\sigma })}{3(2{c}_{\sigma }-{c}_{a})},\hspace{25.0pt}\mathcal{S}\equiv \frac{{m}_{B}{s}_{A}+{m}_{A}{s}_{B}}{{m}_{A}+{m}_{B}}.\end{equation} \tag{ 36 }$$

Note that the dipole flux is proportional to ζ₂ and to ${({s}_{1}-{s}_{2})}^{2}$ (just like in scalar–tensor theories of the Fierz–Jordan–Brans–Dicke type [26, 28, 72]). Therefore, it may dominate over GR's quadrupole emission at low frequencies, depending on the sensitivities and the coupling parameters of the theory [76].

In order to compute the sensitivities, we start from the observation that the metric and æther solutions for a single point particle (equations (22)–(26) with ${~{m}}_{2}=0$) depend on the sensitivity σ already at linear order in the particle's velocity. Moreover, σ regulates the decay of the metric and æther components at large radii and enters already at $\mathcal{O}(1/r)$. The sensitivity is of course a free parameter in the metric and æther solutions for a point particle, but it can be determined by replacing the point particle with a body of finite size. Once a fully non-linear solution for such a body (e.g. in our case, an NS) has been obtained, one can extract its sensitivity from the asymptotic fall-off of the metric and æther fields. Obviously, since σ appears at linear order in velocity, the NS must be moving relative to the preferred foliation. Here, we follow [73] and consider a star in slow motion with respect to the æther, solve the field equations through linear order in the star's velocity, and extract the sensitivities from the asymptotic decay of the fields.

4.1. Metric ansatz

Here, we consider the case of a non-spinning NS at rest with a background æther field moving relative to it. The system is in a stationary regime, i.e. there is no dependence of the metric and the æther field on the time coordinate. For this configuration, letting vⁱ be the velocity of the star relative to the æther, we consider the following ansatz for the metric and the æther

$$\begin{align}\mathrm{d}{s}^{2}& ={\text{e}}^{\nu (r)}\mathrm{d}{t}^{2}-{\left(1-\frac{2M(r)}{r}\right)}^{-1}\mathrm{d}{r}^{2}-{r}^{2}(\mathrm{d}{\theta }^{2}+{\mathrm{sin}}^{2}\enspace \theta \enspace \mathrm{d}{\phi }^{2})\\ & \quad +2vV(r,\theta )\mathrm{d}t\enspace \mathrm{d}r+2vS(r,\theta )\mathrm{d}t\enspace \mathrm{d}\theta +\mathcal{O}({v}^{2}),\end{align} \tag{ 37 }$$

$$\begin{equation}{U}_{\mu }={\text{e}}^{\nu (r)/2}{\delta }_{\mu }^{t}+vW(r,\theta ){\delta }_{\mu }^{r}+vQ(r,\theta ){\delta }_{\mu }^{\theta }+\mathcal{O}({v}^{2}),\end{equation} \tag{ 38 }$$

and we will set G_N = 1 from here on. Note that M(r) has dimensions of length and M(r) → M_⋆ as r → ∞, with M_⋆ thus being the measured mass of the star. Here we have adopted a coordinate system that is comoving with the fluid elements of the NS, by aligning the time coordinate vector to the fluid four-velocity u^μ [73]. In these comoving coordinates, the fluid elements are at rest while the æther is moving. The fluid four-velocity field is

$$\begin{equation}{u}^{\mu }={\text{e}}^{-\nu /2}{\delta }_{t}^{\mu }.\end{equation} \tag{ 39 }$$

The ansatz of equations (37)–(39) depends on M(r) and ν(r) at $\mathcal{O}({v}^{0})$ and on four potentials V(r, θ), S(r, θ), W(r, θ) and Q(r, θ) at $\mathcal{O}(v)$. However, one can perform a coordinate transformation of the form

$$\begin{equation}{t}^{\prime }=t+vH(r,\theta ),\end{equation} \tag{ 40 }$$

which allows for any one of the four potentials to be set to zero while keeping the ansatz valid at $\mathcal{O}(v)$. Here we choose to set Q = 0 without loss of generality.

4.2. Zeroth order in velocity

From equations (6) and (7) let us derive the field equations at zeroth order in velocity, which we will solve to construct the background NS solution. The (t, t), (r, r) and (θ, θ) components of the field equations are the only non-trivial ones and give three independent equations [73]

$$\begin{align}& 16\frac{\mathrm{d}M}{\mathrm{d}r}-4{c}_{a}r(r-2M)\frac{{\mathrm{d}}^{2}\nu }{\mathrm{d}{r}^{2}}-{c}_{a}r(r-2M){\left(\frac{\mathrm{d}\nu }{\mathrm{d}r}\right)}^{2}\\ & \quad +4{c}_{a}\left(r\frac{\mathrm{d}M}{\mathrm{d}r}-2r+3M\right)\frac{\mathrm{d}\nu }{\mathrm{d}r}=64\pi \rho {r}^{2},\end{align} \tag{ 41 }$$

$$\begin{equation}{c}_{a}{r}^{2}(r-2M){\left(\frac{\mathrm{d}\nu }{\mathrm{d}r}\right)}^{2}+8r(r-2M)\frac{\mathrm{d}\nu }{\mathrm{d}r}-16M=64\pi {r}^{3}P,\end{equation} \tag{ 42 }$$

$$\begin{equation*}4{r}^{2}(r-2M)\frac{{\mathrm{d}}^{2}\nu }{\mathrm{d}{r}^{2}}-({c}_{a}-2){r}^{2}(r-2M){\left(\frac{\mathrm{d}\nu }{\mathrm{d}r}\right)}^{2}\end{equation*}$$

$$\begin{equation}\quad -4r\left(r\frac{\mathrm{d}M}{\mathrm{d}r}-r+M\right)\frac{\mathrm{d}\nu }{\mathrm{d}r}-8r\frac{\mathrm{d}M}{\mathrm{d}r}+8M=64\pi {r}^{3}P,\end{equation} \tag{ 43 }$$

where P(r) and ρ(r) are rescaled NS pressure and density respectively (rescaled because ${G}_{N}=\frac{2G}{2-{c}_{a}}$). These can be expressed as

$$\begin{equation}P\equiv \frac{2-{c}_{a}}{2}~{P},\hspace{25.0pt}\rho \equiv \frac{2-{c}_{a}}{2}~{\rho },\end{equation} \tag{ 44 }$$

with $~{P}$ and $~{\rho }$ representing the pressure and energy density that enter directly in the stress–energy tensor for the matter field, which we take to be of a perfect fluid form

$$\begin{equation}{T}_{\mu \nu }^{\text{mat}}=\left(~{\rho }+~{P}\right){u}_{\mu }{u}_{\nu }-~{P}{g}_{\mu \nu }+\mathcal{O}({v}^{2}).\end{equation} \tag{ 45 }$$

Note, that there is a bijective correspondence between the original parametrization (c_a , c_θ , c_ω , c_σ ) and (α₁, α₂, c_ω , c_σ ) that can be derived from equations (13) and (14). As discussed in section 2, c_σ ≪ 10⁻¹⁵ from gravitational wave observations; we thus rewrite equations (41)–(43) in terms of (α₁, α₂, c_σ , c_ω ). This is justified because the bound on c_σ is much stronger than those on the other parameters, which are constrained by solar-system and stability requirements (absence of gradient instabilities and vacuum Cherenkov radiation, and positive energy) to satisfy

$$\begin{equation}{\alpha }_{1}{< }8{\alpha }_{2}{< }0\enspace ,\hspace{25.0pt}{c}_{\omega }{ >}-\frac{{\alpha }_{1}}{2},\end{equation} \tag{ 46 }$$

with |α₁| ≲ 10⁻⁴ and |α₂| ≲ 10⁻⁷, in the limit c_σ → 0. Upon simplification, we get the modified Tolman–Oppenheimer–Volkoff (TOV) equations

$$\begin{align}\frac{\mathrm{d}M}{\mathrm{d}r}& =\frac{1}{{\alpha }_{1}({\alpha }_{1}+8)r}\left\{-4\sqrt{r-2M}({\alpha }_{1}+8)\sqrt{(-{\alpha }_{1}+8)M-4P{\alpha }_{1}\pi {r}^{3}+4r}\right.\\ & \left.\quad -\left({\alpha }_{1}^{2}+24{\alpha }_{1}+128\right)M-16r\right.\left.{\times}\left({\alpha }_{1}\pi {r}^{2}({\alpha }_{1}+2)P-2\pi {r}^{2}{\alpha }_{1}\rho -\frac{{\alpha }_{1}}{2}-4\right)\right\},\end{align} \tag{ 47 }$$

$$\begin{equation}\frac{\mathrm{d}\nu }{\mathrm{d}r}=\frac{1}{{\alpha }_{1}(r-2M)r}\left[-8\sqrt{r-2M}\sqrt{(-{\alpha }_{1}-8)M-4P{\alpha }_{1}\pi {r}^{3}+4r}+16r-32M\right],\end{equation} \tag{ 48 }$$

$$\begin{equation}\frac{\mathrm{d}p}{\mathrm{d}r}=\frac{1}{{\alpha }_{1}(r-2M)r}4(P+\rho )\left[\sqrt{r-2M}\sqrt{(-{\alpha }_{1}-8)M-4P{\alpha }_{1}\pi {r}^{3}+4r}-2r+4M\right].\end{equation} \tag{ 49 }$$

Modifications to the GR TOV equations can be singled out by expanding the above equations (41)–(43) in a small coupling approximation, i.e. c_a ≪ 1 or α₁ ≪ 1 [57, 66, 73].

4.3. First order in velocity

We derive field equations at first order in velocity from equations (6) and (7), which include the potentials as functions of r and θ, at first order in velocity. We can separate variables in r and θ using a Legendre decomposition [73] to obtain

$$\begin{equation}V(r,\theta )=\sum\limits _{n}{K}_{n}(r){P}_{n}(\mathrm{cos}\enspace \theta ),\end{equation} \tag{ 50 }$$

$$\begin{equation}S(r,\theta )=\sum\limits _{n}{S}_{n}(r)\frac{{P}_{n}(\mathrm{cos}\enspace \theta )}{\mathrm{d}\theta },\end{equation} \tag{ 51 }$$

$$\begin{equation}W(r,\theta )=\sum\limits _{n}{W}_{n}(r){P}_{n}(\mathrm{cos}\enspace \theta ),\end{equation} \tag{ 52 }$$

where P_n is the Legendre polynomial of order n. More details on tensor harmonic decomposition can be found in [65]. By separation of variables, we arrive at $\mathcal{O}(v)$ equations, where only the (t, r) and (t, θ) components of the modified Einstein equations and the r and θ components of the æther field equations are non-trivial. We are only interested in n = 1 component of Legendre decomposition since these functions determine sensitivities and consequently the change in orbital period.

Since c_σ ≪ 10⁻¹⁵ (cf section 2), we proceed with calculations in the limit c_σ → 0 to obtain [73]

$$\begin{align}\frac{\mathrm{d}{S}_{1}}{\mathrm{d}r}& =\frac{1}{{\alpha }_{1}r(r-2M)}\left\{-2{S}_{1}\sqrt{r-2M}({\alpha }_{1}+4)\right.\\ & \left.\quad {\times}\sqrt{(-{\alpha }_{1}-8)M-4P{\alpha }_{1}\pi {r}^{3}+4r}-(r-2M)\right.\\ & \quad \left.{\times}\left({J}_{1}{\text{e}}^{\nu /2}{c}_{\omega }{\alpha }_{1}-(3{\alpha }_{1}+16){S}_{1}+{K}_{1}{\alpha }_{1}({c}_{\omega }-1)\right)\right\},\end{align} \tag{ 53 }$$

$$\begin{align}\frac{\mathrm{d}{K}_{1}}{\mathrm{d}r}& =\frac{1}{{c}_{\omega }(r-2M)({\alpha }_{1}^{2}+(2-2{\alpha }_{2}){\alpha }_{1}-16{\alpha }_{2}){\alpha }_{1}({\alpha }_{1}+8)r}\left\{2\left[(({c}_{\omega }+1){\alpha }_{1}+2{c}_{\omega }){J}_{1}{\text{e}}^{\nu /2}\right.\right.\\ & \left.\quad -(6{\alpha }_{1}+32){S}_{1}+{c}_{\omega }{K}_{1}{\alpha }_{1}\right]({\alpha }_{1}^{2}+(2-2{\alpha }_{2}){\alpha }_{1}-16{\alpha }_{2})({\alpha }_{1}+8)\sqrt{r-2M}\\ & \quad {\times}\sqrt{(-{\alpha }_{1}-8)M-4P{\alpha }_{1}\pi {r}^{3}+4r}+\left[-({\alpha }_{1}+8)(r-2M){\alpha }_{1}r\right.\\ & \quad {\times}(({c}_{\omega }+1){\alpha }_{1}^{2}-2{c}_{\omega }({\alpha }_{2}-1){\alpha }_{1}-16{\alpha }_{2}{c}_{\omega })\left(\frac{\partial {J}_{1}}{\partial r}\right)-16\left(-\frac{1}{8}\left(3\left(\left({c}_{\omega }+\frac{5}{3}\right){\alpha }_{1}^{3}\right.\right.\right.\\ & \quad +\left(\left(-2{\alpha }_{2}+\frac{14}{3}\right){c}_{\omega }-\frac{8{\alpha }_{2}}{3}+\frac{8}{3}\right){\alpha }_{1}^{2}+\left(\left(-\frac{64{\alpha }_{2}}{3}+\frac{16}{3}\right){c}_{\omega }-\frac{64{\alpha }_{2}}{3}\right){\alpha }_{1}\\ & \quad -\left.\left.\frac{128{\alpha }_{2}{c}_{\omega }}{3}\right)({\alpha }_{1}+8)M\right)+\left(\pi ({\alpha }_{1}+2){\alpha }_{1}{r}^{2}\left(\right.({c}_{\omega }+1){\alpha }_{1}^{2}-2{c}_{\omega }({\alpha }_{2}-1){\alpha }_{1}\right.\\ & \quad -16{\alpha }_{2}{c}_{\omega }\left.\right)P-2\pi {\alpha }_{1}{r}^{2}(({c}_{\omega }+1){\alpha }_{1}^{2}-2{c}_{\omega }({\alpha }_{2}-1){\alpha }_{1}-16{\alpha }_{2}{c}_{\omega })\rho \\ & \quad +\frac{1}{4}\left(({\alpha }_{1}+8)\left(\left({c}_{\omega }+\frac{3}{2}\right){\alpha }_{1}^{3}+((-2{\alpha }_{2}+4){c}_{\omega }-2{\alpha }_{2}+2)\right.\right.\enspace {\alpha }_{1}^{2}\\ & \quad +\left.\left.\left.\left.((-20{\alpha }_{2}+4){c}_{\omega }-16{\alpha }_{2}){\alpha }_{1}-32{\alpha }_{2}{c}_{\omega }\right)\right)\right)r\right){J}_{1}\left.\right]{\text{e}}^{\nu /2}\\ & \quad +6\left(-\frac{2{({\alpha }_{1}+8)}^{2}{S}_{1}}{3}+{c}_{\omega }{K}_{1}{\alpha }_{1}\right)({\alpha }_{1}+8)\left(\right.{\alpha }_{1}^{2}+(2-2{\alpha }_{2}){\alpha }_{1}-16{\alpha }_{2}\left.\right)M\\ & \quad -16\left[\left(\pi ({\alpha }_{1}^{2}+(2-2{\alpha }_{2}){\alpha }_{1}-16{\alpha }_{2})({\alpha }_{1}+4)\right.\right.\left.\left.{\alpha }_{1}{r}^{2}P-\frac{11{\alpha }_{2}^{3}}{8}\right.\right.\\ & \quad \left.+(3{\alpha }_{2}-11){\alpha }_{1}^{2}+(40{\alpha }_{2}-16){\alpha }_{1}+128{\alpha }_{2}\right)({\alpha }_{1}+8){S}_{1}+{c}_{\omega }{K}_{1}\left(\right.{\alpha }_{1}^{2}\\ & \quad +(2-2{\alpha }_{2}){\alpha }_{1}-16{\alpha }_{2}\left.\right){\alpha }_{1}\left.\left(\pi {r}^{2}({\alpha }_{1}+2)P-2\pi {r}^{2}\rho +{\alpha }_{1}/4+2\right)\right]r\left.\right\},\end{align} \tag{ 54 }$$

$$\begin{align}\frac{{\mathrm{d}}^{2}{J}_{1}}{\mathrm{d}{r}^{2}}& =\frac{1}{\sqrt{(-{\alpha }_{1}-8)M-4P{\alpha }_{1}\pi {r}^{3}+4r}{(r-2M)}^{3/2}{\alpha }_{1}^{3}{r}^{2}({\alpha }_{1}+8)}\left(\right.4\left\{\left[({\alpha }_{1}+8)\right.\right.\\ & \quad {\times}\left(({\alpha }_{1}^{2}-4{\alpha }_{1}{\alpha }_{2}-32{\alpha }_{2}){S}_{1}\right.\left.+\left(-\frac{{\alpha }_{1}^{2}}{2}+{c}_{\omega }({\alpha }_{2}-1){\alpha }_{1}+8{\alpha }_{2}{c}_{\omega }\right){K}_{1}\right){\alpha }_{1}r{\text{e}}^{-\nu /2}\\ & \quad -12{\alpha }_{1}^{2}\left(\left(\frac{1}{24}{\alpha }_{1}^{2}+{\alpha }_{1}+\frac{16}{3}\right)M+\left({\alpha }_{1}\pi {r}^{2}({\alpha }_{1}+2)P-2\pi {r}^{2}{\alpha }_{1}\rho \right.\right.\left.\left.+\frac{{\alpha }_{1}^{2}}{24}-\frac{8}{3}\right)r\right)\\ & \quad {\times}r\left(\frac{\partial {J}_{1}}{\partial r}\right)+\left(8\left(\frac{\partial \rho }{\partial r}\right){\alpha }_{1}^{3}\pi {r}^{4}+\frac{1}{2}\left(({\alpha }_{1}+8)\right.\right.\left.\left.{\times}\left({\alpha }_{1}^{3}+(-8{\alpha }_{2}+56){\alpha }_{1}^{2}\right.\right.\right.\\ & \quad +\left.\left.(-192{\alpha }_{2}+128){\alpha }_{1}-1024{\alpha }_{2}\right)M\right)+\left(\right.8\pi \left({\alpha }_{1}^{3}+(-5{\alpha }_{2}+14)\right.{\alpha }_{1}^{2}\\ & \left. \quad +\left.(-56{\alpha }_{2}+24){\alpha }_{1}-128{\alpha }_{2}\right){\alpha }_{1}{r}^{2}P+4\pi {\alpha }_{1}^{2}{r}^{2}\left({\alpha }_{1}^{2}+(-2{\alpha }_{2}+16){\alpha }_{1}-16{\alpha }_{2}\right.\right.\\ & \quad +\left.16\right)\rho +\left(\frac{{\alpha }_{1}^{3}}{2}+\left({\alpha }_{2}{c}_{\omega }\right.\right.\left.\left.-{c}_{\omega }-12\right){\alpha }_{1}^{2}+(-32+(8{c}_{\omega }+32){\alpha }_{2}){\alpha }_{1}+256{\alpha }_{2}\right)\\ & \quad {\times}\left.\left.({\alpha }_{1}+8)r\right){J}_{1}\right]\sqrt{r-2M}\sqrt{(-{\alpha }_{1}-8)M-4P{\alpha }_{1}\pi {r}^{3}+4r}\\ & \quad +64\left(\left(\frac{{\alpha }_{1}}{4}+2\right)M+P{\alpha }_{1}\pi {r}^{3}-r\right)\\ & \quad {\times}\left[\frac{1}{8}\left(-{S}_{1}{\alpha }_{1}{\alpha }_{2}r{({\alpha }_{1}+8)}^{2}{\text{e}}^{-\nu /2}+{\alpha }_{1}^{2}r({\alpha }_{1}+8)(r-2M)\frac{\partial {J}_{1}}{\partial r}\right)\right.\left.+{J}_{1}\left(\frac{3({\alpha }_{1}+8)}{4}\right.\right.\\ & \quad {\times}\left({\alpha }_{1}^{2}+\left(\frac{-8{\alpha }_{2}}{3}+\frac{8}{3}\right){\alpha }_{1}-\frac{64{\alpha }_{2}}{3}\right)M+\left({r}^{2}\pi {\alpha }_{1}^{2}({\alpha }_{1}+2)P\right.\\ & \quad +\left.\left.\left.\left.\left.{\alpha }_{1}^{2}\pi {r}^{2}({\alpha }_{1}+2)\rho -\frac{3{\alpha }_{1}^{3}}{8}+({\alpha }_{2}-4){\alpha }_{1}^{2}\right.\right.\right.\right.\right.\left.\left.\left.\left.\left.+(16{\alpha }_{2}-8){\alpha }_{1}+64{\alpha }_{2}\right)r\right)\right]\right\}\right),\end{align} \tag{ 55 }$$

where we have defined J_n = W_n + e^−ν/2 K_n [73]. With the above set of equations at hand, the next section describes the methods of solving these equations at each order in velocity.

The sensitivities are calculated by solving the coupled differential equations in equations (47)–(49) and equations (53)–(55), which are obtained from the modified Einstein and the æther field equations in a v ≪ 1 expansion at $\mathcal{O}({v}^{0})$ and $\mathcal{O}(v)$ respectively [73]. In sections 5.1 and 5.2 we describe and apply two methods to solve these equations and find the NS sensitivities. The first method, outlined in section 5.1, was used previously in reference [73], but we will explain how it leads to unstable solutions in particular regions of parameter space. A second method outlined in section 5.2 provides stable results in all regions of parameter space.

The $\mathcal{O}({v}^{0})$ solutions are common to both methods, as they both involve solving $\mathcal{O}({v}^{0})$ differential equations (47)–(49) numerically once in the interior and then in the exterior of the NS. The initial and boundary conditions to it are obtained by imposing regularity at the NS center, while imposing asymptotic flatness at spatial infinity respectively. The differential equations are solved from a core radius (i.e. some small initial radius) to the stellar surface radius, where the pressure goes to zero. These numerical solutions at the NS surface are now used as initial conditions to solve the exterior evolution equations from the stellar surface to an extraction radius r_b . Using continuity and differentiability of the solutions, the asymptotic solutions at spatial infinity are matched to the numerical solutions evaluated at r_b . This gives the observed mass of the NS and the integration constant corresponding to ν(0) (obtained by solving the $\mathcal{O}({v}^{0})$ differential equations [73]) which will be used in solving the $\mathcal{O}(v)$ equations discussed further.

5.1. Method 1: direct numerical solutions

In this method, the aforementioned $\mathcal{O}(v)$ differential equations are solved in two regions, the interior of the star, and the exterior. The initial conditions at $\mathcal{O}(v)$ are obtained by solving the corresponding differential equations asymptotically about a core radius, while imposing regularity at the core, and asymptotically about spatial infinity, while imposing asymptotic flatness [73]. In both cases, the solutions depend on integration constants—$~{C}$ and $~{D}$ in the interior asymptotic solution and $~{A}$ and $~{B}$ in the exterior asymptotic solution—that must be chosen so as to guarantee that the numerical interior and exterior solutions are continuous and differentiable at the stellar surface, where pressure becomes significantly smaller than their core values.

As defined above, the global solution reduces to finding the right constants $(~{A},~{B},~{C},~{D})$, which in turn is a shooting problem. In practice, reference [73] solved this shooting problem by first picking two sets of values for interior constants, ${ \overrightarrow {c}}_{(1)}=({~{C}}_{(1)},{~{D}}_{(1)})$ and ${ \overrightarrow {c}}_{(2)}=({~{C}}_{(2)},{~{D}}_{(2)})$, and then solving the interior equations twice from the core radius r_c to the NS surface R_⋆ to find the solutions ${ \overrightarrow {f}}_{(1)}^{\text{int}}(r)=[{S}_{1}^{(1,\text{int})}(r),{K}_{1}^{(1,\text{int})}(r),{J}_{1}^{(1,\text{int})},{J}_{1}^{\prime (1,\text{int})}(r)]$ and ${ \overrightarrow {f}}_{2}^{\text{int}}(r)=[{S}_{1}^{(2,\text{int})}(r),{K}_{1}^{(2,\text{int})}(r),{J}_{1}^{(2,\text{int})}(r),{J}_{1}^{\prime (2,\text{int})}]$. Then, each interior numerical solution is evaluated at the stellar surface and used as initial conditions for a numerical evolution in the exterior, leading to two exterior solutions ${ \overrightarrow {f}}_{\hspace{-1.5pt}1}^{\text{ext}}(r)=[{S}_{1}^{(1,\text{ext})}(r),{K}_{1}^{(1,\text{ext})}(r),{J}_{1}^{(1,\text{ext})}(r),{J}_{1}^{\prime (1,\text{ext})}(r)]$ and ${ \overrightarrow {f}}_{\hspace{-1.5pt}2}^{\text{ext}}(r)=[{S}_{1}^{(2,\text{ext})}(r),{K}_{1}^{(2,\text{ext})}(r),{J}_{1}^{(2,\text{ext})}(r),{J}_{1}^{\prime (2,\text{ext})}(r)]$.

The global solutions ${ \overrightarrow {f}}_{\hspace{-1.5pt}1,2}^{\text{glo}}(r)={ \overrightarrow {f}}_{\hspace{-1.5pt}1,2}^{\text{int}}(r)\cup { \overrightarrow {f}}_{\hspace{-1.5pt}1,2}^{\text{ext}}(r)$ are then automatically continuous and differentiable at the surface, but in general they will not satisfy the boundary conditions at spatial infinity. Because of the linear and homogeneous structure of the differential system, one can find the correct global solution through linear superposition

$$\begin{equation}{ \overrightarrow {f}}^{\text{glo}}(r;{C}^{\prime },{D}^{\prime })={C}^{\prime }{ \overrightarrow {f}}_{1}^{\text{glo}}(r)+{D}^{\prime }{ \overrightarrow {f}}_{2}^{\text{glo}}(r),\end{equation} \tag{ 56 }$$

where C' and D' are new constants, chosen to guarantee that ${ \overrightarrow {f}}^{\text{glo}}$ satisfies the correct asymptotic conditions near spatial infinity, which in turn depend on $(~{A},~{B})$, i.e.

$$\begin{equation}{ \overrightarrow {f}}^{\text{glo}}({r}_{b};{C}^{\prime },{D}^{\prime })={ \overrightarrow {f}}^{\text{glo},\infty }({r}_{b};~{A},~{B}),\end{equation} \tag{ 57 }$$

where r_b ≫ R_⋆ is the matching radius, ${ \overrightarrow {f}}^{\text{glo}}({r}_{b};{C}^{\prime },{D}^{\prime })$ is given by equation (56) evaluated at r = r_b (which depends on (C', D')) and ${ \overrightarrow {f}}^{\text{glo},\infty }({r}_{b};~{A},~{B})$ is the asymptotic solution to the differential equations near spatial infinity evaluated at the matching radius (which depends on $(~{A},~{B})$).

With this at hand, one can calculate the NS sensitivities via [73]

$$\begin{equation}\sigma =2~{A}\frac{{\alpha }_{1}}{{\alpha }_{1}+8},\end{equation} \tag{ 58 }$$

where $~{A}$ is the coefficient of 1/r in the near-spatial infinity asymptotic solution of ${W}_{1}^{\text{ext}}$ such that

$$\begin{equation}{W}_{1}^{\text{ext}}(r)=~{A}\frac{{M}_{\star }}{r}+\mathcal{O}\left(\frac{{M}_{\star }^{2}}{{r}^{2}}\right),\end{equation} \tag{ 59 }$$

while we recall that α₁ ≲ 10⁻⁴ (cf section 2). Because of the latter constraint, it is obvious that the sensitivities are essentially controlled by $\sigma \approx ~{A}{\alpha }_{1}$, so the numerical stability of its calculation relies entirely on the numerical stability of the calculation of this coefficient. Unfortunately, as we show below, this calculation is not numerically stable in the region of parameter space we are interested in.

Figure 1 shows S₁ as a function of radius, assuming (α₁, α₂, c_ω ) = (10⁻⁴, 4 × 10⁻⁷, −0.1), and setting (r_c , r_b ) = (10², 2 × 10⁷) cm. Observe that both ${S}_{1}^{(1,\text{glo})}$ and ${S}_{1}^{(2,\text{glo})}$ diverge at spatial infinity, so in order to find an ${S}_{1}^{\text{glo}}$ that is finite at spatial infinity, a very delicate cancellation of large numbers needs to take place. This cancellation needs to lead to $~{A}{\alpha }_{1}\approx 0$ but in general $~{A}{\alpha }_{1}\ne 0$, since $\sigma \approx ~{A}{\alpha }_{1}/4\ll 1\ne 0$, and precisely by how much $~{A}{\alpha }_{1}$ deviates from 0 is what determines the value of the sensitivity. We find in practice that σ is highly sensitive to the accuracy of the numerical algorithm used to solve for ${ \overrightarrow {f}}_{(1,2)}^{\text{glo}}$, as well as the choice of r_c , r_b and the value of p(R_⋆) that defines the stellar surface. Figure 1 is in the parameter region that is outside of interest but it indicates how sensitive the calculations are to the aforementioned cancellation, making it difficult to find numerically stable solutions.

Zoom In Zoom Out Reset image size

5.2. Method 2: post-Minkowskian approach

Given that the first method does not allow us to robustly compute the sensitivities in the regime of interest, we developed a new post-Minkowskian method, which we describe here. In this method, the background O(v⁰) equations are solved by direct integration, as done in method 1. The differential equations at O(v), however, are expanded in compactness $\mathcal{C}$ and solved order by order. This is a post-Minkowskian approximation because the compactness always appears multiplied by G/c², so in this sense it is a weak-field expansion. NSs are not weak-field objects, but their compactness is always smaller than ∼1/3 (and usually between [0.1, 0.3]), so provided enough terms are kept in the series, this approximation has the potential to be valid. Moreover, and perhaps more importantly, we will show below that such a perturbative scheme stabilizes the numerical solution for the NS sensitivities.

The procedure presented above is not technically a standard post-Minkowskian series solution because the background equations (or their solutions) are not expanded in powers of $\mathcal{C}$. Had we expanded in powers of $\mathcal{C}$ everywhere, we would have encountered terms in the differential equations with derivatives of the equation of state (EoS). Such derivatives would introduce numerical noise because 'realistic' (tabulated) EoSs are not usually smooth functions, potentially introducing steep jumps, see, e.g. [75]. By not expanding the $\mathcal{O}({v}^{0})$ differential equations, we are implicitly adding higher order terms in compactness, so this procedure could be seen as a resummation technique.

Through this approach, the differential equations at $\mathcal{O}({v}^{1})$ turn into n sets of differential equations for an expansion carried out to $\mathcal{O}({\mathcal{C}}^{n})$, with n therefore labeling the compactness order. In order to derive these equations, however, one must first establish the order of the background solutions, which can be shown to satisfy

$$\begin{equation}M(r)=\mathcal{O}(\mathcal{C}),\hspace{25.0pt}P(r)=\mathcal{O}({\mathcal{C}}^{2}),\end{equation} \tag{ 60 }$$

$$\begin{equation}\rho (r)=\mathcal{O}(\mathcal{C}),\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad \nu (r)=\mathcal{O}(\mathcal{C}),\end{equation} \tag{ 61 }$$

by looking at the differential equations these functions obey at $\mathcal{O}({v}^{0})$. The metric perturbation functions at $\mathcal{O}({v}^{1})$ are then expanded in powers of compactness through

$$\begin{equation}{Y}_{i}(r)=\sum\limits _{j=1}^{\text{n}}{Y}_{ij}{{\epsilon}}^{j},\end{equation} \tag{ 62 }$$

where Y_i ≡ (S₁, K₁, W₁), is a bookkeeping parameter of $\mathcal{O}(\mathcal{C})$ and j indicates the order of $\mathcal{C}$ to be summed over. We work with W₁(r) instead of J₁(r) to avoid introducing numerical error during the conversion between these two functions.

Using these expansions in the differential equations at $\mathcal{O}(v)$, and re-expanding them in powers of compactness, one finds n sets of differential equations. At $\mathcal{O}(\mathcal{C})$, the differential system becomes

$$\begin{equation}\frac{\mathrm{d}{S}_{11}}{\mathrm{d}r}=\frac{2r{K}_{11}-2r\left({S}_{11}+{\alpha }_{1}{W}_{11}\right)-\left(4+{\alpha }_{1}\right)M}{2{r}^{2}},\end{equation} \tag{ 63 }$$

$$\begin{align}\frac{\mathrm{d}{K}_{11}}{\mathrm{d}r}& =\frac{-1}{2{c}_{\omega }{r}^{2}{\alpha }_{1}}\left[8{c}_{\omega }\pi {r}^{3}{\alpha }_{1}\rho +4{c}_{\omega }r{\alpha }_{1}({K}_{11}-{S}_{11})-4r{\alpha }_{1}^{2}{W}_{11}\right.\\ & \quad +{c}_{\omega }(8{\alpha }_{1}+{\alpha }_{1}^{2}-16{\alpha }_{2}-2{\alpha }_{1}{\alpha }_{2})M-4{c}_{\omega }{r}^{2}{\alpha }_{1}{W}_{11}^{\prime }\\ & \left.\quad -2{r}^{2}({\alpha }_{1}^{2}+{c}_{\omega }{\alpha }_{1}^{2}-16{c}_{\omega }{\alpha }_{2}-2{c}_{\omega }{\alpha }_{1}{\alpha }_{2}){W}_{11}^{\prime }\right],\end{align} \tag{ 64 }$$

$$\begin{equation}\frac{{\mathrm{d}}^{2}{W}_{11}}{\mathrm{d}{r}^{2}}=\frac{2{W}_{11}}{{r}^{2}},\end{equation} \tag{ 65 }$$

where S₁₁(r), K₁₁(r) and W₁₁(r) are metric functions at $\mathcal{O}(\mathcal{C})$, and recall that the density ρ(r) is related to pressure through the EoS as ρ(r) = ρ(P(r)). We note that W₁₁(r) is decoupled from S₁₁(r) and K₁₁(r), and so we can solve equation (65) separately and analytically in the regions r ⩽ r_b and r ⩾ r_b . The solutions to them are ${W}_{11}(r)={~{D}}_{1}{r}^{2}$ for r ⩽ r_b and ${W}_{11}^{\infty }(r)={~{A}}_{1}/r$ for r ⩾ r_b , where ${~{A}}_{1}$ and ${~{D}}_{1}$ are integration constants. By requiring continuity and differentiability of metric functions W₁₁(r) and W₁₁'(r), we match the solutions at the extraction radius r_b . This fixes the values of two integration constants ${~{A}}_{1}=0$ and ${~{D}}_{1}=0$.

The remaining two equations, namely equations (63) and (64), are solved numerically with initial conditions obtained using regularity at the NS center and asymptotic flatness at spatial infinity. At the center, we have

$$\begin{align}{S}_{11}(r)& ={~{C}}_{1}-\frac{1}{120{\alpha }_{1}}\pi \left[(240{\alpha }_{1}+40{\alpha }_{1}^{2}-128{\alpha }_{2}-16{\alpha }_{1}{\alpha }_{2}){\rho }_{c}\right.\\ & \left.\quad +(48{\alpha }_{1}{\alpha }_{2}-72{\alpha }_{1}^{2}-15{\alpha }_{1}^{3}+6{\alpha }_{1}^{2}{\alpha }_{2}){p}_{c}\right]{r}^{2}\\ & \quad +\frac{1}{1260}{\pi }^{2}\left[(-160{\alpha }_{1}-35{\alpha }_{1}^{2}+80{\alpha }_{2}+10{\alpha }_{1}{\alpha }_{2}){\rho }_{c}^{2}\right.\\ & \quad +(288{\alpha }_{2}-576{\alpha }_{1}-126{\alpha }_{1}^{2}+36{\alpha }_{1}{\alpha }_{2}){p}_{c}^{2}\left.\right.\\ & \left.\quad +(336{\alpha }_{2}-672{\alpha }_{1}-147{\alpha }_{1}^{2}+42{\alpha }_{1}{\alpha }_{2}){\rho }_{c}{p}_{c}\right]{r}^{4}+\mathcal{O}({r}^{6}),\end{align} \tag{ 66 }$$

$$\begin{align}{K}_{11}(r)& ={~{C}}_{1}-\frac{1}{120{\alpha }_{1}}\pi \left[(400{\alpha }_{1}+40{\alpha }_{1}^{2}-384{\alpha }_{2}-48{\alpha }_{1}{\alpha }_{2}){\rho }_{c}\right.\\ & \left.\quad +(144{\alpha }_{1}{\alpha }_{2}-96{\alpha }_{1}^{2}-15{\alpha }_{1}^{3}+18{\alpha }_{1}^{2}{\alpha }_{2}){p}_{c}\right]{r}^{2}\\ & \quad +\frac{1}{1260}{\pi }^{2}\left[(400{\alpha }_{2}-240{\alpha }_{1}-35{\alpha }_{1}^{2}+50{\alpha }_{1}{\alpha }_{2}){\rho }_{c}^{2}\right.\\ & \quad +(1440{\alpha }_{2}-864{\alpha }_{1}-126{\alpha }_{1}^{2}+180{\alpha }_{1}{\alpha }_{2}){p}_{c}^{2}\\ & \left.\quad +(1680{\alpha }_{2}-1088{\alpha }_{1}-147{\alpha }_{1}^{2}+210{\alpha }_{1}{\alpha }_{2}){\rho }_{c}{p}_{c}\right]{r}^{4}+\mathcal{O}({r}^{6}),\end{align} \tag{ 67 }$$

where ${~{C}}_{1}$ is an integration constant ¹¹ . At spatial infinity, we have

$$\begin{align}{S}_{11}^{\infty }(r)& =-\frac{{~{B}}_{1}}{2{r}^{3}}-\frac{1}{2{c}_{\omega }r{\alpha }_{1}}\left[{~{A}}_{1}({\alpha }_{1}^{2}-2{c}_{\omega }{\alpha }_{1}-2{\alpha }_{1}^{2}{c}_{\omega }+16{c}_{\omega }{\alpha }_{2}+2{c}_{\omega }{\alpha }_{1}{\alpha }_{2})\right.\\ & \left.\quad +{c}_{\omega }(8{\alpha }_{2}-6{\alpha }_{1}-{\alpha }_{1}^{2}+{\alpha }_{1}{\alpha }_{2}){M}_{\star }\right]+(16{\alpha }_{2}-8{\alpha }_{1}-{\alpha }_{1}^{2}+2{\alpha }_{1}{\alpha }_{2})\frac{{M}_{\star }^{2}}{64{r}^{2}}\\ & \quad +(16{\alpha }_{2}-16{\alpha }_{1}-3{\alpha }_{1}^{2}+2{\alpha }_{1}{\alpha }_{2})\frac{{M}_{\star }^{3}}{192{r}^{3}}\\ & \quad +(8{\alpha }_{2}-2{\alpha }_{1}+{\alpha }_{1}{\alpha }_{2})\mathrm{ln}\left(\frac{r}{{M}_{\star }}\right)\frac{{M}_{\star }^{3}}{48{r}^{3}}+\mathcal{O}\left(\frac{{M}_{\star }^{4}}{{r}^{4}}\right),\end{align} \tag{ 68 }$$

$$\begin{align}{K}_{11}^{\infty }(r)& =\frac{{~{B}}_{1}}{{r}^{3}}+\frac{4{M}_{\star }+2{~{A}}_{1}{\alpha }_{1}+{M}_{\star }{\alpha }_{1}}{2r}-({\alpha }_{1}^{2}+16{\alpha }_{2}+2{\alpha }_{1}{\alpha }_{2})\frac{{M}_{\star }^{2}}{64{r}^{2}}\\ & \quad +(2{\alpha }_{1}-8{\alpha }_{2}-{\alpha }_{1}{\alpha }_{2})\mathrm{ln}\left(\frac{r}{{M}_{\star }}\right)\frac{{M}_{\star }^{3}}{24{r}^{3}}+\mathcal{O}\left(\frac{{M}_{\star }^{4}}{{r}^{4}}\right),\enspace \end{align} \tag{ 69 }$$

where ${~{B}}_{1}$ is an integration constant and M_⋆ is the mass of the star.

We next explain how to solve equations (63) and (64) to construct the solution for S₁₁ and K₁₁. First, homogeneous solutions are given by ${S}_{11}^{\mathrm{hom}}={K}_{11}^{\mathrm{hom}}={~{C}}_{1}$. Next, one can construct particular solutions ${S}_{11}^{\text{part}}$ and ${K}_{11}^{\text{part}}$ by setting ${~{C}}_{1}=0$ and numerically integrate the equations from r_c to R_⋆. We then use the numerically calculated interior solutions, evaluated at R_⋆, as the initial conditions to solve the exterior evolution equations with zero pressure and density from R_⋆ to r_b . The true solutions are simply the sum of the homogeneous and particular solutions, namely

$$\begin{equation}{S}_{11}(r)={~{C}}_{1}+{S}_{11}^{\text{part}},\end{equation} \tag{ 70 }$$

$$\begin{equation}{K}_{11}(r)={~{C}}_{1}+{K}_{11}^{\text{part}}.\end{equation} \tag{ 71 }$$

By requiring continuity and differentiability of all metric functions, we match the true numerical solution to the analytic asymptotic solution in equations (68) and (69) at r_b . Applying this matching condition gives the values of ${~{B}}_{1}$ and ${~{C}}_{1}$.

Now let us focus on the solution to $\mathcal{O}({\mathcal{C}}^{2})$ differential equations. The equation for the metric function W₁₂(r) is

$$\begin{align}\frac{{\mathrm{d}}^{2}{W}_{12}}{\mathrm{d}{r}^{2}}& =\frac{2{W}_{12}}{{r}^{2}}+4\pi r{W}_{11}{\rho }^{\prime }-2\pi \rho \left(-2{K}_{11}+2{S}_{11}\right.\\ & \left.\quad +\left(4+{\alpha }_{1}+\frac{2{\alpha }_{1}}{{c}_{\omega }}\right){W}_{11}-6r{W}_{11}^{\prime }\right)\\ & \quad -2\pi \rho (6+{\alpha }_{1})\frac{M}{r}-{W}_{11}^{\prime }\left(4+{\alpha }_{2}+\frac{8{\alpha }_{2}}{{\alpha }_{1}}\right)\frac{M}{{r}^{2}}\\ & \quad -\left(3{K}_{11}-3{S}_{11}-10{W}_{11}-2{\alpha }_{1}{W}_{11}-\frac{2{\alpha }_{1}{W}_{11}}{{c}_{\omega }}\right)\frac{M}{{r}^{3}}\\ & \quad +\left(5+{\alpha }_{1}+\frac{{\alpha }_{2}}{2}+\frac{4{\alpha }_{2}}{{\alpha }_{1}}\right)\frac{{M}^{2}}{{r}^{4}},\end{align} \tag{ 72 }$$

where ρ'(r) is the derivative of ρ(P(r)) obtained from the EoS. This equation is decoupled from the remaining metric functions at $\mathcal{O}({\mathcal{C}}^{2})$, S₁₂ and K₁₂, and can be solved numerically on its own. The initial condition obtained at the center of the NS is

$$\begin{align}{W}_{12}(r)& ={~{D}}_{2}{r}^{2}+\frac{1}{529200{\alpha }_{1}}{\pi }^{2}{r}^{4}\left[560{\rho }_{c}^{3}\pi {r}^{2}{\alpha }_{1}({\alpha }_{1}+8)(5{\alpha }_{1}+4{\alpha }_{2})\right.\\ & \quad -27{p}_{c}^{2}{\alpha }_{1}^{2}\left(\right.70(8\pi {p}_{c}{r}^{2}-7){\alpha }_{1}^{2}-8\pi {p}_{c}{r}^{2}{\alpha }_{1}({\alpha }_{2}-422)\\ & \quad +49{\alpha }_{1}({\alpha }_{2}-62)+8(49-8\pi {p}_{c}{r}^{2}){\alpha }_{2}\left.\right)-252{\rho }_{c}{p}_{c}{\alpha }_{1}\\ & \quad {\times}\left\{70{p}_{c}\pi {r}^{2}{\alpha }_{1}^{3}+{\alpha }_{1}^{2}(70-{p}_{c}\pi {r}^{2}({\alpha }_{2}-382))\right.\\ & \left.\quad +64(7-4\pi {p}_{c}{r}^{2}){\alpha }_{2}-8{\alpha }_{1}(5\pi {p}_{c}{r}^{2}({\alpha }_{2}+8)-7({\alpha }_{2}+10))\right\}\\ & \quad -12{\rho }_{c}^{2}\left(350\pi {p}_{c}{r}^{2}{\alpha }_{1}^{4}-5\pi {p}_{c}{r}^{2}{\alpha }_{1}^{3}(-226+{\alpha }_{2})-31360{\alpha }_{2}\right.\\ & \left.\left.\quad -784{\alpha }_{1}(-40+(5+8\pi {p}_{c}{r}^{2}){\alpha }_{2})\right.\right.\\ & \left.\left.\quad -8{\alpha }_{1}^{2}(-490+\pi {p}_{c}{r}^{2}(980+103{\alpha }_{2}))\right)\right]+\mathcal{O}({r}^{6}),\end{align} \tag{ 73 }$$

where ${~{D}}_{2}$ is an integration constant. The boundary condition to equation (72) at spatial infinity is

$$\begin{align}{W}_{12}^{\infty }(r)& =\frac{{~{A}}_{2}}{r}+\frac{1}{320{r}^{3}{\alpha }_{1}{c}_{\omega }}\left[{~{A}}_{1}{M}_{\star }(40r-{M}_{\star }{\alpha }_{1})\left(\right.{\alpha }_{1}^{2}+{c}_{\omega }\left(\right.4{\alpha }_{1}^{2}\right.\\ & \quad +{\alpha }_{1}(34-4{\alpha }_{2})-32{\alpha }_{2}\left.\right)\left.\right)+{c}_{\omega }\left(\right.80{M}_{\star }^{2}r({\alpha }_{1}+8)({\alpha }_{1}-{\alpha }_{2})\\ & \left.\quad +{M}_{\star }^{3}{\alpha }_{1}(-4{\alpha }_{1}^{2}+56{\alpha }_{2}+{\alpha }_{1}(-38+7{\alpha }_{2}))\left.\right)\right]+\mathcal{O}\left(\frac{{M}_{\star }^{4}}{{r}^{4}}\right),\end{align} \tag{ 74 }$$

where ${~{A}}_{2}$ is an integration constant. To construct the solution, we first note that the homogeneous solution is given by ${W}_{12}^{\mathrm{hom}}={~{D}}_{2}{r}^{2}$. Next, we set ${~{D}}_{2}=0$ and find the particular solution ${W}_{12}^{\text{part}}(r)$ numerically in the interior of the NS by solving equation (72). This interior solution evaluated at the NS surface now serves as initial conditions to solve the differential equations in the exterior up to the boundary radius r_b . The correct solution in the entire numerical domain is then

$$\begin{equation}{W}_{12}(r)={~{D}}_{2}{r}^{2}+{W}_{12}^{\text{part}}(r),\end{equation} \tag{ 75 }$$

where the values of ${~{A}}_{2}$ and ${~{D}}_{2}$ are obtained using the matching condition at r_b . The equations for S₁₂ and K₁₂ are solved similar to the way equations (63) and (64) are solved, so we omit a more detailed description here for brevity. We can use the above method to solve differential equations at higher order in $\mathcal{C}$.

5.3. Comparison between numerical and analytical approaches

5.3.1. Tabulated APR4 EoS

The sensitivity in the æther theory σ for an isolated NS depends on the EoS chosen, and here we perform the calculations of the previous section for the APR4 tabulated EoS [4]. The results are representative of what one finds with other EoSs.

Equation (58) gives the expression of sensitivity in terms of the integration constant $~{A}$ [73] where $~{A}$ can be expressed as

$$\begin{equation}~{A}\equiv \frac{1}{{M}_{\star }}\sum\limits _{j=2}^{n}{~{A}}_{j}{{\epsilon}}^{j},\end{equation} \tag{ 76 }$$

where n is the order of the compactness expansion, with ${~{A}}_{1}=0$. The coefficients ${~{A}}_{j}$ can be calculated numerically as described in the previous subsection. Notice that the leading contribution to $~{A}$ (and hence to the sensitivities) is of $\mathcal{O}({\epsilon})$, since ${M}_{\star }=\mathcal{O}({\epsilon})$.

The calculation of the sensitivity as described above requires one to choose the truncation order n of the post-Minkowskian expansion. We will choose n by the sensitivities computed from methods 1 and 2 in a regime of parameter space where method 1 yields stable results [73]. In particular, we will focus on the choice (α₁, α₂, c_ω , c_σ ) = (10⁻⁴, 4 × 10⁻⁷, 10⁻⁴, 0). Figure 2 compares the sensitivities computed with the two methods with this parameter choice. Observe that as the order of post-Minkowskian approximation increases (i.e. as n increases), the curves approach the method 1 solution, but in an oscillatory manner. The bottom panel of figure 2 shows the stability of post-Minkowskian method at order n = 3.

Zoom In Zoom Out Reset image size

In the weak field limit, the sensitivity can be well-approximated as the ratio of the binding energy to the NS mass (Ω/M_*) through [33]

$$\begin{equation}{s}_{\text{wf}}=\left({\alpha }_{1}-\frac{2}{3}{\alpha }_{2}\right)\frac{{\Omega}}{{M}_{\star }},\end{equation} \tag{ 77 }$$

where the stellar binding energy Ω is [32]

$$\begin{equation}{\Omega}=-\frac{1}{2}\int {\mathrm{d}}^{3}x\enspace \rho (r)\int {\mathrm{d}}^{3}{x}^{\prime }\frac{\rho ({r}^{\prime })}{\vert {\mathbf{x}-\mathbf{x}}^{\prime }\vert },\end{equation} \tag{ 78 }$$

with r = |x| and r' = |x'|. We can use a Legendre expansion of the Green's function to evaluate this integral, and to leading order in $\mathcal{C}$, we find a result that is identical to that computed in the weak field limit by [32]. This can also be seen numerically in figure 2, where the weak field curve coincides with $\mathcal{O}({\mathcal{C}}^{1})$ post-Minkowskian approximation.

5.3.2. Tolman VII EoS

We now focus on the sensitivity of an NS using the Tolman VII EoS. The latter is an analytic model that accurately describes non-rotating NSs [66] by the energy density profile

$$\begin{equation}\rho (r)={\rho }_{c}\left(1-\frac{{r}^{2}}{{R}_{\star }^{2}}\right).\end{equation} \tag{ 79 }$$

The advantage of using the Tolman VII EoS is that the background solution is known analytically in GR [47, 66]. We expand analytically both the $\mathcal{O}({v}^{0})$ and $\mathcal{O}(v)$ equations order by order in compactness. The sensitivity obtained is then

$$\begin{align}s& =\frac{5}{21}\mathcal{C}\left(-3{\alpha }_{1}+2{\alpha }_{2}\right)\\ & \quad +5\left(\frac{573{\alpha }_{1}^{3}+{\alpha }_{1}^{2}(67669-764{\alpha }_{2})+96416{\alpha }_{2}^{2}+68{\alpha }_{1}{\alpha }_{2}(-2632+9{\alpha }_{2})}{252252{\alpha }_{1}}\right){\mathcal{C}}^{2}\\ & \quad +\frac{1}{1801079280{c}_{\omega }{\alpha }_{1}^{2}}\left\{\right.{(4{\alpha }_{1})}^{2}(8+{\alpha }_{1})\left(\right.36773030{\alpha }_{1}^{2}-39543679{\alpha }_{1}{\alpha }_{2}\\ & \quad +11403314{\alpha }_{2}^{2}\left.\right)+{c}_{\omega }\left[\right.-1970100{\alpha }_{1}^{5}+13995878400{\alpha }_{2}^{3}\\ & \quad +640{\alpha }_{1}{\alpha }_{2}^{2}(-49528371+345040{\alpha }_{2})+5{\alpha }_{1}^{4}(-19596941+788040{\alpha }_{2})\\ & \quad +{\alpha }_{1}^{3}(-2699192440+440184934{\alpha }_{2}-5974000{\alpha }_{2}^{2})\\ & \quad {\times}16{\alpha }_{1}^{2}{\alpha }_{2}(1294533212-29152855{\alpha }_{2}+212350{\alpha }_{2}^{2})\left.\right]\left.\right\}{\mathcal{C}}^{3}+\mathcal{O}({\mathcal{C}}^{4}).\end{align} \tag{ 80 }$$

Note that the above expression is not regular in the limit of α₁ → 0 while keeping α₂ finite or c_ω → 0 while keeping α₁ or α₂ finite. This is a known feature of Einstein-æther theory, which recovers GR only when a certain combination of coupling constants is taken to zero at a specific rate.

With this EoS, the compactness can be expressed as a function on Ω/M_⋆ as

$$\begin{equation}\mathcal{C}=-\frac{7{\Omega}}{5{M}_{\star }}+\frac{35819{\alpha }_{1}{{\Omega}}^{3}}{85800{M}_{\star }^{3}}+\mathcal{O}\left(\frac{{{\Omega}}^{4}}{{M}_{\star }^{4}}\right).\end{equation} \tag{ 81 }$$

Here $\mathcal{C}$ and M_⋆ are the observed values with æther corrections included. With this at hand, we can rewrite the sensitivity as a function of Ω to find

$$\begin{align}s& =\frac{(3{\alpha }_{1}+2{\alpha }_{2})}{3}\frac{{\Omega}}{{M}_{\star }}\\ & \quad +\left(\frac{573{\alpha }_{1}^{3}+{\alpha }_{1}^{2}(67669-764{\alpha }_{2})+96416{\alpha }_{2}^{2}+68{\alpha }_{1}{\alpha }_{2}(9{\alpha }_{2}-2632)}{25740{\alpha }_{1}}\right)\frac{{{\Omega}}^{2}}{{M}_{\star }^{2}}\\ & \quad +\frac{1}{656370000{c}_{\omega }{\alpha }_{1}^{2}}\left\{\right.-4{\alpha }_{1}^{2}({\alpha }_{1}+8)\left[36773030{\alpha }_{1}^{2}-39543679{\alpha }_{1}{\alpha }_{2}\right.\\ & \left.\quad +11403314{\alpha }_{2}^{2}\right]+{c}_{\omega }\left[\right.1970100{\alpha }_{1}^{5}-13995878400{\alpha }_{2}^{3}\\ & \quad -640{\alpha }_{1}{\alpha }_{2}^{2}(-49528371+345040{\alpha }_{2})-5{\alpha }_{1}^{4}(19548109+788040{\alpha }_{2})\\ & \quad -16{\alpha }_{1}^{2}{\alpha }_{2}(1294533212-29152855{\alpha }_{2}+212350{\alpha }_{2}^{2})\\ & \quad +{\alpha }_{1}^{3}(2699192440-309701434{\alpha }_{2}+5974000{\alpha }_{2}^{2})\left.\right]\left.\right\}\frac{{{\Omega}}^{3}}{{M}_{\star }^{3}}+\mathcal{O}\left(\frac{{{\Omega}}^{4}}{{M}_{\star }^{4}}\right).\end{align} \tag{ 82 }$$

Note that this expression matches identically to that of [33] when working to leading order in the binding energy.

One may wonder whether the above analytic expression is capable of approximating the sensitivity when using other EoSs. Figure 3 shows the absolute magnitude of the sensitivity as a function of Ω computed analytically with equation (82), as well as numerically with six other EoSs. Here we have chosen to work in a different region of parameter space, namely (α₁, α₂, c_ω ) = (−10⁻⁴, −4 × 10⁻⁷, 10⁻³), where we obtain a stable smoothly varying sensitivity curve. Observe that the sensitivities differ by less than 3%, exhibiting an approximate universality already discovered in [73] as a function of compactness. Given these results, in all future calculations we will use the analytic sensitivities computed with the Tolman VII EoS.

Zoom In Zoom Out Reset image size

The majority of millisecond pulsars are found in binary and triple systems. The orbital dynamics of these systems modulate the time of arrival of radio waves and allow for precise measurements of the orbital parameters [27, 40, 63, 64]. In this section, we discuss the use of precise orbital parameter data to place constraints on the Lorentz-violating Einstein-æther theory. In GR energy is carried away at quadrupolar order due to propagation of tensor modes whereas in this theory (and many of other modified theories of gravity), one usually finds radiation from extra scalar and vector modes which are responsible for energy loss at dipole order, i.e. −1PN order (cf the term proportional to $\mathcal{E}$ in equation (30)) as compared to GR. Hence, energy is radiated faster than what is predicted in GR. This results in a decrease in the orbital separation and orbital period (P_b ) of the binary. The modified orbital period decay rate ($\dot {{P}_{b}}$) relates to the total energy of the binary, i.e. equation (30) via

$$\begin{equation}\frac{\dot {{P}_{b}}}{{P}_{b}}=-\frac{3}{2}\frac{\dot {{E}_{b}}}{{E}_{b}},\end{equation} \tag{ 83 }$$

suggesting a strong dependence of $\dot {{P}_{b}}$ on the sensitivity of the NS [17]. Since the GR predictions agree with the observed value of $\dot {{P}_{b}}$ within observational uncertainty, this allows for stringent constraints to be placed on æther theory.

6.1. Observations

Following from equation (30) (with LV terms set to zero) and (83), one can relate the post-Keplerian parameter $\dot {{P}_{b}}$ in GR to the Keplerian parameter P_b via [73]

$$\begin{equation}{\left(\frac{\dot {{P}_{b}}}{{P}_{b}}\right)}_{\text{GR}}=-\frac{384\pi }{5}{2}^{2/3}{\left(\frac{\pi ({m}_{1}+{m}_{2})}{{P}_{b}}\right)}^{5/3}\frac{{m}_{1}{m}_{2}}{{({m}_{1}+{m}_{2})}^{2}}\frac{1}{{P}_{b}}.\end{equation} \tag{ 84 }$$

Here m₁ and m₂ are the masses in the binary system. In principle, if we can measure the masses, orbital periods and the orbital period decay rates with some uncertainty and find that they are consistent with GR predictions, then we can place constraints on Einstein-æther theory. In this section we focus on data from the measurements of Keplerian and post Keplerian parameters of four different pulsar systems PSR J1738+0333, PSR J0348+0432, PSR J1012+5307 and PSR J0737−3039 (table 1) and a stellar triple system [67]. The first three are pulsar-white dwarf binaries in orbits with $\mathcal{O}(1{0}^{-7})$ eccentricity and 8.5 h period, 0.17 eccentricity and 4.74 h period and $\mathcal{O}(1{0}^{-7})$ eccentricity and 14.5 h period respectively. The fourth is the double pulsar binary system with 0.088 eccentricity and 2.45 h period. Because of the small eccentricity of these systems, we will ignore it in the following, i.e. we will consider quasi-circular binaries.

Table 1. Orbital parameters as measured for the binary systems studied in this paper. Table shows the estimated values of the parameters and the 1-σ uncertainty in the last digits in parentheses. Here, ${\dot {{P}_{b}}}^{\text{obs}}$ is the observed value of $\dot {{P}_{b}}$.

Pulsar System	m1(M⊙)	m2(M⊙)	Pb (days)	${\dot {{P}_{b}}}^{\text{obs}}$
PSR J1738+0333 [36]	$1.4{6}_{-0.05}^{+0.06}$	$0.18{1}_{-0.007}^{+0.008}$	0.354 790 739 8724(13)	−25.9(3.2) × 10−15
PSR J0348+0432 [6]	2.01(4)	0.172(3)	0.102 424 062 722(7)	−0.273(45) × 10−12
PSR J1012+5307 [21, 52]	1.64(0.22)	0.16(0.02)	0.604 672 723 55(3)	−1.5(1.5) × 10−14
PSR J0737−3039 [51]	1.3381(7)	1.2489(7)	0.102 252 562 48(5)	−1.252(17) × 10−12

6.2. Parameter estimation and Bayesian analysis

Our goal is to constrain the theory parameters using measurements of $\dot {{P}_{b}}$. We discuss briefly the Bayesian formalism with Markov-chain Monte-Carlo (MCMC) exploration used to calculate the posteriors on the model parameters (cf section 6.2.1) and derive robust constraints.

For the parameter estimation, we need the expression for the orbital period decay, which depends on both the relative velocity of the binary constituents v₂₁ and the center-of-mass velocity w of the binary's center of mass with respect to the æther field. A natural choice for the æther field direction is provided by the cosmic microwave background (i.e. the æther is expected to be approximately aligned with the cosmological background time direction). In this case, a typical value for the center-of-mass velocity is w ∼ 10⁻³ [33], which for binary pulsar observations is of the same order as v₂₁. If so, the w-dependent corrections in the rate of change of the binding energy (equation (30)) and orbital period (equation (85)) enter at the same PN order (0PN) as the quadrupole emission terms of GR, but multiplied by either (s₁ − s₂) or powers of it. As such, these w-dependent corrections are negligible for both white dwarf-pulsar systems (for which the dominant term is the −1PN dipole emission) and also for the relativistic double pulsar system (for which s₁ − s₂ ≈ 0 as a result of the similar pulsar masses, which kills both the dipole emission and the w-dependent corrections to quadrupole emission). Therefore, our results are independent of the exact value of w as long as that is of order w ∼ 10⁻³ or smaller [73].

From the above assumption and using equations (83) and (30), the orbital period decay rate in Einstein-æther theory is a function of the individual masses (m₁, m₂), ¹² pulsar radii (R_⋆,1, R_⋆,2), orbital period (P_b ) and coupling constants (α₁, α₂, c_ω ) as shown below

$$\begin{align}\frac{\dot {{P}_{b}}}{{P}_{b}}& =\frac{1}{5{({m}_{1}+{m}_{2})}^{4}{P}_{b}\sqrt{\frac{{\alpha }_{1}}{({\alpha }_{1}-8{\alpha }_{2})}}\sqrt{\frac{-{c}_{\omega }}{{\alpha }_{1}}}{\alpha }_{1}^{3}{c}_{\omega }^{2}}\left\{\right.3\enspace {2}^{8/3}\pi {\left(\frac{\pi ({m}_{1}+{m}_{2})}{{P}_{b}}\right)}^{5/3}\\ & \quad {\times}{m}_{1}{m}_{2}\left[\right.-\frac{1}{12}\left(5{\alpha }_{1}{c}_{\omega }{({m}_{1}+{m}_{2})}^{2}{({s}_{1}-{s}_{2})}^{2}\left({2}^{5/6}{\alpha }_{1}^{2}({\alpha }_{1}+8)\right.\right.\\ & \left.\left.\quad {\times}\sqrt{\frac{{\alpha }_{1}}{({\alpha }_{1}-8{\alpha }_{2})}}-{2}^{13/3}\sqrt{\frac{-{c}_{\omega }}{{\alpha }_{1}}}{c}_{\omega }({\alpha }_{1}-8{\alpha }_{2})\right){\left(\frac{{P}_{b}}{\pi m}\right)}^{2/3}\right)\\ & \quad +{\left(({s}_{1}-1)({s}_{2}-1)\right)}^{2/3}\left(\right.({\alpha }_{1}+8){\alpha }_{1}^{3}\left(-4{c}_{\omega }^{2}{({m}_{1}+{m}_{2})}^{2}\sqrt{\frac{-{c}_{\omega }}{{\alpha }_{1}}}\right.\\ & \left.\quad +\sqrt{2}{\alpha }_{1}{({m}_{1}{s}_{2}+{m}_{2}{s}_{1})}^{2}\right)\sqrt{\frac{{\alpha }_{1}}{({\alpha }_{1}-8{\alpha }_{2})}}\\ & \quad +\frac{2\sqrt{\frac{-{c}_{\omega }}{{\alpha }_{1}}}{({\alpha }_{1}-8{\alpha }_{2})}^{2}{c}_{\omega }^{2}{(({m}_{1}+{m}_{2}){\alpha }_{1}+8{m}_{1}{s}_{2}+8{m}_{2}{s}_{1})}^{2}}{3}\left.\right)\left.\right]\left.\right\},\end{align} \tag{ 85 }$$

where s₁ and s₂ are functions of $\mathcal{C}$ and coupling constants. One may worry that the terms inside the square roots in the above expression may be negative, leading to a complex orbital decay rate, but this is not the case because when α₁ > 0, then c_ω ⩽ −α₁/2, while when α₁ < 0, then c_ω ⩾ −α₁/2. As noted in figure 3, sensitivities are independent of the EoS. Here we choose to work with the Tolman VII EoS since it gives stable analytic solutions for the sensitivities.

There are some phenomenological constraints on the æther coupling constants as discussed in section 2, i.e. |α₁| ≲ 10⁻⁴, |α₂| ≲ 10⁻⁷ (Solar System constraints), α₁ < 0, α₁ < 8α₂ < 0 and c_ω > −α₁/2 (positive energy, absence of vacuum Cherenkov radiation and gradient instabilities) and c_σ ≲ 10⁻¹⁵ (GW constraint) ¹³ . Using these pre-existing constraints and by determining if the estimated value of $\dot {{P}_{b}}$ lies within the range ${\dot {{P}_{b}}}^{\text{obs}}{\pm}\delta {\dot {{P}_{b}}}^{\text{obs}}$ (table 1), we determine the consistency of points in the parameter space with observations.

One important point is that we are not using the pulsar timing data directly [5], but instead we are using existing constraints on $\dot {{P}_{b}},{m}_{1}$ and m₂ derived from the primary pulsar timing data as a sufficient statistic. Unfortunately, the published results only quote values for the individual parameters and their uncertainties, so we do not have access to the joint posterior distributions. For simplicity we assume that the parameter correlations are negligible. To check the impact of this assumption, we compared results with zero correlations with a case with 90% correlation between the parameters, and found that it only changed the results by a maximum of 17%.

6.2.1. Bayesian analysis

We are interested in constructing a posterior distribution on a set of model parameters $\overrightarrow {\lambda }$ = (m₁, m₂, P_b , R_⋆,1, R_⋆,2, α₁, α₂, c_ω ) and using an MCMC algorithm to explore the parameter space. According to Bayes' theorem, the probability density for parameters $\overrightarrow {\lambda }$ given data D and hypothesis H (the theory) is

$$\begin{equation}P( \overrightarrow {\lambda }\vert \text{D},\text{H})=\frac{P(\text{D}\vert \overrightarrow {\lambda },\text{H})P( \overrightarrow {\lambda }\vert \text{H})}{P(\text{D}\vert \text{H})},\end{equation} \tag{ 86 }$$

where $P( \overrightarrow {\lambda }\vert \text{H})$ is called the prior which represents the state of knowledge about the parameters before we analyze the data. $P(\text{D}\vert \overrightarrow {\lambda },\text{H})$ is called the likelihood which describes the probability of measuring data D given the model H and a set of parameters $\overrightarrow {\lambda }$. P(D|H) is called the model evidence which represents the overall normalization factor. In practice it is better to work with log probability densities to better cover the dynamic range of the densities.

We assumed uniform priors on α₂ and c_ω such that −4 × 10⁻⁷ ⩽ α₂ ⩽ 4 × 10⁻⁷ and −10⁵ ⩽ c_ω ⩽ 10⁵ and a Gaussian prior for α₁, m₁, m₂, ${\dot {P}}_{b}$ with mean and standard deviation given by the existing bounds listed in table 1. We use Gaussian priors on R_⋆,1 and R_⋆,2 with mean and standard deviation given 12.4 ± 1.1 km based on LIGO and NICER measurements [25]. While these bounds are derived assuming GR, the corrections due to LV effects are sub-dominant compared to those impacting ${\dot {P}}_{b}$ (cf e.g. footnote 11). Using lunar laser ranging experiments, the bounds on α₁ were obtained to be α₁ = (−0.7 ± 0.9) × 10⁻⁴ [55] (cf also section 2).

The log likelihood function is

$$\begin{equation}\mathrm{ln}(P(\text{D}\vert \overrightarrow {\lambda },\text{H}))\propto -\frac{1}{2}\frac{{\left({\left(\dot {{P}_{b}}/{P}_{b}\right)}^{\text{obs}}-{\left(\dot {{P}_{b}}/{P}_{b}\right)}^{\text{th}}\right)}^{2}}{{\sigma }_{(\dot {{P}_{b}}/{P}_{b})}^{2}},\end{equation} \tag{ 87 }$$

where ${(\dot {{P}_{b}}/{P}_{b})}^{\text{th}}$ is the theoretically predicted value of $(\dot {{P}_{b}}/{P}_{b})$ from the model. With the likelihood and the priors in place, we can find the posterior using an MCMC algorithm.

We start the MCMC simulation near the mean values for the model parameters, calculate the posterior and iterate through these steps. Model parameters are allowed to explore the entire range of parameter space and that gives the joint posterior distribution on all parameters $\overrightarrow {\lambda }$. For the proposal distribution we use the prior distribution for a certain set of parameters, and a relative jump from the current position for the remaining. Proposed jumps are accepted or rejected based on the Metropolis–Hastings acceptance probability

$$\begin{equation}\text{H}=\mathrm{min}\left(\frac{P({ \overrightarrow {\lambda }}_{\text{new}}\vert \mathrm{H})P(\mathrm{D}\vert { \overrightarrow {\lambda }}_{\text{new}},\mathrm{H})Q({ \overrightarrow {\lambda }}_{\text{old}}\vert { \overrightarrow {\lambda }}_{\text{new}})}{P({ \overrightarrow {\lambda }}_{\text{old}}\vert \mathrm{H})P(\mathrm{D}\vert { \overrightarrow {\lambda }}_{\text{old}},\mathrm{H})Q({ \overrightarrow {\lambda }}_{\text{new}}\vert { \overrightarrow {\lambda }}_{\text{old}})},1\right).\end{equation} \tag{ 88 }$$

A random number u ∼ U[0, 1] is drawn, and if H > u the proposed jump is accepted, otherwise it is rejected. This process is repeated multiple times to ensure convergence.

We begin by considering a single observation from the pulsar-white dwarf system PSR J1738+0333. Since the sensitivity of a white dwarf (WD) is negligible compared to the NS we can set s₂ = s_WD = 0 (thus R_⋆,2 is excluded from $\overrightarrow {\lambda }$) but for a double pulsar binary we should have s₂ ≠ 0. Figure 4 shows the prior and posterior distribution on the model parameters. We are recovering our priors on the masses and radii, given that these are well constrained as can be noted from table 1. The distribution on α₁ and α₂ is such that α₁ < 0 and α₂ < 0 from existing constraints. This pulsar system further constrains the value of parameter α₁ by approximately a factor of 2 while the coupling constants α₂ and c_ω remain unconstrained. Notice that the posteriors on α₂ and c_ω are flat and very similar to the priors. Therefore, one cannot model them as Gaussian and construct confidence region, as no information is gained for the values of these parameters.

Zoom In Zoom Out Reset image size

We then consider constraints on the coupling parameters by stacking all four different binary systems from table 1 and computing the joint constraints (figure 5). These joint constraints also restrict the region of α₁ by a factor of 2 better than the existing constraints.

Zoom In Zoom Out Reset image size

In the coming years we expect to have more observations, as the sensitivities of radio telescopes will improve as a result of larger collecting areas (e.g. the square kilometre array (SKA) project [20, 68]), which will allow for discovering more pulsars. Moreover the longer observation time (T) will reduce the error in measurements of $\dot {{P}_{b}}$ by T^−5/2 [14], allowing for more precise measurements of the orbital parameters. One may wonder, whether we can get tighter constraints from finding N similar systems or a single system measured with higher signal-to-noise ratio (SNR). The SNR² grows linearly with number of sources N and the observing time T, and quadratically with the effective collecting area of the radio telescope. Significant improvements in the measurement sensitivity are more likely to come from some combination of larger telescopes and additional observing time than from discovering large numbers of systems similar to those known. Figure 6 illustrates the kind of bounds we will get for a PSR J0348+0432 system if ${\dot {{P}_{b}}}^{\text{obs}}$ matches the GR prediction ${\dot {{P}_{b}}}^{\text{GR}}$ and the uncertainties are tightened by a factor of 10. The improved constraints on $\dot {{P}_{b}}$ translate directly into similarly improved bounds on α₁. We also considered an alternative scenario, in which the uncertainties in $\dot {{P}_{b}}$ improved by a factor of 10, but stayed centered on the current observed value. As shown in figure 7, this leads to a value for α₁ bounded away from zero. In other words, in a scenario where the observed period derivative stays at the current value while the uncertainty drops by a factor of ten, we would find that Einstein-æther theory would be favored over GR!

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.3. Constraints from the triple system

Next we have constraints coming from a pulsar in a stellar triple system PSR J0337+1715 consisting of an inner millisecond pulsar-white dwarf binary and a second WD in an outer orbit [7]. Due to the gravitational pull of the outer WD, the pulsar and the inner WD experience accelerations that differ fractionally. If the strong equivalence principle is violated (as a result of the sensitivities), the triple system constrains the fractional acceleration difference parameter δ_a to (+0.5 ± 1.8) × 10⁻⁶ [67]. The relation between δ_a and the sensitivity parameter σ_pulsar (before rescaling) in Einstein-æther theory is [9, 70]

$$\begin{equation}\vert {\delta }_{a}\vert =\left\vert \frac{{\sigma }_{\text{pulsar}}}{1+{\sigma }_{\text{pulsar}}/2}\right\vert \approx \vert {\sigma }_{\text{pulsar}}\vert ,\end{equation} \tag{ 89 }$$

as can be obtained directly from equation (29) (in the Newtonian limit).

We use MCMC simulations in Bayesian analysis similar to that for the $\dot {{P}_{b}}$ constraint and with the likelihood

$$\begin{equation}P(\text{D}\vert \overrightarrow {\lambda },\text{H})\propto \mathrm{exp}\left(-\frac{1}{2}\frac{{\left({\sigma }_{\text{pulsar}}^{\text{obs}}-{\sigma }_{\text{pulsar}}^{\;\text{th}}\right)}^{2}}{{\sigma }_{({\sigma }_{\text{pulsar}})}^{2}}\right),\end{equation} \tag{ 90 }$$

where ${\sigma }_{\text{pulsar}}^{\text{obs}}=(+0.5{\pm}1.8){\times}1{0}^{-6}$ from equation (89) and ${\sigma }_{\text{pulsar}}^{\text{th}}$ is given by equation (58), to constrain the model parameters. Figure 8 shows the joint pulsar and triple system constraints on the model parameters $\overrightarrow {\lambda }$ assuming uniform distribution in α₂ and c_ω . The 95% upper limit on α₁, which was −2.4 × 10⁻⁴ (from the prior constraints) has now shifted to α₁ = −2.4 × 10⁻⁵. It shows that the preferred frame parameter α₁ is constrained by a factor of 10 better than the lunar laser ranging experiments.

Zoom In Zoom Out Reset image size

Table 2 shows bounds on α₁ from binary and triple systems mentioned in this paper. The data from joint binary + triple system allows us to put a stringent constraint on α₁, which is an order of magnitude stronger than the bounds from lunar laser ranging experiments [55, 71].

Table 2. Our bounds on α₁ from different pulsar systems shown in figures 4–8 with 1-σ uncertainty. The first half shows the bounds from existing measurements, while the second half shows projected future bounds assuming that the measurement error on ${\dot {P}}_{b}^{\text{obs}}$ reduces by a factor of 10 with the central value of ${\dot {P}}_{b}^{\text{obs}}$ at the GR predicted value (−27.3 × 10⁻¹⁴) and at the current measured value (−25.3 × 10⁻¹⁴).

Pulsar system	α1
PSR J1738+0333	(−3.975 ± 2.968) × 10−5
Joint binary system	(−4.073 ± 2.936) × 10−5
Joint binary + triple system	(−1.111 ± 0.674) × 10−5
PSR J0348+0432 (${\dot {{P}_{b}}}^{\text{obs}}=-27.3{\times}1{0}^{-14}$)	(−8.119 ± 4.622) × 10−5
PSR J0348+0432 (${\dot {{P}_{b}}}^{\text{obs}}=-25.3{\times}1{0}^{-14}$)	(−1.729 ± 1.805) × 10−5

We have investigated Einstein-æther theory in the context of binary pulsars and NSs. We have recalculated the sensitivities in the regime of coupling parameter space that still survives after the recent measurement of the speed of GWs. This required the development of a new post-Minkowskian approach that allows for stable numerical evaluation of the sensitivities, in addition to the derivation a closed form analytic solution for the Tolman VII EoS. We used these results to place a constraint on certain coupling constants of Einstein-æther theory using Bayesian analysis of binary pulsar observations, including recent observations on the triple system. We find that these data allows for constraints on a certain combination of the coupling constants, α₁, of $\mathcal{O}(1{0}^{-5})$, improving current Solar System constraints by one order of magnitude.

The work carried out here opens the door to several avenues for future research. One such avenue is to use gravitational wave data directly to place constraints on Einstein-æther theory, now that the sensitivities have been analytically calculated. This can be done today to leading PN order in the inspiral, and it remains to be seen whether it is enough to lead to interesting constraints. To include the very late inspiral and merger phase, numerical simulations of coalescing NSs would have to be carried out in Einstein-æther theory. However, since the parameter space of the theory is already quite well constrained, it is not clear whether stronger bounds can be achieved with gravitational wave data.

Another avenue for future research concerns computing sensitivities for BHs. This has been done in khronometric theory [59] but not yet in Einstein-æther theory. Once the BH sensitivities are in hand, and assuming they do not vanish, one could use the existing GW data for binary BH mergers to constrain Einstein-æther theory, including the dipole radiation effect in the gravitational waveform.

One more avenue for future work would be along the lines of improving the analysis in this paper by directly analyzing binary pulsar data and carrying out a parameter estimation and model selection study with a GR and a non-GR timing model. For this, it would be ideal to compute the derivative of the NS sensitivities that enter in the conservative post-Keplerian parameters, such as the periastron precession and Shapiro time delay.

We would like to thank Clifford Will and Ted Jacobson for many insightful discussions, which motivated part of this work. TG and NJC acknowledge the support of NASA EPSCoR Grant MT-80NSSC17M0041 and NSF Grant PHY-1912053. KY acknowledges support from NSF Grant PHY-1806776, NASA Grant No. 80NSSC20K0523, a Sloan Foundation Research Fellowship and the Owens Family Foundation. KY would like to also thank the support by the COST Action GWverse CA16104 and JSPS KAKENHI Grant No. JP17H06358. NY acknowledges support from NASA Grant Nos. NNX16AB98G, 80NSSC17M0041, 80NSSC18K1352 and NSF Grant PHY-1759615. EB and MHV acknowledge financial support provided under the European Union's H2020 ERC Consolidator Grant 'GRavity from Astrophysical to Microscopic Scales' Grant Agreement No. GRAMS-815673.

The data that support the findings of this study are available upon reasonable request from the authors.

In this appendix, we start from the 1PN acceleration (29) and analyze the effect of the 1PN conservative dynamics on the orbital parameters of a binary of compact objects. We will follow the osculating-orbits technique of [70], which will lead us to amend the calculation of the preferred frame parameters ${\hat{\alpha }}_{1}$ and ${\hat{\alpha }}_{2}$ presented in [73]. In doing so, we will also correct a few typos that we found in the expressions of [70]. ¹⁴

The relative acceleration between the two gravitating bodies is obtained by simply letting

$$\begin{equation}\boldsymbol{a}=\frac{\mathrm{d}{\boldsymbol{v}}_{1}}{\mathrm{d}t}-\frac{\mathrm{d}{\boldsymbol{v}}_{2}}{\mathrm{d}t},\end{equation} \tag{ A.1 }$$

while the position of the center of mass is not accelerated. Thus, we can set $\dot {\boldsymbol{X}}=\boldsymbol{X}=0$ at Newtonian order without any loss of generality, getting

$$\begin{equation}{\boldsymbol{x}}_{1}=\left(\frac{{m}_{2}}{m}+\mathcal{O}({\epsilon})\right)\boldsymbol{x},\end{equation} \tag{ A.2 }$$

$$\begin{equation}{\boldsymbol{x}}_{2}=-\left(\frac{{m}_{1}}{m}+\mathcal{O}({\epsilon})\right)\boldsymbol{x},\end{equation} \tag{ A.3 }$$

where ${\epsilon}\sim m/r\sim {v}_{21}^{2}$ is a book-keeping parameter that counts PN order.

Here the acceleration of every individual body is given by (29). Hereinafter we will borrow the notation of [70] and thus we define

$$\begin{equation}m={m}_{1}+{m}_{2},\hspace{25.0pt}\eta =\frac{{m}_{1}{m}_{2}}{{m}^{2}},\hspace{25.0pt}{\Delta}=\frac{{m}_{2}-{m}_{1}}{m},\end{equation} \tag{ A.4 }$$

and the functions of the sensitivities

$$\begin{align}& \mathcal{G}={\mathcal{G}}_{12},\hspace{25.0pt}{\mathcal{B}}_{+}={\mathcal{B}}_{(12)},\hspace{25.0pt}{\mathcal{B}}_{-}={\mathcal{B}}_{[12]},\hspace{25.0pt}\mathcal{D}=\frac{{m}_{2}}{m}{\mathcal{D}}_{122}+\frac{{m}_{1}}{m}{\mathcal{D}}_{211},\\ & \mathcal{E}={\mathcal{E}}_{12},\hspace{25.0pt}{\mathcal{A}}^{(n)}={\left(\frac{{m}_{2}}{{m}_{1}}\right)}^{n}{\mathcal{A}}_{1}-{\left(-\frac{{m}_{1}}{m}\right)}^{n}{\mathcal{A}}_{2}.\end{align} \tag{ A.5 }$$

Using this, the relative acceleration can be written in a compact form

$$\begin{equation}\boldsymbol{a}={\boldsymbol{a}}_{\mathrm{L}}+{\boldsymbol{a}}_{\text{PF}},\end{equation} \tag{ A.6 }$$

where we have separated the purely local contributions and those coming from preferred frame effects. The former reads

$$\begin{equation}{\boldsymbol{a}}_{\mathrm{L}}=\frac{m}{{r}^{2}}\left[\boldsymbol{n}\left({\hat{A}}_{1}{v}_{21}^{2}+{\hat{A}}_{2}{\dot {r}}^{2}+{\hat{A}}_{3}\frac{m}{r}\right)+\dot {r}\hat{B}{\boldsymbol{v}}_{21}\right],\end{equation} \tag{ A.7 }$$

where

$$\begin{equation}{\boldsymbol{v}}_{21}={\dot {\boldsymbol{x}}}_{2}-{\dot {\boldsymbol{x}}}_{1},\end{equation} \tag{ A.8 }$$

$$\begin{equation}{\hat{A}}_{1}=\frac{1}{2}\left[\mathcal{G}(1-6\eta )-3{\mathcal{B}}_{+}-3{\Delta}{\mathcal{B}}_{-}-\eta ({\mathcal{C}}_{12}+2\mathcal{E})+\mathcal{G}{\mathcal{A}}^{(3)}\right],\end{equation} \tag{ A.9 }$$

$$\begin{equation}{\hat{A}}_{2}=\frac{3\eta }{2}(\mathcal{G}+\mathcal{E}),\end{equation} \tag{ A.10 }$$

$$\begin{equation}{\hat{A}}_{3}=\mathcal{D}+\mathcal{G}\left[2\eta \mathcal{G}+3{\mathcal{B}}_{+}+\eta ({\mathcal{C}}_{12}+\mathcal{E})+3{\Delta}{\mathcal{B}}_{-}\right],\end{equation} \tag{ A.11 }$$

$$\begin{equation}\hat{B}=\mathcal{G}(1-2\eta )+3{\mathcal{B}}_{+}+3{\Delta}{\mathcal{B}}_{-}+\eta \mathcal{G}+\mathcal{G}{\mathcal{A}}^{(3)}.\end{equation} \tag{ A.12 }$$

These expressions agree with those of [70]. However, we find a difference in the acceleration due to preferred frame effects

$$\begin{align}{\boldsymbol{a}}_{\text{PF}}& =\frac{m}{{r}^{2}}\left\{-\boldsymbol{n}\left[\left(\frac{{\hat{\alpha }}_{1}}{2}+2\mathcal{G}{\mathcal{A}}^{(2)}\right)(\boldsymbol{\omega }\cdot {\boldsymbol{v}}_{21})+\frac{3}{2}\left({\hat{\alpha }}_{2}+\mathcal{G}{\mathcal{A}}^{(1)}\right){(\boldsymbol{\omega }\cdot \boldsymbol{n})}^{2}\right]\right.\\ & \left.\quad -\boldsymbol{\omega }\left[\frac{{\hat{\alpha }}_{1}}{2}(\boldsymbol{n}\cdot {\boldsymbol{v}}_{21})+{\hat{\alpha }}_{2}(\boldsymbol{n}\cdot \boldsymbol{\omega })\right]+\mathcal{G}{\mathcal{A}}^{(2)}{\boldsymbol{v}}_{21}(\boldsymbol{n}\cdot \boldsymbol{\omega })\right\}-\frac{m{\omega }^{2}}{2{r}^{2}}\left({\mathcal{C}}_{12}+\mathcal{G}{\mathcal{A}}^{(1)}\right)\boldsymbol{n}.\end{align} \tag{ A.13 }$$

where we have already specified the generic boost velocity w in (29) to match the velocity of the preferred frame ω .

This differs from the result of [70] in a sign multiplying the first whole line, as well as in the last term, which is absent in [70]. Here we have defined the following compact-body effective PPN parameters

$$\begin{equation}{\hat{\alpha }}_{1}={\Delta}({\mathcal{C}}_{12}+\mathcal{E})-6{\mathcal{B}}_{-}-2\mathcal{G}{\mathcal{A}}^{(2)},\hspace{25.0pt}{\hat{\alpha }}_{2}=\mathcal{E}-\mathcal{G}{\mathcal{A}}^{(1)}.\end{equation} \tag{ A.14 }$$

These are the strong field versions of the parameters α₁ and α₂, which are contained inside the definition of the calligraphic objects in (A.14), so that ${\hat{\alpha }}_{1}$ and ${\hat{\alpha }}_{2}$ are implicit functions of them. They are directly proportional to each other only in the case in which the sensitivities vanish exactly.

In the absence of PN corrections, the motion of the two-body system describes a Keplerian orbit, parametrized by x = r n with

$$\begin{align}\boldsymbol{n}& =[-\mathrm{cos}\enspace {\Omega}\mathrm{cos}(\omega +f)-\mathrm{cos}\enspace \iota \enspace \mathrm{sin}\enspace {\Omega}\mathrm{sin}(\omega +f)]{\boldsymbol{e}}_{X}+\left[\right.\mathrm{sin}\enspace {\Omega}\enspace \mathrm{cos}(\omega +f)\\ & \quad +\mathrm{cos}\enspace \iota \enspace \mathrm{cos}\enspace {\Omega}\enspace \mathrm{sin}(\omega +f)\left.\right]{\boldsymbol{e}}_{Y}+\mathrm{sin}\enspace \iota \enspace \mathrm{sin}(\omega +f){\boldsymbol{e}}_{Z},\end{align} \tag{ A.15 }$$

where the orbital elements are: inclination ι, longitude of the ascending node Ω and pericenter angle ω. The element f = ω − ϕ is the true anomaly, with ϕ the orbital phase measured from the ascending node. The reference vectors e _i form an orthonormal basis.

When the extra force (A.6) is included, Keplerian orbits are not solutions to the equations of motion anymore. However, provided that the force is small enough relative to the Newtonian force, we can use perturbation theory and translate the dependence on time of the motion to the orbital parameters. This is the method of osculating orbits described in [71], which leads to a secular variation of the orbital elements under the effect of a . In order to parametrize this change in terms of the velocity vector of the preferred frame, we decompose the latter by projecting it onto the orbital plane by defining

$$\begin{equation}{\omega }_{P}=\boldsymbol{\omega }\cdot {\boldsymbol{e}}_{P},\hspace{25.0pt}{\omega }_{Q}=\boldsymbol{\omega }\cdot {\boldsymbol{e}}_{Q},\hspace{25.0pt}{\omega }_{Z}=\boldsymbol{\omega }\cdot \boldsymbol{z},\end{equation} \tag{ A.16 }$$

as well as onto the angular momentum vector

$$\begin{equation}{\omega }_{h}=\boldsymbol{\omega }\cdot \boldsymbol{h}={\omega }_{Z}\sqrt{\mathcal{G}mp}.\end{equation} \tag{ A.17 }$$

Following the computation in [71], we thus find that the local terms in a _L induce a change only on the pericenter angle ω, which in an orbit changes by

$$\begin{equation}{{\Delta}}_{\mathrm{L}}\omega =\frac{6\pi m}{\mathcal{G}p}\left[\mathcal{G}{\mathcal{B}}_{+}+\frac{1}{6}\left({\mathcal{G}}^{2}-\mathcal{D}\right)+\frac{1}{6}\mathcal{G}\left(6{\Delta}{\mathcal{B}}_{-}+\eta (2{\mathcal{C}}_{12}+\mathcal{E})+\mathcal{G}{\mathcal{A}}^{(3)}\right)\right],\end{equation} \tag{ A.18 }$$

and is of course independent of ω . Again, this agrees with [70] up to a typographical error in their result. The rest of secular changes vanish, either because they are identically zero or because they compensate along the orbit.

On the other hand, the force induced by preferred frame effects produces a secular change in all orbital parameters

$$\begin{equation}{{\Delta}}_{\text{PF}}a=\frac{2\pi e{\omega }_{P}}{{(1-{e}^{2})}^{2}}{\left(\frac{mp}{\mathcal{G}}\right)}^{\frac{1}{2}}\left({\hat{\alpha }}_{1}+4{\mathcal{A}}^{(2)}\mathcal{G}\right),\end{equation} \tag{ A.19 }$$

$$\begin{equation}{{\Delta}}_{\text{PF}}\iota =\pi {\hat{\alpha }}_{1}{\left(\frac{m}{\mathcal{G}p}\right)}^{\frac{1}{2}}{\omega }_{h}\enspace \mathrm{sin}(\omega )eF(e)-\frac{2\pi {\hat{\alpha }}_{2}{\omega }_{h}{\omega }_{R}F(e)}{\mathcal{G}\sqrt{1-{e}^{2}}},\end{equation} \tag{ A.20 }$$

$$\begin{equation}{{\Delta}}_{\text{PF}}{\Omega}=-\pi {\hat{\alpha }}_{1}{\left(\frac{m}{\mathcal{G}p}\right)}^{\frac{1}{2}}\frac{{\omega }_{h}}{\mathrm{sin}(\iota )}\enspace \mathrm{cos}(\omega )eF(e)-\left.\right)\frac{2\pi {\hat{\alpha }}_{2}{\omega }_{h}{\omega }_{S}F(e)}{\mathcal{G}\enspace \mathrm{sin}(\iota )\sqrt{1-{e}^{2}}},\end{equation} \tag{ A.21 }$$

$$\begin{align}{{\Delta}}_{\text{PF}}\varpi & =-\pi {\hat{\alpha }}_{1}{\left(\frac{m}{\mathcal{G}p}\right)}^{\frac{1}{2}}{\omega }_{Q}\frac{\sqrt{1-{e}^{2}}F(e)}{e}-\pi {\hat{\alpha }}_{2}({\omega }_{P}^{2}{\omega }_{Q}^{2})F{(e)}^{2}\\ & \quad +\frac{\pi {\omega }_{Q}}{e}{\left(\frac{m}{\mathcal{G}p}\right)}^{\frac{1}{2}}({\hat{\alpha }}_{1}+4{\mathcal{A}}^{(2)}\mathcal{G}),\end{align} \tag{ A.22 }$$

$$\begin{align}{{\Delta}}_{\text{PF}}e& =-\pi {\hat{\alpha }}_{1}{\left(\frac{m}{\mathcal{G}p}\right)}^{\frac{1}{2}}{\omega }_{P}(1-{e}^{2})F(e)+2\pi {\hat{\alpha }}_{2}{\omega }_{P}{\omega }_{Q}e\sqrt{1-{e}^{2}}F{(e)}^{2}\\ & \quad +\pi {\omega }_{P}{\left(\frac{m}{\mathcal{G}p}\right)}^{\frac{1}{2}}({\hat{\alpha }}_{1}+4{\mathcal{A}}^{(2)}\mathcal{G}),\end{align} \tag{ A.23 }$$

where Δ_PF ϖ = Δ_PF ω + cos(ι)Δ_PFΩ and

$$\begin{equation}F(e)=\frac{1}{1+\sqrt{1-{e}^{2}}},\end{equation} \tag{ A.24 }$$

$$\begin{equation}{\omega }_{R}={\omega }_{P}\enspace \mathrm{cos}\enspace \omega -{\omega }_{Q}\sqrt{1-{e}^{2}}\enspace \mathrm{sin}\enspace \omega ,\end{equation} \tag{ A.25 }$$

$$\begin{equation}{\omega }_{S}={\omega }_{P}\enspace \mathrm{sin}\enspace \omega +{\omega }_{Q}\sqrt{1-{e}^{2}}\enspace \mathrm{cos}\enspace \omega .\end{equation} \tag{ A.26 }$$

Out of these deviations, the most relevant one is the variation of the semimajor axis, which can be related to the change in the period of the orbit by using Kepler's third law

$$\begin{equation}\frac{{\Delta}T}{T}=\frac{3}{2}\frac{{\Delta}a}{a}.\end{equation} \tag{ A.27 }$$

Note however that this change is sub-leading with respect to the change expected from emission of gravitational radiation in a binary system like the one considered throughout this paper (cf equation (30)), which is actually the dominant factor.