Gravitational waves in non-local gravity

Theories of physics describing elementary interactions are local, that is, fields are evaluated at the same point, and are governed by point-like Lagrangians from which one derives equations of motion. However, already at classical level, it is possible to observe non-locality in Electrodynamics of continuous media, when spatial or temporal dispersions, due to the non-local constitutive relation between the fields $\left(\mathbf{D},\mathbf{H}\right)$ and $\left(\mathbf{E},\mathbf{B}\right)$, occur [1–4] as

$$\begin{equation}{D}_{i}\left(x\right)=\int \enspace {\mathrm{d}}^{4}{x}^{\prime }{{\epsilon}}_{ij}\left({x}^{\prime }\right){E}^{j}\left(x-{x}^{\prime }\right).\end{equation} \tag{ 1.1 }$$

It is a sort of memory-dependent phenomenon taking into account both the past history of the fields and their values taken in other points of the medium. Also at quantum level, some effective actions show non-local terms and therefore, the associated field equations are integro-differential ones. Recently, non-locality has been considered in cosmology taking into account non-local models to explain early and late-time cosmic acceleration as well as structure formation, without introducing dark energy and dark matter. Specifically, non-locality can play interesting roles to address problems like cosmological constant, Big Bang and black hole singularities, and, in general, coincidence and fine-tuning problems, which affect the ΛCDM model [5–10].

Some non-local field theories are of infinite order because they have an infinite number of derivatives. This feature is due to the presence of operators like $f\left(\square \right)$ which can be expanded in series, assuming that f is an analytic function, as

$$\begin{equation}f(\square )=\sum\limits _{n=0}^{\infty }{a}_{n}{\square }^{n}.\end{equation} \tag{ 1.2 }$$

Here □ is the d'Alembert operator. This procedure is aimed to make the theory ghost free [11–14].

In particular, we know that general relativity describes gravity as a local interaction while quantum mechanics shows non-local aspects. Several approaches have been proposed to achieve a self-consistent quantum gravity as discussed, for example, in [15, 16]. A possibility toward quantum gravity is considering non-local corrections to the Hilbert–Einstein action [17, 18]. It is a natural way to cure ultraviolet and infrared behaviors of general relativity. Introducing non-local terms can work also in alternative theories like teleparallel gravity [19].

In all these approaches, it is important to study the linearized versions of the theories and to derive gravitational waves (GWs). In fact, gravitational radiation allows to detect possible effects of non-local gravity [20] as well as to classify the degrees of freedom of a given theory [21–23].

In this paper, we want to investigate the effect of non-locality in a theory of gravity where the Ricci scalar R of general relativity is corrected by R□⁻¹ R, which is the first interesting non-local curvature term. Specifically, we want to investigate how polarization, helicity, and mass of GWs are affected by this kind of terms and how the E(2) Petrov classification changes.

Section 2 is devoted to a procedure for the localization of the non-local gravitational action by the method of Lagrangian multipliers. Starting from this localized action, it is possible to derive the field equations. Subsequently, in section 3, we linearize the field equations and solve them in harmonic gauge and in a further gauge in view of eliminating ghost modes. Finally, we get the GWs. In sections 4 and 5, polarizations are analyzed by using both geodesic deviation and the Newman–Penrose (NP) formalism. The method consists in expanding the nearly null massive plane waves in term of the exactly null plane waves. The expansion is achieved in terms of a parameter γ_k . In section 6, results are summarized and possible future developments are discussed.

Let us study gravitational interaction governed by the following non-local action

$$\begin{equation}S[g]=\frac{1}{2{\kappa }^{2}}\int {\mathrm{d}}^{4}x\sqrt{-g}\left(R+{a}_{1}R{\square }^{-1}R\right)+\int {\mathrm{d}}^{4}x\sqrt{-g}{\mathcal{L}}_{m}[g],\end{equation} \tag{ 2.1 }$$

where k² = 8πG/c⁴. Its field equations are non-linear integro-differential equations due to the non-local term. We introduce the auxiliary field ϕ(x) defined as

$$\begin{equation}\phi (x)={\square }^{-1}R,\end{equation} \tag{ 2.2 }$$

and then the Ricci scalar is

$$\begin{equation}R=\square \phi (x).\end{equation} \tag{ 2.3 }$$

According to this definition, a Lagrange multiplier can be considered so that the gravitational action becomes [5]

$$\begin{equation}{S}_{\mathrm{g}}[g,\phi ,\lambda ]=\frac{1}{2{\kappa }^{2}}\int {\mathrm{d}}^{4}x\sqrt{-g}\left[R\left(1+{a}_{1}\phi \right)+\lambda \left(\square \phi -R\right)\right],\end{equation} \tag{ 2.4 }$$

where λ(x) is a further scalar field. Using integration by parts and imposing that fields and their derivatives vanish onto the boundary of integration domain, we obtain

$$\begin{equation}{S}_{\mathrm{g}}[g,\phi ,\lambda ]=\frac{1}{2{\kappa }^{2}}\int {\mathrm{d}}^{4}x\sqrt{-g}\left[R\left(1+{a}_{1}\phi -\lambda \right)-{\nabla }^{\nu }\lambda {\nabla }_{\nu }\phi \right].\end{equation} \tag{ 2.5 }$$

Varying with respect to ϕ, we get

$$\begin{equation}\frac{1}{\sqrt{-g}}\frac{\delta S}{\delta \phi }=\frac{1}{2{\kappa }^{2}}\left[{a}_{1}R+\square \lambda \right]=0{\Rightarrow}\square \lambda =-{a}_{1}R,\end{equation} \tag{ 2.6 }$$

while varying with respect to λ, the functional derivative takes the form

$$\begin{equation}\frac{1}{\sqrt{-g}}\frac{\delta S}{\delta \lambda }=\frac{1}{2{\kappa }^{2}}\left[-R+\square \phi \right]=0{\Rightarrow}\square \phi =R.\end{equation} \tag{ 2.7 }$$

Finally, the variation with respect to the metric g^μν gives

$$\begin{align}\frac{2{\kappa }^{2}}{\sqrt{-g}}\frac{\delta {S}_{\mathrm{g}}}{\delta {g}^{\mu \nu }}& =\left({G}_{\mu \nu }+{g}_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu }\right)\left(1+{a}_{1}\phi -\lambda \right)-{\nabla }_{\left(\right.\mu }\phi {\nabla }_{\nu \left.\right)}\lambda +\frac{1}{2}{g}_{\mu \nu }{\nabla }^{\sigma }\phi {\nabla }_{\sigma }\lambda ,\end{align} \tag{ 2.8 }$$

and

$$\begin{equation}\frac{\delta \left(\sqrt{-g}{\mathcal{L}}_{m}\right)}{\delta {g}^{\mu \nu }}=-\frac{\sqrt{-g}}{2}{T}_{\mu \nu }.\end{equation} \tag{ 2.9 }$$

The final field equations are

$$\begin{equation}\begin{array}{|l| }\hline \left({G}_{\mu \nu }+{g}_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu }\right)\left(1+{a}_{1}\phi -\lambda \right)-{\nabla }_{\left(\right.\mu }\phi {\nabla }_{\nu \left.\right)}\lambda +\frac{1}{2}{g}_{\mu \nu }{\nabla }^{\sigma }\phi {\nabla }_{\sigma }\lambda ={\kappa }^{2}{T}_{\mu \nu },\\ \hline\end{array}\end{equation} \tag{ 2.10 }$$

$$\begin{equation}\begin{array}{|l|}\hline \square \phi =R,\\ \hline\end{array}\end{equation} \tag{ 2.11 }$$

$$\begin{equation}\begin{array}{|l|}\hline \square \lambda =-{a}_{1}R,\\ \hline\end{array}\end{equation} \tag{ 2.12 }$$

where G_μν is the Einstein tensor

$$\begin{equation}{G}_{\mu \nu }={R}_{\mu \nu }-\frac{1}{2}{\eta }_{\mu \nu }R.\end{equation} \tag{ 2.13 }$$

According to equations (2.11) and (2.12), the trace of equation (2.10) is

$$\begin{equation}\begin{array}{|l|}\hline \left[1+{a}_{1}\left(\phi -6\right)-\lambda \right]R-{\nabla }^{\rho }\phi {\nabla }_{\rho }\lambda =-{\kappa }^{2}T,\\ \hline\end{array}\end{equation} \tag{ 2.14 }$$

and now we can start our considerations on the weak field behavior of this theory.

In order to analyze gravitational radiation, let us perturb the metric tensor g_μν around the flat metric η_μν and the two scalar fields ϕ and λ around their values in Minkowskian spacetime ϕ₀ and λ₀. It is

$$\begin{equation}{g}_{\mu \nu }\sim {\eta }_{\mu \nu }+{h}_{\mu \nu },\end{equation} \tag{ 3.1 }$$

$$\begin{equation}\phi \sim {\phi }_{0}+\delta \phi \enspace ,\end{equation} \tag{ 3.2 }$$

$$\begin{equation}\lambda \sim {\lambda }_{0}+\delta \lambda .\end{equation} \tag{ 3.3 }$$

At first order in h_μν , the Ricci tensor R_μν and the Ricci scalar R become

$$\begin{equation}{R}_{\mu \nu }^{(1)}=\frac{1}{2}\left({\partial }_{\sigma }{\partial }_{\mu }{h}_{\nu }^{\sigma }+{\partial }_{\sigma }{\partial }_{\nu }{h}_{\mu }^{\sigma }-{\partial }_{\mu }{\partial }_{\nu }h-\square {h}_{\mu \nu }\right),\end{equation} \tag{ 3.4 }$$

$$\begin{equation}{R}^{(1)}={\partial }_{\mu }{\partial }_{\nu }{h}^{\mu \nu }-\square {h}_{\mu \nu },\end{equation} \tag{ 3.5 }$$

where h is the trace of perturbation h_μν . In vacuum and under Lorentz gauge, equations (3.4) and (3.5) and the Einstein tensor G_μν become

$$\begin{equation}{R}_{\mu \nu }^{(1)}=-\frac{1}{2}\square {h}_{\mu \nu },\end{equation} \tag{ 3.6 }$$

$$\begin{equation}{R}^{(1)}=-\frac{1}{2}\square h,\end{equation} \tag{ 3.7 }$$

