Constraints on Einstein-\AE ther theory and Horava gravity from binary pulsar observations

Binary pulsars are ideal to test the foundations of General Relativity, such as Lorentz symmetry, which requires that experiments produce the same results in all free-falling (i.e.inertial) frames. We here break this symmetry in the gravitational sector by specifying a preferred time direction, and thus a preferred frame, at each spacetime point. We then examine the consequences of this gravitational Lorentz symmetry breaking in the orbital evolution of binary pulsars, focusing on the dissipative effects. We find that Lorentz symmetry breaking modifies these effects, and thus the orbital dynamics, in two different ways. First, it generically causes the emission of dipolar radiation, which makes the orbital separation decrease faster than in General Relativity. Second, the quadrupole component of the emission is also modified. The orbital evolution depends critically on the sensitivities of the stars, which measure how their binding energies depend on the motion relative to the preferred frame. We calculate the sensitivities numerically and compute the predicted orbital decay rate of binary pulsars in Lorentz-violating gravity. By testing these predictions against observations, we place very stringent constraints on gravitational Lorentz violation.


I. INTRODUCTION
The social scientist and epistemologist Karl Popper believed that scientific hypotheses could never be proved correct, but rather could only be disproved [1]. In his view, the role of a scientist is to attempt to disprove the canonical set of hypotheses of the day, the status quo, by whatever means possible (experimental, observational or mathematical). These efforts are an important engine for scientific progress -in particular to generate what Kuhn would later call scientific revolutions that leap us forward in our understanding of the Universe [2].
General Relativity (GR) is today the status quo when it comes to describing the gravitational interaction and the motion of large bodies. As Popper suggested, scientists have spent decades putting Einstein's theory to the test, essentially since its conception. Einstein himself would repeatedly look for observational confirmation of the predictions that his theory would make. Today, we have an overwhelming amount of data that confirms GR to incredible precision [3], and one may wonder why we should bother with alternative theories.
Three reasons come immediately to mind. First, almost all the existing tests of GR involve systems that are governed by the weakly-gravitating, mildly-relativistic regime of the Einstein equations. In the Solar System, typical orbital velocities are much below 0.1% of the speed of light, and curvatures are barely measurable. Stronger tests of GR come from binary pulsars, e.g. from observations of their orbital decay rate or time-dilation effects when photons from the secondary graze the surface of the primary. The orbital velocity in observed pulsar binaries, however, is not very different from velocities in the Solar System. In fact, their shortest orbital periods are on the order of hours, and thus orbital velocities are still 1 % the speed of light, although curvatures are much higher inside pulsars than in the Solar System. Thus, the highly-relativistic and dynamical strong-field regime of Einstein's equations, i.e. the regime characterized by velocities comparable to the speed of light, by large, dynamical curvatures, has not yet been tested and could in principle highlight large deviations from GR's predictions (see e.g. Refs. [4][5][6][7] and [8][9][10][11] for two theories of gravity that are indistinguishable from GR with current Solar System and binary pulsar observations, but deviate from it in strong-gravity systems). Tests of this highly-relativistic strong-field regime will probably have to wait for the detection of gravitational waves (GWs) from merging neutron stars (NSs) or black holes [5][6][7][12][13][14]. To understand how much will be learnt about gravitation from these observations, it is important to clarify the corresponding predictions from GR and alternative theories [14].
Another reason why GR cannot be the final word for gravitation in Nature is its intrinsic incompatibility with quantum mechanics and the presence of mathematical pathologies, like singularities in gravitational collapse. One may conjecture that the modifications to GR induced by any approach addressing these issues will be suppressed by a high-energy scale, e.g. the Planck mass M P ∼ 10 19 GeV where GR fails as an effective field the-ory. This would mean that no information about quantum gravity could be obtained by any current or foreseen experiment. However, such arguments hold only on the very thin thread of dimensional analysis. In fact, they are known to fail for existing candidates of quantum gravity, such as models with large extra-dimensions, where the scale of suppression can be much lower [15], or in Hořava gravity [16], where GR modifications can enter at any scale.
Finally, cosmological observations [17,18] reveal a model of the Universe (the ΛCDM model) that is not completely satisfactory from a theoretical standpoint. The most famous puzzle is the conflict between the value of the energy density of the source behind the current cosmological acceleration and the order of magnitude that would be expected based on naïve dimensional arguments. Any explanation addressing this cosmological-constant problem typically requires a modification of gravitation at cosmological (and possibly smaller) scales [19,20].

A. Lorentz-Violating Gravity
Lorentz symmetry has been put to the test in a variety of circumstances. Particle physics experiments have stringently done so in the matter sector [21][22][23][24]. A model-independent formalism, the standard model extension (SME) [25][26][27], has been devised to translate a wide variety of observations into tests of Lorentz invariance. This is very efficient for bounds on violations of Lorentz symmetry primarily in the matter sector [22] or in the sector coupling matter to gravity [28].
Given the previous constraints, one could question the possibility of having observable effects in gravitation from Lorentz violation. Indeed, one could argue that any degree of gravitational Lorentz violation should percolate into particle physics. Therefore, stringent constraints from the matter sector would require that Lorentz violation in gravity also be very small (and thus undetectable with current observations). This argument, however, is not watertight. Different mechanisms have been put forward to justify the possibility of Lorentz violation in gravity to a degree that is orders of magnitude larger than any violation in the matter sector (see e.g. Ref. [29]).
A first possibility would be to assume that this happens due to the finely tuned (small) value of the operators that violate Lorentz invariance in the matter sector, as compared to those in the gravity sector. More interestingly, Lorentz invariance in the matter sector could be an emergent feature at low energies [30], either due to a renormalization group phenomenon [31,32] or to it being an accidental symmetry [33]. Finally, it has been pointed out recently that two sectors with different degrees of Lorentz violation can easily coexist provided that the interaction between them is suppressed by a high energy-scale [34]. This could be the case for the matter and gravity sectors, each with a very different scale of Lorentz violation.
Thus, it is reasonable to seek for independent tests of Lorentz invariance in gravity. These are neither as developed nor as stringent as for the matter sector. Gravitational Lorentz invariance bounds in the context of the SME can be found in Refs. [35][36][37]. More relevant for this work are the bounds related to models where Lorentz invariance is broken by the presence of a preferred timelike vector. Those include constraints placed with Solar System [3,38,39] and cosmological observations [39][40][41][42], but these are rather weak. In the Solar System, observations can only probe certain aspects of gravitational Lorentz violation, namely leading-order post-Newtonian (PN) effects, i.e. leading in the expansion in the ratio of the orbital velocity to the speed of light. In cosmology, only the linear perturbative regime over an expanding background has been studied. Binary pulsar observations have also been used to place constraints on gravitational Lorentz invariance [39,[43][44][45], but these concentrate only on preferred-frame corrections to the conservative orbital dynamics, neglecting dissipative effects.
In this paper, we will focus on theories that modify GR by assuming the existence of a preferred time direction (frame) at every spacetime point. This directly implies the violation of boost symmetry and therefore of Lorentz invariance. The preferred time direction will be described by a timelike unit vector field U µ (the AEther field), whose dynamics may be described by two different theories: Einstein-AEther [46] and khronometric theory [47]. As we will argue in Sec. II, these are generic theories that may arise from more fundamental theories at low energies.
Einstein-AEther theory is a Lorentz-violating, metric theory of gravity where the AEther field is generic. The theory is characterized by the AEther's couplings to gravity. A direct coupling of the AEther field to matter is absent, because that would have consequences that are not observed experimentally, e.g. a fifth force would arise in the matter interactions and the weak equivalence principle (i.e. the universality of free fall for weakly gravitating bodies) would not hold. Einstein-AEther theory can be thought of as representing the low-energy description of some high-energy unknown dynamics, i.e. it can be understood as an effective field theory [48].
Einstein-AEther theory passes all theoretical and phenomenological constraints in a certain region of the parameter space of the AEther's couplings. Of particular importance are the bounds related to Solar System observations, which force the theory to depend on only two combinations of coupling constants, c + and c − , up to corrections of O(10 −4 ) [3,49]. These two combinations, however, are currently weakly constrained, mostly by requiring the linear stability of Minkowski space (i.e. absence of tachyonic instabilities) and the absence of gravitational Cherenkov radiation [50] (i.e. requiring the speeds of the gravitational modes be larger than the speed of light, so that photons and high-energy particles do not lose energy in a Cherenkov-like process). An order of magnitude constraint (of about O(10 −2 )) was also obtained using the orbital decay rate of binary pulsars [51], but this analysis is only valid in the small coupling region (c i 1). Khronometric theory is defined through the same action as that of Einstein-AEther theory, but with the additional requirement that the AEther field be hypersurfaceorthogonal. This choice reduces the dynamical degrees of freedom of the theory, and more importantly it forces the theory to exactly coincide with the low-energy limit of Hořava gravity [16]. Hořava gravity is a proposal for quantum gravity where Lorentz invariance is broken by the existence of a preferred foliation of spacetime into spacelike hypersurfaces, and which provides a powercounting renormalizable completion of GR in the ultraviolet regime. In this work, we will focus on the version of this theory introduced in Ref. [47], which has Minkowski spacetime as a ground-state, and which reduces exactly to khronometric theory at low energies.
As in Einstein-AEther theory, any viable khronometric theory is characterized by two coupling constants, β and λ, after imposing constraints from Solar System observations. As in Einstein-AEther theory, β and λ are weakly constrained by requiring linear stability of Minkowski space and the absence of gravitational Cherenkov radiation [52,53]. However, unlike in Einstein-AEther theory [39][40][41]54], cosmological observations can place stringent bounds on the couplings of khronometric theory [roughly |β, λ| O(10 −1 )]. In particular, strong constraints can be obtained by requiring that big bang nucleosynthesis (BBN) in khronometric theory produces element abundances in agreement with observations [42]. Cosmological constraints are instead much weaker in Einstein-AEther theory because, unlike in khronometric theory, the Newtonian constant regulating the cosmological evolution coincides with the locally measured constant to within O(10 −4 ), once Solar System constraints are imposed.

