1 Introduction

Over the years several alternative theories of gravity have emerged in the attempt to solve some of the problems presented by the general theory of relativity, such as the dark matter and dark energy problems. Besides solving these problems, for an alternative theory of gravity to be considered consistent, it must also reproduce the successful predictions of general relativity. One of these recently confirmed predictions is the existence of gravitational waves [1,2,3,4].

Among the many alternative theories of gravity that have already studied the gravitational waves phenomenology is conformal gravity (CG). It was shown that the plane wave of this theory is composed of the usual plane wave of general relativity plus a plane wave that grows linearly in time [5], which causes the energy carried by the CG plane wave to diverges in momentum space [6].

In this paper, we intend to study the behavior of gravitational waves in another alternative theory of gravity with conformal symmetry called massive conformal gravity (MCG) [7]. In Sect. 2, we present an introduction of the MCG theory. In Sect. 3, we find the plane wave solution of the linearized MCG field equations. In Sect. 4, we discuss the energy–momentum tensor of the MCG plane wave. In Sect. 5, we evaluate the radiated energy from a binary system in MCG. Finally, in Sect. 6, we provide a brief conclusion about the results found in the paper.

2 Massive conformal gravity

Let us consider the total MCG actionFootnote 1 [8]

$$\begin{aligned} S_{\text {tot}}= & {} \frac{1}{\kappa ^2}\int {d^{4}x} \, \sqrt{-g}\bigg [\varphi ^{2} R + 6\partial ^{\mu }\varphi \partial _{\mu }\varphi \nonumber \\&- \frac{1}{2m^2} C^{\alpha \beta \mu \nu }C_{\alpha \beta \mu \nu } \bigg ] + \int {d^{4}x}{\mathcal {L}}_{m}, \end{aligned}$$
(1)

where \(\kappa ^{2} = 32\pi G/3\), \(\varphi \) is a scalar field called dilaton, m is a constant with dimension of mass,

$$\begin{aligned} C^{\alpha \beta \mu \nu }C_{\alpha \beta \mu \nu }= & {} R^{\alpha \beta \mu \nu } R_{\alpha \beta \mu \nu } - 4R^{\mu \nu }R_{\mu \nu } + R^2 \nonumber \\&+ 2\left( R^{\mu \nu }R_{\mu \nu } - \frac{1}{3}R^{2}\right) \end{aligned}$$
(2)

is the Weyl tensor squared, \(R^{\alpha }\,\!\!_{\mu \beta \nu } = \partial _{\beta }\varGamma ^{\alpha }_{\mu \nu } + \cdots \) is the Riemann tensor, \(R_{\mu \nu } = R^{\alpha }\,\!\!_{\mu \alpha \nu }\) is the Ricci tensor, \(R = g^{\mu \nu }R_{\mu \nu }\) is the scalar curvature, and \({\mathcal {L}}_{m} = {\mathcal {L}}_{m}(g_{\mu \nu },\varPsi )\) is the Lagrangian density of the matter field \(\varPsi \).

The variation of the total action (1) with respect to \(g^{\mu \nu }\) and \(\varphi \) gives the MCG field equations

$$\begin{aligned}&\varphi ^{2}G_{\mu \nu } + 6 \partial _{\mu }\varphi \partial _{\nu }\varphi - 3g_{\mu \nu }\partial ^{\rho } \varphi \partial _{\rho }\varphi + g_{\mu \nu } \nabla ^{\rho }\nabla _{\rho } \varphi ^{2} \nonumber \\&\quad - \nabla _{\mu }\nabla _{\nu } \varphi ^{2} - m^{-2}W_{\mu \nu } = \frac{1}{2}\kappa ^2 T_{\mu \nu }, \end{aligned}$$
(3)
$$\begin{aligned}&\left( \nabla ^{\rho }\nabla _{\rho } - \frac{1}{6}R\right) \varphi = 0, \end{aligned}$$
(4)

where

$$\begin{aligned} W_{\mu \nu }= & {} \nabla ^{\rho }\nabla _{\rho }R_{\mu \nu } - \frac{1}{3}\nabla _{\mu }\nabla _{\nu }R -\frac{1}{6}g_{\mu \nu }\nabla ^{\rho } \nabla _{\rho }R \nonumber \\&+ 2R^{\rho \sigma }R_{\mu \rho \nu \sigma } -\frac{1}{2}g_{\mu \nu }R^{\rho \sigma }R_{\rho \sigma } \nonumber \\&- \frac{2}{3}RR_{\mu \nu } + \frac{1}{6}g_{\mu \nu }R^2 \end{aligned}$$
(5)

is the Bach tensor,

$$\begin{aligned} G_{\mu \nu } = R_{\mu \nu } - \frac{1}{2}g_{\mu \nu }R \end{aligned}$$
(6)

is the Einstein tensor, and

$$\begin{aligned} T_{\mu \nu } = \frac{2}{\sqrt{-g}} \frac{\delta {\mathcal {L}}_{m}}{\delta g^{\mu \nu }} \end{aligned}$$
(7)

is the matter energy–momentum tensor.

Besides being invariant under coordinate transformations, the field equations (3) and (4) are also invariant under the conformal transformations

