Constraints on Horndeski Theory Using the Observations of Nordtvedt Effect, Shapiro Time Delay and Binary Pulsars

Alternative theories of gravity not only modify the polarization contents of the gravitational wave, but also affect the motions of the stars and the energy radiated away via the gravitational radiation. These aspects leave imprints in the observational data, which enables the test of General Relativity and its alternatives. In this work, the Nordtvedt effect and the Shapiro time delay are calculated in order to constrain Horndeski theory using the observations of lunar laser ranging experiments and Cassini time-delay data. The effective stress-energy tensor is also obtained using the method of Isaacson. Gravitational wave radiation of a binary system is calculated, and the change of the period of a binary system is deduced for the elliptical orbit. These results can be used to set constraints on Horndeski theory with the observations of binary systems, such as PSR J1738+0333. Constraints have been obtained for some subclasses of Horndeski theory, in particular, those satisfying the gravitational wave speed limits from GW170817 and GRB 170817A.


I. INTRODUCTION
General Relativity (GR) is one of the cornerstones of modern physics. However, it faces several challenges. For example, GR cannot be quantized, and it cannot explain the present accelerating expansion of universe, i.e., the problem of dark energy. These challenges motivate the pursuit of the alternatives to GR, one of which is the scalar-tensor theory. The scalar-tensor theory contains a scalar field φ as well as a metric tensor g µν to describe the gravity. It is the simplest alternative metric theory of gravity. It solves some of GR's problems. For example, the extra degree of freedom of the scalar field might account for the dark energy and explain the accelerating expansion of the universe. Certain scalar-tensor theories can be viewed as the low energy limit of string theory, one of the candidates of quantum gravity [1].
The detection of gravitational waves by the Laser Interferometer Gravitational-Wave Observatory (LIGO) and Virgo confirms GR to an unprecedented precision [2][3][4][5][6][7] and also provides the possibility to test GR in the dynamical, strong field limit. The recent GW170814 detected the polarizations for the first time, and the result showed that the pure tensor polarizations are favored against pure vector and pure scalar polarizations [5]. The newest GW170817 is the first neutron star-neutron star merger event, and the concomitant gamma-ray burst GRB 170817A was later observed by the Fermi Gammaray Burst Monitor and the Anti-Coincidence Shield for the Spectrometer for the International Gamma-Ray Astrophysics Laboratory, independently [6,8,9]. This opens the new era of multi-messenger astrophysics. It is thus interesting to study gravitational waves in alternative metric theories of gravity, especially the scalar-tensor theory. * shou1397@hust.edu.cn † yggong@hust.edu.cn In 1974, Horndeski [10] constructed the most general scalar-tensor theory whose action contains higher derivatives of φ and g µν , but still yields at most the second order differential field equations, and thus has no Ostrogradsky instability [11]. Because of its generality, Horndeski theory includes several important specific theories, such as GR, Brans-Dicke theory [12], and f (R) gravity [13][14][15] etc..
In Refs. [16][17][18], we discussed the gravitational wave solutions in f (R) gravity and Horndeski theory, and their polarization contents. These works showed that in addition to the familiar + and × polarizations in GR, there is a mixed state of the transverse breathing and longitudinal polarizations both excited by a massive scalar field, while a massless scalar field excites the transverse breathing polarization only. In this work, it will be shown that the presence of a dynamical scalar field also changes the amount of energy radiated away by the gravitational wave affecting, for example, the inspiral of binary systems. Gravitational radiation causes the damping of the energy of the binary system, leading to the change in the orbital period. In fact, the first indirect evidence for the existence of gravitational waves is the decay of the orbital period of the Hulse-Taylor pulsar (PSR 1913+16) [19].
Previously, the effective stress energy tensor was obtained by Nutku [20] using the method of Landau and Lifshitz [21]. The damping of a compact binary system due to gravitational radiation in Brans-Dicke theory was calculated in Refs. [22][23][24][25], then Alsing et al. [26] extended the analysis to the massive scalar-tensor theory. Refs. [27,28] surveyed the effective stress-energy tensor for a wide class of alternative theories of gravity using several methods. However, they did not consider Horndeski theory. Refs. [29,30] studied the gravitational radiation in screened modified gravity and f (R) gravity. Hohman [31] developed parameterized post-Newtonian (PPN) formalism for Horndeski theory. In this work, the method of Isaacson is used to obtain the effective stress-energy tensor for Horndeski theory. Then the effective stress-energy tensor is applied to calculate the rate of energy damping and the period change of a binary system, which can be compared with the observations on binary systems to constrain Horndeski theory. Nordtvedt effect and Shapiro time delay effect will also be considered to put further constraints. Ashtekar and Bonga pointed out in Refs. [32,33] a subtle difference between the transverse-traceless part of h µν defined by ∂ ν h µν = 0, η µν h µν = 0 and the one defined by using the spatial transverse projector, but this difference does not affect the energy flux calculated in this work.
There were constraints on Horndeski theory and its subclasses in the past. The observations of GW170817 and GRB 170817A put severe constraints on the speed of gravitational waves [34]. Using this limit, Ref. [35] required that ∂G 5 /∂X = 0 and 2∂G 4 /∂X + ∂G 5 /∂φ = 0, while Ref. [36] required ∂G 4 /∂X ≈ 0 and G 5 ≈ constant. Ref. [37] obtained the similar results as Ref. [36], and also pointed out that the self-accelerating theories should be shift symmetric. Arai and Nishizawa found that Horndeski theory with arbitrary functions G 4 and G 5 needs fine-tuning to account for the cosmic accelerating expansion [38]. For more constraints derived from the gravitational wave speed limit, please refer to Refs. [39][40][41], and for more discussions on the constraints on the subclasses of Horndeski theory, please refer to Refs. [42][43][44][45][46].
In this work, the calculation will be done in the Jordan frame, and the screening mechanisms, such as the chameleon [47,48] and the symmetron [49,50], are not considered. Vainshtein mechanism was first discovered to solve the vDVZ discontinuity problem for massive gravity [51], and later found to also appear in theories containing the derivative self-couplings of the scalar field, such as some subclasses of Horndeski theory [52][53][54][55][56]. When Vainshtein mechanism is in effect, the effect of nonlinearity cannot be ignored within the so-called Vainshtein radius r V from the center of the matter source. Well beyond r V , the linearization can be applied. The radius r V depends on the parameters defining Horndeski theory, and can be much smaller than the size of a celestial object. So in this work, we consider Horndeski theories which predict small r V , if it exists, compared to the sizes of the Sun and neutron stars. The linearization can thus be done even deep inside the stars. In this case, one can safely ignore Vainshtein mechanism.
The paper is organized as follows. In Section II, Horndeski theory is briefly introduced and the equations of motion are derived up to the second order in perturbations around the flat spacetime background. Section III derives the effective stress-energy tensor according to the procedure given by Isaacson. Section IV is devoted to the computation of the metric and scalar perturbations in the near zone up to Newtonian order and the discussion of the motion of self-gravitating objects that source gravitational waves. In particular, Nordtvedt effect and Shapiro time delay are discussed. In Section V, the metric and scalar perturbations are calculated in the far zone up to the quadratic order, and in Section VI, these solutions are applied to a compact binary system to calculate the energy emission rate and the period change. Section VII discusses the constraints on Horndeski theory based on the observations. Finally, Section VIII summarizes the results. Throughout the paper, the speed of light in vacuum is taken to be c = 1.