B. Executive Summary
In this paper, we explain in detail and extend the analysis of Ref. [55] to place very stringent constraints on Lorentz-violating gravity, focusing on Einstein-AEther and khronometric theory. These constraints are obtained by comparing the evolution of binary systems in these theories to binary pulsar observations. Let us first consider modifications to the dissipative sector, which regulates how fast binary systems lose energy and shrink, forcing the orbital period to decay. A generic property of Einstein-AEther and khronometric theory is the excitation of propagating modes that are absent in GR, and which carry energy away from the system at dipole order [38,51,56]. Since dipolar radiation is generally stronger than GR's quadrupole radiation, the binary's separation decreases much faster in Lorentz-violating gravity, leading to a strong modi-fication to the predicted evolution of the orbital period. Furthermore, these extra modes and the modification to the propagation speed of gravitons affect the quadrupolar emission, which is important for systems where dipolar radiation is suppressed. Since binary pulsar observations of the orbital period's decay rate agree with GR's predictions within the observational uncertainties, they allow for stringent constraints to be placed on Lorentz-violating gravity.
The modified orbital decay rate for NS binaries in Lorentz-violating theories, and more in general the motion of these systems, strongly depend on the sensitivities of the stars [51] 1 . These quantities measure how much the binding energy of an isolated star changes with its motion relative to the preferred frame (i.e. relative to the AEther). In order to calculate the sensitivities, we therefore need to find solutions describing NSs moving with respect to the AEther, as suggested in [57]. More specifically, we will show that in order to extract the sensitivities, such solutions are only needed at linear order in a perturbative expansion in the velocity v, i.e. without loss of generality, we can restrict attention to solutions of stationary, non-spinning NS moving slowly with respect to the AEther. We show that two seemingly different definitions of the sensitivities, one applicable in the weakfield regime and another in the strong-field regime, lead to the same result.
To find these slowly-moving solutions, we first write down the most generic ansatz for the metric and AEther field at O(v), ensuring compatibility with the symmetries of the problem (stationarity and rotational invariance around the motion's direction). With this ansatz, we write down the field equations at zeroth-and firstorder in v. The resulting system of partial differential equations is then expanded in tensor spherical harmonics, leading to an ordinary differential system that we solve numerically.
The equations must be solved twice (in the interior and in the exterior of the star) and then matched. The matching ensures that the potentials (defined in terms of the AEther and metric), which enter the equations through derivatives, are continuous everywhere in the spacetime, regular at the center of the star, and produce an asymptotically flat geometry at spatial infinity. The interior solution depends on the equation of state (EoS) for the NS matter. We here investigate a variety of EoSs that are thought to represent realistic NS configurations. We restrict attention to non-rotating, cold (and thus old) NSs, as these are appropriate simplifications for binary pulsar studies.
With this at hand, we use the numerical solutions that we obtained to compute the sensitivities for isolated NSs for the first time in the theories under consideration. We then construct fitting functions that allow us to analytically model the sensitivities as a function of the compactness of the star and the coupling constants of the theory.
Once the sensitivities have been found, we calculate the energy carried away by all propagating degrees of freedom in Einstein-AEther and khronometric theory. We show that dipolar radiation is produced and generically constitutes the dominant modification to the GR predictions [3]. Given the energy flux carried away by the propagating modes, we use conservation of energy and a balance law to compute the rate of change of a binary's binding energy, and from this, the orbital decay rate.
We compare our predictions for the orbital decay rate to observations of binary pulsars PSR J1141-6545 [58], PSR J0348+0432 [59], and PSR J0737-3039 [60]. The first two are pulsars on a 0.17-eccentricity, 4.74-hour orbit and on a O(10 −6 )-eccentricity, 2.46-hour orbit respectively, around a white dwarf companion. The third is the relativistic double pulsar binary, on a 0.088-eccentricity and 2.45-hour orbit. These comparisons allow us to place constraints on the coupling constants of the theory.
Another way to place constraints on Lorentz-violating theories is to consider modifications to the conservative sector , controlled by the Hamiltonian, which for example affects the orbital shape and precession rate. Lorentz-violating corrections to the Hamiltonian induce precession of the spin and orbital angular momentum vectors. Since such non-GR precession is not found in pulsar observations, one can then place constraints on Lorentz-violation. The constraints are cast in a model-independent language by considering strongfield generalizations of the parametrized post-Newtonian (PPN) Hamiltonian. For example, observations of PSR J1738+0333 [61] can be used to constrain the strongfield PPN parameters associated with preferred-frame effects. We here calculate these parameters for Einstein-AEther and khronometric theory, and then use PSR J1738+0333 [61] to place constraints on the couplings.
Combining all of these constraints, we obtain the allowed coupling parameter space shown in Fig. 1 (Einstein-AEther theory in the left panel and khronometric theory in the right panel). The colored regions are those allowed after requiring stability and absence of gravitational Cherenkov radiation [50,52,53] (light blue), BBN constraints [39,40,42,54] (dark orange) and binary pulsar constraints (dark purple). The red dashed line corresponds to the values of the coupling constants for which the orbital decay rate equals the GR prediction, assuming the sensitivities vanish and working at leadingorder in a weak-field expansion [49]. Observe that the new constraints obtained here are much stronger than all other constraints. The constraint found in this paper is consistent with the order of magnitude estimate of Foster's [51] in the small coupling region, but the constraint on c − found here is slightly stronger.
Binary pulsar constraints lead to regions of viable coupling parameter space. This is because in deriving these constraints one has to allow for different values of the coupling constants (within the Solar System constraints) and the sensitivities (because of the different possible EoSs), and account for the observational error in the orbital decay rate and the orbital period, as well as the error in the inferred masses of the binary. The particular shape of these viable regions is a result of the combination of the constraints associated with different binary pulsars. For example, dipole radiation is suppressed for double pulsar binaries relative to pulsar-white dwarf systems, because dipole radiation is proportional to the difference of the sensitivities and NSs have similar sensitivities. Therefore, for double pulsar systems quadrupolar radiation becomes comparable to dipole radiation, which leads to a different shape for the allowed region of coupling constants.

C. Layout and Conventions
The remainder of this paper presents the details of the calculations described above and is organized as follows: Section II presents the basics of Einstein-AEther and khronometric theory, including the action, field equations and current constraints from low-energy phenomenology. Section III describes how Lorentz-symmetry breaking leads to violations of the strong equivalence principle, how this is encoded in the sensitivities, and how it affects the motion of compact objects. Section IV and Section V are devoted to the construction of slowly-moving NS solutions to first order in velocity, in both Einstein-AEther and khronometric theory. Section VI presents results derived by numerically solving the modified field equations, focusing on the sensitivities. Section VII presents constraints on Einstein-AEther and khronometric theory based on binary pulsar observations. Section VIII concludes and points to future research. Appendices A, B and C present further mathematical details.
Our conventions are as follows. We use natural units where c = 1 = , where is the (reduced) Planck constant and c is the speed of light. We restore powers of c when presenting PN expressions. We use Greek letters in index lists to denote spacetime components of tensors, while Latin letters in the middle of the alphabet (i, j, k, . . .) denote purely spatial components. We also employ the metric signature (+, −, −, −). Indices in spatial vectors are raised and lowered with the Kronecker delta (except in Sec. II B, where the 3-metric γ ij is used to raise and lower indices).
We use several gravitational constants, different masses, and different velocities throughout the paper, which we list here for convenience. Regarding gravitational constants: • G AE is the bare gravitational constant appearing in the action [Eq. (1)]; • G N is the "Newtonian" gravitational constant measured locally by Cavendish-type experiments [Eq.  [58], PSR J0348+0432 [59], PSR J0737-3039 [60] and PSR J1738+0333 [61]. The areas outside the (allowed) shaded regions are ruled out by stability/Cherenkov considerations (light blue), BBN (dark orange) and the combined binary pulsar constraints (dark purple). The red dotted line corresponds to the values of the coupling constants required for the orbital decay rate to agree with the GR prediction in the zero-sensitivity/weak-field limit. Observe that the new constraints are much more stringent than all others.
• G C is the "cosmological" gravitational constant that appears in the Friedmann equations [Eq. (28)]; • G is the "effective" gravitational constant in a binary system [Eq. (88)].

Regarding the masses:
•m A is the gravitational mass of the A-th body in a point-particle approximation [Eq. (31)]; • m A is the "active" gravitational mass of the A-th body in a point-particle approximation [Eq. (87)]; • M tot is the total gravitational mass of a star, which includes the gravitational, AEther and baryonic contributions [Eq. (57)]; this mass generalizesm A to regimes where the point-particle approximation does not hold; • M obs is the mass measured by Keplerian experiments, which turns out to coincide with M tot ; • M * ≡ G N M tot = G N M obs is the length scale associated with the total mass M tot = M obs ; • M (r) is a function with dimension of length, defined by Eq. (137) and approaching M * as r → +∞; • m ≡ m 1 + m 2 is the total active mass of a binary system in the point-particle approximation; • µ ≡ m 1 m 2 /m is the active reduced mass of a binary system in the point-particle approximation.
Regarding the velocities, we use • v i or v i A are both the 3-velocity of an object relative to the AEther field; is the relative velocity of the two bodies in a binary; • V i CM is the center-of-mass velocity of the binary relative to the AEther.

II. MODIFIED GRAVITY THEORIES
In this section, we define the theories we focus on. We begin with a description of Einstein-AEther theory and follow with khronometric theory (the low-energy limit of Hořava gravity). In both cases, we first introduce the action that defines the theory and then describe its current experimental constraints.

A. Einstein-AEther Theory
Einstein-AEther theory describes gravity by means of a metric g αβ and a unit-norm timelike dynamical vector field U α (the "AEther field"). The latter locally defines a preferred time direction, which breaks boost-and therefore Lorentz-invariance. Up to total divergences, the most generic covariant action that depends only on the fields and their first derivatives, and is quadratic in the latter, is [39,46,62] where g and R are the metric determinant and the Ricci scalar, U α satisfies the unit constraint U α U α = 1, (c 1 , c 2 , c 3 and c 4 being dimensionless coupling constants), and the "bare" gravitational constant G AE is related to the "Newtonian" constant measured with Cavendishtype experiments through [54] .
In order to enforce the weak equivalence principle 2 and the absence of a fifth force in matter interactions, one has to (minimally) couple the matter fields to the metric but not to the AEther. Collectively denoting the matter fields by ψ, one can then write the total action as The matter action is diffeomorphism invariant, which implies the covariant conservation of the matter stressenergy tensor T µν mat , i.e. ∇ µ T µν mat = 0 [63], and as a consequence test particles follow geodesics [64], thus satisfying the weak equivalence principle. However, as we will show in this paper, in Einstein-AEther (and khronometric) theory, strongly gravitating bodies follow trajectories that depend on the ratio between the body's binding energy and its total mass. This effect (known as the Nördvedt effect [65,66]) generically appears when one introduces extra gravitational degrees of freedom coupled non-minimally to the metric (e.g. it is present also in scalar tensor theories [3]), and amounts to a violation of the strong equivalence principle. Its physical origin lies in the fact that the matter fields couple to the extra gravitational fields through the metric, and this coupling becomes important in strong-gravity regimes.
The absence of a direct coupling of the matter fields to U α in Eq. (4) is also enforced by tests of Lorentz invariance in particle physics. Indeed, if the matter fields were directly coupled to U α , this would generically produce Lorentz-violating effects in the standard model of particle physics. The bounds on these effects are very tight [22,23], which means that for gravitational experiments it is reasonable to put all the couplings between AEther and matter to zero 3 . Naturalness arguments seem to be at odds with this choice, since they suggest that Lorentz violation may percolate from the gravitational sector into the standard model, thus producing big (and experimentally forbidden) effects. However, this is not necessarily the case. As already mentioned, different mechanisms have been put forth leading to Lorentz-violating effects in gravity that can be much larger than those in particle physics, see e.g. Ref. [29] for a review, Ref. [32] for recent developments, and Sec. I A for further details.
Variation of the action given in Eq. (4) with respect to the metric and AEther field [imposing the unit constraint U α U α = 1 directly or by means of a Lagrange multiplier (U µ U µ − 1) in the action] provides a set of modified Einstein equations and the AEther equations Here, G αβ = R αβ − R g αβ /2 is the Einstein tensor, is the AEther stress-energy tensor, where Finally, the matter stress-energy tensor is defined as usual by Experimental constraints and stability requirements greatly reduce the viable parameter space for the four couplings c 1 , c 2 , c 3 and c 4 . At 1PN order 4 , the dynamics of Einstein-AEther theory matches that of GR, with the exception of two "preferred-frame" parameters α 1 and α 2 , which are exactly zero in GR but not in Einstein-AEther theory [49]. In the weak-field, Solar System observations constrain these parameters to very small values |α 1 | 10 −4 and |α 2 | 10 −7 [3], using lunar laser ranging and solar alignment with the ecliptic. The strongfield counterparts of these parameters (denoted here with an overhead hat) are constrained through binary and isolated pulsar observations also to very small values, |α 1 | 10 −5 and |α 2 | 10 −9 [44,45]. We will discuss the Einstein-AEther form of (α 1 , α 2 ) in Sec. III A 2 and the relation to their strong-field generalizations in Sec. VII C.
Since current constraint force (α 1 , α 2 ) to be small dimensionless numbers, we can expand the theory in these quantities and shrink the parameter space to just two parameters. In Einstein-AEther theory, this implies [39,49] to leading order in α 1 and α 2 . The only free parameters in Einstein-AEther theory are then c ± ≡ c 1 ± c 3 . Further constraints on c ± come from requiring that perturbations about a Minkowski background are stable and have positive energy [68], and that matter does not emit gravitational Cherenkov radiation [50]. Such radiation would be emitted if the speed of the AEther's and metric's propagating modes were smaller than the speed of light. These requirements result in the bounds to leading order in α 1 and α 2 . In addition to these, Foster [51] estimated that roughly c ± < O(10 −2 ), assuming a small-coupling approximation (c i 1) and using order-of-magnitude information about the orbital decay of binary pulsars.
For arbitrary c i , the dependence of the gravitational constant G C appearing in the Friedmann equation on the coupling constants is different than that of the gravitational constant G N measurable via Cavendish-type experiments. For generic values of the couplings, therefore, G N and G C differ. This difference leads to deviations from the metal abundances predicted by Big Bang Nucleosynthesis (BBN) [54] and to a modified growth of cosmological perturbations that affects the Cosmic Microwave Background (CMB) and the distribution of matter at large scales [42]. However, once Eqs. (10) and (11) are satisfied, G N ≈ G C up to terms of O(10 −4 ), and thus, cosmological observations do not significantly reduce the viable region defined by Eqs. (12) and (13) [39,40]. Similarly, although black hole solutions in Einstein-AEther theory differ from GR [53,69,70], current electromagnetic observations of black hole candidates are still not accurate enough to provide constraints competitive with current bounds (see e.g. Ref. [71]).

B. Khronometric Theory and Hořava Gravity
As mentioned in the previous section, Einstein-AEther theory breaks Lorentz invariance by locally specifying a preferred time direction through an AEther vector field U α . This vector field may appear from the existence of a global preferred time. In this case, it will be orthogonal to the hypersurfaces of constant preferred time, i.e. proportional to the gradient of a foliation-defining scalar field T : which satisfies normalization condition U α U α = 1. The AEther vector field is timelike as a consequence of the scalar defining a preferred time coordinate. The scalar T is often referred to as the "khronon" 5 , and a metric theory in which Lorentz invariance is broken globally by such a scalar is called "khronometric theory." The most generic action (quadratic in derivatives) for khronometric theory is given by Eq. (1) with the definition in Eq. (14) [52,74]. The condition in Eq. (14) allows one to express one of the AEther terms in the action in terms of the other ones, and thus there are only three free independent AEther terms. In particular, we can absorb the c 1 term into the other three terms, by multiplying the second, third and forth term respectively by the new couplings Another form for the action can be derived by choosing the time coordinate to coincide with T = constant hypersurfaces. In this gauge, Eq. (14) becomes where N = (g T T ) −1/2 is the lapse. The action in Eq. (1) with the definition in Eq. (14) then becomes [75] where K ij , (3) R and h ij are respectively the extrinsic curvature, the 3-Ricci curvature and the 3-metric of the T = constant hypersurfaces, and where we have defined the "acceleration" of the AEther flow 6 , a i ≡ ∂ i ln N . Latin indices are manipulated with the 3-metric of the T = constant hypersurfaces. An additional motivation for khronometric theory is that it coincides with the low-energy limit of a theory of gravity that has remarkable properties at high energies. 5 From the Greek χρóνoς (khronos) -time. The romanization was chosen to avoid confusion with previous uses of the prefix 'chrono', e.g. [72,73], which are not related to the theories under study here. 6 This vector is related to the acceleration of AEther congruence by U µ ∇µU ν = a i δ ν i This theory is known as Hořava gravity [16], and achieves power-counting renormalizability by adding higher-order derivative terms to the khronometric action [16,47]: where the L 2 term corresponds to the gravitational Lagrangian of khronometric theory [cf. Eq. (17)], while L 4 and L 6 are terms that are suppressed by a suitable mass scale M . These terms are respectively of fourth and sixth order in the spatial derivatives, but contain no derivatives with respect to T .
Let us stress that the experimentally viable range for the mass scale M is rather broad. On the one hand, M is bound from above (M 10 16 GeV) to allow the theory to remain perturbative at all scales [75][76][77], so that the power-counting renormalizability arguments of Ref. [16] can be applied. On the other hand, M is also bound from below from tests of Lorentz invariance. That lower bound, however, depends on the details of the percolation of Lorentz violations from the gravity to the matter sector, which is not yet completely understood [29,34,52,[78][79][80]. Taking instead into account only the constraints from Lorentz violation in the gravity sector, one gets the relatively weak bound M 10 −2 eV [3,52]. In this paper, we will only consider khronometric theory, i.e. we will focus on the low-energy limit of Hořava gravity and neglect the higher-order terms L 4 and L 6 . Neglecting these terms is, however, an excellent approximation as far as astrophysical studies are concerned. In fact, on purely dimensional grounds, the error introduced by the L 4 and L 6 terms on a NS solution of mass . Taking the lowest conceivable value for M , i.e M ∼ 10 −2 eV, one gets an error of roughly 10 −16 (M /M tot ) 2 or smaller, when neglecting L 4 and L 6 in binary pulsar systems, which is clearly negligible.
The field equations for khronometric theory can be obtained by varying the action in Eq. (1), with the definition in Eq. (14) replaced in it before the variation. Varying with respect to g αβ and T one finds where E αβ and AE α are defined as in Einstein-AEther theory [cf. Eqs. (5)-(6)]. As can be seen, any hypersurface-orthogonal solution to Einstein-AEther theory is also a solution to khronometric theory, but the converse is not necessarily true. It can be shown that spherically symmetric, static, asymptotically flat solutions are indeed the same in the two theories [52,74,81,82]. However, the same is not true in more general cases; for example, the sets of slowly rotating black hole solutions of the two theories do not over-lap [70,82,83]. Note that because of the Bianchi identity, Eq. (20) is actually implied by Eq. (19) and by the equations of motion of matter (which imply ∇ ν T µν mat = 0), i.e. the only independent equations of khronometric theory are actually the modified Einstein equations and the equations of motion of matter [74].
Finally, let us discuss the current experimental bounds on the coupling constants of khronometric theory. Like in the case of Einstein-AEther theory, Solar System tests require α 1 10 −4 and α 2 10 −7 , and similar stringent bounds are imposed on (α 1 ,α 2 ) by pulsar observations (see Sec. VII for a detailed discussion). The fact that for khronometric theory [38] yields two possible ways to enforce the previous bounds. One way is to choose α 1 to be small enough so that both bounds are simultaneously satisfied. This can be enforced by only one condition of the form which leaves the coupling parameters λ and β unconstrained. Another way is to saturate both bounds for α 1 and α 2 , which leads to a smaller parameter space, i.e. one can require For our purposes, Eqs. (22) and (23) are equivalent, since our bounds will certainly not be sensitive to differences below O(10 −4 ). Thus, the choice given by Eqs. (23) and (24) is more restrictive, and without loss of generality we can neglect it [as it merely selects a very small subset of the parameter space allowed by Eq. (22) at leading order in α 1 and α 2 ]. We stress, however, that the formalism presented in this paper is perfectly suited to the constraints of Eqs. (23) and (24) as well, and we do investigate this choice in more detail in a later section (see e.g. Fig 6).
Requiring that the propagating modes are stable and have positive energy in flat space, and not allowing gravitational Cherenkov emission by matter, one finds the conditions [50,52,53] Unlike in Einstein-AEther theory, cosmological constraints are rather strong in khronometric theory. In fact, one has to leading order in α 1 . The agreement between observations and the metal abundances predicted by BBN requires |G C /G N − 1| 1/8. Once this bound is combined with the stability/Cherenkov constraints of Eqs. (25)- (27), it selects a rather limited region of the (λ, β) plane (cf. the orange region in Fig. 1). Finally, when Eqs. (23) and (24) are imposed, the ratio G C /G N is almost one, and BBN observations do not impose a new constraint, as in the case in Einstein-AEther theory. Further constraints may come from observations of the CMB and large scale structure of the Universe. As for BBN, those constraints are expected to be efficient except for the case β = −λ. This has been demonstrated explicitly in Ref. [42], which extended the action in Eq. (1) by including a dynamical dark energy component. However, because dynamical dark-energy models affect the evolution of cosmological perturbations, we will not consider the bounds of Ref. [42] here. Also, like in the Einstein-AEther case, black hole solutions in khronometric theory differ from the GR ones [53,69,70,[81][82][83], and may in principle allow to test the theory in the future.
Finally, Ref. [84] recently found a relation between khronometric theory and Einstein-AEther theory at the level of the solutions to the field equations. More precisely, Ref. [84] showed that the solutions of khronometric theory can be obtained from the Einstein-AEther theory solutions in the limit c 1 − c 3 → ∞ (with c 1 + c 3 , c 2 and c 1 + c 4 held fixed). While we have not used this correspondence in this paper to derive the solutions to the khronometric theory field equations, we have used it to test the correctness of our results.

III. VIOLATIONS OF THE STRONG EQUIVALENCE PRINCIPLE AND THE SENSITIVITY PARAMETERS
The prediction of observables in physical theories requires knowledge of the solution to the field equations for the system under consideration. However, when studying compact binary systems, e.g. of pulsars, exact solutions are not available in modified gravity theories, or for that matter, in GR. One is then forced to rely on approximations, such as the PN scheme, where the system's dynamics is expanded in the ratio between the characteristic velocity of the system and the speed of light. In this scheme, compact objects are modeled effectively with point particles. Clearly, this approximation breaks down in the strong-field regions near the compact objects, but this will be irrelevant in this paper; the PN scheme is ideal to study the orbital decay rate of binary pulsars.
As already mentioned in Sec. II, the strong equivalence principle is violated in Einstein-AEther and khronometric theory because of the Nördvedt effect [65,66], i.e. because an effective coupling between matter and the AEther appears in the strong-gravity regime. To model this effect, the point-particle action used to describe compact objects must depend on "AEther charges" or "sensitivity parameters" [85,86] coupling the particles to the AEther. In this section, we will explain how this comes about, what these parameters mean physically, and how they can be computed from the PN metric tensor. We then proceed with a strong-field definition of the sensitivities and conclude with a description of how the sensitivities percolate into observables, with a particular focus on GW fluxes.

Point-Particle Action for Compact Objects
In PN theory, one models the motion of compact objects sufficiently far from each other effectively through point particles. The way these particles couple to the different fields of the theory encapsulates finite-size or strong-field effects. For weakly gravitating objects (i.e. ones with negligible binding energy compared to their total gravitational mass), the point particles effectively describing them in the PN scheme can only couple to the metric [c.f. Eq. (4) and discussion in Sec. II]. However, even though the matter fields do not couple directly to the AEther, the metric does so non-minimally. Therefore, because matter couples to the metric (although weakly) through gravity, an effective coupling appears between the AEther and the matter fields when gravity is strong (i.e. when the metric perturbations produced on the background by the presence of the object are large, as in the case of NSs). The existence of this effective coupling has long been known in the context of scalar-tensor theories (where it is called the Nördvedt effect [65,66]), and introduces deviations away from geodesic motion for strongly gravitating objects, thus violating the strong equivalence principle. The sensitivities will parametrize this effective coupling.
In order to clarify how violations of the strong equivalence principle come about, let us briefly review the physical meaning of the sensitivities in scalar-tensor theories [85]. In those theories, the gravitational interaction is mediated by the metric and by a gravitational scalar field φ, coupled minimally to the matter but non-minimally to the metric. The gravitational constant, and therefore the binding energy of a compact object, depend on the local value of the scalar field φ. Because the binding energy contributes to the gravitational mass, one has to include an explicit dependence of the gravitational massm on the scalar field φ when using a point-particle model to describe compact objects, i.e. S pp A = − dτ AmA (φ), with A the object's label. From this, it is clear that the strong equivalence principle will be violated, because the scalar field depends on position, the mass becomes position-dependent, and the variation of S pp A does not yield the geodesic equation.
We can parametrize deviations from geodesic motion by exploiting the fact that in the PN regime all motion is slow compared to the speed of light, and one expects that changes in scalar field are also generally small and slow. At leading PN order, it is sufficient to consider a leading-order Taylor expansion of the massm A (φ), and parametrize deviations from geodesic motion through the sensitivity [85] where φ 0 is the constant value of the scalar field far from the object. The partial derivative in Eq. (29) is to be taken along a reversible transformation of the object [85]. This is because under a change of the local scalar field (and therefore of the local value of the gravitational constant), the binding energy also changes, resulting in the body expanding or shrinking. The extra kinetic energy gained by the volume elements of the body is then assumed to be transformed into potential energy without production of heat (i.e. without production of entropy).
In other words, the body is thought to gradually adapt its structure to the change in the scalar field. This is a good approximation if the scalar field changes slowly enough, which is the case in the PN regime. In Einstein-AEther and khronometric theory the situation is similar. In general, because of the effective coupling between the AEther vector field and matter in the strong field regime, the compact object's structure, its binding energy and its gravitational mass will be a function of the motion relative to the AEther. This is exactly why Lorentz symmetry is violated in these theories, i.e. if Lorentz symmetry is broken, then the motion of any compact body with respect to the AEther should be experimentally detectable. Drawing inspiration from scalartensor theories, let us model systems of strongly gravitating objects with a point-particle action of the form [51] where A labels the object, γ A ≡ U µ u µ A (with u µ A denoting the particle's four-velocity and U µ , as usual, the AEther vector field) is the Lorentz factor of each particle relative to the AEther, and dτ A is the proper length along the particle's trajectory. Because PN theory is a perturbative expansion in the system's characteristic velocity, and because γ A ∼ 1 corresponds to an object moving slowly relative to the AEther, we Taylor expand the action of Eq. (30) as wherem A ≡m A (1) is a constant and we have defined the sensitivity parameters σ A and σ A via in analogy with scalar tensor theories. In what follows, it will sometimes be convenient to use the rescaled sensitivity parameter Clearly, for weakly gravitating objects, σ A ≈ s A ≈ 0 and σ A ≈ 0. Let us stress that the assumption that is the velocity of the A-th compact object relative to the AEther, is implicit in Eq. (31). The PN expansion, on the other hand, is in terms of the ratio v 12 /c, where v 12 is the the binary's relative velocity (i.e. the binary system's characteristic speed). In our analysis we will consider both v 12 /c 1 and v A /c 1. This follows from the combination of these two facts: (i) cosmologically, the AEther must be almost aligned with the CMB frame (otherwise there would be strong and non-viable effects on the cosmological evolution; see also Ref. [87] for a dynamical study of the alignment); (ii) the peculiar velocity of our Galaxy relative to the CMB is ∼ 10 −3 c, thus the center of mass of binary systems moves slowly relative to the AEther.
For an action of the form of Eq. (30), there is a subtlety in the derivation of the field equations. In fact, the field equations presented in Sec. II assume that the AEther does not couple directly to matter [c.f. Eq. (4)], but as mentioned above, this is not the case in a pointparticle model aiming at effectively describing systems of compact objects. Because of the presence of this effective coupling, the field equations get modified by terms proportional to the sensitivities. In particular, for Einstein-AEther theory, the equations for the AEther [Eq. (6)] are modified to include a source term where x i A is the worldline of the A-th point particle. Also, the Einstein equations [Eq. (5)] remain valid, but the matter stress-energy tensor picks up additional terms with respect to the definition Eq. (9). Mathematically, these extra terms appear because the modified AEther equation [Eq. (35)] changes the value of the Lagrange multiplier enforcing the unit constraint. Doing the calculation explicitly, it turns out that the matter stress energy tensor becomes with the stress-tensor for the AEther field is still given by Eq. (7).
Similarly, in khronometric theory, Eqs. (19) and (20) are still valid if AE µ is replaced byAE µ as defined by Eq. (35), and if the point-particle stress-energy tensor is given by This implies, in particular, that hypersurface-orthogonal solutions of Einstein-AEther theory are also solutions of khronometric theory, although the converse is not always true.

PN Solution to the Modified Field Equations
In this section, we review the solution to the field equations for a binary system of compact objects. In particular, we will solve the equations within a PN approximation (v 12 /c 1), for binaries moving slowly relative to the AEther (v A /c 1) and described by the effective action of Eqs. (30) and (31). Note that the solution will intrinsically depend on the parameterm A and the sensitivities σ A and σ A , but σ A will enter at higher order than σ A in the perturbative expansion [c.f. Eq. (31)]. The solution presented here for the Einstein-AEther theory was derived in Ref. [51], while the solution for khronometric theory is novel (but see Ref. [38] for a restricted solution with vanishing couplings between the point particles and the AEther). The derivations are standard and we will simply present the results here for completeness.
Let us work in the standard PN gauge and PN coordinates (t , x , y , z ) [88]. In both Einstein-AEther and khronometric theory, the 1PN metric takes the form wherem A is the mass of the A-th point particle [cf. Eq. (31)], v i A its velocity, r 12 the binary's separation, r A the distance from the A-th particle to the field point, n i A the unit vector associated with r A and the symbol 1 ↔ 2 means that one is to add all terms on the righthand side of the equality with the exchange 1 ↔ 2. We have restored here factors of 1/c to make the PN order counting clearer. A term proportional to (m A /r A ) N and where the weak-field PPN parameters are and we have defined In khronometric theory, the values of the constants B ± where now the weak-field PPN parameters are The AEther field takes the form where for Einstein-AEther while for khronometric theory

The Sensitivities from the PN Metric
The PN solution presented in the previous section provides a way to calculate the sensitivities in practice. Consider in particular the solution for the metric tensor, but specialize it to a single object (rather than a binary system) moving relative to the AEther. Mathematically, this is achieved by setting one of the masses (e.g.m 2 ) to zero (or the separation r 12 to infinity), but physically it corresponds to placing oneself at a distance r 1 from object 1, such that R * r 1 r 12 (with R * being the radius of the object). The existence of this "buffer" region is possible because the PN scheme requires r 12 R * , in order to model the stars as point particles.
The advantage of considering an (effectively) singleobject system is that exact strong-field solutions that do not rely on the point-particle approximation might exist for such a system. Recall that this is not the case for binary systems, for which only approximate PN solutions exist, even in GR. If such a solution for a single object moving relative to the AEther can be found, the metric near spatial infinity will be given by the PN solution of Eqs. (38)- (40), with one of the masses set to zero. Because this solution depends on σ A , one can in principle read off this quantity from the behavior of the exact solution near spatial infinity.
What components of the metric near spatial infinity should we use to extract σ A ? One cannot use the spatial part of the metric g i j , because this does not depend on σ A . We could use the (t , t ) piece of the metric, but σ A enters multiplied by a term of O(1/c 4 ). One can find σ A more easily by using the gravitomagnetic sector of the metric [Eqs. (39), (41) and (47)], and in particular the O(v) terms. By using Eq. (41) and Eq. (47), one easily finds in Einstein-AEther theory and in khronometric theory.
The coefficients B ± A that are needed to compute the sensitivities are to be extracted from the exact solution for an isolated NS moving relative to the AEther, as measured by an observer near spatial infinity. The exterior solution will depend on the interior solution, since both must be properly matched, as we will see in Sec. IV. Therefore, the coefficients B ± A will depend on the strong-field behavior of the solution in the stellar interior. Note also that these coefficients appear in the metric [Eqs. (39), (41) and (47)] at O(v), and thus, to extract the sensitivities we only need the strong-field solution at linear order in velocities.
Clearly, this procedure can in principle be pushed to next order in v, i.e. σ A can be extracted from the O(v 2 ) terms of the strong-field solution near spatial infinity. However, as already mentioned, σ A enters the dynamics of compact binaries at higher order than σ A in a small-velocity expansion. Because the orbital velocities of observed binary pulsars are at most v 12 /c ∼ 10 −3 (a value reached by the double binary pulsar PSR J0737-3039A [60,[89][90][91]), and because the center-of-mass velocities relative to the AEther are on the order of the peculiar velocity of our Galaxy (i.e. V CM /c ∼ 10 −3 ), effects proportional to σ A will be negligible relative to effects proportional to σ A .

B. A Strong-Field Route
As shown in the previous section, the sensitivity parameters can be calculated from the metric describing a body moving with velocity v relative to the AEther. One unsatisfactory aspect of that derivation, however, is that it uses a "weak-field" PN approach that models bodies with point-particles. Here, we will relax the weak-field and point-particles assumptions, and show that Eqs. (54) and (55) also follow from a strong-field definition of the sensitivities (i.e. one that relates the sensitivities to the interior structure of the star, where the PN/pointparticle treatment of the previous section breaks down). In particular, we will confirm that the calculation of the sensitivities only requires knowledge of the metric near spatial infinity, and at linear order in the velocity.
While the former is not that surprising (e.g. in GR, genuinely strong-field quantities, such as the ADM mass, can be defined through the asymptotic behavior of the metric or in terms of integrals over the extent of the body), the latter is. In fact, one can generalize the weakfield definition of the sensitivities [Eqs. (31) and (32)] to the strong field by replacing the gravitational mass of the point-particle,m(γ), with its strong-field counterpart, the total mass-energy of the body, M tot , and thus define It can be shown [92,93] that M tot is indeed the strongfield generalization of the point-particle massm(γ), because far away from the star g tt = 1 − 2G N M tot /r + ..., which agrees with the point-particle solution of Eq. (38) form(γ) = M tot . As for v (and γ), as mentioned above, they are the velocity (and Lorentz factor) of the body relative to the AEther, so in a coordinate system comoving with the star (that is, one where the star is at rest), they are defined in terms of the behavior of the AEther far from the star. In other words, in a system of asymptotically Cartesian coordinates where the motion is along the z-axis, one has What is confusing about Eq. (56), however, is that it would seem to imply that the mass M tot (and therefore the metric) need to be calculated at order O(v 2 ) to extract the sensitivities, in contrast with the results of the previous section. In what follows, we will show that Eq. (56) does however lead to the same expressions [Eqs. (54) and (55)] for the sensitivities as the PN treatment of the previous section, essentially thanks to Gauss' theorem. This will confirm that one can extract the sensitivity from the O(v) pieces of the metric alone. The procedure we outline below is similar to that used in scalar-tensor theories to extract the sensitivities of NSs from the asymptotic behavior of the scalar field at spatial infinity [86] (see also Ref. [94] for an approach similar to ours.)

Einstein-AEther Theory
As mentioned above, we consider a star at rest, and the AEther moving relative to it, so that (1/r) and g µν = η µν +O(1/r) in asymptotically Cartesian coordinates. We also assume that the system is in a stationary regime (i.e. that the metric and AEther do not depend on the time coordinate), which is expected to be the case after an initial transient. For such a system, we can write the total mass as [86] where L g is the Einstein Lagrangian density 7 , L AE is the AEther Lagrangian density [the Lagrange multiplier term (U µ U ν −1) is not included here because Eq. (57) is to be evaluated on shell], L mat is the matter's Lagrangian density, and Σ is a hypersurface of constant time. We stress that one can write the total mass in terms of the Lagrangian densities alone because the fields have no explicit time dependence (i.e. in the absence of time derivatives, the Hamiltonian is simply H =qp − L = −L), and because the Einstein Lagrangian density L g only depends on the metric and (quadratically) on its first derivatives. We also note that the definition given by Eq. (57) was also used in Ref. [92] for the study of black hole mechanics, and later in Ref. [56] in the context of GW emission. For khronometric theory, it was used in Ref. [38], where a brief summary with references to related work can be found. The sensitivities can then be obtained by varying the mass of Eq. (57). When doing so, the bulk terms evaluate to zero because of the Euler-Lagrange equations and one is left with surface terms alone [38,86,92]. For a star, the hypersurface Σ can be extended to its center, and we thus have only surface terms at the outer boundary (spatial infinity).
More in detail, because we are taking the difference between two neighboring solutions of the modified field equations, the AEther and metric variations δU µ and δg µν preserve the unit constraint U µ U µ = 1. Thus, the bulk terms lead to the modified field equations and thus vanish, even though Eq. (57) does not contain the Lagrange multiplier term. Similarly, the surface terms coming from the variation of L g vanish, because we are working in a gauge where g µν = η µν + O(1/r) [86]. Also, the surface terms coming from the variation of the matter Lagrangian vanish because we assume the matter fields to be confined within the star, while the boundary of the hypersurface Σ can be pushed to spatial infinity.
As a result, the only surface terms that contribute in Einstein-AEther theory are those coming from the AEther Lagrangian, Let us then evaluate Eq. (59) for a moving stationary configuration. Clearly, δU µ is the difference between the AEther 4-velocities of two neighboring moving solutions with velocities v and v + δv, while the derivative 7 The Einstein Lagrangian density, unlike the Einstein-Hilbert one (from which it differs by a total divergence), depends only on first derivatives of the metric.
∂L U /∂(∂ i U µ ) must be calculated for a moving solution with velocity v. It is more convenient to use asymptotically spherical coordinates, in which Eq. (59) becomes At lowest order and in Cartesian coordinates, the 4velocity variation is . In spherical coordinates, this becomes To calculate J r µ , we need to solve the modified field equations for a system comprised of a star at rest and an AEther moving slowly relative to it. In particular, expanding the equations at linear order in the AEther's velocity and near spatial infinity, one obtains the solution with where the length scale M * can be shown to be related to the total mass of the star, M tot , by M * = G N M tot [92,93], whileσ AE is a free parameter characterizing the asymptotic behavior of the solution, and which a priori might have no relation to the sensitivity of Eq. (56). To compute this parameter, one needs to solve the modified field equations in the stellar interior and match to an exterior solution, as we do in the next section. Not surprisingly, Eqs. (62) and (63) can also be obtained directly by transforming the PN solution in Eqs. (38), (39), (40), (50), (51) to the appropriate gauge and identifying σ →σ.
Using these expressions, one gets with which Eq. (60) can be evaluated explicitly to find where we have used the relation M * = G N M tot mentioned above. We therefore have and the sensitivity is which reduces exactly to Eq. (54).

Khronometric Theory
As discussed in Sec. II, khronometric theory can be thought of as Einstein-AEther theory with hypersurface orthogonality [Eq. (14)] being imposed in the action before variation. Equivalently, khronometric theory can be derived from the same action as Einstein-AEther theory, plus four Lagrange multipliers α ω α that enforce that the vorticity vector Such a formulation of the action has the advantage of containing only first (and not second) derivatives of the fields (and in particular of the AEther field U α ). This allows one to define Lagrangian densities in the usual way and proceed as in the previous section for Einstein-AEther theory to define the sensitivities in the strong field regime. Because the Lagrange multiplier terms vanish on shell, the mass of a time-independent configuration is given again by Eq. (57), exactly as in the Einstein-AEther case. When taking the difference between the mass of two neighboring solutions, the bulk terms vanish as a result of the modified field equations, and one is left with surface terms alone. Exactly as in the Einstein-AEther case, the surface terms coming from L g vanish in a gauge where g µν = η µν +O(1/r), and so do the terms coming from the variation of the matter Lagrangian, because the matter field are confined within the star.
As for the surface terms coming from the Einstein-AEther Lagrangian, using Eq. (14) we can express the AEther variation δU µ in terms of the variation δT of the foliation as A simple calculation then yields From Gauss' theorem, we then have (73) This expression must be evaluated with boundary con- In a suitable gauge, the asymptotic solution for a slowly moving star is given by with where M * = G N M tot [92,93] andσ kh is again a constant (a priori unrelated to the sensitivity) that characterizes the solutions and which needs to be calculated from the full slowly moving solution (i.e. from the solution describing both the interior and exterior of the star). As in the Einstein-AEther case, Eqs. (74) and (75) can be obtained by solving the modified field equations directly at O(v) near spatial infinity, or also by transforming Eqs. (38), (39), (40), (50), (51) after identifying σ →σ. From Eqs. (74) and (75), one obtainsǢ µ = O(v) and T = t + vr cos θ, so one can rewrite Eq. (73) as Finally, we can write and use Eqs. (74) and (75) to explicitly compute Inserting these expressions in Eq. (78), one easily obtains where we have used again M * = G N M tot . We then have and the sensitivity The strong equivalence principle is defined as the universality of free fall for strongly gravitating bodies. GR satisfies this principle, but this is clearly not the case for theories in which the sensitivities are not zero. This is because the sensitivities, as described previously, characterize how the structure of a compact object (i.e. a NS) changes with the motion relative to the ambient field in which the object is immersed (i.e. the AEther field in Einstein-AEther or khronometric theory). As we will show later, the sensitivities depend on the particular object under consideration, and in particular, on its compactness. Therefore, unless the sensitivities are exactly zero as in GR, different bodies respond differently to motion relative to the ambient field, and thus move along different trajectories, violating the strong-equivalence principle.
In the case of Einstein-AEther and khronometric theory, the sensitivities affect both the conservative and dissipative sectors. As for the former, the sensitivities modify Newton's universal gravitation law, i.e. at Newtonian order the motion of a binary is described by [51] v where Let us rewrite this expression as [51] v where we define the active gravitational masses as and the 2-body coupling constant Similarly, the sensitivities enter the equations of motion also at higher PN order in the conservative sector [51].
When it comes to the dissipative sector, and in particular to binary pulsars, the most well-known observable is the rate of change of the orbital period. In fact, it was the monitoring of this quantity that led to the first indirect detection of GWs by Hulse and Taylor [96][97][98]. For orbits satisfying Eq. (86), this observable can be written asṖ where P b is the orbital period and E b is the binary's binding energy. In deriving this relation, we have used the fact that the conservative sector is corrected only as in Eq. (85) to leading PN order. Thus, the binding energy is given by E b = −Gµm/(2a) and the orbital period by P b = 2πa 3/2 /(Gm) 1/2 , as in Newtonian orbital mechanics, but in terms of the 2-body coupling constant G and the active gravitational masses m A . Here, a is the semi-major axis and are the reduced (active) mass and the total (active) mass respectively. Equation (89) can be further manipulated by relating the rate of change of the binding energy to the total flux of energy F carried away from the system, In GR, F is only due to the propagation of tensor modes (i.e. GWs), but in modified theories one normally finds radiation from scalar and vector modes. The balance law in Eq. (91) should hold because the total energy of the system is a conserved quantity [38,92,93]. The energy flux carried by all propagating degrees of freedom can be computed by investigating the Noether charges and currents in the theory under consideration. These charges and currents depend on the solution to the evolution equations for the metric perturbations on a Minkowski background in the far-zone, i.e. at distances much larger than the GW wavelength as measured in the system's center of mass frame. In what follows, we will summarize how to compute this flux in Einstein-AEther and khronometric theory, which then leads to a prediction for the orbital decay rate via Eq. (89).

Dissipative PN Dynamics in Einstein-AEther Theory
The energy flux in Einstein-AEther theory and for a matter action of the form of Eq. (31) was first calculated by Foster [51]. After re-deriving these results we have found important discrepancies with Foster's work, which justifies showing the calculation in some detail here. In the following, we present the correct formulas that correct the miscalculations in Refs. [51].
The energy flux is carried by degrees of freedom that propagate away from the source. In this region, it is convenient to treat the metric and AEther fields as perturbations around a Minkowski background η µν with the AEther pointing in the time direction and decompose their spatial components as where ∆ is the Laplacian differential operator and Once the previous decomposition is used in the field equations, Eqs. (5) and (35), one can follow the calculation in Ref. [51] to find the waveform at a distance r from the source (note the opposite sign due to different conventions) where is the speed of propagation of the φ ij modes and with I ij given by for a system of A bodies. Recall that here x i A (t) are the trajectories of the A-th point particle, while v i A (t) =ẋ i A are their 3-velocities.
The transverse-traceless projector is built with the unit-norm vectorn i ≡ r i /r, where r i are the coordinates of the point where the fields are being computed [38]. To derive the previous expressions we used the conservation properties of an improved energy-momentum tensor as defined in Ref. [51]. For the vector modes, we choose the gauge φ i = 0. The field γ h i can be solved as up to terms with vanishing time derivative. Finally, the equation for ν i yields Here, the different tensors on the right-hand side are to be evaluated at a retarded time t − r/w 1 , where w 1 is the speed of propagation of the vector modes, Also, Q ij is the trace-free part of the rescaled mass quadrupole moment and we have also introduced the tensors 8 8 The tensor V ij is absent in the calculations of Ref. [51], which is a mistake. For the bounds derived later, this term plays a subleading role.
Finally, for the scalar sector, once we impose the gauge ν = γ h = 0, the equations for the fields φ h and h 00 in the far-zone read 9 while where we have defined and the speed of the F mode is  43)]. Note that although w 0 is the same as that defined in Ref. [49,51,56], Z and the structure of Eq. (106) differ from Ref. [51,56]. We have checked 10 that our expressions are consistent when taking the weak-field limit [56], which is not the case for Ref. [51].
Using the previous equations and the expressions in Ref. [51,56] we can compute the total flux F, and via Eq. (89) the rate of change of the orbital period: 9 At the PN order we work to, we do not need the O(1/c 4 ) in Eq. (104) since it corresponds to a conserved quantity, see Ref. [51]. 10 Another way to check the consistency of the approach is to use the limit of Ref. [84] to compare with the khronometric results to be shortly derived. Also in this case, our results are consistent.
In these equations, the angled-brackets stand for an average over several wavelengths and we have defined the shorthands Let us now evaluate the different terms in Eq. (109) for a compact binary system in a generic orbit in the center of mass frame. Using basic results from Newtonian and PN orbital mechanics, we can evaluate the moments in Eqs. (96)-(97) (cf. Ref. [51]). The only expression not found in Ref. [51] iṡ /m are the center-of-mass coordinates [see also Eq. (90) and the definitions around Eq. (85)]. These moments are to be used in Eq. (109). As explained in Ref. [51], the secular terms that depend on X i CM can be neglected, as they are canceled by an opposite contribution in the action. The final result iṡ where recall that v i 12 ≡ṙ i 12 is the relative 3-velocity, V i CM ≡Ẋ i CM is the center-of-mass velocity of the system relative to the AEther field and S ≡ s 1 m 2 /m + s 2 m 1 /m.
At this juncture, let us stop to discuss some qualitative features of the orbital decay rate predicted in Einstein-AEther theory [Eq. (116)]. First, note that in the limit in which all the coupling parameters go to zero, i.e. c i → 0, Eq. (116) reduces to the GR result. Second, the Einstein-AEther terms can in principle dominate over the GR ones, under appropriate conditions (i.e. unequal and sufficiently large sensitivities, and sufficiently large couplings). This is because the leading-order Einstein-AEther term scales with fewer powers ofm/r 12 than the GR result. Indeed, the Einstein-AEther terms enter first at absolute order O(1/c 8 ), whereas the GR prediction is of absolute order O(1/c 10 ). Such scalings are typical of dipolar radiation, and those terms are proportional to the difference of the sensitivity parameters squared. The existence of a preferred frame is clear in the previous expression from the dependence on the velocity of the center of mass with respect to the AEther, V i CM . Finally, note that for circular orbits,ṙ 12 = 0 andn i 12 v i 12 = 0, and the above expressions simplify greatly.

Dissipative PN Dynamics in Khronometric Theory
The energy flux in khronometric theory was first calculated in Ref. [38], which focused on weak-field sources. We will here present the main features of their analysis, together with the generalization to the strong field case.
The total energy flux carried by propagating degrees of freedom can be computed in khronometric theory, using the decomposition in Eq. (92), to be [38] is the speed of propagation of the scalar modes. The termȮ corresponds to a boundary term given by a time derivative. This term has two components: one that cancels when averaged over the orbit and another one that cancels against secular terms (cf. Ref. [51]). The far-zone fields appearing in Eq. (117) can be obtained by solving the modified field equations from Sec. III A 1, and using the method described in Ref. [38] (see also Ref. [51]). The results are In the expression for F , the different moments should be evaluated at the retarded time t − r/c s , and we have defined After performing the angular integration in Eq. (117), the orbital decay rate is still given by Eq. (109), but now with the coefficients Clearly then, the same observations made after Eq. (109) are still valid here. In particular, the orbital decay rate is of dipolar structure and depends on the sensitivities.

IV. NEUTRON STAR SOLUTIONS IN EINSTEIN-AETHER THEORY
In this section, we construct Einstein-AEther theory solutions that describe isolated non-spinning NSs moving slowly relative to the AEther. We begin by presenting the metric and AEther ansatz, and expand the field equations at zeroth-and first-order in the velocity v relative to the AEther. We then explain how to solve these equations numerically.

A. Metric Ansatz, AEether Field Ansatz and Neutron Star Model
Let us consider a generic ansatz for the metric and AEther fields, for a stationary configuration describing a non-spinning NS in slow motion with velocity v i relative to the AEther.
One can adopt a coordinate system comoving with the fluid elements of the NS by aligning the time coordinate vector to the fluid's 4-velocity u µ . More specifically, we can choose a spacelike hypersurface Σ and assign spatial coordinates on it. One can then define the spatial coordinates of an event p as those of the intersection between Σ and the worldline of the fluid element passing through p. 11 In these comoving coordinates, the fluid elements are at rest, while the AEther is moving. In particular, the fluid 4-velocity field is Because the system is invariant under spatial rotations around the direction of the motion relative to the AEther, it is convenient to adopt asymptotically spherical coordinates (t, r, θ, ϕ) in which the AEther's velocity direction corresponds to the polar axis. Clearly, in these coordinates both the metric and AEther will not depend on the azimuthal coordinate ϕ. Noting then that the system's configuration must be invariant under the simultaneous reflection of it is clear that at first order in the velocity, the only nonzero components of the metric are the cross-terms g tr and g tθ , and the only non-vanishing components of the AEther field are U r and U θ . (Note that a contribution to U t at O(v) is forbidden by the normalization condition U µ U µ = 1.) More explicitly, the most generic ansatz at linear order in v for the metric and AEther is given by ds 2 = e ν(r) dt 2 − e µ(r) dr 2 − r 2 dθ 2 + sin 2 θdϕ 2 + 2vV (r, θ)dtdr + 2vrS(r, θ)dtdθ + O(v 2 ) , (128) and the fluid 4-velocity Eq. (125) becomes Note that for a star, the AEther cannot have any U r component at O(v 0 ), as shown explicitly in Ref. [57] using the field equations. The ansatz of Eqs. (128)-(130) therefore depends on two potentials µ(r) and ν(r) at order O(v 0 ), and on four potentials V (r, θ), S(r, θ), W (r, θ), Q(r, θ) at O(v). This ansatz, however, has been derived by choosing spatial coordinates comoving with the fluid without specifying the time coordinate. In particular, one is free to perform a coordinate transformation of the form 11 Note that this construction yields a bona fide coordinate chart in the absence of caustics for the fluid flow (which is a reasonable assumption for the system that we are considering) and in the absence of closed timelike curves (which would in any case violate causality and make the spacetime pathological).
which can be used to set any one of the potentials V (r, θ), S(r, θ), W (r, θ), Q(r, θ) to zero while keeping the ansatz of Eqs. (128)-(130) valid at O(v) (modulo a redefinition of the remaining three potentials). In what follows, we will set therefore Q = 0 without loss of generality. We use a perfect fluid stress-energy tensor to describe the NS matter: whereρ(r) andp(r) are the fluid's energy density and pressure at O(v 0 ), respectively. The density and pressure do not have any O(v) contributions because they must be invariant under the transformations in Eqs. (126) and (127). Moreover, one can explicitly show that any O(v) pieces in the density and pressure must vanish from the O(v) pieces of the (t, t) and (r, r) components of the modified field equations. Note also that, in spite of these facts, the (t, r) and (t, θ) components of T mat µν have O(v) contributions due to the O(v) terms in the metric.
The relation between internal pressure and energy density is parametrized by the EoS. We here investigate four different, realistic EoSs: APR [99], SLy [100], Shen [101,102] and LS220 [103]. For the Shen and LS220 EoSs, we use a temperature of 0.1 MeV and consider neutrino-less, β-equilibrium. These EoSs are not the only ones that are observationally viable, but they represent a sufficiently large sample of the EoSs to allow us to investigate how the different observables depend on this choice.

O(v 0 ) Equations
Let us first look at the field equations at O(v 0 ). The (t, t), (r, r), and (θ, θ) components of the field equations give three independent equations respectively, where p(r) and ρ(r) are the rescaled NS pressure and density with the length scale M (r) defined by One may find it instructive to expand these equations further in a small coupling approximation, i.e. c 14 1. Doing so, the above equations become The first terms in the above equations agree with the GR result for the TOV equation, see e.g. Ref. [104,105].

O(v) Equations
Let us now look at the field equations at O(v). The modified field equations and the AEther field equations become a system of partial differential equations for the unknown functions V (r, θ), S(r, θ) and W (r, θ). One can separate the variables r and θ by using the following Legendre decomposition: S(r, θ) = n s n (r) dP n (cos θ) dθ , W (r, θ) = n w n (r)P n (cos θ) , where P n is the n-th Legendre polynomial. This structure can be inferred from symmetry principles by looking at how each component of the metric behaves under spatial rotations and reflections with respect to the equatorial plane, i.e. this is nothing but a tensor spherical harmonic decomposition (see e.g. Ref. [106]).
With this decomposition at hand, the (t, r) component of the modified field equations and the r and θ components of the AEther field equations become Note that the n-th mode is independent of all other modes, i.e. there is no mode coupling (see Appendix A for an explanation). Note also that, as expected, the new unknown functions at O(v) depend on the solutions to the O(v 0 ) equations of the previous subsection.

C. Asymptotic Solutions, Matching Conditions and Numerical Techniques
The above system of differential equations must be solved order by order in velocity twice: once in the interior of the NS with a given EoS, and once in the exterior where the density and pressure vanish. We solve these equations as an initial value problem. The solutions will generically depend on integration constants. These are determined by imposing that the solutions be continuous and differentiable at some matching radius, usually the NS surface.
Let us first concentrate on the initial conditions for the O(v 0 ) equations. By imposing regularity at the NS center, the behavior of the solution to the O(v 0 ) equations near the center is where ρ c , p c and ν c are the values of the density, pressure and ν at the NS center. The quantity ρ 2 is a constant that can be expressed in terms of ρ c and p c via Eqs. (152), (153) and the EoS. Similarly, by requiring asymptotic flatness, the asymptotic behavior of the solution near spatial infinity is where M * ≡ M (∞) = G N M obs , with M obs the NS mass observed by a Keplerian experiment.
Given these asymptotic solutions, we numerically solve the O(v 0 ) differential equations as follows. First, we choose a specific value of the NS central density ρ c , and with it, we obtain p c from the EoS p c = p(ρ c ). We obtain the interior solution by solving Eqs. (139) and (141), together with the EoS, from a core radius r R * up to the NS surface, defined as the radius R * at which the pressure vanishes. We then use the values obtained from this interior solution evaluated at the NS surface as initial conditions for a new evolution of the equations in the exterior. The exterior solution for M (r) is then evaluated at a large boundary radius r b R * and set equal to the asymptotic expression in Eq. (156), evaluated at the same radius r b , to solve for the observed NS mass M obs .
Next, we solve Eq. (140) for ν(r) in the interior (from r = r to the NS surface) using a trial value for ν c . As before, we then use the value of ν(r = R * ) as an initial condition to solve for ν(r) in the exterior. Since Eq. (140) is shift-invariant, the trial solution plus a constant ν tr (r) + C ν is also a solution. We choose C ν such that ν(r) satisfies the asymptotic behavior of Eq. (157) at r = r b [11,107] Note that in practice r b < ∞, but we choose it to be sufficiently large such that the error incurred is smaller than other numerical errors, e.g. those due to the discretization of the differential equations.
Let us now discuss the initial conditions for the O(v) equations. Equations (39) and (62) imply that we are only interested in the n = 1 component of the Legendre decomposition, since these functions determine the sensitivities, and thus, the orbital period decay rate. Imposing regularity at the NS center, the solution to Eqs. (148)-(150) near the center is where C and D are constants of integration. Similarly, we can impose asymptotic flatness at spatial infinity and obtain the asymptotic solution where A and B are constants of integration, s A1 and s c1 are given by and the other coefficients, k A2 , k c2 , k A31 , k c31 , s A2 , s c2 , s A30 , s c30 , s A31 , s c31 , w A2 , w c2 , w A30 and w c30 , similarly depend only on c i (i = 1, 2, 3, 4). We do not present explicit expressions for these quantities here because they are lengthy and unilluminating. One might be worried that the first terms in Eqs. (162) and (163) seem to be inconsistent with the asymptotic flatness condition, because they lead to terms in the metric of the form g ti ∼ vdx i dt. However, such terms are coordinate artefacts; for example, the Schwarzschild metric in Eddington-Finkelstein coordinates has a component that goes as g vr ∼ dvdr. An inverse coordinate transformation t ≈ t − vz with z = r cos θ would then eliminate such terms in Eqs. (162) and (163), but we choose not to do so here because the equations of structure are simpler in these coordinates. The asymptotic behavior above is consistent with the one found by Foster [51], as explained in Appendix B.
The numerical solution to the modified field equations at O(v), i.e. Eqs. (148)-(150), can be obtained by exploiting the linearity and homogeneity of the differential system. First, let us arbitrarily choose two sets of values for the constants (C, D) in Eqs. (159)-(161) and solve Eqs. (148)-(150) from r = r to the NS surface. Then, we evaluate these interior solutions at the NS surface and use them as initial conditions to numerically find the exterior solution. Doing so, we obtain two homogeneous solutions A (1) i and A (2) i , where A i ≡ (k 1 , s 1 , w 1 ), everywhere in the numerical domain. With this at hand, the general solutions A i can be obtained by linear superposition The new constants C and D must be determined, together with A and B in Eqs. (162)-(164), by requiring the following matching conditions at spatial infinity: As before, we evaluate the matching at a fixed radius r b R * that is large enough for the errors to be smaller than those due to the discretization of the equations. We will later check that our results are insensitive to the choice of core and boundary radii.
The numerical techniques used to solve the O(v 0 ) and O(v) differential system are the following. As already mentioned, we treat the system as an initial value problem, with initial conditions given by the asymptotic behavior of the solution about the NS center. We then numerically solve the equations with an adaptive stepsize, 4th-order Runge-Kutta method [108]. All numerical solutions presented in this paper are obtained with this numerical algorithm. In all cases, we have checked that increasing the numerical resolution does not affect the results presented.
Finally, to obtain the sensitivities we transform the asymptotic solution for the metric [Eqs. (162) and (163)] to the gauge of Eqs. (38)-(40) (see Appendix B for the explicit calculation). This relates the integration constant A, which we can determine from our numerical solution, to the constants B ± appearing in Eqs. (54) and (69) through With this relation at hand, we can then write the sensitivities for an isolated star in Einstein-AEther theory in terms of the constant A as (173)

V. NEUTRON STAR SOLUTIONS IN KHRONOMETRIC THEORY
In this section, we construct solutions of khronometric theory that describe non-spinning NSs moving slowly relative to the AEther. We follow the same organizational structure as in Sec. IV, and thus begin by introducing the metric and khronon field ansatz for a generic nonspinning stationary and slowly moving system, and then present the field equations and explain how to solve them numerically.

C. Asymptotic Solutions, Matching Conditions and Numerical Techniques
The initial conditions for the O(v 0 ) equations are exactly the same as in Einstein-AEther theory, as described in Sec. IV C, since the differential equations are identical. The initial conditions for the O(v) equations, on the other hand, will be different. Requiring that the metric be asymptotically flat, the asymptotic behavior of the solution to Eqs. (177) and (178) for the n = 1 modes is k ∞ 1 (r) = −1 + A M * r + (7 + 11β + 18λ)A − 6(6 + α)λ − 2(2α + 11)β − 2(α + 7) 8(β + λ) Imposing regularity at the NS center, the solution near the center is The matching conditions and procedure are identical to those described in Sec. IV C for Einstein-AEther theory. Similarly, the numerical techniques are the same, and the extraction of the sensitivity is done by transforming the asymptotic solution for the metric [Eqs. (183)

VI. NUMERICAL NEUTRON STAR SOLUTIONS AND SENSITIVITIES
In this section we present the results obtained by numerically solving the modified field equations. In particular, we concentrate on deriving numerical results for the sensitivities and developing an analytic fitting formula. We first tackle the Einstein-AEther case, and then move on to khronometric theory.

A. Einstein-AEther Theory
Let us first focus on the numerical solutions at O(v 0 ). At this order, the main observable is the relation between the NS mass and its radius, for a sequence of NSs in a given EoS family. As can be seen from the equations of Sec. IV B 1 (and as already noted in Ref. [57]), these solutions depend only on the c 14 combination of the coupling constants [see Eq. (45)]. Figure 2 shows the mass-radius relation in Einstein-AEther theory for different values of c 14 . The horizontal line at M = 1.97M is the lower mass bound derived from observations of PSR J0348+0432 [59]. Observe that, as one increases c 14 , the NS mass decreases for a fixed radius, which is consistent with the conclusions in Ref. [57].  [59]. Observe that as c14 is increased, the NS mass decreases for a fixed radius.
This figure is a perfect example of the strong degeneracy between the EoS and modified gravity effects, which in turn prevents us from constraining modified theories with the mass-radius relation alone. For example, if we knew that the LS220 EoS was the correct one, then we could argue that the observation of PSR J0348+0432 [59] requires that c 14 < 0.1. However, we do not know the correct EoS and, for instance, the APR EoS could be the correct one. If that were the case, we would not be able to place competitive constraints on c 14 . We then conclude that the observation of PSR J0348+0432 [59] (or any other system for that matter) is ineffective at constraining Einstein-AEther theory through the mass-radius relation.
Let us now consider the O(v) solutions, and in particular, the sensitivities in Einstein-AEther theory. In the weak-field limit, i.e. expanding in the ratio of the binding energy Ω to the NS mass M obs , one can show that the sensitivity scales as [51] with α AE 1 and α AE 2 given by Eqs. (42) and (43), while [49,56] with r = |x| and r = |x |. When plotting the weakfield sensitivity using Eq. (184), we evaluate Ω by using the Legendre expansion of the Green's function of the Laplacian operator [109]. Of course, this integral depends on the EoS through ρ(r). Figure 3 shows the absolute magnitude of the sensitivity in Einstein-AEther theory, calculated from the numerical solution to the O(v) modified field equations, for different EoS as a function of NS compactness C * = M * /R * = G N M obs /R * . For comparison, we also plot the weak-field expression for the sensitivity [Eq. (184)] with the APR EoS. The bottom panel shows the fractional difference between the actual sensitivity and its weak field approximation. Observe that as C * increases, Eq. (184) becomes highly inaccurate, with errors of roughly 15-30% for realistic NS compactnesses, i.e. for compactnesses 0.1 C * 0.3. Figure 4 shows the absolute magnitude of the sensitivity as a function of c − for c + = 1 and c + = 0.1, saturating the Solar System constraint for α AE 1 and α AE 2 . Observe that the sensitivity increases monotonically with c − and c + .
One may wonder whether the results presented here are robust to the method used to obtain the numerical solutions of the modified field equations. Two key features of this method are the choice of core radius (where the interior integration is numerically started) and the value of the boundary radius (where the matching is carried out). Figure 5 shows the sensitivity as a function of the matching radius (top panel) and the core radius (bottom panel) for the APR EoS. Observe that the fractional error on the sensitivities is always smaller than 0.1%. . We saturate the PPN parameters α AE 1 and α AE 2 such that they satisfy the Solar System constraint, and c+ = 10 −4 = c−. Observe that, as expected, the weak-field result is highly accurate for small compactness, but it becomes inaccurate for realistic NS compactnesses.  Figures 3 and 4 show that the sensitivities in Einstein-AEther theory do not depend strongly on the EoS. Because of this, it is natural to develop an EoS-independent fitting formula. First, we choose the LS220 EoS as a representative EoS with which to compute sensitivities and we set α AE 1 = 10 −4 and α AE 2 = 4 × 10 −7 , thus saturat- ing Solar System constraints. With this, we then numerically compute the sensitivity in the range c + ∈ 3 × (10 −4 , 10 −1 ), c − ∈ 3 × (10 −5 , 10 −3 ) and C * ∈ (0.11, 0.28) for a total of 10 3 points. We choose this region in the parameters (c + , c − ) as it will still be allowed after imposing binary pulsar constraints (see Fig. 1).
With the data in hand, we choose a functional form for the fitting function. By plotting the sensitivity as a function of c + , c − and C * , we empirically find that s AE is well-described by a polynomial in these quantities (at least in the region of parameter space described above). We could use a fitting function that corrects Eq. (184) with higher powers of C * , but we find that this is not as good as We fit the data described above with Eq. (186), using the fact that each numerical point is known to at least a fractional numerical accuracy of 1%. The fit returns an r 2 value of 0.959, with fitted coefficients given in Table I. We also carried out a more accurate fit by keeping more terms in the polynomial expansion of Eq. (186); the r 2 is then 0.9994 and the coefficients are given in Appendix C. Note that c 0,0,0 = O(α AE 1,2 ) = 0, because we have chosen α AE 1,2 such that Solar System constraints are saturated. Thus, the fitting function does not go to zero as c + → 0 and c − → 0.

B. Khronometric Theory
Spherically symmetric solutions in khronometric theory are identical to those of Einstein-AEther theory [after the substitution in Eq. (15)] [52,74,81,82] . Thus, the mass-radius relation for NSs will also be identical to what we presented in Fig. 2. We recall, however, that stability/Cherenkov requirements, Solar System experiments and BBN bounds restrict the couplings λ, β and α to very small values (cf. Fig. 1), so the deviations away from the GR mass-radius relation will be even smaller than in Einstein-AEther theory. For this reason, and because of the degeneracies between the mass-radius relation and the EoS, khronometric theory cannot be constrained by measurements of NS pulsar masses alone.
Let us then focus on the numerical solutions to the modified field equations at O(v). In the weak-field limit, one can show that the sensitivities in khronometric theory scale in the same exact way as in Einstein-AEther theory [Eq. (184)], except that now α kh 1 and α kh 2 are given by Eqs. (48) and (49). Figure 6 shows the absolute magnitude of the sensitivity in khronometric theory, calculated from the numerical solution to the O(v) modified field equations, for different EoSs as a function of β (left panel) and NS compactness C * (right panel). Here we have chosen α and λ by saturating Eqs. (23) and (24). The bottom right panel shows the fractional difference between the actual sensitivity and its weak field approximation, where we recall that the latter is given by Eq. (184) with α AE 1 and α AE 2 replaced by α kh 1 and α kh 2 . Observe first that the behavior of the sensitivity in khronometric theory as a function of β is rather different from that of the sensitivity in Einstein-AEther theory as a function of c − . Observe also that as C * increases, the weak-field result becomes highly inaccurate, with errors of roughly 50-100% for realistic NS compactnesses.
As in the Einstein-AEther case, Fig. 6 shows that the sensitivities in khronometric theory are also rather insensitive to the EoS. Because of this, it is natural to develop an EoS-independent fitting formula. Again, we choose a fitting function of the form and fit this to sensitivities computed numerically using the LS220 EoS and choosing α = 2β due to current Solar System constraints [see Eq. (22)]. We carry out this fit with data in the region β ∈ (5 × 10 −5 , 5 × 10 −3 ), λ ∈ (3 × 10 −3 , 10 −1 ) and C * ∈ (0.11, 0.26) with a total of 10 3 points. This range in β and λ for the fits is chosen in this way because these values remain viable after placing binary pulsar constraints. As in the Einstein-AEther case, when carrying out the fits, we use the fact that each numerical point is known to at least a fractional accuracy of 1%. The fit returns an r 2 value of 0.99996, with fitted coefficients given in Table I. Notice that this fit is signif-AE c0,0,0 c0,0,1 c0,0,2 c0,1,0 c0,1,1 c0,1,2 c0,2,0 c0,2,1 c0,2,2  In both cases, we set the PPN parameters α kh 1 and α kh 2 by saturating the Solar System constraint [we saturate Eqs. (23) and (24)]. The bottom right panel shows the fractional difference between the numerical sensitivity and its weak-field value. Observe that for large compactnesses, the weak-field result can be very inaccurate.
icantly better than the Einstein-AEther one, so we do not need to re-do the fit with a higher-order polynomial.
We will not repeat Fig. 5 for khronometric theory, but we have checked that indeed our results are accurate to better than 0.1%, and thus, unaffected by the way the numerical solution is obtained.

VII. PULSAR CONSTRAINTS
In this section, we discuss how to place constraints on Lorentz-violating theories with binary pulsar observations. We begin with an overview of how binary pulsars can be used to test GR. We then specialize our discussion to Lorentz-violating theories. We conclude with an implementation of these ideas to constrain Einstein-AEther and khronometric theory.

A. Binary Pulsars As Laboratories for Fundamental Physics
Einstein-AEther and khronometric theory predict that the orbital period decay of binary systems should generically be faster than that predicted by GR because of the emission of dipolar radiation, i.e. due to a correction to the dissipative sector of GR. All observations of binary pulsars to date, however, have verified the GR prediction. Therefore, binary pulsars can be used to place constraints on these theories, by essentially requiring that the non-GR effects be smaller than observational uncertainties [110].
Binary pulsar astrophysics is carried out by fitting a suitably-averaged set of observed pulses to a given pulsar model. The output of this procedure is then a set of best-fit model parameters that describe the observation (e.g. the rate of change of pericenter, the orbital decay rate, the rate of change of the line of nodes, etc.) with a given observational error estimate for each of these parameters. This observational error, of course, decreases with observation time, as more data is accumulated.
Some parameters of the binary system, however, are not directly measurable, but are rather inferred from other observables. Let us for example consider the individual masses of the binary. These are inferred by noting that other observables are functions of the individual masses, once one chooses a gravitational theory. Since these observables are measured up to some uncertainty, this error also propagates into the inferred masses.
More specifically, the inferred individual masses can be determined from two binary observables that depend only on the conservative sector of the theory, e.g. the rate of change of the pericenter and the Shapiro time-delay. This inference can be carried out with the leading-order, Newtonian expressions for these post-Keplerian observables, since 1PN order corrections and higher will be greatly subdominant. In Einstein-AEther and in khronometric theory, the conservative sector is modified to leading PN order only through the substitutions [cf. Eq. (85)]. This leads to corrections of O(σ A ) to the masses measured assuming GR as the fiducial gravitational theory.
One may wonder whether one has to account for corrections of O(σ A ) in the masses, as well as 1PN order corrections to the conservative sector, when placing constraints on Einstein-AEther and khronometric theory with observations of the orbital decay rate. This is not the case because in these theoriesṖ b /P b differs from the GR prediction at leading PN order due to corrections to the dissipative sector. More precisely, deviations away from GR appear at -1PN order for NS-white dwarf systems, and at 0PN when the -1PN terms are subdominant (such as in the case of binary NSs; see discussion after Eq. (193) for more details). Thus, corrections of O(σ A ) and those of 1PN to the conservative sector are always smaller than observational errors.
Once one has determined the inferred individual masses and measured the orbital period of the system, one can calculate the predicted rate of change of the orbital period and compare this to the observed value. To date, the GR prediction has always agreed with the ob-served value up to the observational uncertainty δ obs inṖ b and the inferred error δ inf induced in the predicted value ofṖ b /P b due to the error that propagates from the individual masses.
In Einstein-AEther and in khronometric theory, the predicted orbital decay rate differs from the GR one: where δṖ AE,kh is the fractional difference in the prediction. In fact, because of the presence of dipolar radiation, δṖ AE,kh is much larger than unity for small velocities, unless the coupling constants are very small. As we will see later in this section, 10 12 /c 10 ), with v 12 the magnitude of the binary's orbital velocity.
Comparing Eqs. (188) to (189), one then finds that the most conservative constraint comes from demanding In what follows, we will instead use the equivalent relation to place constraints on the coupling constants κ i AE,kh , where κ i AE = (c + , c − ) in Einstein-AEther theory and κ i kh = (λ, β) in khronometric theory. Here m A and P b are the best-fit parameter for the (inferred) individual masses and the orbital period, while δm A and δP b are their associated statistical errors. Recall that, in these theories, the non-GR correction to the prediction of the orbital decay rate not only depends on the individual masses and the orbital period, but also on the EoS through the sensitivities. Thus, we will search for the region in the twodimensional coupling parameter space where Eq. (191) is satisfied, with δṖ AE,kh evaluated within m A ± δm A , P b ± δP b and various values of the sensitivities, determined by the EoS and the individual masses. Let us emphasize that Eq. (191) does not lead to a line in the two-dimensional coupling parameter space, but rather a two-dimensional region.

B. Binary and Isolated Pulsars As Laboratories for Lorentz-Violation
To date, all robust pulsar constraints on Lorentzviolation come from studies of its effects on the conservative dynamics of the system that then propagate to the observed pulse sequence. Such effects are modeled with the parameterized post-Newtonian (PPN) formalism [65,[111][112][113][114][115], where the point-particle Lagrangian is modified to include preferred-frame effects. The latter are proportional to certain PPN parameters (the relevant ones here are α 1 and α 2 ), as well as contractions of preferred-frame-related velocity vectors and PN vector potentials [3,88]. Such modifications to the Lagrangian induce effects on pulsar observables, which can then be used to place constraints on α 1 and α 2 . Let us classify these in terms of the sources used: isolated pulsars and binary pulsars.
Let us first consider the isolated pulsar case. In Lorentz-violating theories of gravity, the presence of a preferred frame induces precession of the pulsar's spin axis relative to this frame, pushing the pulsar beam in and out of the line of sight [116,117]. This effect has not been observed in isolated pulsars, which then allows for constraints [45] on the strong-field generalization (α 1 ,α 2 ) [117] of the PPN parameters (α 1 , α 2 ). The angular precession velocity of the spin angular momentum vector due to preferred frame effects, and thus, the rate of change of the angle between the spin axis and the line of sight, is proportional toα 2 [116,117], but also depends on the magnitude and direction of the relative 3-velocity of the isolated pulsar with respect to the preferred frame. Using data from PSR B1937+21 [118] and PSR J1744-1134 [119] that shows no such precession, Shao, et al [45] placed the constraintα 2 < 1.6 × 10 −9 at 95% confidence. Here, they associate the preferred-frame 3-velocity with our Solar System barycenter velocity relative to the CMB, randomizing over the pulsar's radial velocity with respect to us and the spin orientation.
Let us now consider the binary pulsar case. For systems with small eccentricity, preferred-frame effects parametrized byα 1 andα 2 decouple. The former induce a correction to the rate of change of the eccentricity, which in turn affects the precession of the pericenter. The latter induce precession of the orbital angular momentum around the center of mass' 3-velocity vector in the preferred frame, leading to a change in the inclination angle with respect to the line of sight, and thus, precession of the projected semi-major axis [110]. The orbital plane and pericenter precession produces noticeable effects in the pulse profiles, all of which depend on the magnitude and direction of the 3-velocity of the binary system in the preferred frame.
All current pulsar observations match GR predictions without these extra precession effects, which then allows for constraints onα 1 andα 2 . Using data from pulsarwhite dwarf binaries with known three-dimensional velocity, PSRs J1012+5307 [120] and J1738+0333 [61], Shao and Wex [44] placed the constraintsα 2 < 1.8×10 −4 andα 1 10 −5 at 95% confidence. We see that the constraint onα 2 is not as strong as that from isolated pulsars. These constraints depend on the magnitude and direction of the 3-velocity of the binary system in the preferred frame, which Shao and Wex [44] took to be its velocity relative to the CMB.
We will show later that the isolated pulsar constraint is not suitable to bound Einstein-AEther or khronometric theory. This is for one main reason. PSR B1937+21 [118] and PSR J1744-1134 [119] are isolated, millisecond pulsars, and thus their masses have not been measured, precisely because they are isolated. As we will see later, α 1 andα 2 in Einstein-AEther and khronometric theory depend on the sensitivities, which cannot be determined without knowing the mass, even if the EoS were known. Without a priori knowledge of the mass, one cannot map the isolated pulsar constraint onα 1 andα 2 to a constraint on (c + , c − ) or (λ, β). Note that this will be the case for any Lorentz-violating theory of gravity.

C. Constraining Lorentz-Violating Gravity with Binary Pulsars
Let us first concentrate on observations of the orbital decay rate of PSR J1141-6545 [58], PSR J0348+0432 [59] and PSR J0737-3039 [60]. The first two are binary systems composed of a NS and white dwarf, in a 0.17 eccentricity, 4.74-hour orbit and in a O(10 −6 ) eccentricity, 2.46-hour orbit respectively. Therefore, these systems have an orbital velocity of v 12 /c = O(10 −3 ). The last system is a double binary pulsar with 0.088 eccentricity and 2.45-hour orbit. In all cases, the eccentricity is small and will thus be neglected. Our constraints are robust to this assumption, because eccentricity corrections to δṖ AE,kh are of O(e 2 ), which is negligible in Eq. (191) relative to the uncertainties due to the dependence on the individual masses, orbital period, and EoS. Of course, one would not be able to neglect eccentricity if considering certain binary pulsar systems, such as the very wellstudied PSR 1913+16 [96][97][98], which has an eccentricity of roughly 0.6.
One can neglect the spin of the pulsar when modeling the orbital dynamics of binary systems to test GR for the following reason. In the small-spin or slow-rotation approximation [121][122][123][124], all observables can be expanded in a series inā 1. To zeroth-order in this approximation, one recovers the non-spinning results used in this paper. To first-order in this approximation, only observ-ables that depend on the gravitomagnetic sector of the metric, i.e. the "cross" time-space components, are modified. One of these observables is precisely the orbital decay rate, which acquires a spin-orbit coupling correction. Such a correction, however, appears at 1.5PN order), i.e. the leading-order spin correction to the orbital decay rate is actually of O(ā v 3 12 /c 3 ) smaller than the leading-order terms, which makes this effect negligible for our purposes. Similarly, there will be spin corrections to observables associated with the conservative dynamics, i.e. spin corrections to the Hamiltonian. As in the orbital decay rate, however, they are suppressed in the PN sense, and thus, can be neglected. Given this, one can model the pulsars in the binary systems that we consider as if they were non-spinning. This allows us to use the results of the previous sections, which do not account for NS spins.
With all of this at hand, consider then a binary system with non-spinning components in a circular orbit in either Einstein-AEther or khronometric theory. Because of the circularity condition, all terms proportional toṙ 12 vanish in Eq. (116) to leading PN order. After orbit-averaging, the orbital decay rate to leading PN order reduces tȯ where we recall one more time that m A are the active masses, m = m 1 + m 2 is the total active mass, P b is the orbital period and we have defined There are many interesting features in Eqs. (192) and (193) that deserve further discussion. First, notice that Eq. (192) reduces to the GR prediction when A = 1, which is the case when the constants c i → 0 [recall that C also depends on the c i via Eq. (114)]. Second, notice that the first term in Eq. (193) (the one proportional to C) leads to dipolar radiation, which scales as (v/c) −2 relative to the GR quadrupole term. Third, notice that in Eq. (193) we have included the order of uncontrolled remainders explicitly, which can be read off from Eq. (116). In particular, the dipole term is corrected by standard 1PN terms (proportional to v 2 ), as well as terms proportional to V 2 CM ; the latter can be safely assumed to be small because the preferred frame is to be identified with that in which the CMB is isotropic. The dipole term also possesses uncontrolled remainders proportional to the vV CM product, which scale linearly with the difference in sensitivities, as opposed to quadratically as the leading-order dipole term thus.
The rate of orbital decay is dominated by the term with the least powers of velocity, or equivalently the least powers of Gm/P b , provided the binary is widely separated. Indeed, this is the case for all observed binary pulsars, with Gm/P b = O(10 −10 ) a typical value for a binary with a 1-hour orbital period. Clearly then, dipolar radiation dominates the orbital decay rate, unless s 1 − s 2 ≈ 0. On the other hand, the sensitivities of a binary system will be similar to each other if their EoSs and their masses are similar. This is the case for the double pulsar PSR J0737-3039 [60], for which the quadrupole term [the second term in Eq. (193)] becomes comparable to the dipole (the first term in this equation). Of course, the uncontrolled remainders proportional to V CM that correct the dipole term at 1PN order are irrelevant here, as they are also multiplied by the difference in the sensitivities. As we will see, the inclusion of this system allows for better constraints on Lorentz-violating theories.
Let us now concentrate on mapping the constraints onα 1 andα 2 to constraints on (c + , c − ) and (λ, β) from observations of PSRs J1012+5307 [120] and J1738+0333 [61]. The mapping can be obtained by replacing α 1 →α 1 , α 2 →α 2 , s 1 → 0 and s 2 → 0 in Eq. (41) [86,117]. By equating the latter to Eq. (41) without these replacements, we can read off the relation between the strong-field PPN coefficients and the standard PPN coefficients and sensitivities. Notice that this is done here for the first time, since previous work (e.g. [43][44][45]117]) did not calculate the sensitivities. This procedure is completely generic and does not depend on the parameters of the orbit. Doing so in Einstein-AEther theory, we find while in khronometric theory we find where we recall that the c i and (α, β) are coupling constants in Einstein-AEther and khronometric theory. Observe that both of these quantities depend not only on the weak-field values of the PPN parameters (α 1 , α 2 ), but also on the sensitivities of the NSs, which depend on the masses; this is why we use observations of PSRs J1012+5307 [120] and J1738+0333 [61], since we know the NS masses for these systems. Saturating the Solar System constraints on the weak-field parameters (α 1 , α 2 ), and using the strong-field constraints on (α 1 ,α 2 ), one can place constraints on (c − , c + ) and (λ, β).

