Post-Newtonian binary dynamics in the effective field theory of Horndeski gravity

General relativity (GR) has been the most successful theory of gravitation since its proposal by Einstein. It has been tested by various precision experiments and provides an essential theoretical framework for cosmology. However, there are still several cosmological puzzles that GR cannot provide satisfactory explanations for, including the observed small cosmological constant, dark energy, the nature of dark matter, the beginning of our universe, etc. Partly motivated by these observational puzzles and theoretical issues, various modified gravity theories have been proposed and offered as alternative frameworks for the currently well-accepted ΛCDM paradigm.

One of the simplest ways to modify GR is to add a new degree of freedom [1−6], for instance, the Jordan-Brans-Dicke scalar-tensor theory [7], k-essence theory [8−10], Galileon theory [11−13], gauge theories of gravity [14−17], Starobinsky model [18] and its extension with Weyl symmetry [19−21], etc. To be theoretically consistent, such models should satisfy well-definedness conditions. For example, the equation of motion should be of second order in derivatives to avoid instabilities or so-called "ghosts" [22]. With these considerations, Horndeski gravity stands out as a well-formulated and highly general theory [23−25]. It is constructed by simply adding one scalar degree of freedom to GR. Therefore, it serves as a good starting point for exploring the parameter space of modified gravity theories [25−30].

Since 2015 the observation of gravitational waves (GW) from compact binaries has provided a new platform for testing gravitation theories [31]. For example, the merger of neutron stars GW170817 has shown that the speed of a GW is very close to the speed of light [32]. In general, compact binaries are canonical objects for gravitating systems. Besides the GW signals, the emission of electromagnetic pulse signals has provided useful information about their orbital dynamics as well, and such binary pulsars can be used to test gravity theories [33, 34]. For our purpose, we take the precession of the orbital motion to test Horndeski gravity due to the relatively high precision. We compare our theoretical predictions with the values observed for two binary pulsar systems, PSR B 1534+12 [35] and PSR J0737-3039 [36]. Formulating the post-Newtonian (PN) dynamics of a binary system in Horndeski gravity, we provide a self-contained approach to constrain the parameter space of such a theory.

Concretely, we investigate the dynamics of an inspiral compact binary system in Horndeski gravity. We extend a PN effective field theoretical (EFT) approach in GR [37−39] in the perturbative regime. By taking the orbital velocity as an expansion parameter, we calculate the leading contributions and integrate out the potential mode to get an effective action for the binary system. This formulation gives us a systematic way to analytically evaluate the effective dynamics of any order of orbital velocity. In particular, we can include the contribution from the self-interaction of the scalar field and give a bound for the coupling constant from the observation of binary pulsars. This approach can be systematically extended to include higher-order dynamics.

This paper is organized as follows. In section II, we introduce the formalism of Horndeski gravity for a gravitating compact binary system and build up the full Lagrangian. Then, in section III, we introduce the post-Newtonian EFT approach for calculating the binary dynamics and give the detailed Feynman rules from the expansion. In section IV, we obtain the orbital precession in conservative dynamics of a binary system at first PN (1PN) order by calculating a series of the Feynman diagrams needed for the 2-body dynamics in Horndeski gravity. In section V, we compare our results with the observations of binary pulsars and constrain the model's parameters. Finally, we summarize the results and give our conclusion.

Throughtout our discussion, we use mostly-minus signature $(+,-,-,-)$ for the metric $\eta_{\mu\nu}$.

Horndeski theory is the class of modified theories of gravity that introduce an extra scalar degree of freedom ϕ while having a second-order equation of motion for the fields. Thus, it avoids the Ostrogradski instability caused by higher derivatives. Its action has the general form