$$\begin{aligned} {\tilde{\varPhi }} = \varOmega (x)^{-\varDelta _{\varPhi }}\varPhi , \end{aligned}$$
(8)

where \(\varOmega (x)\) is an arbitrary function of the spacetime coordinates, and \(\varDelta _{\varPhi }\) is the scaling dimension of the field \(\varPhi \), whose values are \(-2\) for the metric field, 0 for gauge bosons, 1 for scalar fields, and 3/2 for fermions. Using the conformal invariance of the theory, we can impose the unitary gauge \(\varphi = \varphi _{0} = 1\). In this case, the field equations (3) and (4) becomes

$$\begin{aligned}&G_{\mu \nu } - m^{-2}W_{\mu \nu } = \frac{1}{2}\kappa ^2 T_{\mu \nu }, \end{aligned}$$
(9)
$$\begin{aligned}&R = 0. \end{aligned}$$
(10)

Taking the trace of (9) and comparing with (10), we find that MCG couples only with matter whose energy–momentum tensor is traceless, which is a feature of conformally invariant theories. This is the case of the standard model of particle physics without a Higgs mass term, for instance. For simplicity, we consider the conformally invariant matter Lagrangian density [9]

$$\begin{aligned} {\mathcal {L}}_{m}= & {} -\sqrt{-g}\Bigg [S^{2}R + 6\partial ^{\mu }S\partial _{\mu }S + \lambda S^{4} \nonumber \\&+ \frac{i}{2}\left( \, {\overline{\psi }} \gamma ^{\mu }D_{\mu }\psi - D_{\mu }{\overline{\psi }}\gamma ^{\mu }\psi \right) + \mu S{\overline{\psi }}\psi \Bigg ], \end{aligned}$$
(11)

where S is a scalar Higgs field, \(\lambda \) and \(\mu \) are dimensionless coupling constants, \({\overline{\psi }} = \psi ^{\dagger } \gamma ^{0}\) is the adjoint fermion field, \(D_{\mu } = \partial _{\mu } + [\gamma ^{\nu },\partial _{\mu }\gamma _{\nu }]/8 - [\gamma ^{\nu },\gamma _{\lambda }] \varGamma ^{\lambda }_{\mu \nu }/8\), and \(\gamma ^{\mu }\) are the general relativistic Dirac matrices, which satisfy the anticommutation relation \(\{\gamma ^{\mu },\gamma ^{\nu }\} = 2g^{\mu \nu }\).

Before proceeding, it is worth noting that both the coordinate and the conformal symmetries of the theory allow the introduction of a quartic self-interaction term of the dilaton field in the gravitational part of the total MCG action (1). The reason why we do not consider such a term is that its inclusion makes the flat metric no longer a solution of the vacuum field equations, which invalidates the usual S-matrix formulation. Additionally, we can include a coupling between the dilaton and the Higgs fields in the matter Lagrangian density (11). However, we neglect this coupling because it leads to a nonvanishing trace of the matter energy–momentum tensor, as we can see by considering the variation of the modified \({\mathcal {L}}_{m}\) with respect to \(\varphi \) on the right side of (4) and comparing the resulting field equation with the trace of (3).

The variation of (11) with respect to \({\overline{\psi }}\) and \(\psi \) gives the field equations

$$\begin{aligned}&i\gamma ^{\mu }D_{\mu }\psi + \mu S \psi = 0, \end{aligned}$$
(12)
$$\begin{aligned}&iD_{\mu }{\overline{\psi }}\gamma ^{\mu } - \mu S {\overline{\psi }} = 0. \end{aligned}$$
(13)

Substituting (11) into (7), and using (12) and (13), we obtain the matter energy–momentum tensor

$$\begin{aligned} T_{\mu \nu }= & {} 2g_{\mu \nu }\nabla ^{\rho }S\nabla _{\rho }S -8\nabla _{\mu }S\nabla _{\nu }S +4 S\nabla _{\mu }\nabla _{\nu } S \nonumber \\&- 4g_{\mu \nu }S\nabla ^{\rho }\nabla _{\rho } S + 2S^{2}G_{\mu \nu } + T^{f}_{\mu \nu } - g_{\mu \nu }\lambda S^4, \end{aligned}$$
(14)

where

$$\begin{aligned} T^{f}_{\mu \nu } = \frac{i}{4}\big (\, {\overline{\psi }} \gamma _{\mu }D_{\nu }\psi - D_{\nu }{\overline{\psi }}\gamma _{\mu }\psi + {\overline{\psi }}\gamma _{\nu }D_{\mu }\psi - D_{\mu }{\overline{\psi }}\gamma _{\nu } \psi \big )\nonumber \\ \end{aligned}$$
(15)

is the fermion energy–momentum tensor.

Considering that, at scales below the electroweak scale, the Higgs field acquires a spontaneously broken constant vacumm expectation value \(S_{0}\), we find that (14) reduces to

$$\begin{aligned} T_{\mu \nu } = 2S_{0}^{2}G_{\mu \nu } + T^{f}_{\mu \nu } - g_{\mu \nu }\lambda S_{0}^4. \end{aligned}$$
(16)

Taking the trace of (16) and substituting into the trace of (9), we obtain