Einstein-AEther Theory
Let us now concentrate on Einstein-AEther theory. The predicted rate of orbital decay is given by Eq. (193) with (C, A 1 , A 2 , A 3 ) functions of the coupling constants (c + , c − ), given in Eqs. (110)-(112), with Eqs. (107), (108), (100) and (95). The predicted strong-field generalization of the PPN parameters (α AE 1 ,α AE 2 ) are given by Eqs. (194) and (195) in terms of their weak-field versions in Eqs. (42) and (43), as well as the sensitivities. Comparing the prediction of the orbital decay rate to observations of PSR J1141-6545 [58], PSR J0348+0432 [59] and PSR J0737-3039 [60], we can then place constraints on the (c + , c − ) coupling parameter space of this theory. Comparing the prediction of (α AE 1 ,α AE 2 ) to the constraint on these quantities derived from observations of PSR J1738+0333 [61], we can also place constraints on this parameter space. We do not consider constraints derived from the orbital decay rate of PSR J1738+0333 [61] and from the (α AE 1 ,α AE 2 ) constraint of PSR J1012+5307 [120] because they are weaker than constraints derived with the systems described above. Figure 7 shows these constraints. The area below and to the right of the solid black line is the allowed (c + , c − ) region after imposing Cherenkov and stability constraints [49]. The green-shaded area with red dotted borders is the allowed (c + , c − ) region after imposing constraints from observations of PSR J1141-6545 [58] (top, left panel), PSR J0348+0432 [59] (top, right panel), PSR J0737-3039 [60] (bottom left panel) and PSR J1738+0333 [61] (bottom right panel). As discussed earlier, Eq. (191) must be evaluated varying over (c + , c − ), but also allowing for the inferred error in the masses, the observational error in the orbital decay rate, the allowed range on (α AE 1 , α AE 2 ) given Solar System constraints, and various NS sensitivities associated with different EoSs. Similarly, Eqs. (194) and (195) must be evaluated by trying different values of α AE 1 and α AE 2 within the Solar System constraints, as well as different values of the sensitivities for different EoSs. The result, thus, is not a line in (c + , c − ) space, but rather a two-dimensional region, as shown in the figures.
Observe that the allowed (c + , c − ) regions in Fig. 7 associated with PSR J1141-6545 [58] and PSR J0348+0432 [59] are rather similar, while they are both quite different from that associated with PSR J0737-3039 [59]. The first two systems are NS-white dwarf binaries, and since white dwarfs have much smaller binding energies than NSs, s NS − s WD ∼ s NS = 0, which implies that the orbital decay rate is dominated by dipolar radiation, i.e. by the first term in Eq. (193). On the other hand, PSR J0737-3039 [60] is a double pulsar binary with similar masses and thus the difference between the sensitivities is small, making the dipolar and quadrupolar terms in the orbital decay rate comparable, i.e. making the first and second terms in Eq. (193) of the same order. Observe also that constraints on (c + , c − ) derived from bounds on (α AE 1 ,α AE 2 ) in the bottom right panel of Fig. 7 also possess a characteristically different shape. This is because these constraints derive from comparing the effect of the conservative sector of the theory to observations, i.e. the Hamiltonian, instead of the effect of the dissipative sector, i.e. the radiation-reaction force. For comparison, Fig. 7 also shows the values of (c + , c − ) (dashed black curve) that would be required for the orbital decay rate to be identical to that of GR in the limit of zero sensitivity [56]. Weakly-gravitating objects have small sensitivities, as the latter scales with the binding energy, which is why we label this curve "weak-field". This curve was obtained in Ref. [56] by requiring that α AE 1 and α AE 2 be identically zero, which results in vanishing sensitivities in the weak-field limit [c.f. Eq. (184)]. Observe that the strong-field binary pulsar constraints derived here overlap with this weak-field constraint only in a small (c + , c − ) region. Observe also that, although the slope of this weak-field curve is similar to the average slope of the allowed region with PSR J0737-3039 [60] (bottom left panel), it disagrees with the slope of the other allowed regions. This is because setting the sensitivities to zero automatically discards the dipolar term in the orbital decay rate, which in fact dominates over the quadrupolar one for mixed binary pulsar systems.
Reference [56] argues that the width of the allowed region given binary pulsar constraints should be approximately O(0.1%), but this is not the case for the systems studied here. The main difference between our work and that of Ref. [56] is that the latter only considered constraints from PSR 1913+16 [96][97][98], which consists of two NSs, in the weak-field/vanishing sensitivity limit, while we here considered other more relativistic double pulsar and mixed pulsar binaries in the strong-field/nonvanishing sensitivity limit. Moreover, the width of the allowed region is not controlled just by the observational error in the orbital decay rate, but also by the inferred error in the individual masses, the observational error in the orbital decay rate, the allowed range of (α AE 1 , α AE 2 ) given Solar System constraints, and also the different NS sensitivities that depend on the EoS, all of which depend strongly on the system that is observed. As we discussed earlier, one must vary over all such parameters to obtain the correct allowed region, which as shown in Fig. 7 is between O(0.1%) and O(1%).
Reference [51] argues that current tests will be satisfied if the weak field conditions are imposed and the coupling constants satisfy (c + , c − ) ≤ O(10 −2 ). For any given system, we find that the allowed region is not a box with sides of O(10 −2 ) but rather a band with width O(10 −3 ). However, when we combine all of these allowed regions in the left panel of Fig. 1, their intersection is indeed almost a box with width of O(10 −3 ) in the c − direction and O(10 −2 ) in the c + direction. The change in shape of the combined allowed region is because different observations lead to allowed bands with different average slopes, which in turn is because for some systems dipolar radiation dominates, while for others it does not. The width of the box is determined by the smallest width of The first three panels use observations of the orbital decay rate through Eq. (191), allowing for inferred error in the masses, observational error in the orbital period, the allowed range of the PPN parameters (α AE 1 , α AE 2 ) given Solar System constraints, and different sensitivities given different EoSs. The last panel uses constraints on the strong-field PPN parameters and Eqs. (194) and (195). For comparison, we also plot the allowed (c+, c−) region after imposing Cherenkov and stability constraints [49] (below and to the right of the solid black line), as well as the values of (c+, c−) required for the orbital decay rate to agree exactly with the GR prediction in the zero-sensitivity limit. Observe that the new, strong-field constraints are much more stringent than the Cherenkov/stability bounds. the bands, which is of O(10 −3 ). All together, these constraints constitute the first, accurate strong-field test of Einstein-AEther theory.