where the spacetime integration is defined as $\int _x\equiv \int \mathrm{d}^4 x \sqrt{-g}$, g is the determinant of the metric field $g_{\mu\nu}$, $S_{\rm EH}=-2 m_p^2 \int_x R$ is the Einstein-Hilbert action of the Ricci scalar R determined by $g_{\mu\nu}$, $m^{-2}_{p} = 32\pi G$, G is related to Newton's constant, $G^{\mu\nu}$ is the Einstein tensor, and $G_n(\phi,X)$ are functions of ϕ and X with notations

We impose that the speed of the GW propagation is the same as the speed of light, then only the $G_2(\phi,X)$, $G_3(\phi,X)$ , and $G_4(\phi)$ terms remain [25]. To calculate the leading PN effective Lagrangian, we only need to expand the full Lagrangian for 3-point vertices. In addition, the linear ϕ or X terms in $G_2$ can be shifted away by redefinition of ϕ. Furthermore, we consider massless scalars so that $\phi^2$ terms are not included. The 3-point vertex $g_2 \phi X$ term can be combined with the kinematic term of X and redefine $\sqrt{1+g_2 \phi}\,d\phi \to d\Phi$ to canonically normalize the kinetic term. Therefore, the leading term in $G_2(\phi,X)$ , for our interest, is $\phi^3$. For a detailed discussion about the removal of the $g_2 \phi X$ term, see Appendix A.

Let us consider the following action

In particular, we consider the following case to illustrate the dynamics

where $g_2, g_3$, and $g_4$ are normalized dimensionless constants that are to be determined or constrained by observations. Note that $g_3$ is parameterized to be dimensionless by introducing the inverse of Planck mass which is large compared to the energy scales of the typical process in a binary pulsar system. This suggests that the viable $g_3$ can be a big number, still allowed by the bound of binary pulsars, as we shall see in later discussions.

In addition, we should add the gauge-fixing (GF) term associated with the gauge invariant action $S_{\rm EH}$, $S_{\rm GF}$, which in harmonic gauge is given by

where $\Gamma^{\mu}_{\; \rho\sigma}$ is the connection determined by the metric $g_{\mu\nu}$.

As we only consider the inspiral phase of the compact binary in which the stars are moving much slower than light, we can treat the binary system classically. Each star of the binary system is regarded as a point particle with mass $m_a$, located at $x^{\mu}_a(t)$, and interacting with the gravitational field which is described by a worldline Lagrangian $S_{pp}$:

To summarize, the full action is

We can observe that the scalar field does not couple with the worldline Lagrangian (that is, matter) directly, but only indirectly through gravity, which shall simplify the calculations, unlike previous investigations on scalar-tensor theories in which the scalar field couples directly to matter [40−43].

To work out Horndeski theory for a compact binary system in the inspiral phase, we use an effective field theoretical approach based on PN expansion [37], which is often also called the non-relativistic GR (NRGR) method. In the inspiral phase of a compact binary system, the orbital velocity v is in a non-relativistic regime and much lower than the speed of light. Hence, it is meaningful to separate the scales of energy with respect to the post-Newtonian expansion parameter v. In particular, we have 3 relevant physical scales in the system: the Schwarzschild radius of the system, $r_s = 2 G m$, the radius of the orbit, r, and the wavelength of the gravitational wave, λ.

In a non-relativistic regime, the system satisfies the virial theorem $v^2 \sim Gm/r$, thus we have $r_s/r \sim v^2 \ll 1$. While from $\lambda \sim r/v$, we have $r/\lambda \sim v \ll 1$. So these scales are separated in the non-relativistic limit $v \ll 1$. The physics at the scale $r_s$ is only relevant in strong field regime and contained in the effective mass of point particles. The main physics in the system is carried by gravitons with 4 wave-vector of order $(\omega, |\mathbf{k}|) \sim (v/r,1/r)$ in the potential region and of order $(\omega, |\mathbf{k}|) \sim (v/r, v/r)$ in the radiation region. Furthermore, by integrating the potential gravitons, we will get the effective potential and the effective coupling of the radiation of the systems.