$$\begin{aligned} -R = \frac{1}{2}\left( -2S_{0}^{2}R + T^{f} - 4\lambda S_{0}^4 \right) , \end{aligned}$$
(17)

where \(T^{f} = g^{\mu \nu }T^{f}_{\mu \nu }\). The additional use of (10) then gives the relation

$$\begin{aligned} \lambda S_{0}^4 = \frac{1}{4}T^{f}. \end{aligned}$$
(18)

Finally, substituting this relation back into (16), we arrive at

$$\begin{aligned} T_{\mu \nu } = 2S_{0}^{2}G_{\mu \nu } + T^{T}_{\mu \nu }, \end{aligned}$$
(19)

where

$$\begin{aligned} T^{T}_{\mu \nu } = T^{f}_{\mu \nu } - \frac{1}{4}g_{\mu \nu }T^{f} \end{aligned}$$
(20)

is the traceless part of the fermion energy–momentum tensor.

3 Plane gravitational waves

In order to find the MCG gravitational wave equations, we must perturb the metric according to

$$\begin{aligned} g_{\mu \nu } = \eta _{\mu \nu } + h_{\mu \nu }, \end{aligned}$$
(21)

where \(\eta _{\mu \nu } = \mathrm {diag}( -1, +1, +1, +1)\) and \(|h_{\mu \nu }| \ll 1\). Then, using (19), we find that up to first order in the perturbation \(h_{\mu \nu }\), the field equations (9) and (10) becomeFootnote 2

$$\begin{aligned}&G^{(1)}_{\mu \nu } - m^{-2}W^{(1)}_{\mu \nu } = \frac{1}{2}\kappa ^2 T^{(0)}_{\mu \nu }, \end{aligned}$$
(22)
$$\begin{aligned}&R^{(1)} = 0, \end{aligned}$$
(23)

where

$$\begin{aligned} T^{(0)}_{\mu \nu } = 2S_{0}^{2}G^{(0)}_{\mu \nu } + T^{T(0)}_{\mu \nu } = T^{f(0)}_{\mu \nu } - \frac{1}{4}\eta _{\mu \nu }T^{f(0)} \end{aligned}$$
(24)

is the zero-order matter energy–momentum tensor,

$$\begin{aligned} W^{(1)}_{\mu \nu } = \Box R^{(1)}_{\mu \nu } - \frac{1}{3}\partial _{\mu }\partial _{\nu }R^{(1)} -\frac{1}{6}\eta _{\mu \nu } \Box R^{(1)} \end{aligned}$$
(25)

is the first-order Bach tensor,

$$\begin{aligned} G^{(1)}_{\mu \nu } = R^{(1)}_{\mu \nu } - \frac{1}{2}\eta _{\mu \nu }R^{(1)} \end{aligned}$$
(26)

is the first-order Einstein tensor,

$$\begin{aligned} R^{(1)}_{\mu \nu } = \frac{1}{2} \left( \partial _{\mu }\partial ^{\rho } h_{\rho \nu } + \partial _{\nu }\partial ^{\rho }h_{\rho \mu } - \Box h_{\mu \nu } - \partial _{\mu }\partial _{\nu }h \right) \end{aligned}$$
(27)

is the first-order Ricci tensor, and

$$\begin{aligned} R^{(1)} = \partial ^{\mu }\partial ^{\nu }h_{\mu \nu } - \Box h \end{aligned}$$
(28)

is the first-order scalar curvature, with \(\Box = \partial ^{\mu }\partial _{\mu }\) and \(h = \eta ^{\mu \nu }h_{\mu \nu }\).

It is not difficult to see that both (22) and (23) are invariant under the coordinate gauge transformation

$$\begin{aligned} h'_{\mu \nu } = h_{\mu \nu } + \partial _{\mu }\xi _{\nu } + \partial _{\nu }\xi _{\mu }, \end{aligned}$$
(29)

where \(\xi _{\mu }\) is an arbitrary spacetime dependent vector field. By imposing the gauge condition

$$\begin{aligned} \partial ^{\mu }h_{\mu \nu } - \frac{1}{2}\partial _{\nu }h = 0, \end{aligned}$$
(30)

which fix the coordinate gauge freedom up to residual gauge parameter satisfying the subsidiary condition

$$\begin{aligned} \Box \xi _{\mu } = 0, \end{aligned}$$
(31)

and using (23), we find that (22) reduces to

$$\begin{aligned} \left( m^{-2}\Box - 1 \right) \Box h_{\mu \nu } = \kappa ^2 T^{(0)}_{\mu \nu }. \end{aligned}$$
(32)

The simplest physical solution to (32) in vacuum (\(T^{(0)}_{\mu \nu }=0\)) is the plane wave

$$\begin{aligned} h_{\mu \nu } = a_{\mu \nu }\cos (k_{\rho }x^{\rho }) + b_{\mu \nu }\cos (q_{\rho }x^{\rho }), \end{aligned}$$
(33)

where \(a_{\mu \nu }\) and \(b_{\mu \nu }\) are symmetric wave polarization tensors, and \(k_{\mu }\) and \(q_{\mu }\) are wave vectors, which satisfies