$$\begin{equation}{G}_{\mu \nu }^{(1)}=-\frac{1}{2}\square {h}_{\mu \nu }+\frac{1}{4}{\eta }_{\mu \nu }\square h.\end{equation} \tag{ 3.8 }$$

Later, according to equations (2.10)–(2.12), the linearized field equations in vacuum are

$$\begin{equation}\left(1+{a}_{1}{\phi }_{0}-{\lambda }_{0}\right)\left(2\square {h}_{\mu \nu }-{\eta }_{\mu \nu }\square h\right)+4{\partial }_{\mu }{\partial }_{\nu }\left({a}_{1}\delta \phi -\delta \lambda \right)+4{a}_{1}{\eta }_{\mu \nu }\square h=0,\end{equation} \tag{ 3.9 }$$

$$\begin{equation}\square \left(2\delta \phi +h\right)=0,\end{equation} \tag{ 3.10 }$$

$$\begin{equation}\square \left(2\delta \lambda -{a}_{1}h\right)=0,\end{equation} \tag{ 3.11 }$$

and, from equations (3.10) and (3.11), the linearized trace equation is

$$\begin{equation}\left[1+{a}_{1}\left({\phi }_{0}-6\right)-{\lambda }_{0}\right]\square h=0.\end{equation} \tag{ 3.12 }$$

The trace equation (3.12) admits solutions if

$$\begin{equation}1+{a}_{1}{\phi }_{0}-{\lambda }_{0}\ne 6{a}_{1}\quad \to \quad \square h=0,\end{equation} \tag{ 3.13 }$$

or

$$\begin{equation}1+{a}_{1}{\phi }_{0}-{\lambda }_{0}=6{a}_{1}\quad \to \quad \square h\ne 0\quad \text{and}\quad \square h=0.\end{equation} \tag{ 3.14 }$$

In the case (3.13), equations (3.9)–(3.11) become

$$\begin{equation}\left(1+{a}_{1}{\phi }_{0}-{\lambda }_{0}\right)\square {h}_{\mu \nu }+2{\partial }_{\mu }{\partial }_{\nu }\left({a}_{1}\delta \phi -\delta \lambda \right)=0,\end{equation} \tag{ 3.15 }$$

$$\begin{equation}\square \delta \phi =0,\end{equation} \tag{ 3.16 }$$

$$\begin{equation}\square \delta \lambda =0.\end{equation} \tag{ 3.17 }$$

In k-space, considering the Fourier transform, equation (3.15) takes the following form

$$\begin{equation}\left(1+{a}_{1}{\phi }_{0}-{\lambda }_{0}\right){k}^{2}{\tilde {h}}_{\mu \nu }\left(k\right)+2{k}_{\mu }{k}_{\nu }\left({a}_{1}\tilde {\delta \phi }\left(k\right)-\tilde {\delta \lambda }\left(k\right)\right)=0.\end{equation} \tag{ 3.18 }$$

It implies

$$\begin{equation}\square h=0{\Rightarrow}{k}^{2}\tilde {h}\left(k\right)=0{\Rightarrow}{k}^{2}=0,\end{equation} \tag{ 3.19 }$$

and, putting (3.19) into (3.18), we get a solution if

$$\begin{equation}{a}_{1}\tilde {\delta \phi }\left(k\right)=\tilde {\delta \lambda }\left(k\right),\end{equation} \tag{ 3.20 }$$

which, in x-space, becomes

$$\begin{equation}{a}_{1}\delta \phi \left(x\right)=\delta \lambda \left(x\right).\end{equation} \tag{ 3.21 }$$

Inserting (3.21) into (3.15), for 1 + a₁ ϕ₀ − λ₀ ≠ 0, we obtain

$$\begin{equation}\square {h}_{\mu \nu }=0,\end{equation} \tag{ 3.22 }$$

which, together with equations (3.16) and (3.17), yield massless, two-helicity transverse waves solutions for ${k}_{1}^{2}=0$, namely the standard GWs of general relativity. In the case (3.14), excluding □h = 0 because the previous massless case returns, equations (3.9)–(3.11) become a coupled partial differential equations system, i.e.

$$\begin{equation}6{a}_{1}\square {h}_{\mu \nu }-{a}_{1}{\eta }_{\mu \nu }\square h+2{\partial }_{\mu }{\partial }_{\nu }\left({a}_{1}\delta \phi -\delta \lambda \right)=0,\end{equation} \tag{ 3.23 }$$

$$\begin{equation}\square \delta \phi =-\frac{1}{2}\square h,\end{equation} \tag{ 3.24 }$$

$$\begin{equation}\square \delta \lambda =\frac{{a}_{1}}{2}\square h,\end{equation} \tag{ 3.25 }$$

with the additional condition

$$\begin{equation}\square h\ne 0.\end{equation} \tag{ 3.26 }$$

From equation (3.26), in the momentum space, we have

$$\begin{equation}{k}_{2}^{2}\tilde {h}\left(k\right)\ne 0{\Rightarrow}{k}_{2}^{2}\ne 0\quad \text{and}\quad \tilde {h}\left(k\right)\ne 0,\end{equation} \tag{ 3.27 }$$

where ${\left({k}_{2}\right)}^{\mu }=({\omega }_{2},\mathbf{k})$ is the wave four-vector and ${k}_{2}^{2}={M}^{2}$. We assume that equation (3.26) is

$$\begin{equation}\square h\left(x\right)=g\left(x\right),\end{equation} \tag{ 3.28 }$$

choosing g(x) with any k² ≠ 0 such as