II. HORNDESKI THEORY
The action of Horndeski theory is given by [57], where ψ m represents matter fields, S m is the action for ψ m , and the terms in the integrand are In these expressions, ;ρ for simplicity. G i (i = 2, 3, 4, 5) are arbitrary functions of φ and X [58]. For notational simplicity and clarity, we define the following symbol for the function f (φ, X), so in particular, f (0,0) = f (φ 0 , 0) with φ 0 the value of φ in the flat spacetime background. Suitable choices of G i reproduce interesting subclasses of Horndeski theory. For instance, one obtains GR by choosing G 4 = (16πG N ) −1 and the remaining G i = 0, with G N Newton's constant. Brans-Dicke theory is recovered with G 2 = 2ω BD X/φ, G 4 = φ, G 3 = G 5 = 0, while the massive scalar-tensor theory with a potential U (φ) [26] is obtained with

A. Matter action
Although there are no coupling terms between matter fields ψ m and φ, matter fields ψ m indirectly interact with φ via the metric tensor. For example, in Brans-Dicke theory, φ acts effectively like the gravitational constant, which influences the internal structure and motion of a gravitating object, so the binding energy of the object depends on φ. Since the total energy E is related to the inertial mass m, then m depends on φ, too. When their spins and multipole moments can be ignored, the gravitating objects can be described by point like particles, and the effect of φ can be taken into account by the following matter action according to Eardley's prescription [59], whose stress-energy tensor is where x λ a (τ ) describes the worldline of particle a and u µ = dx µ (τ )/dτ . Therefore, if there is no force other than gravity acting on a self-gravitating object, this object will not follow the geodesic. This causes the violation of the strong equivalence principle (SEP).
In this work, the gravitational wave is studied in the flat spacetime background with g µν = η µν and φ = φ 0 , so we expand the masses around the value φ 0 in the following way, Here, ϕ = φ − φ 0 is the perturbation, and m a = m a (φ 0 ) for simplicity. This expansion also requires that φ 0 = 0, so the present discussion does not apply to f (R) gravity. s a and s a are the first and second sensitivities of the mass m a , The sensitivities measure the violation of SEP.