To perform the calculations in effective field theory, we have to expand a field in different modes for different energy scales. For the Horndeski theory, apart from the metric field whose expansion yields two tensor modes of gravitons, we will have an extra scalar field. Both the scalar and tensor degrees of freedom have their potential and radiation modes. Expanding the metric tensor and scalar fields in different regions, we have

where $H_{\mu\nu}$ and Φ are referred to as radiation modes, while $h_{\mu\nu}$ and $\phi_{\rm pot}$ are referred to as potential modes. In the following discussion, we will use ϕ without subscription to refer to the potential part of the scalar field in order to avoid confusion.

In the gravitating binary system, our goal is to calculate the effective Lagrangian for the two stars and to give observable quantities. We consider the amplitude of two classically moving point particles scattering by exchanging potential modes. After integrating these potential modes, the amplitude will yield the effective Lagrangian of the classical motion of the 2 binary stars. Thus we can use the agreement of the calculated quantities of the effective motion with the observed values to determine a bound for our model parameters.

To perform the post-Newtonian EFT calculations, we need to derive the Feynman rules from an expansion of the full Lagrangian and collect the relevant Feynman diagrams for a binary system. To calculate the effective Lagrangian of the compact stars up to 1PN order in an EFT approach, we derive the propagators of tensor and scalar modes, and relevant three-point vertices between these modes and the worldline. Here we only show the propagators, and leave the vertices in Appendix B.

To obtain the propagators, we expand the action up to $\phi^2, \phi h$ and $h^2$ order, and turn to momentum space representation,

where $\int_p$ stands for integration over momentum space $\int {\rm d}^4 p/{(2\pi)^4}$. Then the relevant action can be written as follows

where we have defined the following quantities

The inverse of the kinetic terms then gives the tensor and scalar propagators as displayed in Fig. 1,

Figure 1.  Tensor and scalar propagators. The solid line represents the scalar mode, and the dashed line represents the tensor mode.

where we can calculate

Then we can read the propagator for the tensor modes,

The scalar-tensor mixing term is given by the off-diagonal term

For off-shell potential modes, $p=(\omega, \mathbf{p}),\; \omega \sim |\mathbf{p}| v$. Hence at the leading order, we can estimate that

For the off-diagonal term, the components of $q_{\mu\nu}$ are

at the leading order. To perform the integration in potential modes, we turn to three-dimensional Fourier-transformed configuration space in which the propagators are written as

Note that we could diagonalize this propagator matrix, which then will induce coupling between matter and scalar field. However, in the basis without diagonalization, the contributing diagrams would be simpler due to the absence of direct coupling between matter and scalar field, although the two bases are equivalent.

As we need the classical dynamics of the system, we have to extract the pure classical part (leading in $\hbar$) from the EFT calculations. In other words, this means that we should limit the exchanged graviton modes to the potential modes. Such contributions involve the diagrams satisfying 2 criteria:

1. Every loop must include an edge on one of the worldlines, otherwise, the integration of loop momenta at a large energy scale will introduce radiative contributions.

2. No internal edge must start and end on the same worldline, otherwise there will be internal energy contribution for the particle on the worldline.

With these restrictions, the number of relevant diagrams can be considerably reduced.

When integrating internal momenta in loop diagrams, we have to regularize and renormalize the divergences in the divergent integrals. We choose dimensional regularization and minimal subtraction scheme in our EFT approach [37]. So the dimensionless integrals vanish and can be dropped. For 1PN calculations, the most general integrals can be written as follows

We set the dimension $d=3$ and drop the scaleless integrals. Among integrals with different $n_i$, there are linear relations generated from the integration of full derivatives or integration by part. Using these relations, all the integrals are reduced to $f(1,1,0,0,0;r)$ multiplied by a factor. After integration, $f(1,1,0,0,0;r)$ yields simply

For example, we have

In this section, we discuss the classical binary dynamics in Horndeski gravity up to 1PN and show the leading contributions in diagrammatic pictures. By calculating the scattering amplitude, we can obtain the corresponding Lagrangian and Hamiltonian, which determine the classical dynamics.