$$\begin{aligned} k^{\rho }k_{\rho } = 0,\quad q^{\rho }q_{\rho } = -m^2. \end{aligned}$$
(34)

By substituting (33) into (23) and (30), and using (34), we obtain

$$\begin{aligned} k^{\mu }a_{\mu \nu } = \frac{1}{2}k_{\nu }a, \quad q^{\mu }b_{\mu \nu } = 0, \quad b = 0, \end{aligned}$$
(35)

where \(a = \eta ^{\mu \nu }a_{\mu \nu }\) and \(b = \eta ^{\mu \nu }b_{\mu \nu }\).

The symmetry of the polarization tensors \(a_{\mu \nu }\) and \(b_{\mu \nu }\) means that each of them has ten independent components. The conditions (35) reduce the independent components of \(a_{\mu \nu }\) to six and of \(b_{\mu \nu }\) to five. In addition, we can choose a solution of (31) to impose four more conditions on \(a_{\mu \nu }\). For instance, choosing

$$\begin{aligned} \xi _{\mu } = \epsilon _{\mu }\sin (k_{\rho }x^{\rho }) , \end{aligned}$$
(36)

and substituting into (29) together with (33), we arrive at

$$\begin{aligned} a'_{\mu \nu } = a_{\mu \nu } + k_{\mu }\epsilon _{\nu } + k_{\nu }\epsilon _{\mu }. \end{aligned}$$
(37)

Since \(\epsilon _{\mu }\) is arbitrary, we can select it to impose four more conditions on \(a_{\mu \nu }\). In particular, we can choose \(\epsilon _{\mu }\) such that

$$\begin{aligned} a_{0i} = 0, \ \ \ \ \ \ \delta ^{ij}a_{ij} = 0, \end{aligned}$$
(38)

which reduce the independent components of \(a_{\mu \nu }\) to just two. Thus, we conclude that the MCG plane gravitational wave (33) has seven propagating degrees of freedom, represented by the two independent components of \(a_{\mu \nu }\) and the five independent components of \(b_{\mu \nu }\), which is consistent with recent results obtained in the literature [8, 10].

4 Gravitational energy–momentum tensor

We can see from (32) that the total linearized MCG Lagrangian density is dynamically equivalent to

$$\begin{aligned} {\mathcal {L}}_{\text {tot}}= & {} -\frac{1}{4\kappa ^2}\left( m^{-2}\Box h^{\mu \nu }\Box h_{\mu \nu } + \partial ^{\rho }h^{\mu \nu }\partial _{\rho }h_{\mu \nu } \right) \nonumber \\&+ \frac{1}{2}h^{\mu \nu }T^{(0)}_{\mu \nu }. \end{aligned}$$
(39)

Inserting the gravitational part of (39) into the canonical energy–momentum tensor [11]

$$\begin{aligned} t^{\mu \nu }= & {} \bigg \langle \left[ \partial _{\rho }\frac{\partial {\mathcal {L}}}{\partial (\partial _{\mu }\partial _{\rho } h_{\alpha \beta })} - \frac{\partial {\mathcal {L}}}{\partial (\partial _{\mu }h_{\alpha \beta })} \right] \partial ^{\nu }h_{\alpha \beta } \nonumber \\&- \frac{\partial {\mathcal {L}}}{\partial (\partial _{\mu }\partial _{\rho } h_{\alpha \beta })}\partial ^{\nu }\partial _{\rho }h_{\alpha \beta } + \eta ^{\mu \nu }{\mathcal {L}} \bigg \rangle , \end{aligned}$$
(40)

integrating by parts, and using (32) in vacuum, we obtain

$$\begin{aligned} t_{\mu \nu } = \frac{1}{2\kappa ^2 }\left\langle 2m^{-2}\partial _{\mu }\partial _{\nu }h^{\alpha \beta }\Box h_{\alpha \beta } + \partial _{\mu }h^{\alpha \beta } \partial _{\nu }h_{\alpha \beta } \right\rangle , \end{aligned}$$
(41)

where the angle brackets denote the average over a macroscopic region.

The substitution of the plane wave (33) into (41) gives

$$\begin{aligned} t_{\mu \nu } = \frac{1}{2\kappa ^2} \left[ \left( a^{\alpha \beta }a_{\alpha \beta }\right) k_{\mu }k_{\nu } - \left( b^{\alpha \beta }b_{\alpha \beta }\right) q_{\mu }q_{\nu } \right] , \end{aligned}$$
(42)

where we used (34). We can see from (42) that the energies of the waves with the two physical polarizations \(a_{\mu \nu }\) are positive, while those of the ones with the five physical polarizations \(b_{\mu \nu }\) are negative. Since these waves do not interact with each other, energy cannot flow between them so that there is no violation of energy conservation.

5 Gravitational waves from a binary system

In order to analyze the gravitational waves created from binary systems, we need to solve the fourth-order differential equation (32) in the presence of matter (\(T^{(0)}_{\mu \nu } \ne 0\)). For simplicity, we can split \(h_{\mu \nu }\) according to

$$\begin{aligned} h_{\mu \nu } = A_{\mu \nu } + B_{\mu \nu }, \end{aligned}$$
(43)

where \(A_{\mu \nu }\) and \(B_{\mu \nu }\) obey the second-order differential equations