Khronometric Theory
Let us now concentrate on khronometric theory. The predicted rate of orbital decay is given by Eq. (193) with (C, A 1 , A 2 , A 3 ) functions of the coupling constants (λ, β), given in Eqs. (121)- (122). The predicted strong-field generalization of the PPN parameters (α kh 1 ,α kh 2 ) is given by Eqs. (196) and (197) in terms of their weak-field versions in Eqs. (48) and (49), as well as in terms of the sensitivities. Comparing the prediction of the orbital decay rate to observations of PSR J1141-6545 [58], PSR J0348+0432 [59] and PSR J0737-3039 [60], we can then place constraints on the (λ, β) coupling parameter space of this theory. Comparing the prediction of (α kh 1 ,α kh 2 ) to the constraint on these quantities derived from observations of PSR J1738+0333 [61], we also place constraints on this parameter space. Figure 8 shows these constraints. The area above and to the left of the solid black line is the allowed (λ, β) region after imposing Cherenkov and stability constraints. The green-shaded area with red dotted borders is the allowed (λ, β) region after imposing constraints from observations of PSR J1141-6545 [58] (top, left panel), PSR J0348+0432 [59] (top, right panel), PSR J0737-3039 [60] (bottom left panel) and PSR J1738+0333 [61] (bottom right panel). As discussed earlier, Eq. (191) must be evaluated varying over (λ, β), but also allowing for the inferred error in the masses, the observational error in the orbital decay rate, the allowed range on (α kh 1 , α kh 2 ) given Solar System constraints, and various NS sensitivities associated with different EoSs. Similarly, Eqs. (196) and (197) must be evaluated by choosing different values of α kh 1 and α kh 2 within the Solar System constraints, as well as different values of the sensitivities for different EoSs. As in the Einstein-AEther case, the result is not a line in the (λ, β) plane, but rather a two-dimensional region.
Observe that, as in the Einstein-AEther case, the allowed (λ, β) regions associated with PSR J1141-6545 [58] and PSR J0348+0432 [59] are similar to each other, but different from that associated with PSR J0737-3039 [59]. As before, this is because the first two systems are mixed binaries (NS-white dwarf), and thus, s NS −s WD ∼ s NS = 0, since s WD s NS . In turn, this implies that the orbital decay rate is dominated by dipolar radiation, i.e. by the first term in Eq. (193). PSR J0737-3039 [60] is a double pulsar binary, i.e. two NSs with similar masses, and thus the difference in sensitivities is small. This makes dipolar and quadrupolar terms in the orbital decay rate comparable, and thus, changes the shape of the allowed coupling region. The allowed region in the bottom right panel of this figure is also different from the other three, because the former derives from constraints on (α kh 1 , hatα kh 2 ), and thus from the conservative sector of the theory, instead of the dissipative one. Observe also that the first three binary systems leave λ unconstrained in the β → 0 limit. This is to be expected since in this limit the constant dominating the dipole radiation, C, vanishes when α kh 1 = 0 = α kh 2 . Furthermore, the constants in the quadrupole contribution also approach their GR value, A 1 → 1, A 2 → 0, A 3 → 0 in this limit, and the only effect comes from the sensitivities which are too small to produce a strong constraint. When we combine all the panels of Fig. 8, one obtains the right panel of Fig. 1. Observe that binary pulsar observations further constrain the parameter region previously allowed by stability and BBN considerations.