For classical binary dynamics, at the Newtonian order, we have only one diagram, as shown in Fig. 3. The amplitude of the diagram in Fig. 3 is calculated as

Figure 3.  Leading-order diagram for binary dynamics.

At the first post-Newtonian (1PN) order, we have the following diagrams to calculate. The diagrams in Fig. 4 have the same topology as the leading order one but with higher order $O(v^2)$ vertices, thus they scale as $O(Gv^2)$. We calculate the amplitude as

Figure 4.  Next-to-leading diagrams at tree level with tensor modes only. (a)

The diagrams in Fig. 5 are one-loop ones with leading-order vertices and scale as $O(G^2)$. We obtain

Figure 5.  Next-to-leading diagrams at one-loop level with tensor modes only. (b)

Because of the $\phi^3$ vertex (spatial derivatives are allowed), we also have the diagram shown in Fig. 6. As we are working in the base with $\phi-h$ mixing-type propagators, only h gravitons are interacting with matter worldline directly. Therefore, this diagram does not vanish when the mixing is present. We obtain

Figure 6.  The next-to-leading diagram at one-loop level with scalar self-interaction. (c)

This contribution is proportional to $1/r^4$ and comes from the scalar self-coupling with the derivatives $g_3 X\square\phi$. Note that the presence of the Planck scale $l_{p}^2$ does not mean it is a quantum correction. This term appears because our parameterization of $g_3$ in Eq. (3) introduces such a length scale for such a coupling. Moreover, the coupling $g_2\phi^3$ does not contribute at this leading order, since it is proportional to the integral

This is scaling as $p^0$ , hence the result is a constant in r and does not affect the dynamics.

Additionaly, other contributing diagrams are shown in Fig. 7 and Fig. 8. Their results are respectively

Figure 7.  Next-to-leading diagrams at one-loop level with graviton-scalar mixing. (d)

Figure 8.  Next-to-leading diagrams at one-loop level with graviton-scalar mixing. (e)

and

In the amplitudes we have calculated, we should make the substitution $G(1+\kappa) \to G$, so that the Newtonian-order effective Lagrangian is the same as Einstein's effective Lagrangian for two particles. The effective Lagrangian is extracted from the matching relation

Exploiting the matching relation at each PN order, we have

where we have defined the following quantities

After a Legendre transformation of the 1PN two-body effective Lagrangian, we can get the Hamiltonian. In the meantime, the two-body motion can be reduced to the motion of a reduced mass $\mu = m_1m_2/(m_1+m_2)$. After rescaling the following quantities

we have the Hamiltonian for the motion of the reduced mass

where

According to Damour and Schafer [44], we can use Hamilton-Jacobi Formalism to derive the periastron advance $\Delta \Phi$ per period from the reduced Hamiltonian, $k \equiv {\Delta \Phi}/({2\pi})$,

where

Rewriting these with observable parameters, eccentricity e, orbital period $P_b$, and inferred total mass M, and considering the following relations

we obtain the periastron advance in the discussed Horndeski theory

Now we can use the observation data from several binary systems to constrain the model's parameters. Here, we choose two typical pulsar binary systems, PSR B 1534+12 [35] and PSR J0737-3039 A/B [36], which have a rather high observation precision of periastron advance. The observed data for the two systems are shown in Table 1. We use the observed values of periastron advance and their error to obtain a bound for physical values in the $g_3-g_4$ plane. In principle, the total radiation power is another observable that can test the 1PN dynamics of the model. However, the precision is much lower than that of periastron advance, therefore here we only use information from the latter.