$$\begin{equation}g(x)=\frac{1}{{(2\pi )}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\enspace \tilde {q}\left(\mathbf{k}\right){\text{e}}^{\text{i}\varphi }\enspace {\text{e}}^{\text{i}{\omega }_{2}t}\enspace {\text{e}}^{-\text{i}\mathbf{k}\cdot \mathbf{x}}+\mathrm{c}.\mathrm{c}.\end{equation} \tag{ 3.29 }$$

$$\begin{equation}=\frac{1}{{(2\pi )}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\enspace 2\tilde {q}\left(\mathbf{k}\right)\mathrm{cos}\left({\omega }_{2}t-\mathbf{k}\cdot \mathbf{x}+\varphi \right),\end{equation} \tag{ 3.30 }$$

where $\tilde {q}\left(\mathbf{k}\right)\in {L}^{2}\left({\mathbf{R}}^{3}\right)$ is a square integrable function in R³ and φ is the phase. By performing the Fourier transform with respect to the spatial coordinates, the non-homogeneous wave equation (3.28) becomes, in k-space,

$$\begin{equation}\left({\partial }_{0}^{2}+\vert \mathbf{k}{\vert }^{2}\right)\tilde {h}\left(t,\mathbf{k}\right)=\tilde {q}\left(\mathbf{k}\right){\text{e}}^{\text{i}\varphi }\enspace {\text{e}}^{\text{i}{\omega }_{2}t},\end{equation} \tag{ 3.31 }$$

which has a particular solution

$$\begin{equation}\tilde {h}\left(t,\mathbf{k}\right)=-\frac{\tilde {q}\left(\mathbf{k}\right){\text{e}}^{\text{i}\varphi }}{{k}_{2}^{2}}\enspace {\text{e}}^{\text{i}{\omega }_{2}t}=\tilde {A}\left(\mathbf{k}\right){\text{e}}^{\text{i}{\omega }_{2}t},\end{equation} \tag{ 3.32 }$$

where ${\omega }_{2}=\sqrt{{M}^{2}+\vert \mathbf{k}{\vert }^{2}}$. The value M² corresponds to the squared mass of a new effective scalar field Ψ which can be defined as the combination

$$\begin{equation}{\Psi}\left(x\right)=\frac{\lambda \left(x\right)}{{a}_{1}}-\phi \left(x\right),\end{equation} \tag{ 3.33 }$$

whose linear perturbation δΨ satisfies the Klein–Gordon equation

$$\begin{equation}\left(\square +{M}^{2}\right)\delta {\Psi}=0,\end{equation} \tag{ 3.34 }$$

under the assumptions that g(x) behaves as a wave packet with k² exactly equal to M². It is worth noticing that the linear perturbation of Ψ in k-space, $\delta \tilde {{\Psi}}(\mathbf{k})$, corresponds to $\tilde {h}(\mathbf{k})=\tilde {A}(\mathbf{k})$ see equation (3.32). It is the trace of metric perturbation in k-space. From equations (3.24), (3.25) and (3.28), it is

$$\begin{equation}\square {\Psi}\left(x\right)=\square \left(\frac{\lambda \left(x\right)}{{a}_{1}}-\phi \left(x\right)\right)=\square h,\end{equation} \tag{ 3.35 }$$

which, in k-space with non-null k², becomes

$$\begin{equation}\delta \tilde {{\Psi}}\left(\mathbf{k}\right)=\frac{\delta \tilde {\lambda }\left(\mathbf{k}\right)}{{a}_{1}}-\delta \tilde {\phi }\left(\mathbf{k}\right)=\tilde {A}\left(\mathbf{k}\right).\end{equation} \tag{ 3.36 }$$

With these considerations in mind, a solution of equation (3.28) is

$$\begin{equation}h(x)=\frac{1}{{(2\pi )}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\enspace \tilde {A}\left(\mathbf{k}\right){\mathrm{e}}^{\mathrm{i}{k}_{2}\cdot x}+\mathrm{c}.\mathrm{c}.,\end{equation} \tag{ 3.37 }$$

and, according to the two wave equations (3.24) and (3.25) in k-space, it is

$$\begin{equation}\left({\partial }_{0}^{2}+\vert \mathbf{k}{\vert }^{2}\right)\left[2\enspace \tilde {\delta \phi }\left(t,\mathbf{k}\right)+\tilde {h}\left(t,\mathbf{k}\right)\right]=0,\end{equation} \tag{ 3.38 }$$

$$\begin{equation}\left({\partial }_{0}^{2}+\vert \mathbf{k}{\vert }^{2}\right)\left[2\enspace \tilde {\delta \lambda }\left(t,\mathbf{k}\right)-{a}_{1}\tilde {h}\left(t,\mathbf{k}\right)\right]=0.\end{equation} \tag{ 3.39 }$$

From equation (3.32), we obtain two particular solutions for fixed k

$$\begin{equation}\tilde {\delta \phi }\left(t,\mathbf{k}\right)=-\frac{1}{2}\tilde {A}\left(\mathbf{k}\right){\text{e}}^{\text{i}{\omega }_{2}t},\end{equation} \tag{ 3.40 }$$

$$\begin{equation}\tilde {\delta \lambda }\left(t,\mathbf{k}\right)=\frac{{a}_{1}}{2}\tilde {A}\left(\mathbf{k}\right){\text{e}}^{\text{i}{\omega }_{2}t},\end{equation} \tag{ 3.41 }$$

that allows us to derive the waves ${k}_{1}^{2}=0$ and ${k}_{2}^{2}={M}^{2}$ for the two scalar linear perturbations δϕ and δλ, i.e. we have

$$\begin{equation}\delta \phi (x)=\frac{1}{{\left(2\pi \right)}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\enspace \tilde {D}\left(\mathbf{k}\right){\text{e}}^{\text{i}{k}_{1}\cdot x}+\frac{1}{{(2\pi )}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\left(-\frac{1}{2}\right)\tilde {A}\left(\mathbf{k}\right){\mathrm{e}}^{\mathrm{i}{k}_{2}\cdot x}+\mathrm{c}.\mathrm{c}.,\end{equation} \tag{ 3.42 }$$

and

$$\begin{align}\delta \lambda (x)& =\frac{1}{{\left(2\pi \right)}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\enspace \left(-{a}_{1}\tilde {D}\left(\mathbf{k}\right)\right){\text{e}}^{\text{i}\enspace {k}_{1}\cdot x}+\frac{1}{{(2\pi )}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\left(\frac{{a}_{1}}{2}\right)\tilde {A}\left(\mathbf{k}\right){\mathrm{e}}^{\mathrm{i}{k}_{2}\cdot x}+\mathrm{c}.\mathrm{c}..\end{align} \tag{ 3.43 }$$

Now, carrying out the spatial Fourier transform of equation (3.23) and by means of equations (3.32), (3.40) and (3.41), we get

$$\begin{equation}6{a}_{1}\left({\partial }_{0}^{2}+\vert \mathbf{k}{\vert }^{2}\right){\tilde {h}}_{\mu \nu }\left(t,\mathbf{k}\right)+{a}_{1}{\eta }_{\mu \nu }{k}_{2}^{2}\tilde {A}\left(\mathbf{k}\right){\text{e}}^{\text{i}{\omega }_{2}t}+2{a}_{1}{\left({k}_{2}\right)}_{\mu }{\left({k}_{2}\right)}_{\nu }\tilde {A}\left(\mathbf{k}\right){\text{e}}^{\text{i}{\omega }_{2}t}=0,\end{equation} \tag{ 3.44 }$$

that simplifying gives

$$\begin{equation}\left({\partial }_{0}^{2}+\vert \mathbf{k}{\vert }^{2}\right){\tilde {h}}_{\mu \nu }\left(t,\mathbf{k}\right)=-\left(\frac{{\eta }_{\mu \nu }{k}_{2}^{2}+2{\left({k}_{2}\right)}_{\mu }{\left({k}_{2}\right)}_{\nu }}{6}\right)\tilde {A}\left(\mathbf{k}\right){\text{e}}^{\text{i}{\omega }_{2}t}.\end{equation} \tag{ 3.45 }$$

To solve the non-homogeneous equation (3.45), first we get the related homogeneous solution linked to the massless wave ${k}_{1}^{2}=0$,

$$\begin{equation}{\tilde {h}}_{\mu \nu }^{(o)}\left(t,\mathbf{k}\right)={C}_{\mu \nu }^{(1)}\left(\mathbf{k}\right){\text{e}}^{\text{i}{\omega }_{1}t}+{C}_{\mu \nu }^{(2)}\left(\mathbf{k}\right){\text{e}}^{-\text{i}{\omega }_{1}t},\end{equation} \tag{ 3.46 }$$

and then a particular solution, linked to the massive wave ${k}_{2}^{2}={M}^{2}$, is

$$\begin{equation}{\tilde {h}}_{\mu \nu }^{(p)}\left(t,\mathbf{k}\right)=\frac{{\eta }_{\mu \nu }{k}_{2}^{2}+2{\left({k}_{2}\right)}_{\mu }{\left({k}_{2}\right)}_{\nu }}{6{k}_{2}^{2}}\tilde {A}\left(\mathbf{k}\right){\text{e}}^{\text{i}{\omega }_{2}t}.\end{equation} \tag{ 3.47 }$$

Remembering that the solutions must be real, from the following decomposition of spatial coordinates, we have

$$\begin{equation}{h}_{\mu \nu }\left(x\right)=\frac{1}{{\left(2\pi \right)}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\enspace {\tilde {h}}_{\mu \nu }\left(t,\mathbf{k}\right){\text{e}}^{-\text{i}\mathbf{k}\cdot \mathbf{x}},\end{equation} \tag{ 3.48 }$$

and then we can reconstruct ${h}_{\mu \nu }\left(x\right)$ as

$$\begin{equation}\begin{array}{|l|}\hline \begin{aligned}{h}_{\mu \nu }\left(x\right)=\enspace & \frac{1}{{\left(2\pi \right)}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\enspace {C}_{\mu \nu }\left(\mathbf{k}\right){\text{e}}^{\text{i}{k}_{1}\cdot x}\\ \quad & +\frac{1}{{\left(2\pi \right)}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\left[\frac{1}{3}\left(\frac{{\eta }_{\mu \nu }}{2}+\frac{{\left({k}_{2}\right)}_{\mu }{\left({k}_{2}\right)}_{\nu }}{{k}_{2}^{2}}\right)\right]\tilde {A}\left(\mathbf{k}\right){\text{e}}^{\text{i}\enspace {k}_{2}\cdot x}+\mathrm{c}.\mathrm{c}.\end{aligned},\\ \hline\end{array}\end{equation} \tag{ 3.49 }$$

that is, the GWs in non-local linear gravity. For a similar approach in higher order gravity theories see [24]. It is worth noticing that the part of solution due to non-locality, related to ${k}_{2}^{2}\ne 0$, appears only if the constraint (3.14) is verified.

Equation (3.49) has the disadvantage of presenting ghost modes, so it is more useful to choose a suitable gauge in order to suppress waves without physical meaning. In this perspective, perturbing field equations (2.10)–(2.12) to first order in h_μν , δϕ and δλ, we get

$$\begin{equation}{\left(1+{a}_{1}\phi -\lambda \right)}^{(0)}{G}_{\mu \nu }^{(1)}+{\left({g}_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu }\right)}^{(0)}{\left(1+{a}_{1}\phi -\lambda \right)}^{(1)}={\kappa }^{2}{T}_{\mu \nu }^{(0)},\end{equation} \tag{ 3.50 }$$

that is

$$\begin{equation}\left(1+{a}_{1}{\phi }_{0}-{\lambda }_{0}\right){G}_{\mu \nu }^{(1)}+\left({\eta }_{\mu \nu }\square -{\partial }_{\mu }{\partial }_{\nu }\right)\left({a}_{1}\delta \phi -\delta \lambda \right)={\kappa }^{2}{T}_{\mu \nu }^{(0)},\end{equation} \tag{ 3.51 }$$

where □ = η^μν ∂_μ ∂_ν , and

$$\begin{equation}\square \delta \phi ={R}^{(1)},\end{equation} \tag{ 3.52 }$$

$$\begin{equation}\square \delta \lambda =-{a}_{1}{R}^{(1)}.\end{equation} \tag{ 3.53 }$$

In a particular coordinates frame $\left\{{x}^{\mu }\right\}$, being our equations gauge invariant, we define a new gauge as

$$\begin{equation}{\overline{h}}_{\mu \nu }={h}_{\mu \nu }-\frac{1}{2}{\eta }_{\mu \nu }h-\frac{{\eta }_{\mu \nu }}{1+{a}_{1}{\phi }_{0}-{\lambda }_{0}}\left({a}_{1}\delta \phi -\delta \lambda \right),\end{equation} \tag{ 3.54 }$$

such that, in our reference frame, it is

$$\begin{equation}{\partial }_{\mu }{\overline{h}}^{\mu \nu }=0,\end{equation} \tag{ 3.55 }$$

with 1 + a₁ ϕ₀ − λ₀ ≠ 0. The trace of equation (3.54) yields

$$\begin{equation}\overline{h}=-h-\frac{4}{1+{a}_{1}{\phi }_{0}-{\lambda }_{0}}\left({a}_{1}\delta \phi -\delta \lambda \right).\end{equation} \tag{ 3.56 }$$

In this gauge, we get the following expressions to first order for the Ricci tensor R_μν , the Ricci scalar R and the Einstein tensor G_μν

$$\begin{equation}{R}_{\mu \nu }^{(1)}=-\frac{1}{2}\square {h}_{\mu \nu }+\frac{1}{1+{a}_{1}{\phi }_{0}-{\lambda }_{0}}{\partial }_{\mu }{\partial }_{\nu }\left({a}_{1}\delta \phi -\delta \lambda \right),\end{equation} \tag{ 3.57 }$$

$$\begin{equation}{R}^{(1)}=-\frac{1}{2}\square h+\frac{1}{1+{a}_{1}{\phi }_{0}-{\lambda }_{0}}\square \left({a}_{1}\delta \phi -\delta \lambda \right),\end{equation} \tag{ 3.58 }$$

$$\begin{equation}{G}_{\mu \nu }^{(1)}=-\frac{1}{2}\square {h}_{\mu \nu }+\frac{1}{4}{\eta }_{\mu \nu }\square h+\frac{1}{1+{a}_{1}{\phi }_{0}-{\lambda }_{0}}\left({\partial }_{\mu }{\partial }_{\nu }-\frac{1}{2}{\eta }_{\mu \nu }\square \right)\left({a}_{1}\delta \phi -\delta \lambda \right).\end{equation} \tag{ 3.59 }$$

In terms of barred quantities ${\overline{h}}_{\mu \nu }$ and $\overline{h}$, we obtain

$$\begin{equation}{R}^{(1)}=\frac{1}{2}\square \overline{h}+\frac{3}{1+{a}_{1}{\phi }_{0}-{\lambda }_{0}}\square \left({a}_{1}\delta \phi -\delta \lambda \right),\end{equation} \tag{ 3.60 }$$

$$\begin{equation}{G}_{\mu \nu }^{(1)}=-\frac{1}{2}\square {\overline{h}}_{\mu \nu }-\frac{1}{1+{a}_{1}{\phi }_{0}-{\lambda }_{0}}\left({\eta }_{\mu \nu }\square -{\partial }_{\mu }{\partial }_{\nu }\right)\left({a}_{1}\delta \phi -\delta \lambda \right),\end{equation} \tag{ 3.61 }$$

which replaced in equations (3.51)–(3.53) give, in vacuum,

$$\begin{equation}\square {\overline{h}}_{\mu \nu }=0,\end{equation} \tag{ 3.62 }$$

$$\begin{equation}\square \delta \phi =\frac{3}{1+{a}_{1}{\phi }_{0}-{\lambda }_{0}}\left({a}_{1}\square \delta \phi -\square \delta \lambda \right),\end{equation} \tag{ 3.63 }$$

and

$$\begin{equation}\square \delta \lambda =-\frac{3{a}_{1}}{1+{a}_{1}{\phi }_{0}-{\lambda }_{0}}\left({a}_{1}\square \delta \phi -\square \delta \lambda \right),\end{equation} \tag{ 3.64 }$$

because equation (3.62) implies $\square \overline{h}=0$. Equations (3.63) and (3.64) can be rewritten as

$$\begin{equation}\left(\frac{1+{a}_{1}{\phi }_{0}-{\lambda }_{0}}{3}-{a}_{1}\right)\square \delta \phi +\square \delta \lambda =0,\end{equation} \tag{ 3.65 }$$

and

$$\begin{equation}{a}_{1}\square \delta \phi +\left(\frac{1+{a}_{1}{\phi }_{0}-{\lambda }_{0}}{3{a}_{1}}-1\right)\square \delta \lambda =0.\end{equation} \tag{ 3.66 }$$

The homogeneous linear system (3.65) and (3.66), in □δϕ and □δλ variables, admits trivial and non-trivial solutions depending on the determinant of the following matrix

$$\begin{align}& \mathrm{det}\left(\begin{matrix}\hfill \frac{1+{a}_{1}{\phi }_{0}-{\lambda }_{0}}{3}-{a}_{1}\hfill & \hfill 1\hfill \\ \hfill {a}_{1}\hfill & \hfill \frac{1+{a}_{1}{\phi }_{0}-{\lambda }_{0}}{3{a}_{1}}-1\hfill \end{matrix}\right)=\frac{1}{{a}_{1}}\left(1+{a}_{1}{\phi }_{0}-{\lambda }_{0}-6{a}_{1}\right)\\ & {\times}\left(1+{a}_{1}{\phi }_{0}-{\lambda }_{0}\right),\end{align} \tag{ 3.67 }$$

which can be equal to zero or non-equal to zero, depending on 1 + a₁ ϕ₀ − λ₀ ≠ 0. So, if the determinant does not vanish, under the constraint

$$\begin{equation}1+{a}_{1}{\phi }_{0}-{\lambda }_{0}\ne 6{a}_{1},\end{equation} \tag{ 3.68 }$$

we get the trivial solution

$$\begin{equation}\square \delta \phi =\square \delta \lambda =0,\end{equation} \tag{ 3.69 }$$

from which, it follows

$$\begin{equation}\square \left({a}_{1}\delta \phi -\delta \lambda \right)=0.\end{equation} \tag{ 3.70 }$$

A particular solution of equation (3.70) is

$$\begin{equation}{a}_{1}\delta \phi \left(x\right)=\delta \lambda \left(x\right).\end{equation} \tag{ 3.71 }$$

Putting this solution together with equations (3.62) and (3.54), we get massless transverse GWs with helicity 2, associated to ${k}_{1}^{2}=0$ under the constraint (3.68)

$$\begin{equation}\begin{array}{|l|}\hline {h}_{\mu \nu }^{({k}_{1})}\left(x\right)=\frac{1}{{\left(2\pi \right)}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\enspace {\tilde {C}}_{\mu \nu }\left(\mathbf{k}\right){\text{e}}^{\text{i}{k}_{1}\cdot x}+\mathrm{c}.\mathrm{c}.,\\ \hline\end{array}\end{equation} \tag{ 3.72 }$$

where ${\tilde {C}}_{\mu \nu }$ is the traceless polarization tensor, i.e. ${\tilde {C}}_{\enspace \sigma }^{\sigma }=0$, in a suitable gauge that leaves ${\partial }_{\mu }{\overline{h}}^{\mu \nu }=0$ invariant, as in the case (3.13).

The non-trivial solution is obtained by imposing that the determinant of the above matrix is equal to zero, which means

$$\begin{equation}1+{a}_{1}{\phi }_{0}-{\lambda }_{0}=6{a}_{1},\end{equation} \tag{ 3.73 }$$

that gives again

$$\begin{equation}{a}_{1}\square \delta \phi -\square \delta \lambda \ne 0,\end{equation} \tag{ 3.74 }$$

or equivalently

$$\begin{equation}\square h\ne 0,\end{equation} \tag{ 3.75 }$$

as in the case (3.14). Then, from the inverse gauge under the constraint (3.73)

$$\begin{equation}{h}_{\mu \nu }={\overline{h}}_{\mu \nu }-\frac{1}{2}{\eta }_{\mu \nu }\overline{h}-\frac{{\eta }_{\mu \nu }}{6{a}_{1}}\left({a}_{1}\delta \phi -\delta \lambda \right),\end{equation} \tag{ 3.76 }$$

we choose the non-homogeneous wave equation

$$\begin{equation}\frac{2}{3}\square \left({a}_{1}\delta \phi -\delta \lambda \right)=-{a}_{1}\square h=-{a}_{1}g\left(x\right),\end{equation} \tag{ 3.77 }$$

with $g\left(x\right)\in L({\mathbf{R}}^{4})$, as already defined in (3.29). The particular solution of equation (3.77) is

$$\begin{equation}{a}_{1}\delta \phi (x)-\delta \lambda (x)=\frac{1}{{(2\pi )}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\enspace \left(-\frac{3}{2}{a}_{1}\right)\tilde {A}\left(\mathbf{k}\right){\mathrm{e}}^{\mathrm{i}{k}_{2}\cdot x}+\mathrm{c}.\mathrm{c}.,\end{equation} \tag{ 3.78 }$$

with any ${k}_{2}^{2}={M}^{2}$. Inserting equation (3.78) into equation (3.76), we obtain a massive wave

$$\begin{equation}\begin{array}{|l|}\hline {h}_{\mu \nu }^{({k}_{2})}\left(x\right)=\frac{1}{{(2\pi )}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\enspace \frac{{\eta }_{\mu \nu }}{4}\tilde {A}\left(\mathbf{k}\right){\mathrm{e}}^{\mathrm{i}{k}_{2}\cdot x}+\mathrm{c}.\mathrm{c}..\\ \hline\end{array}\end{equation} \tag{ 3.79 }$$

In general under the constraint (3.73), the GWs for non-local f(R, □⁻¹ R) gravity are both massless and massive waves, that is

$$\begin{equation}\begin{array}{|l|}\hline \begin{aligned}{h}_{\mu \nu }\left(x\right)=\enspace & \frac{1}{{\left(2\pi \right)}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\enspace {\tilde {C}}_{\mu \nu }\left(\mathbf{k}\right){\text{e}}^{\text{i}{k}_{1}\cdot x} & +\enspace \frac{1}{{(2\pi )}^{3/2}}\int {\mathrm{d}}^{3}\mathbf{k}\enspace \frac{{\eta }_{\mu \nu }}{4}\tilde {A}\left(\mathbf{k}\right){\mathrm{e}}^{\mathrm{i}{k}_{2}\cdot x}+\mathrm{c}.\mathrm{c}.\end{aligned}.\\ \hline\end{array}\end{equation} \tag{ 3.80 }$$

Also here the solution component, related to ${k}_{2}^{2}\ne 0$, appears if the constraint (3.73) is satisfied. With these results in mind, let us investigate gravitational wave polarizations in non-local gravity.

A useful tool to study the polarization of gravitational radiation is the use of geodetic deviation produced by the wave when it invests a small region of spacetime, as the relative acceleration measured between nearby geodesics. Hence, we start from a wave ${h}_{\mu \nu }\left(t-{v}_{\mathrm{g}}z\right)$ propagating in $+\hat{z}$ direction in a local proper reference frame, where v_g is the group velocity in units where c = 1 defined as

$$\begin{equation}{v}_{\mathrm{g}}=\frac{\mathrm{d}\omega }{\mathrm{d}{k}_{z}}=\frac{{c}^{2}{k}_{z}}{\omega },\end{equation} \tag{ 4.1 }$$

and let us consider the equation for geodesic deviation

$$\begin{equation}{\ddot {x}}^{i}=-{R}_{\enspace 0k0}^{i}{x}^{k},\end{equation} \tag{ 4.2 }$$

where the Latin index range over the set $\left\{1,2,3\right\}$ and ${R}_{\enspace 0k0}^{i}$ are the only measurable components called the electric ones [25]. After replacing the linearized electric components of the Riemann tensor ${R}_{\enspace i0j0}^{\left(1\right)}$, expressed in terms of the metric perturbation h_μν

$$\begin{equation}{R}_{\enspace i0j0}^{\left(1\right)}=\frac{1}{2}\left({h}_{i0,j0}+{h}_{j0,i0}-{h}_{ij,00}-{h}_{00,ij}\right),\end{equation} \tag{ 4.3 }$$

in equation (4.2), we get a linear non-homogeneous system of differential equations

$$\begin{equation}\begin{cases}\ddot {x}(t)=-\frac{1}{2}{h}_{11,00}x-\frac{1}{2}{h}_{12,00}y+\frac{1}{2}\left({h}_{10,03}-{h}_{13,00}\right)z\\ \ddot {y}(t)=-\frac{1}{2}{h}_{12,00}x-\frac{1}{2}{h}_{22,00}y+\frac{1}{2}\left({h}_{02,03}-{h}_{23,00}\right)z\\ \ddot {z}(t)=\frac{1}{2}\left({h}_{01,03}-{h}_{13,00}\right)x+\frac{1}{2}\left({h}_{02,03}-{h}_{23,00}\right)y\\ \quad \quad +\enspace \frac{1}{2}\left(2{h}_{03,03}-{h}_{33,00}-{h}_{00,33}\right)z.\end{cases}\end{equation} \tag{ 4.4 }$$

For a massless plane wave traveling in $+\hat{z}$ direction, that is ${k}_{1}^{2}=0$, which propagates at speed c, if we keep k fixed and ${k}_{1}^{\mu }=\left({\omega }_{1},0,0,{k}_{z}\right)$, equation (3.72) yields

$$\begin{equation}{h}_{\mu \nu }^{({k}_{1})}\left(t,z\right)=\sqrt{2}\left[{\tilde {{\epsilon}}}^{(+)}\left({\omega }_{1}\right){{\epsilon}}_{\mu \nu }^{(+)}+{\tilde {{\epsilon}}}^{({\times})}\left({\omega }_{1}\right){{\epsilon}}_{\mu \nu }^{({\times})}\right]{\text{e}}^{\text{i}{\omega }_{1}\left(t-z\right)}+\mathrm{c}.\mathrm{c}.,\end{equation} \tag{ 4.5 }$$

where

$$\begin{equation}{{\epsilon}}_{\mu \nu }^{(+)}=\frac{1}{\sqrt{2}}\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),\end{equation} \tag{ 4.6 }$$

$$\begin{equation}{{\epsilon}}_{\mu \nu }^{({\times})}=\frac{1}{\sqrt{2}}\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),\end{equation} \tag{ 4.7 }$$

and ω₁ = k_z . Furthermore, for a massive plane wave propagating in $+\hat{z}$ direction, that is ${k}_{2}^{2}\ne 0$, always keeping k fixed and instead with ${k}_{2}^{\mu }=\left({\omega }_{2},0,0,{k}_{z}\right)$, equation (3.79) becomes

$$\begin{equation}{h}_{\mu \nu }^{({k}_{2})}\left(t,z\right)=\frac{\tilde {A}\left({k}_{z}\right)}{4}{\eta }_{\mu \nu }\enspace {\text{e}}^{\text{i}\left({\omega }_{2}t-{k}_{z}z\right)}+\mathrm{c}.\mathrm{c}.,\end{equation} \tag{ 4.8 }$$

where here the propagation speed is less than c. In a more compact form, the metric linear perturbation h_μν , traveling in the $+\hat{z}$ direction, assuming k fixed, is

$$\begin{align}{h}_{\mu \nu }\left(t,z\right)& =\sqrt{2}\left[{\tilde {{\epsilon}}}^{(+)}\left({\omega }_{1}\right){{\epsilon}}_{\mu \nu }^{(+)}+{\tilde {{\epsilon}}}^{({\times})}\left({\omega }_{1}\right){{\epsilon}}_{\mu \nu }^{({\times})}\right]{\text{e}}^{\text{i}{\omega }_{1}\left(t-z\right)}+{\tilde {{\epsilon}}}_{\mu \nu }^{\left(s\right)}\left({k}_{z}\right){\text{e}}^{\text{i}\left({\omega }_{2}t-{k}_{z}z\right)}+\mathrm{c}.\mathrm{c}.,\end{align} \tag{ 4.9 }$$

where ${\tilde {{\epsilon}}}_{\mu \nu }^{\left(s\right)}$ is the polarization tensor associated to the mixed scalar mode

$$\begin{equation}{\tilde {{\epsilon}}}_{\mu \nu }^{\left(s\right)}\left({k}_{z}\right)=\left({{\epsilon}}_{\mu \nu }^{(TT)}-\sqrt{2}{{\epsilon}}_{\mu \nu }^{(b)}-{{\epsilon}}_{\mu \nu }^{(l)}\right)\frac{\tilde {A}\left({k}_{z}\right)}{4}.\end{equation} \tag{ 4.10 }$$

and the polarization tensors are explicitly given by

$$\begin{align}{{\epsilon}}_{\mu \nu }^{(TT)}& =\left(\begin{matrix}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),\qquad {{\epsilon}}_{\mu \nu }^{(b)}=\frac{1}{\sqrt{2}}\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right)\enspace ,\\ {{\epsilon}}_{\mu \nu }^{(l)}& =\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right),\end{align} \tag{ 4.11 }$$

where the set of polarization tensors $\left\{{{\epsilon}}_{\mu \nu }^{\left(+\right)},{{\epsilon}}_{\mu \nu }^{\left({\times}\right)},{{\epsilon}}_{\mu \nu }^{\left(TT\right)},{{\epsilon}}_{\mu \nu }^{\left(b\right)},{{\epsilon}}_{\mu \nu }^{\left(l\right)}\right\}$ satisfy the orthonormality relations

$$\begin{equation}\mathrm{Tr}\left\{{{\epsilon}}^{(a)}{{\epsilon}}^{(b)}\right\}={{\epsilon}}_{\mu \nu }^{(a)}{{\epsilon}}^{(b)\mu \nu }={\delta }^{a,b}\quad \text{with}\quad a,b\in \left\{+,{\times},TT,b,l\right\}.\end{equation} \tag{ 4.12 }$$

The scalar mode ${\tilde {{\epsilon}}}_{\mu \nu }^{\left(s\right)}$ is a mixed state obtained from a combination of longitudinal and transverse scalar modes which is produced by the single degree of freedom $\tilde {A}$ as for f(R) gravity which has three d.o.f.: ${\tilde {{\epsilon}}}^{(+)}$, ${\tilde {{\epsilon}}}^{({\times})}$ and $\tilde {A}$ [26, 27]. However, as we will see, the transverse component weighs more than the longitudinal one. In fact, the polarization tensor ${\tilde {{\epsilon}}}_{\mu \nu }^{\left(s\right)}$, restricted to spatial components ${\tilde {{\epsilon}}}_{i,j}^{\left(s\right)}$, is provided by

$$\begin{equation}{\tilde {{\epsilon}}}_{i,j}^{\left(s\right)}=\left(\sqrt{2}{{\epsilon}}_{i,j}^{(b)}+{{\epsilon}}_{i,j}^{(l)}\right)\frac{\tilde {A}\left({k}_{z}\right)}{4},\end{equation} \tag{ 4.13 }$$

where (i, j) range over (1, 2, 3). Hence, from equation (4.4) in the case of massless plane waves ${h}_{\mu \nu }^{({k}_{1})}$, we have

$$\begin{equation}\begin{cases}\ddot {x}(t)=\frac{1}{2}{\omega }_{1}^{2}\left[{\tilde {{\epsilon}}}^{\left(+\right)}\left({\omega }_{1}\right)x+{\tilde {{\epsilon}}}^{\left({\times}\right)}\left({\omega }_{1}\right)y\right]{\text{e}}^{\text{i}{\omega }_{1}\left(t-z\right)}+\mathrm{c}.\mathrm{c}.\quad \\ \ddot {y}(t)=\frac{1}{2}{\omega }_{1}^{2}\left[{\tilde {{\epsilon}}}^{\left({\times}\right)}\left({\omega }_{1}\right)x-{\tilde {{\epsilon}}}^{\left(+\right)}\left({\omega }_{1}\right)y\right]{\text{e}}^{\text{i}{\omega }_{1}\left(t-z\right)}+\mathrm{c}.\mathrm{c}.\quad \\ \ddot {z}(t)=0\quad \end{cases},\end{equation} \tag{ 4.14 }$$

that give us the two standard transverse tensor polarizations predicted by general relativity, conventionally called plus and cross modes.

Otherwise, in the case of a massive plane wave ${h}_{\mu \nu }^{({k}_{2})}$ with ${M}^{2}={\omega }_{2}^{2}-{k}_{z}^{2}$, the linearized geodesic deviation equation (4.4) becomes

$$\begin{equation}\begin{cases}\ddot {x}(t)=-\frac{1}{8}{\omega }_{2}^{2}\tilde {A}\left({k}_{z}\right)x\enspace {\text{e}}^{\text{i}\left({\omega }_{2}t-{k}_{z}z\right)}+\mathrm{c}.\mathrm{c}.\quad \\ \ddot {y}(t)=-\frac{1}{8}{\omega }_{2}^{2}\tilde {A}\left({k}_{z}\right)y\enspace {\text{e}}^{\text{i}\left({\omega }_{2}t-{k}_{z}z\right)}+\mathrm{c}.\mathrm{c}.\quad \\ \ddot {z}(t)=-\frac{1}{8}{M}^{2}\tilde {A}\left({k}_{z}\right)z\enspace {\mathrm{e}}^{\mathrm{i}\left({\omega }_{2}t-{k}_{z}z\right)}+\mathrm{c}.\mathrm{c}.\quad \end{cases},\end{equation} \tag{ 4.15 }$$

that can be integrated, assuming that ${h}_{\mu \nu }\left(t,z\right)$ is small, as

$$\begin{equation}\begin{cases}x(t)=x(0)+\frac{1}{8}\tilde {A}\left({k}_{z}\right)x(0){\text{e}}^{\text{i}\left({\omega }_{2}t-{k}_{z}z\right)}+\mathrm{c}.\mathrm{c}.\quad \\ y(t)=y(0)+\frac{1}{8}\tilde {A}\left({k}_{z}\right)y(0){\text{e}}^{\text{i}\left({\omega }_{2}t-{k}_{z}z\right)}+\mathrm{c}.\mathrm{c}.\quad \\ z(t)=z(0)+\frac{1}{8{\omega }_{2}^{2}}{M}^{2}\tilde {A}\left({k}_{z}\right)z(0){\text{e}}^{\text{i}\left({\omega }_{2}t-{k}_{z}z\right)}+\mathrm{c}.\mathrm{c}.\quad \end{cases}.\end{equation} \tag{ 4.16 }$$

If we suppose that M² is very small, which happens when we are sufficiently far from the radiation source, and that becomes zero for exactly null plane waves, k² = 0, then we can expand our results with respect to a parameter γ defined as [28]

$$\begin{equation}\gamma ={\left(\frac{c}{{v}_{\mathrm{g}}}\right)}^{2}-1,\end{equation} \tag{ 4.17 }$$

which takes into account the difference in speed between nearly null waves with speed v_g and null ones with speed c. Thus, keeping k_z fixed and by using Landau symbols, namely little-o and big-$\mathcal{O}$ notation, we can expand in terms of our parameter γ the following quantities

$$\begin{equation}\frac{{M}^{2}}{{k}_{z}^{2}}=\frac{1}{{k}_{z}^{2}}\left[\frac{{\omega }_{2}^{2}}{{c}^{2}}-{k}_{z}^{2}\right]={\left(\frac{{\omega }_{2}}{c{k}_{z}}\right)}^{2}-1=\gamma ,\end{equation} \tag{ 4.18 }$$

and in units where c = 1

$$\begin{equation}\frac{{\omega }_{2}}{{k}_{z}}=\sqrt{1+\gamma }=\left(1+\frac{1}{2}\gamma \right)+o\left(\gamma \right),\end{equation} \tag{ 4.19 }$$

$$\begin{equation}{\text{e}}^{\text{i}\left({\omega }_{2}t-{k}_{z}z\right)}={\mathrm{e}}^{\mathrm{i}{k}_{z}\left(t-z\right)}+\mathcal{O}\left(\gamma \right),\end{equation} \tag{ 4.20 }$$

that implies

$$\begin{equation}\frac{{M}^{2}}{{\omega }_{2}^{2}}=\gamma +o\left(\gamma \right).\end{equation} \tag{ 4.21 }$$

Within the first order in γ, the solution (4.16) for the mode ω₂ give us

$$\begin{equation}\begin{cases}x(t)=x(0)+\frac{1}{8}\tilde {A}\left({k}_{z}\right)x(0){\text{e}}^{\text{i}\enspace {k}_{z}\left(t-z\right)}+\mathcal{O}\left(\gamma \right)+\mathrm{c}.\mathrm{c}.\quad \\ y(t)=y(0)+\frac{1}{8}\tilde {A}\left({k}_{z}\right)y(0){\text{e}}^{\text{i}\enspace {k}_{z}\left(t-z\right)}+\mathcal{O}\left(\gamma \right)+\mathrm{c}.\mathrm{c}.\quad \\ z(t)=z(0)+\frac{1}{8}\gamma \tilde {A}\left({k}_{z}\right)z(0){\text{e}}^{\text{i}\enspace {k}_{z}\left(t-z\right)}+\mathcal{O}\left({\gamma }^{2}\right)+\mathrm{c}.\mathrm{c}.\quad \end{cases},\end{equation} \tag{ 4.22 }$$

that is

$$\begin{equation}{\Delta}z(t)=o\left({\Delta}x(t)\right)\quad \mathrm{f}\mathrm{o}\mathrm{r}\enspace \gamma \to 0,\end{equation} \tag{ 4.23 }$$

and then only the breathing tensor polarization ^(b) survives because longitudinal modes are infinitesimal in higher order than transverse modes when γ tends to zero. Thus when a GW strikes a sphere of particles of radius $r=\sqrt{{x}^{2}(0)+{y}^{2}(0)+{z}^{2}(0)}$, this will be distorted into an ellipsoid described by

$$\begin{equation}{\left(\frac{x}{{\rho }_{1}(t)}\right)}^{2}+{\left(\frac{y}{{\rho }_{1}(t)}\right)}^{2}+{\left(\frac{z}{{\rho }_{2}(t)}\right)}^{2}={r}^{2},\end{equation} \tag{ 4.24 }$$

where ${\rho }_{1}(t)=1+\frac{1}{4}\tilde {A}\left({k}_{z}\right)\mathrm{cos}\left[{k}_{z}\left(t-z\right)+\phi \right]+\mathcal{O}\left(\gamma \right)$ and ${\rho }_{2}(t)=1+\mathcal{O}\left(\gamma \right)$, that, at zero order in γ only ρ₁, is varying between their maximum and minimum values. This ellipsoid swings only on xy-plane between two circumferences of minimum and maximum radius and represents an additional transverse scalar polarization which has zero helicity within the lowest order in γ [29].

According to these considerations, the d.o.f. of non-local f(R, □⁻¹ R) gravity are three: two of these, ${\tilde {{\epsilon}}}^{\left(+\right)}$ and ${\tilde {{\epsilon}}}^{\left({\times}\right)}$, give rise to the standard tensor modes of general relativity while the degree of freedom $\tilde {A}$ generates a further breathing scalar mode. In summary, f(R, □⁻¹ R) gravity has three polarizations, namely two massless two-helicity tensor modes and one massive zero-helicity scalar mode, all purely transverse within the lowest order in γ, exactly like f(R) gravity (see for a discussion [27, 30–35]).

A further approach to study polarizations can be obtained by adopting the NP formalism for low-mass GWs [28]. Even if it is not directly applicable to massive waves because it was, in origin, worked out for massless waves, it is possible to generalize it to low-mass waves propagating along nearly null geodesics [36]. It is worth noticing that the little group $E\left(2\right)$ classification fails for massive waves but can be recovered within the first order in γ.

Let us introduce a new basis, namely a local quasi-orthonormal null tetrad basis $\left(k,l,m,\overline{m}\right)$ defined as [29, 37, 38]

$$\begin{equation}k={\partial }_{t}+{\partial }_{z},\qquad l=\frac{1}{2}\left({\partial }_{t}-{\partial }_{z}\right),\end{equation} \tag{ 5.1 }$$

$$\begin{equation}m=\frac{1}{\sqrt{2}}\left({\partial }_{x}+\mathrm{i}{\partial }_{y}\right),\qquad \overline{m}=\frac{1}{\sqrt{2}}\left({\partial }_{x}-\mathrm{i}{\partial }_{y}\right),\end{equation} \tag{ 5.2 }$$

which, adopting the Minkowski metric tensor η_μν of signature −2, satisfies the relations

$$\begin{equation}\begin{gathered}k\cdot l=-m\cdot \overline{m}=1\enspace ,\\ k\cdot k=l\cdot l=m\cdot m=\overline{m}\cdot \overline{m}=0\enspace ,\\ k\cdot m=k\cdot \overline{m}=l\cdot m=l\cdot \overline{m}=0\enspace ,\end{gathered}\end{equation} \tag{ 5.3 }$$

that is

$$\begin{equation}{\eta }^{\mu \nu }=2{k}^{\left(\right.\mu }{l}^{\nu \left.\right)}-2{m}^{\left(\right.\mu }{\overline{m}}^{\nu \left.\right)}.\end{equation} \tag{ 5.4 }$$

Therefore we can raise and lower the tetrad indices by the metric of the tetrad η_ab

$$\begin{equation}{\eta }_{ab}={\eta }^{ab}=\left(\begin{matrix}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \end{matrix}\right),\end{equation} \tag{ 5.5 }$$

where (a, b) run over $(k,l,m,\overline{m})$. We now split the Riemann tensor into three irreducible parts, namely: Weyl tensor, traceless Ricci tensor and Ricci scalar, known as NP quantities. The four-dimensional Weyl tensor C_μνρσ is defined as

$$\begin{equation}{C}_{\mu \nu \rho \sigma }={R}_{\mu \nu \rho \sigma }-2{g}_{\left[\right.\mu \vert \left[\right.\rho }{R}_{\sigma \left.\right]\vert \nu \left.\right]}+\frac{1}{3}{g}_{\mu \left[\right.\rho }{g}_{\sigma \left.\right]\nu }R,\end{equation} \tag{ 5.6 }$$

and in tetrad form becomes

$$\begin{equation}{C}_{abcd}={R}_{abcd}-2{\eta }_{\left[\right.a\vert \left[\right.c}{R}_{d\left.\right]\vert b\left.\right]}+\frac{1}{3}{\eta }_{a\left[\right.c}{\eta }_{d\left.\right]b}R,\end{equation} \tag{ 5.7 }$$

taking into account that the tetrad components of the generic tensor P_abcd... express in terms of the local coordinate basis are

$$\begin{equation}{P}_{abcd\dots }={P}_{\alpha ,\beta \gamma \delta \dots }{a}^{\alpha }{b}^{\beta }{c}^{\gamma }{d}^{\delta }\dots ,\end{equation} \tag{ 5.8 }$$

where (a, b, c, d, ...) run over $(k,l,m,\overline{m})$. The fifteen NP-amplitudes are specifically the five complex Weyl-NP Ψ scalars, expressed in tetrad components of the Weyl tensor as

$$\begin{align}{{\Psi}}_{0}\equiv -{C}_{kmkm},\\ {{\Psi}}_{1}\equiv -{C}_{klkm},\\ {{\Psi}}_{2}\equiv -{C}_{km\overline{m}l},\\ {{\Psi}}_{3}\equiv -{C}_{kl\overline{m}l},\\ {{\Psi}}_{4}\equiv -{C}_{\overline{m}l\overline{m}l},\end{align} \tag{ 5.9 }$$

and the ten Ricci-NP scalars Φ, Λ, expressed in tetrad components of Ricci tensor as

$$\begin{equation*}{{\Phi}}_{02}\equiv -\frac{1}{2}{R}_{mm,}\end{equation*}$$

$$\begin{equation*}\begin{cases}{{\Phi}}_{01}\equiv -\frac{1}{2}{R}_{km}\quad \\ {{\Phi}}_{12}\equiv -\frac{1}{2}{R}_{lm}\quad \end{cases},\end{equation*}$$

$$\begin{equation*}\begin{cases}{{\Phi}}_{00}\equiv -\frac{1}{2}{R}_{kk}\quad \\ {{\Phi}}_{11}\equiv -\frac{1}{4}\left({R}_{kl}+{R}_{m\overline{m}}\right)\quad \\ {{\Phi}}_{22}\equiv -\frac{1}{2}{R}_{ll}\quad \end{cases},\end{equation*}$$

$$\begin{equation*}\begin{cases}{{\Phi}}_{10}\equiv -\frac{1}{2}{R}_{k\overline{m}}={{\Phi}}_{01}^{{\ast}}\quad \\ {{\Phi}}_{21}\equiv -\frac{1}{2}{R}_{l\overline{m}}={{\Phi}}_{21}^{{\ast}}\quad \end{cases},\end{equation*}$$

$$\begin{align}{{\Phi}}_{20}\equiv -\frac{1}{2}{R}_{\overline{m}\overline{m}}={{\Phi}}_{02}^{{\ast}}\enspace ,\\ {\Lambda}=\frac{R}{24}.\end{align} \tag{ 5.10 }$$

In order to expand the low-mass gravitational radiation in terms of null plane waves, we first define the 'wave' four-vector ${\left({k}_{2}^{\prime }\right)}^{\mu }$ in units where c = 1 associated to the nearly null plane wave propagating in positive z direction as

$$\begin{equation}{k}_{2}^{\mu }=\left({\omega }_{2},0,0,{k}_{z}\right)={\omega }_{2}{k}_{2}^{\prime \mu },\end{equation} \tag{ 5.11 }$$

with

$$\begin{equation}{k}_{2}^{\prime \mu }=\left(1,0,0,{v}_{\mathrm{g}}\right),\end{equation} \tag{ 5.12 }$$

and we set the time retarded $\tilde {u}$ as

$$\begin{equation}\tilde {u}=t-{v}_{\mathrm{g}}z,\end{equation} \tag{ 5.13 }$$

that gives us

$$\begin{equation}{\nabla }_{\mu }\tilde {u}={\left({k}_{2}^{\prime }\right)}_{\mu }.\end{equation} \tag{ 5.14 }$$

Now, we expand k₂' with respect to our null tetrads basis

$$\begin{equation}{k}_{2}^{\prime \mu }=\left(1+{\gamma }_{k}\right){k}^{\mu }+{\gamma }_{l}{l}^{\mu }+{\gamma }_{m}{m}^{\mu }+{\overline{\gamma }}_{m}{\overline{m}}^{\mu },\end{equation} \tag{ 5.15 }$$

where the expansion coefficients γ_k , γ_l , γ_m are of same order of the γ of previous section. Given the arbitrariness of the observer to orient its reference system, it is possible to perform orientation in such a way that ${k}_{2}^{0}={k}^{0}$ and ${k}_{2}^{3}\propto {k}^{3}$, where k⁰ and k³ are the angular frequency and the third component of vector wave of its null wave, respectively. Therefore we obtain, from equation (5.15), γ_l = −2γ_k and γ_m = 0 that yields

$$\begin{equation}{k}_{2}^{\prime \mu }={k}^{\mu }+{\gamma }_{k}\left({k}^{\mu }-2{l}^{\mu }\right),\end{equation} \tag{ 5.16 }$$

or

$$\begin{equation}{k}_{2}^{\prime \mu }=\left(1,0,0,1+2{\gamma }_{k}\right),\end{equation} \tag{ 5.17 }$$

and, as already observed according to

$$\begin{equation}{\gamma }_{k}=-\frac{1}{4}\gamma +o\left(\gamma \right),\end{equation} \tag{ 5.18 }$$

the parameters γ and γ_k are of the same order. The derivatives of Riemann tensor can be expressed as

$$\begin{equation}{R}_{\alpha \beta \gamma \delta ,\mu }=\frac{\partial {R}_{\alpha \beta \gamma \delta }}{\partial \tilde {u}}{\nabla }_{\mu }\tilde {u}={\left({k}_{2}^{\prime }\right)}_{\mu }{\dot {R}}_{\alpha \beta \gamma \delta },\end{equation} \tag{ 5.19 }$$

where the superscripted dot means the durative with respect to $\tilde {u}$. The following identities

$$\begin{equation}{R}_{\alpha \beta \gamma \delta ,k}=-2{\gamma }_{k}{\dot {R}}_{\alpha \beta \gamma \delta },\end{equation} \tag{ 5.20 }$$

$$\begin{equation}{R}_{\alpha \beta \gamma \delta ,l}=\left(1+{\gamma }_{k}\right){\dot {R}}_{\alpha \beta \gamma \delta },\end{equation} \tag{ 5.21 }$$

$$\begin{equation}{R}_{\alpha \beta \gamma \delta ,m}=0,\end{equation} \tag{ 5.22 }$$

combined with differential Bianchi identity

$$\begin{equation}{R}_{ab[cd,e]}=0,\end{equation} \tag{ 5.23 }$$

involve that the only non-zero tetrad components of Riemann tensor R_αβγδ to zero order in γ_k are terms of form R_lplq with (p, q) range over $(k,m,\overline{m})$. In the nearly null plane waves framework, only four complex NP tetrad components are independent and non-vanishing within the first order in γ, that is, from equation (5.7), they are

$$\begin{equation}{{\Psi}}_{2}\left(\tilde {u}\right)=-\frac{1}{6}{R}_{lklk}+\mathcal{O}\left({\gamma }_{k}\right),\end{equation} \tag{ 5.24 }$$

$$\begin{equation}{{\Psi}}_{3}\left(\tilde {u}\right)=-\frac{1}{2}{R}_{lkl\overline{m}}+\mathcal{O}\left({\gamma }_{k}\right),\end{equation} \tag{ 5.25 }$$

$$\begin{equation}{{\Psi}}_{4}\left(\tilde {u}\right)=-{R}_{l\overline{m}l\overline{m}},\end{equation} \tag{ 5.26 }$$

$$\begin{equation}{{\Phi}}_{22}\left(\tilde {u}\right)=-{R}_{lml\overline{m}}.\end{equation} \tag{ 5.27 }$$

In terms of tetrad components of metric perturbation h_ab , they become

$$\begin{equation}{{\Psi}}_{2}\left(\tilde {u}\right)=\frac{1}{12}{\ddot {h}}_{kk}+\mathcal{O}\left({\gamma }_{k}\right),\end{equation} \tag{ 5.28 }$$

$$\begin{equation}{{\Psi}}_{3}\left(\tilde {u}\right)=\frac{1}{4}{\ddot {h}}_{k\overline{m}}+\mathcal{O}\left({\gamma }_{k}\right),\end{equation} \tag{ 5.29 }$$

$$\begin{equation}{{\Psi}}_{4}\left(\tilde {u}\right)=\frac{1}{2}{\ddot {h}}_{\overline{m}\overline{m}}+\mathcal{O}\left({\gamma }_{k}\right),\end{equation} \tag{ 5.30 }$$

$$\begin{equation}{{\Phi}}_{22}\left(\tilde {u}\right)=\frac{1}{2}{\ddot {h}}_{m\overline{m}}+\mathcal{O}\left({\gamma }_{k}\right).\end{equation} \tag{ 5.31 }$$

These real amplitudes, under the subgroup of Lorentz transformations which leaves k₂ unchanged, namely the little group E(2), show the following helicity values s:

$$\begin{equation}{{\Psi}}_{2}\quad s=0,\qquad {{\Psi}}_{4\quad }s=2,\end{equation} \tag{ 5.32 }$$

$$\begin{equation}{{\Psi}}_{3}\quad s=1,\qquad {{\Phi}}_{22}\quad s=0.\end{equation} \tag{ 5.33 }$$

They allow the Petrov classification for a set of quasi-Lorentz invariant GWs. Remembering that our gravitational radiation travels along the positive $+\hat{z}$ axis, we obtain the following four NP amplitudes expressed both in terms of the electric components of the Riemann tensor R_i0j0 and its linearized components [39–42]. That is, taking into account the identities

$$\begin{equation}{R}_{kml\overline{m}}=0\to \quad {R}_{1313}+{R}_{2323}={R}_{0101}+{R}_{0202}\end{equation} \tag{ 5.34 }$$

$$\begin{equation}{R}_{kmk\overline{m}}=0\to \quad -2{R}_{0131}-2{R}_{0232}=2\left({R}_{0101}+{R}_{0202}\right)\end{equation} \tag{ 5.35 }$$

$$\begin{equation}{R}_{klkm}=0\to \quad {R}_{0301}=-{R}_{0331}\quad \text{and}\quad {R}_{0302}=-{R}_{0332}\end{equation} \tag{ 5.36 }$$

$$\begin{equation}{R}_{kmlm}=0\to \quad {R}_{3132}={R}_{0102}\quad \text{and}\quad {R}_{3131}-{R}_{3232}={R}_{0101}-{R}_{0202}\end{equation} \tag{ 5.37 }$$

$$\begin{equation}{R}_{kmkm}=0\to \quad \begin{cases}_{3101}-2{R}_{3202}=-{R}_{3131}+{R}_{3232}-{R}_{0101}+{R}_{0202}\quad \\ {R}_{0102}=-{R}_{0132}-{R}_{0231}-{R}_{3132},\quad \end{cases}\end{equation} \tag{ 5.38 }$$

we get, in terms of R_i0j0,

$$\begin{equation}{{\Psi}}_{2}\left(\tilde {u}\right)=-\frac{1}{6}{R}_{0303}+\mathcal{O}\left({\gamma }_{k}\right),\end{equation} \tag{ 5.39 }$$

$$\begin{equation}{{\Psi}}_{3}\left(\tilde {u}\right)=-\frac{1}{2\sqrt{2}}{R}_{0301}+\frac{1}{2\sqrt{2}}\mathrm{i}{R}_{0302}+\mathcal{O}\left({\gamma }_{k}\right),\end{equation} \tag{ 5.40 }$$

$$\begin{equation}{{\Psi}}_{4}\left(\tilde {u}\right)=-\frac{1}{2}{R}_{0101}+\frac{1}{2}{R}_{0202}+\mathrm{i}{R}_{0102}+\mathcal{O}\left({\gamma }_{k}\right),\end{equation} \tag{ 5.41 }$$

$$\begin{equation}{{\Phi}}_{22}\left(\tilde {u}\right)=-\frac{1}{2}{R}_{0101}-\frac{1}{2}{R}_{0202}+\mathcal{O}\left({\gamma }_{k}\right),\end{equation} \tag{ 5.42 }$$

while, in term of metric perturbation, the NP scalars take the form

$$\begin{equation}{{\Psi}}_{2}\left(\tilde {u}\right)=-\frac{1}{12}\left(2{h}_{03,03}-{h}_{00,33}-{h}_{33,00}\right)+\mathcal{O}\left({\gamma }_{k}\right),\end{equation} \tag{ 5.43 }$$

$$\begin{equation}{{\Psi}}_{3}\left(\tilde {u}\right)=-\frac{1}{4\sqrt{2}}\left({h}_{01,03}-{h}_{13,00}\right)+\frac{1}{4\sqrt{2}}\mathrm{i}\left({h}_{02,03}-{h}_{23,00}\right)+\mathcal{O}\left({\gamma }_{k}\right),\end{equation} \tag{ 5.44 }$$

$$\begin{equation}{{\Psi}}_{4}\left(\tilde {u}\right)=\frac{1}{4}\left({h}_{11,00}-{h}_{22,00}\right)-2\mathrm{i}{h}_{12,00}+\mathcal{O}\left({\gamma }_{k}\right),\end{equation} \tag{ 5.45 }$$

$$\begin{equation}{{\Phi}}_{22}\left(\tilde {u}\right)=\frac{1}{4}\left({h}_{11,00}+{h}_{22,00}\right)+\mathcal{O}\left({\gamma }_{k}\right).\end{equation} \tag{ 5.46 }$$

From equations (4.5) and (4.8) for non-null and null geodesic congruences of GWs, we get, for a massless mode ω₁ and a massive mode ω₂, at k fixed at the lowest order in γ_k , the expressions

$$\begin{align*}{{\Psi}}_{2}\left(\tilde {u}\right)& =\mathcal{O}\left({\gamma }_{k}\right),\\ {{\Psi}}_{3}\left(\tilde {u}\right)& =\mathcal{O}\left({\gamma }_{k}\right)\\ {{\Psi}}_{4}\left(\tilde {u}\right)& =-{\omega }_{1}^{2}\left[\left({{\epsilon}}^{(+)}\left({\omega }_{1}\right){\text{e}}^{\text{i}{\omega }_{1}\left(t-z\right)}+\mathrm{c}.\mathrm{c}.\right)-\mathrm{i}\left({{\epsilon}}^{({\times})}\left({\omega }_{1}\right){\text{e}}^{\text{i}{\omega }_{1}\left(t-z\right)}+\mathrm{c}.\mathrm{c}.\right)\right]+\mathcal{O}\left({\gamma }_{k}\right)\\ {{\Phi}}_{22}\left(\tilde {u}\right)& =\frac{{k}_{z}^{2}}{8}A\left({k}_{z}\right){\text{e}}^{\text{i}{k}_{z}\left(t-z\right)}+\mathcal{O}\left({\gamma }_{k}\right)+\mathrm{c}.\mathrm{c}.,\end{align*}$$

considering

$$\begin{equation}{\omega }_{2}={k}_{z}+\mathcal{O}\left({\gamma }_{k}\right),\end{equation} \tag{ 5.47 }$$

in units c = 1. Hence, being Ψ₂ = Ψ₃ = 0 and Φ₂₂ ≠ 0 at zero order in γ_k , the quasi-Lorentz invariant E(2) class of non-local gravitational theory R + a₁ R□⁻¹ R is N₃, according to the Petrov classification. Here the presence or absence of all modes is observer-independent. The driving-force matrix S(t) can be expressed in terms of the six new basis polarization matrices W_A (e_z ) along the wave direction $\hat{\mathbf{k}}={\mathbf{e}}_{z}$ as

$$\begin{equation}S\left(t\right)=\sum\limits _{A}{p}_{A}\left({\mathbf{e}}_{z},t\right){W}_{A}\left({\mathbf{e}}_{z}\right),\end{equation} \tag{ 5.48 }$$

where the index A ranges over {1, 2, 3, 4, 5, 6} [37, 38] and explicitly

$$\begin{align}{W}_{1}\left({\mathbf{e}}_{z}\right)=& -6\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right),& {W}_{2}\left({\mathbf{e}}_{z}\right)=& -2\sqrt{2}\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),\\ {W}_{3}\left({\mathbf{e}}_{z}\right)=\enspace & 2\sqrt{2}\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \end{matrix}\right),& {W}_{4}\left({\mathbf{e}}_{z}\right)=& -\left(\begin{matrix}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),\\ {W}_{5}\left({\mathbf{e}}_{z}\right)=& \left(\begin{matrix}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),& {W}_{6}\left({\mathbf{e}}_{z}\right)=& -\left(\begin{matrix}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right).\end{align} \tag{ 5.49 }$$

Consequently, there are six polarizations modes: the longitudinal mode ${p}_{1}^{\left(l\right)}$, the vector-x mode ${p}_{2}^{\left(x\right)}$, the vector-y mode ${p}_{3}^{\left(y\right)}$, the plus mode ${p}_{4}^{\left(+\right)}$, the cross mode ${p}_{5}^{\left({\times}\right)}$, and the breathing mode ${p}_{6}^{\left(b\right)}$. Here ${p}_{A}\left({\mathbf{e}}_{z},t\right)$ are the amplitudes of the wave at the detector in the frame origin [30–32, 43, 44]. Taking into account that the spatial components of matrix S(t) are the electric components of Riemann tensor, that is

$$\begin{equation}{S}_{ij}(t)={R}_{i0j0},\end{equation} \tag{ 5.50 }$$

we can get the six polarization amplitudes ${p}_{A}\left({\mathbf{e}}_{z},t\right)$ in terms of the NP scalars for our non-local waves

$$\begin{align}{p}_{1}^{\left(l\right)}\left({\mathbf{e}}_{z},t\right)& ={{\Psi}}_{2}=\mathcal{O}\left({\gamma }_{k}\right),\\ {p}_{2}^{\left(x\right)}\left({\mathbf{e}}_{z},t\right)& =\text{Re}\enspace {{\Psi}}_{3}={p}_{3}^{\left(y\right)}\left({\mathbf{e}}_{z},t\right)=\text{Im}\enspace {{\Psi}}_{3}=\mathcal{O}\left({\gamma }_{k}\right),\\ {p}_{4}^{\left(+\right)}\left({\mathbf{e}}_{z},t\right)& =\text{Re}\enspace {{\Psi}}_{4}=-{\omega }_{1}^{2}{\tilde {{\epsilon}}}^{\left(+\right)}\left({\omega }_{1}\right){\text{e}}^{\text{i}{\omega }_{1}t}+\mathrm{c}.\mathrm{c}.+\mathcal{O}\left({\gamma }_{k}\right),\\ {p}_{5}^{\left({\times}\right)}\left({\mathbf{e}}_{z},t\right)& =\text{Im}\enspace {{\Psi}}_{4}={\omega }_{1}^{2}{\tilde {{\epsilon}}}^{\left({\times}\right)}\left({\omega }_{1}\right){\text{e}}^{\text{i}{\omega }_{1}t}+\mathrm{c}.\mathrm{c}.+\mathcal{O}\left({\gamma }_{k}\right),\\ {p}_{6}^{\left(b\right)}\left({\mathbf{e}}_{z},t\right)& ={{\Phi}}_{22}=\frac{{k}_{z}^{2}}{8}A\left({k}_{z}\right){\text{e}}^{\text{i}{k}_{z}t}+\mathcal{O}\left({\gamma }_{k}\right)+\mathrm{c}.\mathrm{c}..\end{align} \tag{ 5.51 }$$

It is clear, from equation (5.51), that the two vector modes ${p}_{2}^{\left(x\right)}$, and ${p}_{3}^{\left(y\right)}$, together with the longitudinal scalar mode ${p}_{1}^{\left(l\right)}$, are suppressed at the first order in γ_k , while the two standard plus and cross transverse tensor polarization modes, ${p}_{4}^{\left(+\right)}$ and ${p}_{5}^{\left({\times}\right)}$ of frequency ω₁, survive together with the transverse breathing scalar mode ${p}_{6}^{\left(b\right)}$ of frequency ω₂ at zero order in γ_k .

In summary, the gravitational radiation for non-local gravity R + a₁ R□⁻¹ R shows three polarizations: (+) and (×) massless two-spin transverse tensor modes of frequency ω₁, governed by two degrees of freedom ⁽⁺⁾(ω₁) and ⁽⁺⁾(ω₂) and a massive zero-spin transverse scalar mode of frequency ω₂, namely the breathing mode, governed by one d.o.f. A(k_z ). The main results of this paper are summarized in the table below (table 1).

Table 1. Polarizations and modes for GWs in a theory of gravity with non-local corrections.

Constraints	Order	Frequency	Polarization	Type	d.o.f.	Modes Petrov class	Helicity	Mass
$1+{a}_{1}{\phi }_{0}-{\lambda }_{0}=6{a}_{1}$	2nd	${\omega }_{1}$	2, transverse	Tensor	2	(+), (×)	2	0
						N2
$1+{a}_{1}{\phi }_{0}-{\lambda }_{0}=6{a}_{1}$	2nd	${\omega }_{1}$	3, transverse	Tensor	3	(+), (×), b	2	0
		${\omega }_{2}$		Scalar		N2	0	M

It is worth noticing that the same approach can be used also for higher-order theories of gravity both in metric and in teleparallel formalisms. For details see [45, 46].

The main result of this paper is the presence of a massive scalar gravitational mode in addition to the standard massless tensor ones in a non-local gravity theory of the form R + a₁ R□⁻¹ R. This model can be considered as a straightforward extension of general relativity where a non-local correction is taken into account. In this sense, it can be considered as an extended theory of gravity where the Einstein theory is a particular case [47–49].

The further scalar mode is achieved in the limit of plane waves, that is, the observer is assumed far from the wave source. It exhibits only a transverse polarization and not a mixed one because the longitudinal polarization is suppressed if we retain only the lowest order terms in the small parameters γ and γ_k . Such parameters take into account deviations of the waves from exactly massless ones propagating at the light speed. In addition, this further massive transverse scalar mode has helicity two and it is governed by one degree of freedom A(k_z ). Via the NP formalism, we have found that the GWs belong to the N₃, E(2) classes of Petrov classification. These results have been obtained using both the geodesic deviation and the NP formalism.

The approach can be easily extended to models containing more general terms like R□^−k R which appear as effective non-local quantum corrections. Beside the renormalization and regularization of gravitational field at ultraviolet scales [17, 18], these models can be ghost-free [11] and their infrared counterparts can be interesting at astrophysical [50] and cosmological scales [19] to address the dark side issues.

Finally, detecting further modes as the scalar massive one derived here is a major signature to break the degeneracy of modified theories of gravity which could be discriminated at fundamental level [51]. In a forthcoming paper, we will match these non-local theories with GW observations.

SC is supported by the INFN sezione di Napoli, initiative specifiche MOONLIGHT2 and QGSKY.

No new data were created or analysed in this study.