B. Linearized equations of motion
The equations of motion can be obtained and simplified using xAct package [60][61][62][63][64]. Because of their tremendous complexity, the full equations of motion will not be presented. Interested readers are referred to Refs. [57,65]. As we checked, xAct package gives the same equations of motion as Refs. [57,65]. For the purpose of this work, the equations of motion are expanded up to the second order in perturbations defined as These equations are given in A.
The gravitational wave solutions are investigated in the flat spacetime background, which requires that This can be easily checked by a quick inspection of Eqs.
(A1) and (A2). Then dropping higher order terms in Eqs. (A1) and (A2), the linearized equations of motion are thus given by where T = g µν T µν is the trace, 2 = η µν ∂ µ ∂ ν from now on, and the superscript (1) implies the leading order part of the quantity. The equations of motion can be decoupled by introducing an auxiliary fieldh µν defined as following, where h = η µν h µν is the trace, and the original metric tensor perturbation is, withh = η µνh µν . The equations of motion are gauge invariant under the the following infinitesimal coordinate transformation, with x µ = x µ + ξ µ . Therefore, one can choose the transverse gauge ∂ νh µν = 0, and after some algebraic manipulations, the equations of motion become where T µν , and the mass of the scalar field is Of course, ζ = 0, otherwise ϕ is non-dynamical. From the equations of motion (17) and (18)), one concludes that the scalar field is generally massive unless G 2(2,0) is zero, and the auxiliary fieldh µν resembles the spin-2 graviton fieldh µν = h µν − η µν h/2 in GR.h µν is sourced by the matter stress-energy tensor, while the source of the scalar perturbation ϕ is a linear combination of the trace of the matter stress-energy tensor and the partial derivative of the trace with respect to φ. This is because of the indirect interaction between the scalar field and the matter field via the metric tensor.

III. EFFECTIVE STRESS-ENERGY TENSOR
The method of Isaacson [67,68] will be used to obtain the effective stress-energy tensor for gravitational waves in Horndeski theory in the short-wavelength approximation, i.e., the wavelength λ 1/ √ R with R representing the typical value of the background Riemann tensor components. This approximation is trivially satisfied in our case, as the background is flat and R = 0. In averaging over several wavelengths, the following rules are utilized [69]: 1. The average of a gradient is zero, e.g., where implies averaging. These rules apply to not only terms involvingh but also those involving ϕ. In the case of a curved background, these rules are supplemented by the one that covariant derivatives commute, which always holds in the flat background case. With this method, the effective stress-energy tensor in an arbitrary gauge can be calculated straightforwardly using xAct and given by, It can be checked that this expression is gauge invariant under Eq. (16). In fact, the terms in the first around brackets take exactly the same forms as in GR excerpt for a different factor. The fourth line remains invariant, as ϕ = ϕ in the gauge transformation. To show that the remaining lines are also gauge invariant, making the Far away from the matter, ∂ ρ ∂ ρ ϕ = m 2 s ϕ according to Eq. (17). Substituting this into the fourth line of Eq. (21), one immediately finds total derivatives of the forms ∂ µ (ϕ∂ ρ ∂ ρ ξ ν ) and ∂ σ (ϕ∂ ρ ∂ ρ ξ σ ). So the first averaging rule implies that the last three lines of Eq. (21) vanish. Therefore, the effective stress-energy tensor (20) is indeed gauge invariant.
In vacuum, the transverse-traceless (TT) gauge (∂ νh µν = 0 andh = 0) can be taken, and the effective stress-energy simplifies, whereh TT µν denotes the transverse-traceless part. In the limit that G 4 = (16πG N ) −1 and the remaining arbitrary functions G i vanish, Eq. (20) recovers the effective stress-energy tensor of GR [69]. One can also check that Eq. (20) reduces to the one given in Ref. [25] for Brans-Dicke theory in the gauge of ∂ νh µν = 0 andh = −2ϕ/φ 0 . In order to calculate the energy carried away by gravitational waves, one has to first study the motion of the source. This is the topic of the next section.