Quantity	PSR B 1534+12	PSR J0737-3039A/B
Period	$P_b = 0.420737299122(10) \,d$	$P_b = 0.1022515592973(10) \,d$
Eccentricity	$e = 0.2736775(3)$	$e = 0.087777023(61)$
Pulsar Mass	$m_1 = 1.3332(10) \,M_{\odot}$	$m_A = 1.338185(^{+12}_{-14}) \,M_{\odot}$
Companion Mass	$m_2 = 1.3452(10) \,M_{\odot}$	$m_B = 1.248868(^{+13}_{-11}) \,M_{\odot}$
Total Mass	$M = 2.678428(18) \,M_{\odot}$	$M = 2.587052(^{+9}_{-7}) \,M_{\odot}$
Precession Rate	$\dot{\omega} = 1.755789(9)\rm \;deg/yr$	$\dot{\omega} = 16.899323(13) \rm\; deg/yr$

Table 1.  Observables fitted from pulsar timing data of PSR B 1534+12 [35] and PSR J0737-3039 [36].

From the observed pulsar timing data of the system, people usually extract a set of Kepler parameters from the data

to represent the dynamics at the leading order of relativistic expansion[45]. Here ω is the angle between the line from the orbital center to the point of the periastron and the nodal line, $x=a\sin i$ is the apparent semi-axis, and $T_0$ is the observation time. The higher-order dynamics which is relevant in our model can be expressed as the rate of these experimental parameters. This so-called post-Keplerian includes the periastron advance $\dot{\omega}_0$ as well as the changing rate of the orbital period $\dot{T}_0$, which is related to the radiated power and other parameters. we will make use of the observation data of $\dot{\omega}_0$, which has a high precision, to constrain our model parameters..

We implement the constraints from the observed values for fractional periastron advance per period k:

where $k_{\rm model}$ is given by Eq. (39), and should lie within $2\sigma$ around the observed value $k_{\rm obs}$.

The bound of $g_3$ and $g_4$ is showed in Fig. 9(a) and Fig. 9(b) for the two systems, respectively. Here we use $\Sigma = \sqrt{\sigma_{\rm model}^2 + \sigma_{\rm obs}^2}$ to represent the total standard error. As the observed values can be greater or smaller than the predicted ones, we see the dotted lines have different shapes in these two systems. The reason is that in Fig. 9(a) we have $k_{\rm model,GR} < k_{\rm obs}$, while in Fig. 9(b) we have $k_{\rm model,GR} > k_{\rm obs}$. Here $k_{\rm model,GR}$ is the predicted value if there are no corrections for general relativity. Despite this difference, the predicted and observed values in the two examples both lie within the $2\sigma$ range. Note that the large value of $g_3$ is again due to our parametrization of the $g_3$ term in Eq. (3).

Figure 9.  Constraint of $g_3$ and $g_4$ for the two systems. The solid lines represent the trajectory of $(g_4,g_3)$ values which give $k_{\rm model}=k_{\rm obs}-2\Sigma$, the dotted lines represent the trajectory of $k_{\rm model}=k_{\rm obs}$, and the dashed lines represent the trajectory of $k_{\rm model}=k_{\rm obs}+2\Sigma$.

We have investigated the conservative dynamics of a binary system up to the first-order post-Newtonian in Horndeski gravity by formulating its effective field theory. From the effective Lagrangian or Hamiltonian, we have computed the periastron advance of a binary system in its inspiral phase. Comparing the theoretical predictions of periastron advance for two binary systems with their observations, PSR B 1534+12 and PSR J0737-3039, we have obtained the constraints on the model parameters: $g_3$ and $g_4$ in Eq. (3). The contours of the available parameter space are shown in Fig. 9. As the precision of these observations will improve in the future, we expect the bound of these parameters to become more restrictive. Also, the formalism developed here can also be systematically extended to higher orders by summing higher-loop diagrams and applying it to other binary systems, including binary black holes. Effects on the gravitational waves of such an inspiral system would also be interesting [46−50], which we shall explore for future work.

The $g_2 \phi X$ term can be removed by field redefinition $\phi \to \Phi$ with

Here we use the notation $X_A=\frac12 g^{\mu\nu}\nabla_{\mu} A \nabla_{\nu} A$ and unit $m_p = 1$. In differential form, we demand that