$$\begin{aligned}&\Box A_{\mu \nu } = -\kappa ^2 T^{(0)}_{\mu \nu }, \end{aligned}$$
(44)
$$\begin{aligned}&\left( \Box - m^{2}\right) B_{\mu \nu } = \kappa ^2 T^{(0)}_{\mu \nu }. \end{aligned}$$
(45)

In the frequency domain, the field equations (44) and (45) become

$$\begin{aligned}&\left( \nabla ^{2} + \omega ^2 \right) {\tilde{A}}_{\mu \nu }(\omega ,\mathbf{x }) = -\kappa ^2 {\tilde{T}}^{(0)}_{\mu \nu }(\omega ,\mathbf{x }), \end{aligned}$$
(46)
$$\begin{aligned}&\left( \nabla ^{2} + \omega ^2 - m^{2}\right) {\tilde{B}}_{\mu \nu }(\omega , \mathbf{x }) = \kappa ^2 {\tilde{T}}^{(0)}_{\mu \nu }(\omega ,\mathbf{x }), \end{aligned}$$
(47)

where \(\nabla ^{2}\) is the Laplacian, \(\omega \) is the frequency of the wave, and the tilde denotes the Fourier transform. The general solutions to (46) and (47) are given by

$$\begin{aligned}&{\tilde{A}}_{\mu \nu }(\omega , \mathbf{x }) = - \kappa ^2 \int { d^{3}\mathbf{x }'}{\tilde{G}}_{A}(\omega ,\mathbf{r }){\tilde{T}}^{(0)}_{\mu \nu } (\omega , \mathbf{x }'), \end{aligned}$$
(48)
$$\begin{aligned}&{\tilde{B}}_{\mu \nu }(\omega , \mathbf{x }) = \kappa ^2 \int { d^{3}\mathbf{x }'}{\tilde{G}}_{B}(\omega ,\mathbf{r }){\tilde{T}}^{(0)}_{\mu \nu } (\omega , \mathbf{x }'), \end{aligned}$$
(49)

where the frequency domain Green functions \({\tilde{G}}_{A}(\omega ,\mathbf{r })\) and \({\tilde{G}}_{B}(\omega ,\mathbf{r })\) are defined by

$$\begin{aligned}&\left( \nabla ^{2} + \omega ^2 \right) {\tilde{G}}_{A}(\omega ,\mathbf{r }) = \sqrt{8\pi }\delta ^{(3)}(\mathbf{r }), \end{aligned}$$
(50)
$$\begin{aligned}&\left( \nabla ^{2} + \omega ^2 - m^{2}\right) {\tilde{G}}_{B}(\omega ,\mathbf{r }) = \sqrt{8\pi }\delta ^{(3)}(\mathbf{r }), \end{aligned}$$
(51)

with \(\mathbf{r } = \mathbf{x } - \mathbf{x }'\) being the difference between the positions of the observer (\(\mathbf{x }\)) and the source (\(\mathbf{x }'\)).

It follows from (50) that

$$\begin{aligned} {\tilde{G}}_{A}(\omega ,\mathbf{r }) = -\frac{ e^{i\omega |\mathbf{r }|}}{{4\pi |\mathbf{r }|}}, \end{aligned}$$
(52)

and from (51) that

$$\begin{aligned} {\tilde{G}}_{B}(\omega ,\mathbf{r }) = -\frac{e^{ik_{\omega }|\mathbf{r }|} \varTheta (\omega - m)+ \text{ c.c. }\varTheta (-\omega - m)}{4\pi |\mathbf{r }|} \end{aligned}$$
(53)

for \(m^2 < \omega ^2\), and

$$\begin{aligned} {\tilde{G}}_{B}(\omega ,\mathbf{r }) = -\frac{e^{-k_{m}|\mathbf{r }|} \varTheta (m -|\omega |)}{4\pi |\mathbf{r }|} \end{aligned}$$
(54)

for \(m^2 > \omega ^2\), where c.c. is the complex conjugate of the exponential function, \(\varTheta \) is the Heaviside step function, \(k_{\omega } = \sqrt{\omega ^2-m^2}\), and \(k_{m} = \sqrt{m^2-\omega ^2}\).