IV. THE MOTION OF GRAVITATING OBJECTS IN THE NEWTONIAN LIMIT
The motion of the source will be calculated in the Newtonian limit. The source is modeled as a collection of gravitating objects with the action given by Eq. (6). In the slow motion, weak field limit, there exists a nearly global inertial reference frame. In this frame, a Cartesian coordinate system is established whose origin is chosen to be the center of mass of the matter source. Let x represent the field point whose length is denoted by r = | x|.
In the near zone [70], the metric and the scalar perturbations will be calculated at the Newtonian order. The stress-energy tensor of the matter source is given by [71], and one obtains, In these expressions, the 4-velocity of particle a is u µ a = u 0 a (1, v a ) and v 2 a = v 2 a . With these results, the leading order of the source for the scalar field is with S a = G 4(1,0) − 2G 4(0,0) φ0 s a . Now, the linearized equations (17,18) take the following forms (27) and the leading order contributions to the perturbations are easily obtained, andh 0j =h jk = 0 at this order, where r a = | x − x a | and the scalar field is given by a sum of Yukawa potentials. The leading order metric perturbation can be determined by Eq. (15), with h 0j = 0.

A. Static, spherically symmetric solutions
For the static, spherically symmetric solution with a single point mass M at rest at the origin as the source, the time-time component of the metric tensor is where S M = G 4(1,0) − 2G 4(0,0) s M /φ 0 and s M is the sensitivity of the point mass M . From this, the "Newton's constant" can be read off which actually depends on the distance r because the scalar field is massive. The measured Newtonian constant at the earth is G N (r ⊗ ) with r ⊗ the radius of the Earth. The "post-Newtonian parameter" γ(r) can also be read off by examining g jk , which is In the PPN formalism, the space-space components of the metric take the following form, where the parameter γ is a constant. So The above result can recover the results for f (R) gravity and general scalar-tensor theory [31,[72][73][74] if we keep the equivalence principle. In the massless case (G 2(2,0) = 0), we get Note that G N (r) and γ(r) both depend on S M which reflects the internal structure and motion of the gravitating object in question. Even if the scalar field is massless, this dependence still persists. Therefore, neither of them is universal due to the violation of SEP caused by the scalar field. It is obvious that G N (r ⊗ ) should take the same value as G N .

B. Equations of motion of the matter
With the near zone solutions (28), (30) and (31) one obtains the total matter Lagrangian up to the linear order, where r ab = | x a − x b | is the distance between the particles a and b. The equation of motion for the mass m a can thus be obtained using the Euler-Lagrange equation, yielding its acceleration, withr ab = ( x a − x b )/r ab . In particular, for a binary system, the relative acceleration a j = a j 1 − a j 2 is where m = m 1 + m 2 is the total mass. The first term in the square brackets gives the result that resembles the familiar Newtonian gravitational acceleration, while the second one reflects the effect of the scalar field. In the massless case, the second term no longer depends on r 12 and can be absorbed into the first one, so the binary system moves in a similar way as in Newtonian gravity with a modified Newton's constant.
The Hamiltonian of the matter is where p j a = ∂L m /∂x j a is the j-th component of the canonical momentum of particle a, and the total rest mass has been dropped. In particular, the Hamiltonian of a binary system is given by where v = v 1 − v 2 , and µ = m 1 m 2 /m is the reduced mass. This will be useful for calculating the total mechanical energy of a binary system and the ratio of energy loss due to the gravitational radiation.

C. Nordtvedt effect
The presence of the scalar field modifies the trajectories of self-gravitating bodies. They will no longer follow geodesics. Therefore, SEP is violated in Horndeski theory. This effect is called the Nordtvedt effect [75,76]. It results in measurable effects in the solar system, one of which is the polarization of the Moon's orbit around the Earth [77,78].
To study the Nordtvedt effect, one considers a system of three self-gravitating objects a, b and c and studies the relative acceleration of a and b in the field of c. With Eq. (40) and assuming r ab r ac ≈ r bc , the relative acceleration is where the first term presents the Newtonian acceleration modified by the presence of the scalar field, the second is the tidal force caused by the gravitational gradient due to the object c, and the last one describes the Nordtvedt effect. The effective Nordtvedt parameter is This parameter depends on S c = G 4(1,0) − 2G 4(0,0) s c /φ 0 , so this effect is indeed caused by the violation of SEP.

D. Shapiro time delay effect
Another effect useful for constraining Horndeski theory is the Shapiro time delay [79]. In order to calculate this effect, one considers the photon propagation time in a static (or nearly static) gravitational field produced by a single mass M at the origin. Due to the presence of gravitational potential, the 3-velocity of the photon in the nearly inertial coordinate system is no longer 1 and varies. The propagation time is thus different from that when the spacetime is flat. Let the 4 velocity of the photon be u µ = u 0 (1, v), then u µ u µ = 0 gives where h 00 and h jk are given by Eqs. (30) and (31) specialized to a single mass M case. In the flat spacetime, the trajectory for a photon emitted from position x e at time t e is a straight line x(t) = x e +N (t − t e ), whereN is the direction of the photon. The presence of the gravitational potential introduces a small perturbation δ x(t) so that x(t) = x e +N (t − t e ) + δ x(t). Substituting Eqs. (30) and (31) into Eq. (46), one obtainŝ where r(t) = | x(t)|. Suppose the photon emitted from position x e is bounced back at position x p and finally returns to x e . The total propagation time is where δt is caused by the Shapiro time delay effect, where r e = | x e |, r p = | x p | and r b = |N × x e | is the impact parameter of the photon relative to the source. Since M in Eq. (49) is not measurable, one replaces it with the Keplerian mass with S M = G 4(1,0) − 2G 4(0,0) s M /φ 0 and s M the sensitivity of the source. In terms of M K , the Shapiro time delay is (51) For the Shapiro time delay occurring near the Sun, r in the above equation should be 1 AU, as this is approximately the distance where the Keplerian mass M K of the Sun is measured.

V. GRAVITATIONAL WAVE SOLUTIONS
In the far zone, only the space-space components of the metric perturbation are needed to calculate the effective stress-energy tensor. Since the equation of motion (18) forh µν takes the similar form as in GR, the leading order contribution toh jk is given by, where I jk = a m a x j a x k a is the mass quadrupole moment. As in GR, the TT part ofh jk is also related to the reduced quadrupole moment J jk = I jk − δ jk δ il I il /3, The leading order term for the scalar field ϕ is the mass monopole which does not contribute to the effective stress-energy tensor, so it is necessary to take higher order terms into account. To do so, the scalar equation (A2) is rewritten with the linearized equations substi-tuted in, which is given by In the following discussion, it is assumed that the scalar field is massless for simplicity. The details to obtain the following results can be found in B. The leading order contribution to ϕ comes from the first term on the right hand side of Eq. (54), which is the mass monopole moment, From now on, the superscript [n] indicates the order of a quantity in terms of the speed v, i.e., ϕ [n] is of the order O(v 2n ). ϕ [1] is independent of time, so it does not contribute to the effective stress-energy tensor. The next leading order term is the mass dipole moment, in whichn = x/r. This gives the leading contribution to the effective stress-energy tensor. At the next next leading order, there are more contributions from the remaining terms on the right hand side of Eq. (54). First, there is the mass quadruple moment contribution, And the remaining contribution to the scalar wave is where a,b means summation over a and b with a = b, Note that the penultimate line of Eq. (58) is a sum of terms proportional to r ab , which grows as r ab increases and potentially dominates over other terms. Since matters are confined within the source zone, this line never blows up. The scalar field up to the fourth order in velocity is given by It is easy to check that this result agrees with Eq. (86) in Ref. [26] with m s = 0.

VI. GRAVITATIONAL RADIATION FOR A COMPACT BINARY SYSTEM
This section is devoted to calculating the gravitational radiation for a compact binary system in the case with massless scalar field . According to Eq. (22), the energy carried away by the gravitational wave is at a rate oḟ where the integration is carried out on a 2-sphere in the far zone and in the final step, higher order terms have been dropped. The first term gives the contribution of the spin-2 gravitational wave, while the second one gives the contribution of the scalar field. Next, one has to calculate the motion of the binary system explicitly. By Eq. (41), the relative acceleration is given by where As in GR, one can orient the coordinate system such that the orbit lies in the xOy plane. In the polar coordinate system (r, θ, z), the relative distance is thus given by where with l the angular momentum per unit mass and e the eccentricity. The orbital period is All these above results can be obtained by suitably modifying those in GR as found in Ref. [70]. Using Eq. (43) with m s set to 0, the total mechanical energy of the binary system is where a = p/(1 − e 2 ) is the semi-major axis. Following Ref. [26], the rate of energy loss due to the spin-2 gravitational wave iṡ (67) which reproduces the radiation damping of GR in the appropriate limit [70].
Ignoring the leading order contribution to ϕ, the higher order correction is given by where and The first term at the right hand side of Eq. (68) is a dipolar contribution and oscillates at the orbital frequency. This term is of order v −1 1 relative to the remaining terms. However, it also depends on the difference in the sensitivities (s 1 − s 2 ) of the objects in the binary system, which might be small or even vanish. For example, in the Shift-Symmetric Horndeski theory (SSHT) with G i functions of X only, the stellar sensitivity s a vanishes [80], and in Brans-Dicke theory, the sensitivity of a black hole is 1/2 [26,81,82]. So if the binary system consists of, e.g., two neutron stars in SSHT or if the two stars are black holes in Brans-Dicke theory, the dipolar radiation vanishes.
In the generic case, (s 1 − s 2 ) might not be zero, and the dipolar contribution should be taken into account. So the contribution of the scalar field to the energy flux A straightforward but tedious calculation shows that Eq. (76) reduces to Eq. (3.24) in Ref. [25] for Brans-Dicke theory with sensitivities set to zero and the Hadamard regularization imposed [83][84][85]. The period changeṪ can be measured experimentally, and the fractional period changeṪ /T is given bẏ The first term is caused by the spin-2 gravitational wave, while the remaining ones by the scalar field.
Given the sensitivities (s a , s a ) of all kinds of celestial objects, Eq. (77) can be compared with the observed period change to set bounds on some of parameters characterizing a particular scalar-tensor theory (e.g., φ 0 , G 4(0,0) , G 4(1,0) , ζ etc.) as done in Ref. [26].

VII. OBSERVATIONAL CONSTRAINTS
In this section, constraints on Horndeski theory are obtained using observations from lunar laser ranging experiments, Cassini time-delay measurement and binary pulsars. Since Horndeski theory contains many parameters, the following discussions start with generic constraints on the full Horndeski theory, and then specify to some concrete subclasses of Horndeski theory.

A. Constraints from lunar laser ranging experiments
The lunar laser ranging experiment gave the most precise measurement of the Nordtvedt effect, and the Nordtvedt parameter was determined to be [86] η obs. N = (0.6 ± 5.2) × 10 −4 = δ 1 ± 1 .
To get the constraints, one requires that |η N − δ 1 | < 2 1 at 95% confidential level. Using Eq. (45), one obtains where r = 1 AU and the sensitivity of the Sun is ignored as its sensitivity is expected to be smaller than 10 −4 , which is the white dwarf's sensitivity [26,82].

B. Constraints from Cassini time-delay data
In 2002, the Cassini spacecraft measured the Shapiro time delay effect in the solar system by radio tracking [87]. The PPN parameter γ was given by At 95% confidential level, one requires that |γ(r) − γ meas. | < 2 2 , which leads to in which the Sun's sensitivity is also ignored, and r = 1 AU. In the massless case, this constraint can be trans- [31], which reduces to ω BD when (the massless) Brans-Dicke theory is considered.

C. Constraints from period change for circular motion
Now, one obtains the constraints on Horndeski theory using the data of pulsars. For this end, one considers the circular motion of a binary system, not only for simplicity but also because the first sensitivities s a are known at least in some subclasses of Horndeski theory, such as Brans-Dicke theory [26,81,82] and SSHT [80], while the second sensitivities s a are unknown. In the case of the circular motion (e = 0), one assumes that ω is the orbital angular frequency so that r 12 = a and θ = ωt. The orbital angular frequency can be obtained using Eq. (64), which is The total mechanical energy of the binary system is The rates of radiation damping are greatly simplified, anḋ where the first term comes from the mass dipole moment. The fractional period change iṡ The first two terms are caused by the scalar field, while the last one by the spin-2 gravitational wave.
Provided that the sensitivities (s 1 , s 2 ) of celestial objects are given, Eq. (86) can be compared with the observed period change to set bounds on some parameters in Horndeski theory, using the observational data of the binary system PSR J1738+0333 [88]. This is a 5.85-ms pulsar with a white dwarf companion, orbiting around each other every 8.51 hours. Some of the orbit parameters are tabulated in Table I. The eccentricity of PSR is the observed period change and σ is the uncertainty forṪ obs. . The expression forṪ pred. −Ṫ obs. is too complicated and will not be presented here.
Note that since the Newton's constant G N is measured in the vicinity of the Earth, the Earth's sensitivity s ⊗ is ignored in Eq. (37), and so ζ does not depend on s ⊗ . Plug ζ into Eq. (81), and the Shapiro time delay effect constrains G 4(0,0) , Plug ζ into Eq. (79), and one gets which shows a nice property that the product χ = φ 0 G 4(1,0) appears in the above expression. In fact, after one substitutes ζ into Eq. (86),Ṫ can also be expressed as a function of G 4(0,0) and χ, which is too complicated to be presented. Note that the sensitivities for the pulsar and the white dwarf are taken to be approximately 0.2 and 10 −4 , respectively. So the constraints from the Nordtvedt effect and the period change of the binary pulsar can be represented by the constraints on G 4(0,0) and χ. The result is given in Fig. 1. The shaded area is the commonly allowed parameter space (G 4(0,0) , χ). Finally, since ζ is given in Eq. (II B), one knows that Note that the above constraints cannot be applied to the special case where G 4 ∝ φ, as in this case, G 4(1,0) ∝ G 4(0,0) /φ 0 , i.e., G 4(1,0) and G 4(0,0) are not independent of each other. Example 2: Now, consider a second subclass of Horndeski theory whose G 4 = G 4 (φ) and G 5 = 0. The scalar field is still assumed to be massless. This subclass satisfies the constraints set by the gravitational wave speed limit [35][36][37]. One can introduce a new scalar field φ such that G 4 (φ) = φ /16π, and the form of action (1) remains the same after replacing φ by φ in it. So let us simply call the new scalar field φ, and thus G 4 (φ) = φ/16π and G 4(1,0) = 1/16π. Using all the constraints discussed in the previous subsections, one obtains and this leads to (92) Example 3: One may also consider the constraints set on a massive Horndeski theory. In this case, one can only use the constraints from the Nordtvedt effect and the Shapiro time delay. The mass m s of the scalar field is expected to be very small. As suggested in Ref. [26], if 10 −21 eV < m s < 10 −15 eV, the constraints can also be set on G 4(0,0) and χ, provided that they are independent of each other. The allowed parameter space (G 4(0,0) , χ) is approximately given by the area enclosed by the two vertical dashed curves, and the dot dashed one in Fig. 1. The constraint on the combination G 2(0,1) − 2G 3(1,0) is also approximately given by Eq. (90). If G 4 ∝ φ, the constraints are approximately given by Eqs. (91) and (92).

VIII. CONCLUSION
In this work, the observational constraints on Horndeski theory are obtained based on the observations from the Nordtvedt effect, Shapiro time delay and binary pulsars. For this purpose, the near zone metric and scalar perturbations are first calculated in order to obtain the equations of motion for the stars. These solutions are thus used to study the Nordtvedt effect and the Shapiro time delay. Then, the effective stress-energy tensor of Horndeski theory is derived using the method of Isaacson. It is then used to calculate the rate of energy radiated away by the gravitational wave and the period change of a binary system. For this end, in the far zone, the auxiliary metric perturbation is calculated using the familiar quadratic formula, and the scalar field is calculated with the monopole moment contribution dominating, although it does not contribute to the effective stressenergy tensor. The leading contribution of the scalar field to the energy damping is the dipolar radiation, which is related to the difference in the sensitivities of the stars in the binary system, so the dipolar radiation vanishes if the two stars have the same sensitivity. The energy damping is finally calculated with the far zone field perturbations, and the period change is derived. Finally, the observational constraints are discussed based on the data from lunar laser ranging experiments, the observations made by the Cassini spacecraft, and the observation on the PSR J1738+0333. Explicit constraints have been obtained for both the massless and massive Horndeski theory, and in particular, for the one satisfying the recent gravitational wave speed limits [6].
in perturbations are where 2 = η µν ∂ µ ∂ ν , and the superscript (1) implies the leading order piece of the quantity while the superscript (2) represents the second order piece.

Appendix B: Post-Newtonian Expansion of the Scalar Field
In this appendix, the procedure to derive the post-Newtonian expansion of the scalar field is briefly presented. The basic idea is the following. Suppose a scalar field ψ satisfies the massless Klein-Gordon equation with a source S, where 2 = ∂ µ ∂ µ . In the far zone, the scalar field is given by Here, the integration is over the near zone N , as ψ will be calculated only up to the quadratic order in perturbations. Since r = | x| > | x |, one can expand the integrand in powers of x in the following way, where u = t − r is the retarded time, Q is a multi-index, namely, ∂ Q = ∂ j1 ∂ j2 · · · ∂ jq and I Q = I j1j2···jq , and the repeated indices imply summation. The symbol I Q (u) is in which the integration is over M, the intersection of the near-zone worldtube with the constant retarded time hypersurface u = C. Since ∂ j u = −x j /r = −n j , Eq. (B3) is approximately given by (B5) For the purpose of the present work, one identifies ψ with ϕ and −16πS with the right hand side of Eq. (54) up to the quadratic order. One should further truncate the series in the above expression at an appropriate order in the following discussion.
The leading contribution to ϕ comes from the first term on the right hand side of Eq. (54), which is the mass monopole moment, This gives the leading contribution to the effective stressenergy tensor.
At the next next leading order, there are more contributions from the right hand side of Eq. (54). First, there is the mass quadruple moment, (B8) The above three contributions (B6), (B7) and (B8) all come from the first term in the source (the right hand side of Eq. (54)).
Other contributions to the scalar quadruple moment come from the remaining terms in the source. Firstly, there are the following three contributions, where a,b means summation over a and b with a = b, and in the second step of Eq. (B10), the contribution from η jkh jk T (1) * is dropped since it is of order O(v 2 ) relative toh 00 T (1) * . Secondly, the term containing T (1) µν ∂ µ ∂ ν ϕ in Eq. (54) does not contribute as where each term on the right hand side in the above expression indicates the relative order of that term to T (1) 00 ϕ, and ρ = T (1) 00 . Note that the action of ∂ 0 increases the order by one since ∂ 0 is actually −∂/c∂t. Therefore, these terms are of higher order than those considered in Eqs. (B9), (B10) and (B11), and will be ignored. Similarly, the term containingh µν ∂ µ ∂ ν ϕ is also of higher order and dropped.
Thirdly, the following integral will be useful, (B13) The next useful integral is To compute it, we first consider the terms with a = b, Remember that R defines the boundary separating the near zone from the far zone. However, the scalar field should not depend on R, as shown in Ref. [89]. So this result will be discarded. Second, consider the contributions from terms with a = b. Define y = r a = x − x a , then r b = x − x b = y + r ab . Since the source is located deep inside the near zone, | x a | R. For x ∈ N , | x| 2 = | y + x a | 2 = y 2 + 2 y · x a + x 2 a < R 2 , and one knows that, where y = | y| andŷ = y/y. So