To integrate the equation, we set the zero point of Φ so that $\Phi=0$ when $\phi=0$, then

Assuming that this correction to kinematic term is small $g_2 \phi \ll 1$,

Then the correction to $G_3$ is of higher order

Similarly to the derivation of graviton propagators, we expand the worldline action $S_{pp}$ to find the couplings of gravitons and the point mass source at each order. To calculate an effective Lagrangian up to 1PN order, we need worldline vertices with at most two gravitons and up to the $O(v^2)$ order. We expand the Lagrangian to get

The Lagrangian of particle worldlines in Eq. (B1) are expanded up to the $O(v^2)$ terms. Terms without $h_{\mu\nu}$ contribute to the kinetic energy of the binary. At 1PN order, this yields the kinetic energy

Up to the $v^2 h^2$ order, we have the worldline-graviton vertices as in Fig. 2, these vertices are summarized in Table A1.

Figure 2.  Worldline-graviton vertices. Double lines represent the particle worldline while the dashed lines represent (tensor mode) gravitons.

	$O(v^0)$	$O(v^1)$	$O(v^2)$
$h^1$	$\displaystyle h_{00}:\; -{\rm i}\dfrac{m_a}{2m_p} {\rm e}^{{\rm i} \mathbf{k}\cdot \mathbf{x}_a(t)}$	$\displaystyle h_{0i}:\; -{\rm i}\dfrac{m_a v_a^i}{m_p} {\rm e}^{ {\rm i} \mathbf{k}\cdot \mathbf{x}_a(t)}$	$\begin{array}{*{20}{l} }\displaystyle h_{00}: \; -{\rm i}\dfrac{m_a v_a^2}{4m_p} {\rm e}^{ {\rm i} \mathbf{k}\cdot \mathbf{x}_a(t)} \\\displaystyle h_{ij}: \; -{\rm i}\dfrac{m_a v_a^i v_a^j}{2m_p} {\rm e}^{ {\rm i} \mathbf{k}\cdot \mathbf{x}_a(t)}\end{array}$
$h^2$	$\displaystyle h_{00}^2:\; {\rm i}\dfrac{m_a}{4m_p^2}{\rm e}^{ {\rm i}(\mathbf{k_1}+\mathbf{k}_2)\cdot \mathbf{x}_a(t)}$	$\displaystyle h_{00}h_{0i}:\; {\rm i}\dfrac{m_a v_a^i}{2m_p^2} {\rm e}^{ {\rm i} (\mathbf{k_1} + \mathbf{k_2})\cdot \mathbf{x}_a(t)}$	$\begin{array}{*{20}{l} }\displaystyle h_{00}^2:\; {\rm i}\dfrac{3m_a v_a^2}{8m_p^2} {\rm e}^{ {\rm i} (\mathbf{k_1} + \mathbf{k_2})\cdot \mathbf{x}_a(t)} \\\displaystyle h_{0i}h_{0j}:\; {\rm i}\dfrac{m_a v_a^i v_a^j}{m_p^2} {\rm e}^{ {\rm i} (\mathbf{k_1} + \mathbf{k_2})\cdot \mathbf{x}_a(t)} \\\displaystyle h_{00}h_{ij}:\; {\rm i}\dfrac{m_a v_a^i v_a^j}{4m_p^2} {\rm e}^{ {\rm i} (\mathbf{k_1} + \mathbf{k_2})\cdot \mathbf{x}_a(t)} \end{array}$

Table A1.  Worldline-graviton couplings up to $O(h^2)$, $O(v^2)$.

For the three-point graviton self-interaction vertices, we need to expand the action up to the $O(h^3)$ order. The action is then Fourier transformed and symmetrized in $p_1,p_2, p_3$:

where h without indices stands for traces $h(p)=h^{\; \mu}_{\mu}(p)$ and

Similarly, we have symmetrized the 3-point vertices in the expanded action at $O(\phi^3)$. For instance,