Substituting (52) into (48), transforming back to real space, and using the far zone approximation (\(|\mathbf{x }| \approx |\mathbf{x } - \mathbf{x }'|\)), we can write the spatial components of \(A_{\mu \nu }\) in the form

$$\begin{aligned} A_{ij}(t,\mathbf{r }) = \frac{\kappa ^2}{8\pi r} \left[ \int _{-\infty }^{+\infty }\frac{d\omega }{2\pi }e^{-i\omega (t-r)}\right] \omega ^{2}{\tilde{Q}}_{ij}(\omega ), \end{aligned}$$
(55)

where \(r = |\mathbf{x }|\) is the distance between the observer and the source, and \({\tilde{Q}}_{ij}(\omega )\) is the Fourier transform of the reduced quadrupole moment

$$\begin{aligned} Q_{ij}(t) = \int {d^{3}\mathbf{r }}\left( x_{i}x_{j} - \frac{1}{3}\delta _{ij}r^2\right) T^{(0)}_{00}(t, \mathbf{r }). \end{aligned}$$
(56)

Following the same steps for \(B_{\mu \nu }\), but now with the substitution of (53) and (54) into (49) , we obtain

$$\begin{aligned} B_{ij}(t,\mathbf{r })= & {} -\frac{\kappa ^2}{8\pi r} \bigg [\int _{m}^{\infty }\frac{d\omega }{2\pi }e^{-i\omega t}e^{ik_{\omega }r} \nonumber \\&+\int _{-\infty }^{-m}\frac{d\omega }{2\pi }e^{-i\omega t}e^{-ik_{\omega }r}\bigg ] \omega ^{2}{\tilde{Q}}_{ij}(\omega ) \end{aligned}$$
(57)

for \(m^2 < \omega ^2\), and

$$\begin{aligned} B_{ij}(t,\mathbf{r }) = -\frac{\kappa ^2}{8\pi r}\left[ \int _{-\infty }^{m}\frac{d\omega }{2\pi }e^{-i\omega t}e^{-k_{m}r}\right] \omega ^{2}{\tilde{Q}}_{ij}(\omega ) \end{aligned}$$
(58)

for \(m^2 > \omega ^2\).

In the case of a circular binary system formed by a pair of masses \(m_{1}\) and \(m_{2}\), separated by a distance d, orbiting each other in the xy-plane with frequency \(\omega _{s} = \omega /2\), we have \(T^{f(0)}_{\mu \nu }(t,\mathbf{x }) = \mu \delta ^{0}_{\mu }\delta ^{0}_{\nu } \delta ^3(\mathbf{x })\). The substitution of this value into (24) gives

$$\begin{aligned} T^{(0)}_{00}(t,\mathbf{r }) = \frac{3}{4}\mu \delta ^3(\mathbf{r }), \end{aligned}$$
(59)

where \(\mu = m_{1}m_{2}/(m_{1}+m_{2})\) is the reduced mass. By inserting (59) and the relative coordinates

$$\begin{aligned} x_{1} = -d\sin (\omega _{s}t), \ \ \ \ \ x_{2} = d\cos (\omega _{s}t), \ \ \ \ \ x_{3} = 0, \end{aligned}$$
(60)

into (56), and taking the Fourier transform, we find

$$\begin{aligned} {\tilde{Q}}_{11}(\omega )= & {} \frac{3\pi \mu d^{2}}{8}\left[ \delta (\omega ) - \delta (\omega + 2\omega _{s}) - \delta (\omega - 2\omega _{s})\right] ,\nonumber \\ \end{aligned}$$
(61)
$$\begin{aligned} {\tilde{Q}}_{22}(\omega )= & {} \frac{3\pi \mu d^{2}}{8}\left[ \delta (\omega ) + \delta (\omega + 2\omega _{s}) + \delta (\omega - 2\omega _{s})\right] ,\nonumber \\ \end{aligned}$$
(62)
$$\begin{aligned} {\tilde{Q}}_{12}(\omega )= & {} \frac{3\pi \mu d^{2}}{8i}\left[ \delta (\omega - 2\omega _{s}) - \delta (\omega + 2\omega _{s})\right] , \end{aligned}$$
(63)

where we omitted the term proportional to \(\delta _{ij}\) in (56) because it does not contribute to the radiated energy.

The insertion of (61)–(63) into (55) gives

$$\begin{aligned} A_{11}(t,r)= & {} - A_{22}(t,r) = \frac{2\mu d^{2}\omega _{s}^2}{r} \cos {(2\omega _{s}t_{\text {ret}})}, \end{aligned}$$
(64)
$$\begin{aligned} A_{12}(t,r)= & {} A_{21}(t,r) = \frac{2\mu d^{2}\omega _{s}^2}{r} \sin {(2\omega _{s}t_{\text {ret}})}, \end{aligned}$$
(65)

where \(t_{\text {ret}} = t -r\) is the retarded time. In the same way, substituting (61)–(63) into (57) and (58), we arrive at

$$\begin{aligned} B_{11}(t,r)= & {} - B_{22}(t,r) = -\frac{2\mu d^{2}\omega _{s}^2}{r} \cos {(2\omega _{s}t_{m})}, \end{aligned}$$
(66)
$$\begin{aligned} B_{12}(t,r)= & {} B_{21}(t,r) = -\frac{2\mu d^{2}\omega _{s}^2}{r} \sin {(2\omega _{s}t_{m})}, \end{aligned}$$
(67)

for \(m^2 < 4\omega _{s}^2\), and

$$\begin{aligned} B_{11}(t,r)= & {} - B_{22}(t,r) = -\frac{2\mu d^{2}\omega _{s}^2}{r} e^{-k_{m}r}\cos {(2\omega _{s}t)},\nonumber \\ \end{aligned}$$
(68)
$$\begin{aligned} B_{12}(t,r)= & {} B_{21}(t,r) = -\frac{2\mu d^{2}\omega _{s}^2}{r} e^{-k_{m}r}\sin {(2\omega _{s}t)}, \end{aligned}$$
(69)

for \(m^2 > 4\omega _{s}^2\), where \(t_{m}= t - v_{m}r\) is the travel time, with \(v_{m}=\sqrt{1-m^{2}/(4\omega _{s}^2)}\) being the speed of the massive gravitational wave.

The rate of energy loss from a source, in the far field limit, is given by

$$\begin{aligned} {\dot{E}} = - r^2\int _{\partial V}{d\varOmega } \, t^{0i}n_{i}, \end{aligned}$$
(70)

where the dot is the derivative with respect to time, \(d\varOmega = \sin \theta d\theta d\phi \) is the differential solid angle, \(\partial V\) is the surface of a spherical shell with volume V centered around the source, and \(n_{i}\) are the components of the spatial unit vector pointing from the source to the observer.

By substituting (43) into (41), integrating by parts, and using (44) and (45) in vacuum, we obtain

$$\begin{aligned} t^{0i}n_{i} = \frac{1}{2\kappa ^2}n_{i}\big \langle \partial ^{0} A^{\alpha \beta }\partial ^{i} A_{\alpha \beta } - \partial ^{0}B^{\alpha \beta } \partial ^{i} B_{\alpha \beta } \big \rangle . \end{aligned}$$
(71)

We can see from (64) and (65) that

$$\begin{aligned} \partial _{i}A_{\alpha \beta } \approx - n_{i}\partial _{0}A_{\alpha \beta }, \end{aligned}$$
(72)

and from (66)–(69) that

$$\begin{aligned} \partial _{i}B_{\alpha \beta } \approx - n_{i}v_{m}\partial _{0}B_{\alpha \beta } \end{aligned}$$
(73)

for \(m^2 < 4\omega _{s}^2\), and

$$\begin{aligned} \partial _{i}B_{\alpha \beta } \approx - n_{i}k_{m}B_{\alpha \beta } \end{aligned}$$
(74)

for \(m^2 > 4\omega _{s}^2\), where we neglected terms of order \(1/r^2\). In addition, considering the conservation and the traceless condition of \(T^{(0)}_{\mu \nu }\), and the traceless condition of \(Q_{ij}\), it follows from (48), (49), (55), (57) and (58) that \(A_{\mu \nu }\) and \(B_{\mu \nu }\) obey the traceless-transverse conditions

$$\begin{aligned}&\partial ^{i}A_{ij} = 0, \ \ \ \ \ \ A_{0i} = 0, \ \ \ \ \ \ \delta ^{ij}A_{ij} = 0, \end{aligned}$$
(75)
$$\begin{aligned}&\partial ^{i}B_{ij} = 0, \ \ \ \ \ \ B_{0i} = 0, \ \ \ \ \ \ \delta ^{ij}B_{ij} = 0. \end{aligned}$$
(76)

By using the combination of (72)–(76), and \(n^{i}n_{i} = 1\), we can write (71) as

$$\begin{aligned} t^{0i}n_{i} \approx \frac{1}{2\kappa ^2}\varLambda _{ijkl}\big \langle \partial _{0}A^{ij}\partial _{0} A^{kl} - v_{m}\partial _{0}B^{ij}\partial _{0} B^{kl} \big \rangle \end{aligned}$$
(77)

for \(m^2 < 4\omega _{s}^2\), and

$$\begin{aligned} t^{0i}n_{i} \approx \frac{1}{2\kappa ^2}\varLambda _{ijkl}\big \langle \partial _{0}A^{ij}\partial _{0} A^{kl} - k_{m}B^{ij}\partial _{0} B^{kl} \big \rangle \end{aligned}$$
(78)

for \(m^2 > 4\omega _{s}^2\), where

$$\begin{aligned} \varLambda _{ijkl} = P_{ik}P_{jl} -\frac{1}{2}P_{ij}P_{kl} \end{aligned}$$
(79)

is the Lambda tensor, with

$$\begin{aligned} P_{ij} = \delta _{ij} - n_{i}n_{j} \end{aligned}$$
(80)

being the traceless-transverse projection operator.

Inserting (77) and (78) into (70), and using the surface integral

$$\begin{aligned} \int _{\partial V}{d\varOmega }\varLambda _{ijkl} = \frac{2\pi }{15}\left( 11\delta _{ik}\delta _{jl} -4 \delta _{ij}\delta _{kl} + \delta _{il}\delta _{jk} \right) , \end{aligned}$$
(81)

we find

$$\begin{aligned} {\dot{E}} \approx - \frac{3 r^2}{40}\big \langle \partial _{0}A^{ij}\partial _{0} A_{ij} - v_{m}\partial _{0}B^{ij} \partial _{0} B_{ij} \big \rangle \end{aligned}$$
(82)

for \(m^2 < 4\omega _{s}^2\), and

$$\begin{aligned} {\dot{E}} \approx - \frac{3 r^2}{40}\big \langle \partial _{0}A^{ij}\partial _{0} A_{ij} - k_{m}B^{ij}\partial _{0} B_{ij} \big \rangle \end{aligned}$$
(83)

for \(m^2 > 4\omega _{s}^2\).

Finally, substituting (64)–(69) into (82) and (83), we arrive at

$$\begin{aligned} {\dot{E}} \approx \frac{3}{8} \left( 1-\sqrt{1-\frac{m^{2}}{4\omega _{s}^2}} \,\right) {\dot{E}}_{\text {GR}} \end{aligned}$$
(84)

for \(m^2 < 4\omega _{s}^2\), and

$$\begin{aligned} {\dot{E}} \approx \frac{3}{8}\, {\dot{E}}_{\text {GR}} \end{aligned}$$
(85)

for \(m^2 > 4\omega _{s}^2\), where

$$\begin{aligned} {\dot{E}}_{\text {GR}} = -\frac{32G\mu ^{2}d^{4}\omega _{s}^6}{5} \end{aligned}$$
(86)

is the standard energy loss of general relativity.

Experiments on the inverse square law of the MCG gravitational potential of a point particle with a mass M, which is given by [8]

$$\begin{aligned} \phi (r) = -\frac{GM}{r}\left( 1 - e^{-mr}\right) , \end{aligned}$$
(87)

constrain the graviton mass to the ranges \(m< 10^{-22}\, \text {eV}\) for \(m^2 < 4\omega _{s}^2\) and \(m> 10^{-2}\, \text {eV}\) for \(m^2 > 4\omega _{s}^2\) [12]. In particular, considering the orbital frequency \(\omega _{s} \approx 1.3\times 10^{-20}\, \text {eV}\) of the binary system PSR J1012 + 5307 formed by a neutron star and a white dwarf in quasi-circular motion [13,14,15], we have that \(m^2/4\omega _{s}^2 < 10^{-5}\) for \(m^2 < 4\omega _{s}^2\). By substituting this value into (84), we find

$$\begin{aligned} {\dot{E}} \lessapprox 10^{-6} {\dot{E}}_{\text {GR}} \end{aligned}$$
(88)

for \(m^2 < 4\omega _{s}^2\).

The energy loss of the binary system results in a decay of its orbital period P, which can be written as

$$\begin{aligned} \frac{{\dot{P}}}{P} = \frac{{\dot{d}}}{2d} - \frac{{\dot{\phi }}'}{2\phi '}, \end{aligned}$$
(89)

where \(\phi '(d)=\mu ^{-1}\partial _{d}U(d)\) is the derivative of the gravitational potential \(\phi \) with respect to d, with U being the gravitational potential energy. It follows from (87) that

$$\begin{aligned} U(d) = -\frac{Gm_{1}m_{2}}{d}\left( 1 - e^{-md}\right) . \end{aligned}$$
(90)

Substituting this result into (89), we find

$$\begin{aligned} \frac{{\dot{P}}}{P} \approx -\frac{3}{2}\frac{|{\dot{E}}|}{|E_{\text {GR}}|} \end{aligned}$$
(91)

for both \(m^2 < 4\omega _{s}^2\) and \(m^2 > 4\omega _{s}^2\), where \(|E_{\text {GR}}| = G m_{1}m_{2}/2d\) and we considered \(m^2d^2 \ll 1\) for \(m^2 < 4\omega _{s}^2\).

The insertion of (88) and (85) into (91) then gives

$$\begin{aligned} \frac{{\dot{P}}}{P} \lessapprox 10^{-6}\frac{{\dot{P}}_{\text {GR}}}{P_{\text {GR}}} \end{aligned}$$
(92)

for \(m^2 < 4\omega _{s}^2\), and

$$\begin{aligned} \frac{{\dot{P}}}{P} \approx \frac{3}{8}\frac{{\dot{P}}_{\text {GR}}}{P_{\text {GR}}} \end{aligned}$$
(93)

for \(m^2 > 4\omega _{s}^2\), where \({\dot{P}}_{\text {GR}}/P_{\text {GR}} = -3|{\dot{E}}_{\text {GR}}|/2|E_{\text {GR}}|\). We can see from (92) that the MCG decay of the orbital period for \(m^2 < 4\omega _{s}^2\) is several orders of magnitude smaller than in general relativity, which rules out the theory with a small graviton mass (\(m< 10^{-22}\, \text {eV}\)). On the other hand, the discrepancy seen in (93) between the MCG decay of the orbital period for \(m^2 > 4\omega _{s}^2\) and the general relativity result disappears if we include the first-order term \(2S_{0}^{2}G^{(1)}_{\mu \nu }\) in the matter energy–momentum tensor (24) as done in Ref. [16] to CG.Footnote 3 Taking this into account, the theory with a large graviton mass (\(m> 10^{-2}\, \text {eV}\)) can explain the decrease of the orbital period of binary systems.

6 Final remarks

Here we have shown that the MCG plane wave has seven propagating degrees of freedom, two of which are massless and carry positive energies and the other five are massive and carry negative energies. Despite the presence of the waves with negative energies, the lack of interaction between them and the waves with positive energies means that the energies of all the seven MCG plane waves do not diverge.

The study on the radiated energy from a binary system restricts the MCG to large graviton mass, which give rise to a modification of general relativity only at high energies and small distances. Although this modification is the cause of MCG being renormalizable [18, 19], it makes the theory unable to explain galaxy rotation curves without dark matter. However, the conformal symmetry of the matter part of MCG allows us to consider that the Higgs mass is generated by the symmetry breaking of an extra scalar field, which may be a good candidate for dark matter. Further studies are needed to figure this out.