VIII. CONCLUSIONS AND FUTURE DIRECTIONS
We have investigated Lorentz-violating gravity (focusing in particular on violations of boost invariance) in the light of binary pulsar observations. We focused on Einstein-AEther theory and khronometric theory, which are generic theories that break boost invariance at low energies, and showed that they predict emission of dipolar radiation from binary systems, as first reported in Ref. [38,51,56]. This greatly modifies the orbital evolution of these systems, and in particular the decay rate of their orbital period, as compared to the GR prediction. Furthermore, the emission of quadrupolar radiation is also modified. All these modifications, however, depend on the NS sensitivities, which measure how the binding energy of a NS responds to its motion relative to the additional fields of the theory (the AEther and khronon). Thus, constraints on Lorentz-violating gravity can only be placed once the sensitivities have been computed.
Here, we showed in two independent ways that the sensitivities can be extracted, without loss of generality, from solutions describing slowly-moving NSs to first order in velocity, and we then found these solutions numerically for the first time. We began with a generic ansatz for the metric and for the AEther and khronon fields. We then expanded the modified field equations in tensor spherical harmonics, and succeeded at reducing the field equations to a system of ordinary differential equations, which we solved numerically. These slowly-moving solutions then allowed us to extract the NS sensitivities, for which we also present fits in terms of the coupling constants and NS compactness. We checked that the sensitivities approach the weak-field expression derived by Foster [51] as one decreases the NS compactness, and showed that the two differ rather significantly in the large compactness regime.
With these solutions at hand, we then computed the predictions of the orbital decay rate in these theories, compared them to observations of PSR J1141-6545 [58], PSR J0348+0432 [59], and PSR J0737-3039 [60], and obtained very stringent bounds on the coupling constants of Einstein-AEther and khronometric theory. We also complemented these constraints by computing the bounds on the coupling parameters induced by constraints on the precession of PSR J1738+0333 [61]. All together, the constraints derived here on the full set of parameters that control gravitational Lorentz violation are much stronger than previously-obtained bounds.
Future work could extend the results obtained here in various directions. One possibility would be to verify the results of this paper by carrying a Bayesian hypothesistesting or model-selection analysis [125,126]. Such an analysis would take into account covariances between the system parameters and the coupling constants of the theories, which were here neglected. However, since the modifications to the orbital decay rate enter at lower PN order than the leading-order GR prediction, we expect such covariances to have a small effect.
Another possible avenue for future work is to study other binary pulsar observables, such as rate of advance of perihelion [110]. In order to compute this observable, one would have to first compute the acceleration of each component of the binary system, which depends not only on the sensitivities, but also on their first derivative [51]. To calculate the latter, one would have to obtain moving NS solutions to higher order in the slow-velocity expansion. The framework we developed here should suffice to obtain such solutions. For example, the modified field equations should still separate if one carries out a tensor spherical harmonic decomposition. (This can be seen e.g. by extending the procedure of Appendix A to next order in the velocity.) Yet another possibility would be to investigate the The first three panels use observations of the orbital decay rate through Eq. (191), allowing for inferred error in the masses, observational error in the orbital period, the allowed range of the PPN parameters (α kh 1 , α kh 2 ) given Solar System constraints, and different sensitivities given different EoSs. The last panel uses constraints on the strong-field PPN parameters and Eqs. (194) and (195). For comparison, we also plot the allowed (λ, β) region after imposing Cherenkov and stability constraints (above and to the left of the solid black line). Observe how stringent the new, strong-field constraints are.
constraints that future GW observations could place on Einstein-AEther and khronometric theory, given a future detection. GW observations should be particularly sensitive to the frequency-evolution of the GW phase, which will be modified in these theories due to dipolar emission. Moreover, additional polarizations would be excited in binary NS inspirals, beyond the two transverse-traceless ones of GR, which could lead to even stronger constraints if a signal is detected by multiple interferometers [127].
Finally, it would be interesting to compute the sensitivities in Lorentz-violating gravity for isolated black holes [53,69,70,[81][82][83]. Once this is accomplished, one could perform an analysis similar to the one carried out here to place constraints on Lorentz violation in gravity from the orbital decay rates of low-mass X-ray binaries [128][129][130]. Indeed, the orbital decay rates of such systems have already been used to place constraints on certain modified gravity theories [128,[131][132][133][134]. More-over, black hole binaries are expected to be a primary target for future GW detectors. Constraints on the coupling constants of Einstein-AEther and khronometric theory derived from such observations would require knowledge of the sensitivity parameters. The latter could be obtained by following the analysis carried out here, but generalized to black hole solutions.
sub-award 00001944. NY and KY would like to thank, respectively, the Institute d'Astrophysique de Paris and the Yukawa Institute for Theoretical Physics for their hospitality, while some of this work was being carried out. EB acknowledges support from the European Union's Seventh Framework Programme (FP7/PEOPLE-2011-CIG) through the Marie Curie Career Integration Grant GALFORMBHS PCIG11-GA-2012-321608. Some calculations used the computer algebra-systems MAPLE, in combination with the GRTENSORII package [135].

Appendix A: Separability of Field Equations
In Secs. IV and V, we have shown explicitly that the field equations of Einstein-AEther and khronometric theory for an isolated non-spinning star moving slowly relative to the AEther reduce to a system of ordinary differential equations (ODEs), if the "perturbation potentials" [i.e. the functions of r and θ that describe, at O(v), the effect of the motion relative to the AEther on the system] are decomposed in Legendre polynomials. This is non-trivial and does not simply follow from the fact that the Legendre polynomials are a basis for the space of the regular functions of the polar angle θ. In principle, one may indeed imagine that the equations for the various multipole moments may not decouple, leaving one with a hierarchy of infinite coupled ODEs. In fact, the reason why such a situation does not happen and the equations for the various multipole moments decouple lies in the symmetries of the "unperturbed" [i.e. O(v 0 )] solution, which we assumed to be static and spherically symmetric. In GR, it has long been known that whenever the background has those symmetries, the equations that govern (linear) perturbations decouple and therefore reduce to a (finite) system of ODEs [136,137]. As far as we are aware, a general proof of this fact is not yet available in Einstein-AEther or khronometric theory, but we will present another route to this result in the particular case considered here. Relative to the "brute-force" decomposition of Secs. IV and V, this approach has the advantage of showing more explicitly the role of the symmetries of the background.
We begin by noting that the spacetime outside a spherically symmetric static star can be described, in cylindrical isotropic coordinates, by the following metric with r = ρ 2 + z 2 . Note that the use of isotropic coordinates (rather than the areal coordinates used elsewhere in this paper) has the advantage that the spatial part of the metric is conformally flat. The AEther field is instead simply Consider now a star slowly moving in the z-direction relative to the AEther. Based on the transformation properties of the metric and AEther under a coordinate change z → −z and t → −t, it is clear that only the g tz and g tρ components of the metric and the U t , U z and U ρ components of the AEther can be affected. The normalization condition for the AEther immediately requires that the perturbation to U t be zero, because δ(U µ U µ ) = 2g µν U µ δU ν = 2δU t / f (r) = 0. One is then left with the perturbations δg tz , δg tρ , δU z and δU ρ . We only have two vectors at our disposal to construct the perturbations to these fields, namely n = (ρ, z)/r (clearly, |n| = 1) and v = (0, v), so we can write Clearly, the functions S, V, Q, W can only depend on r, because in order to introduce a dependence on ρ and z singularly one would have to use the vector v one more time, i.e. such a dependence only appears at O(v 2 ). The most generic ansatz describing the system at linear order in the velocity is therefore ds 2 = f (r)dt 2 − h(r)(dρ 2 + ρ 2 dθ 2 + dz 2 ) + 2v S(r) + V (r) z 2 r 2 dzdt + 2vV (r) while that for the AEther is Q(ρ, z) = Q(r) + W (r) z 2 r 2 ,W (ρ, z) = W (r) zρ r 2 . (A7) Because this decomposition is completely general and simply based on the symmetries of the problem, and because the metric and AEther now only depend on potentials S, V, Q, W that are functions of r only, the field equations must reduce to a system of differential equations for S, V, Q, W . (Clearly, if they did not, it would mean that the above general decomposition is inconsistent with the field equations, which would in turn mean that slowly moving stars do not exist in the theory under consideration). In fact, inserting this ansatz into the field equations for Einstein-AEther theory, one finds that the non-trivial equations follow from the r and θ components of the AEther equation and from the (t, r) and (t, θ) components of the Einstein equations, and have the following structure where v = (S(r), S (r), S (r), V (r), V (r), V (r), Q(r), Q (r), Q (r), W (r), W (r), W (r)); (A12) v = (S(r), S (r), S (r), V (r), V (r), Q(r), Q (r), Q (r), W (r), W (r)); (A13) w = (S(r), S (r), V (r), Q(r), Q (r), W (r)); (A14) and where the functions A i (r), B i (r), C i (r), D i (r) depend on the coupling constants as well as on the O(v 0 ) potentials f (r) and h(r). Their explicit form is not particularly illuminating and we will not write it down here, but for our purposes it is sufficient to mention that they are functions of r only, so the system of Eqs. (A8)-(A11) does indeed become a system of differential equations in the radial coordinate, as expected.
To show that these differential equations are indeed ordinary, i.e. that they can be put in the form dy/dr = F (y), where y is a suitable array of variables, a more detailed understanding of the gauge degrees of freedom is required. As already shown in Sec. IV A, it is clear that one can set one of the potentials S, V, Q, W exactly to zero in Einstein-AEther theory with a gauge choice.
[This can be checked explicitly by looking at the transformation properties of the metric and AEther under a change of coordinates t = t + vH(r)z.] Therefore, one of the four equations in Eqs. (A8)-(A11) must be algebraically related to the other three, i.e. one should have just three independent equations for the three physical degrees of freedom. In fact, it is possible to show that the radial derivative de 3 /dr can be expressed as a linear combination of e 1 , e 2 , e 3 and e 4 , once the O(v 0 ) equations are satisfied. This fact can be checked explicitly, but follows elegantly from the existence of a generalized Bianchi identity in Einstein-AEther theory [138] (see also Refs. [53,70]), which relates the Einstein equations to the AEther equations [Eq. (6)] and is given by As a result, one can set e.g. W = 0 with a gauge transformation, and rewrite the system e 1 = e 2 = e 4 = 0 as an ODE system S (r) = S 2 (S, S , Q, Q , V, V ) , (A16) Q (r) = Q 2 (S, S , Q, Q , V, V ) , (A17) V (r) = V 2 (S, S , Q, Q , V, V ) , which can be evolved in the radial coordinate (starting e.g. from an initial radius r outwards). The additional equation e 3 = 0 is an initial value constraint for this evolution system, i.e. e 3 = 0 only involves initial data for the system (A16)-(A18) [as can be seen from the derivative structure of Eq. (A10)], and by virtue of the generalized Bianchi identity (A15), it is satisfied everywhere if imposed at the initial radius r (c.f. Refs. [53,70] for more details). Similarly, in the case of khronometric theory, the only non-trivial equations are the (t, r) and (t, θ) components of the Einstein equations (the equation following from the variation of the khronon field is actually implied by the Einstein equations, and does not need to be imposed explicitly [74]). These equations take the form where F i and H i are again functions of r only [depending on the couplings and on the background potentials f (r) and h(r)]. Also, the requirement that the AEther be hypersurface orthogonal simply imposes a relation W (r) = rQ (r) between the AEther potentials at O(v). Therefore, the field equations reduce to a system of differential equations in the radial isotropic coordinate r, exactly as expected. More precisely, as we have shown in Sec. V, a gauge transformation allows to set both AEther perturbations Q and W to zero. (This can be checked explicitly by considering a gauge transformation t = t + vH(r)z.) By doing so, the system takes the form S (r) = S 2 (S, S , V, V ) , (A21) V (r) = V 2 (S, S , V, V ) , which is indeed a set of ODEs.
The function Eq. (186) used to fit the numerical results for the NS sensitivities in Einstein-AEther theory contained only 27 terms, which thus led to an r 2 value of approximately 0.96. In this Appendix, we repeat the fit but with a function that contains 120 terms, thus yielding an r 2 value closer to 0.999. The function that we fit is s AE = We experimented with different fitting functions (e.g. polynomials of different orders in c + , c − and C * ) and found that Eq. (C1) leads to the highest r 2 value with less than 130 coefficients. After the fit, the latter are