# Gravitational radiation in
higher order non-local gravity

Maurizio Capriolo<EMAIL_ADDRESS>_Dipartimento di Fisica “E. R.
Caianiello”, Università degli Studi di Salerno, via Giovanni Paolo II, 132,
Fisciano, SA I-84084, Italy._

###### Abstract

In this paper we examine gravitational radiation in higher order non-local
gravity described by the non-local gravitational Lagrangian density
$\mathcal{L}_{g}=R+\sum_{h=1}^{n}a_{h}R\Box^{-h}R$. This non-local theory of
gravitation always exhibits the tensor transverse gravitational radiation for
$k_{1}^{2}=0$, corresponding to the angular frequency $\omega_{1}$, composed
of two standard $(+)$ and $(\times)$ polarization modes, massless and of
helicity 2. Furthermore, it shows, under suitable constraint and $n\geq 2$, an
additional massive transverse scalar gravitational radiation with helicity 0.
It is composed of $n-1$ modes associated to $n-1$ angular frequencies
$\omega_{2},\ldots,\omega_{n}$, each of which of breathing polarization $(b)$
to lowest order in $\gamma$, a parameter that takes into account the
difference in speed between the slightly massive wave and the massless one.
Thanks to NP formalism, we find that the $E(2)$ class of non-local
gravitational waves is $N_{3}$, according Petrov classification, where the
presence or absence of all modes are observer independent. Also, the scalar
radiation is forbidden for $n=1$ and $n=2$ cases, when some conditions are
satisfied. Finally, in $\Box^{-1}$ gravity where $n=1$, a possible degenerate
case with a continuous infinity of transverse massive scalar breathing modes
appears under a particular constraint, which reproduces in two-dimensional
spacetime the Polyakov effective action.

Keywords: Non-local gravity; modified gravity; gravitational waves.

Mathematics Subject Classification 2020: 83C40, 83C05, 83C10, 83D99

##  1 Introduction

Local theories of gravity are described by Lagrangians that are formed by a
finite sum of products between fields and their derivatives evaluated at the
same point $x$ of the spacetime and dynamical variables are governed by
differential equations, as in general relativity. Instead, non-local theories
of gravity are described by Lagrangians composed by a finite sum of products
between fields and their derivatives evaluated at different points $x$ and
$x^{\prime}$ of the spacetime and dynamical variables are governed by integro-
differential equations. That is, the value of the field at one point depends
on its value at another point of the spacetime, weighted by a function called
nucleus or kernel. The non-locality in theories of gravity can essentially
enter through three different ways. The first manner is by means of a
convergent series expansion with real coefficients of an analytic non-
polynomial function $\mathcal{F}$ of operator $\Box$, known as Infinite
Derivative of Gravity (IDG) [1, 2, 3, 4, 5, 6, 7, 8]

$\mathcal{F}(\Box)=\sum_{n=1}^{\infty}c_{n}\Box^{n}\ .$ (1)

In the second way, the non-linearity intervenes in the non-analytic operators
as

$\Box^{-n}R(x)=\Phi_{0}(x)+\int\,d^{4}x^{\prime}\,\sqrt{-g(x^{\prime})}G(x,x^{\prime})R(x^{\prime})\
,$ (2)

where $G(x,x^{\prime})$ is the retarded Green function of operator $\Box^{n}$,
that is

$\sqrt{-g(x)}\Box_{x}^{n}G(x,x^{\prime})=\delta(x-x^{\prime})\ ,$ (3)

subject to retarded boundary conditions for the principle of causality

$G(x,x^{\prime})=0\quad\forall\,t<t^{\prime}\ ,$ (4)

and $\Phi_{0}$ is the homogeneous solution of following partial differential
equation

$\Box^{n}\Phi_{0}(x)=0\ .$ (5)

Finally, in a different approach where non-locality enters through a non-local
constitutive law as in the electrodynamics of media (for details see [9, 10]),
where the constitutive relations between
$H_{\mu\nu}\rightarrow(\mathbf{D},\mathbf{H})$ and
$F_{\mu\nu}\rightarrow(\mathbf{E},\mathbf{B})$ involving both the memory of
the electromagnetic field, that is, the temporal dispersion and the anisotropy
and non-homogeneity of the medium, that is, the spatial dispersion, take the
non-local form [11, 12]

$H_{\alpha\beta}(x)=\int\,d^{4}x^{\prime}K_{\alpha\beta}^{\phantom{\alpha\beta}\mu\nu}(x^{\prime})F_{\mu\nu}(x-x^{\prime})\
.$ (6)

Non-locality could be a manifestation of quantum nature of gravity. Non-local
models of gravity, built by adding to Einstein-Hilbert action, non-local terms
like $\Box^{-n}$ or $\mathcal{F}(\Box)$ represent quantum corrections to
classical general relativity (loop corrections), can be regarded as a quantum
effective field theories. These models, today, are widely used in cosmology
and astrophysics, for example, to explain the current cosmic expansion or the
early times acceleration as well as structure formation, without introducing
dark energy and dark matter. Specifically, non-locality could play a crucial
role to address problems like cosmological constant, Big Bang and black hole
singularities, and, in general, coincidence and fine-tuning problems, which
affect the $\Lambda$CDM model [13, 14, 15, 16]. Also, these effective theories
are often ghost-free and stable, consistent with Solar system constraints,
super-renormalizability at quantum level and then less affected by infrared
and ultraviolet divergences [17].

In this paper we analyze gravitational waves (GWs) in non-local theory of
gravity described by Lagrangian density of higher order, linear in $R$ and
$\Box^{-h}$, i.e. $\mathcal{L}_{g}=R+\sum_{h=1}^{n}a_{h}R\Box^{-h}R$. In
particular, we want to investigate how non-local terms can contribute to
adding a further modes in gravitational radiation, different from the
transverse tensor ones, predicted in general relativity.

This paper is organized as follows. In section 2, we localize the non-local
action through the introduction of suitable scalar fields, i.e. via Lagrangian
multipliers approach, deriving the related field equations. Then, in section 3
weakly perturbing the metric tensor and scalar fields we obtain the linearized
equations and in suitable gauge we find their wavelike solutions. Sections 4
and 5 are devoted to study of polarizations and helicities via geodesic
deviation and the Newman–Penrose (NP) formalism. Below, the section 6 provides
some tables while results and final remarks are summarized in section 7.

##  2 Localization of higher order non-local action via Lagrangian multipliers
approach

We want to study the gravity governed by the following higher order non-local
action

$S[g]=\frac{1}{2\kappa^{2}}\int
d^{4}x\sqrt{-g}\biggl{(}R+\sum_{h=1}^{n}a_{h}R\Box^{-h}R\biggr{)}+\int
d^{4}x\sqrt{-g}\mathcal{L}_{m}[g]\ ,$ (7)

where $\kappa^{2}=8\pi G/c^{4}$ and the D’Alembert operator $\Box$, is defined
as $\Box=g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$. Its field equations are non-
linear integro-differential equations hard to solve and therefore a different
procedure is used to analyze this model of non-local gravity, starting from an
action equivalent to that of Eq. (7). Thus, we introduce $n$ auxiliary fields
$\phi_{1}(x),\dots,\phi_{n}(x)$ as

$\displaystyle\phi_{1}(x)$ $\displaystyle=\Box^{-1}R(x)\ ,$ (8)
$\displaystyle\phi_{2}(x)$ $\displaystyle=\Box^{-2}R(x)\ ,$ (9)
$\displaystyle\vdots$ $\displaystyle\phi_{n}(x)$ $\displaystyle=\Box^{-n}R(x)\
,$ (10)

involving

$\displaystyle R(x)$ $\displaystyle=\Box\phi_{1}(x)\ ,$ (11)
$\displaystyle\phi_{1}(x)$ $\displaystyle=\Box\phi_{2}\left(x\right)\ ,$ (12)
$\displaystyle\vdots$ (13) $\displaystyle\phi_{n-1}(x)$
$\displaystyle=\Box\phi_{n}\left(x\right)\ ,$ (14)

and $n$ Lagrange multipliers, new scalar fields
$\lambda_{1}(x),\dots,\lambda_{n}(x)$, so that the free gravitational action
Eq. (7) is equivalent to [18, 19, 20, 21]

$S_{g}[g,\phi_{1},\dots\phi_{n},\lambda_{1},\dots,\lambda_{n}]=\frac{1}{2\kappa^{2}}\int
d^{4}x\sqrt{-g}\Bigg{[}R\biggl{(}1+\sum_{h=1}^{n}a_{h}\phi_{h}\biggr{)}\\\
+\lambda_{1}\left(R-\Box\phi_{1}\right)+\sum_{h=1}^{n-1}\lambda_{h+1}\left(\phi_{h}-\Box\phi_{h+1}\right)\Bigg{]}\
.$ (15)

Using the integration by parts, the Gauss theorem and imposing that fields and
their derivatives vanish on the boundary of integration domain, we get

$S_{g}[g,\phi_{1},\dots\phi_{n},\lambda_{1},\dots,\lambda_{n}]=\frac{1}{2\kappa^{2}}\int
d^{4}x\sqrt{-g}\Bigg{[}\biggl{(}1+\sum_{h=1}^{n}a_{h}\phi_{h}+\lambda_{1}\biggr{)}R\\\
+\sum_{h=1}^{n-1}\lambda_{h+1}\phi_{h}+\sum_{h=1}^{n}\nabla^{\nu}\lambda_{h}\nabla_{\nu}\phi_{h}\Bigg{]}\
.$ (16)

The variational derivative of the gravitational action Eq. (15) with respect
to $\phi_{h}$ and $\phi_{n}$ with $1\leq h\leq n-1$, yields

$\frac{2\kappa^{2}}{\sqrt{-g}}\frac{\delta
S_{g}}{\delta\phi_{h}}=a_{h}R+\lambda_{h+1}-\Box\lambda_{h},$ (17)
$\frac{2\kappa^{2}}{\sqrt{-g}}\frac{\delta
S_{g}}{\delta\phi_{n}}=a_{n}R-\Box\lambda_{n},$ (18)

while the functional derivative with respect to $\lambda_{1}$ and
$\lambda_{l}$, with $2\leq l\leq n$, takes the following form

$\frac{2\kappa^{2}}{\sqrt{-g}}\frac{\delta
S_{g}}{\delta\lambda_{1}}=R-\Box\phi_{1}\ ,$ (19)
$\frac{2\kappa^{2}}{\sqrt{-g}}\frac{\delta
S_{g}}{\delta\lambda_{l}}=\phi_{l-1}-\Box\phi_{l}\ .$ (20)

Finally, varying with respect to $g^{\mu\nu}$ the gravitational and the
material part of the action, Eq. (15) and Eq. (7), we have, respectively, the
following functional derivatives

$\displaystyle\frac{2\kappa^{2}}{\sqrt{-g}}\frac{\delta S_{g}}{\delta
g^{\mu\nu}}$
$\displaystyle=\left(G_{\mu\nu}+g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu}\right)\biggl{(}1+\sum_{h=1}^{n}a_{h}\phi_{h}+\lambda_{1}\biggr{)}-\frac{1}{2}g_{\mu\nu}\sum_{h=1}^{n-1}\phi_{h}\lambda_{h+1}$
$\displaystyle+\sum_{h=1}^{n}\nabla_{(\mu}\phi_{h}\nabla_{\nu)}\lambda_{h}-\frac{1}{2}g_{\mu\nu}\sum_{h=1}^{n}\nabla^{\sigma}\phi_{h}\nabla_{\sigma}\lambda_{h}\
,$ (21)

and

$\frac{\delta S_{m}}{\delta g^{\mu\nu}}=-\frac{\sqrt{-g}}{2}T_{\mu\nu}\ .$
(22)

Thanks to Eqs. (17), (19), (2) and (22) we derive a second order non-linear
field equation in matter, that is

$\displaystyle\left(G_{\mu\nu}+g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu}\right)$
$\displaystyle\biggl{(}1+\sum_{h=1}^{n}a_{h}\phi_{h}+\lambda_{1}\biggr{)}-\frac{1}{2}g_{\mu\nu}\sum_{h=1}^{n-1}\phi_{h}\lambda_{h+1}$
$\displaystyle+\sum_{h=1}^{n}\nabla_{(\mu}\phi_{h}\nabla_{\nu)}\lambda_{h}-\frac{1}{2}g_{\mu\nu}\sum_{h=1}^{n}\nabla^{\sigma}\phi_{h}\nabla_{\sigma}\lambda_{h}=\kappa^{2}T_{\mu\nu}\
,$ (23)

plus $2n$ non-linear differential equations each of them of $2$-th degree

$\displaystyle\Box\phi_{1}$ $\displaystyle=R\ ,$ $\displaystyle\vdots$
$\displaystyle\Box\phi_{l}$ $\displaystyle=\phi_{l-1}\ ,$
$\displaystyle\vdots$ (24) $\displaystyle\Box\lambda_{h}$
$\displaystyle=a_{h}R+\lambda_{h+1}\ ,$ $\displaystyle\vdots$
$\displaystyle\Box\lambda_{n}$ $\displaystyle=a_{n}R\ ,$

with $2\leq l\leq n$ and $1\leq h\leq n-1$, where $G_{\mu\nu}$ is the Einstein
tensor

$G_{\mu\nu}=R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R\ .$ (25)

Thus, the order of theory is $4n+2$ for $\Box^{-n}$ gravity. Thereafter, we
perform the trace of Eq. (2) that yields the trace equation

$\left(R-3\Box\right)\biggl{(}1+\sum_{h=1}^{n}a_{h}\phi_{h}+\lambda_{1}\biggr{)}+2\sum_{h=1}^{n-1}\phi_{h}\lambda_{h+1}+\sum_{h=1}^{n}\nabla^{\alpha}\lambda_{h}\nabla_{\alpha}\phi_{h}=-\kappa^{2}T\
.$ (26)

##  3 Linearization and gravitational waves

We investigate the gravitational waves at very large distances from their
sources, where the weak field approximation is applicable. So, in order to
analyze radiation in non-local gravitational theories, we first perturb the
metric tensor $g_{\mu\nu}$ around the flat metric $\eta_{\mu\nu}$ and the $2n$
scalar fields $\phi_{1},\ldots,\phi_{n}$ and $\lambda_{1},\ldots,\lambda_{n}$
around their constant values assumed in the Minkowskian spacetime,
$\phi^{(0)}$ and $\lambda^{(0)}$ [22], respectively as

$\displaystyle g_{\mu\nu}$ $\displaystyle\sim\eta_{\mu\nu}+h_{\mu\nu}\ ,$ (27)
$\displaystyle\phi_{1}$ $\displaystyle\sim\phi_{1}^{(0)}+\delta\phi_{1}\ ,$
$\displaystyle\vdots$ $\displaystyle\phi_{n}$
$\displaystyle\sim\phi_{n}^{(0)}+\delta\phi_{n}\ ,$ $\displaystyle\lambda_{1}$
$\displaystyle\sim\lambda_{1}^{(0)}+\delta\lambda_{1}\ ,$
$\displaystyle\vdots$ $\displaystyle\lambda_{n}$
$\displaystyle\sim\lambda_{n}^{(0)}+\delta\lambda_{n}\ .$

Then, to first order in $h_{\mu\nu}$, the Riemann tensor $R_{\mu\nu}$, the
Ricci scalar $R$ and the Einstein tensor $G_{\mu\nu}$ become

$R_{\mu\nu}^{(1)}=\frac{1}{2}\left(\partial_{\sigma}\partial_{\mu}h^{\sigma}_{\nu}+\partial_{\sigma}\partial_{\nu}h^{\sigma}_{\mu}-\partial_{\mu}\partial_{\nu}h-\Box
h_{\mu\nu}\right)\ ,$ (28)
$R^{(1)}=\partial_{\mu}\partial_{\nu}h^{\mu\nu}-\Box h\ ,$ (29)
$G_{\mu\nu}^{(1)}=\frac{1}{2}\left(\partial_{\sigma}\partial_{\mu}h^{\sigma}_{\nu}+\partial_{\sigma}\partial_{\nu}h^{\sigma}_{\mu}-\partial_{\mu}\partial_{\nu}h-\Box
h_{\mu\nu}-\eta_{\mu\nu}\partial_{\alpha}\partial_{\beta}h^{\alpha\beta}+\eta_{\mu\nu}\Box
h\right)\ ,$ (30)

where $h$ is the trace of perturbation $h_{\mu\nu}$ and
$\Box=\eta^{\mu\nu}\partial_{\mu}\partial_{\nu}$ is the D’Alembert operator in
the Minkowskian space. Later, according to Eqs. (2), (2), (26), (27), (28),
(29) and (30), the linearized field equations up to first order are

$\biggl{(}1+\sum_{h=1}^{n}a_{h}\phi_{h}+\lambda_{1}\biggr{)}^{(0)}G_{\mu\nu}^{(1)}+\Bigl{(}g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu}\Bigr{)}^{(0)}\biggl{(}1+\sum_{h=1}^{n}a_{h}\phi_{h}+\lambda_{1}\biggr{)}^{(1)}\\\
-\frac{1}{2}g_{\mu\nu}^{(1)}\sum_{h=1}^{n-1}\phi_{h}^{(0)}\lambda_{h+1}^{(0)}-\frac{1}{2}g_{\mu\nu}^{(0)}\sum_{h=1}^{n-1}\phi_{h}^{(1)}\lambda_{h+1}^{(0)}-\frac{1}{2}g_{\mu\nu}^{(0)}\sum_{h=1}^{n-1}\phi_{h}^{(0)}\lambda_{h+1}^{(1)}\\\
-\frac{1}{2}g_{\mu\nu}^{(0)}\sum_{h=1}^{n-1}\phi_{h}^{(0)}\lambda_{h+1}^{(0)}=\kappa^{2}T_{\mu\nu}^{(0)}\
,$ (31) $\begin{gathered}\Box\delta\phi_{1}=R^{(1)}\ ,\\\ \vdots\\\
\Box\delta\phi_{l}=\phi_{l-1}^{(0)}+\delta\phi_{l-1}\ ,\\\ \vdots\\\
\Box\delta\lambda_{h}=a_{h}R^{(1)}+\lambda_{h+1}^{(0)}+\delta\lambda_{h+1}\
,\\\ \vdots\\\ \Box\delta\lambda_{n}=a_{n}R^{(1)}\ ,\end{gathered}$ (32)
$\biggl{(}1+\sum_{h=1}^{n}a_{h}\phi_{h}+\lambda_{1}\biggr{)}^{(0)}R^{(1)}-3\Box^{(0)}\biggl{(}1+\sum_{h=1}^{n}a_{h}\phi_{h}+\lambda_{1}\biggr{)}^{(1)}\\\
+2\sum_{h=1}^{n-1}\phi_{h}^{(0)}\lambda_{h+1}^{(1)}+2\sum_{h=1}^{n-1}\phi_{h}^{(1)}\lambda_{h+1}^{(0)}+2\sum_{h=1}^{n-1}\phi_{h}^{(0)}\lambda_{h+1}^{(0)}=-\kappa^{2}T^{(0)}\
,$ (33)

where $2\leq l\leq n$ and $1\leq h\leq n-1$. Equaling terms to zero-th order
on the left and right of Eqs. (32), we obtain

$\phi_{l-1}^{(0)}=0\quad\text{and}\quad\lambda_{h+1}^{(0)}=0\ .$ (34)

Thus, field equations and the trace become

$\Bigl{(}1+a_{n}\phi_{n}^{(0)}+\lambda_{1}^{(0)}\Bigr{)}G_{\mu\nu}^{(1)}+\Bigl{(}\eta_{\mu\nu}\Box-\partial_{\mu}\partial_{\nu}\Bigr{)}\biggl{(}\sum_{h=1}^{n}a_{h}\delta\phi_{h}+\delta\lambda_{1}\biggr{)}=\kappa^{2}T_{\mu\nu}^{(0)}\
,$ (35) $\displaystyle\Box\delta\phi_{1}$ $\displaystyle=R^{(1)}\ ,$ (36)
$\displaystyle\vdots$ $\displaystyle\Box\delta\phi_{l}$
$\displaystyle=\delta\phi_{l-1}\ ,$ $\displaystyle\vdots$
$\displaystyle\Box\delta\lambda_{h}$
$\displaystyle=a_{h}R^{(1)}+\delta\lambda_{h+1}\ ,$ $\displaystyle\vdots$
$\displaystyle\Box\delta\lambda_{n}$ $\displaystyle=a_{n}R^{(1)}\ ,$
$\Bigl{(}1+a_{n}\phi_{n}^{(0)}+\lambda_{1}^{(0)}\Bigr{)}R^{(1)}-3\Box\biggl{(}\,\sum_{h=1}^{n}a_{h}\delta\phi_{h}+\delta\lambda_{1}\biggr{)}=-\kappa^{2}T^{(0)}\
.$ (37)

Now, we define a new gauge in a coordinate frame $\left\\{x^{\mu}\right\\}$
with $1+a_{n}\phi_{n}^{(0)}+\lambda_{1}^{(0)}\neq 0$ as

$\bar{h}_{\mu\nu}=h_{\mu\nu}-\frac{1}{2}\eta_{\mu\nu}h-\frac{\eta_{\mu\nu}}{1+a_{n}\phi_{n}^{(0)}+\lambda_{1}^{(0)}}\biggl{(}\,\sum_{h=1}^{n}a_{h}\delta\phi_{h}+\delta\lambda_{1}\biggr{)}\
,$ (38)

such that

$\partial_{\mu}\bar{h}^{\mu\nu}=0\ .$ (39)

The trace of Eq. (38) is

$\bar{h}=-h-\frac{4}{1+a_{n}\phi_{n}^{(0)}+\lambda_{1}^{(0)}}\biggl{(}\,\sum_{h=1}^{n}a_{h}\delta\phi_{h}+\delta\lambda_{1}\biggr{)}\
.$ (40)

In terms of barred quantities $\bar{h}_{\mu\nu}$ and $\bar{h}$, we obtain from
Eqs. (29),(30), (38) and (40)

$R^{(1)}=\frac{1}{2}\Box\bar{h}+\frac{3}{1+a_{n}\phi_{n}^{(0)}+\lambda_{1}^{(0)}}\Box\biggl{(}\,\sum_{h=1}^{n}a_{h}\delta\phi_{h}+\delta\lambda_{1}\biggr{)}\
,$ (41)
$G_{\mu\nu}^{(1)}=-\frac{1}{2}\Box\bar{h}_{\mu\nu}-\frac{1}{1+a_{n}\phi_{n}^{(0)}+\lambda_{1}^{(0)}}\left(\eta_{\mu\nu}\Box-\partial_{\mu}\partial_{\nu}\right)\biggl{(}\,\sum_{h=1}^{n}a_{h}\delta\phi_{h}+\delta\lambda_{1}\biggr{)}\
.$ (42)

If we set $B=1+a_{n}\phi_{n}^{(0)}+\lambda_{1}^{(0)}\neq 0$, field equations
can be rewritten as

$\Box\bar{h}_{\mu\nu}=-\frac{2\kappa^{2}}{B}T_{\mu\nu}^{(0)}\ ,$ (43)
$\displaystyle\Box\delta\phi_{1}$
$\displaystyle=\frac{1}{2}\Box\bar{h}+\frac{3}{B}\biggl{(}\,\sum_{h=1}^{n}a_{h}\Box\delta\phi_{h}+\Box\delta\lambda_{1}\biggr{)}\
,$ (44) $\displaystyle\vdots$ $\displaystyle\Box\delta\phi_{l}$
$\displaystyle=\delta\phi_{l-1}\ ,$ $\displaystyle\vdots$
$\displaystyle\Box\delta\lambda_{h}$
$\displaystyle=\frac{a_{h}}{2}\Box\bar{h}+\frac{3a_{h}}{B}\Box\biggl{(}\,\sum_{l=1}^{n}a_{l}\delta\phi_{l}+\delta\lambda_{1}\biggr{)}+\delta\lambda_{h+1}\
,$ $\displaystyle\vdots$ $\displaystyle\Box\delta\lambda_{n}$
$\displaystyle=\frac{a_{n}}{2}\Box\bar{h}+\frac{3a_{n}}{B}\Box\biggl{(}\,\sum_{h=1}^{n}a_{h}\delta\phi_{h}+\delta\lambda_{1}\biggr{)}\
,$

that is, ten linear differential equations in $\bar{h}_{\mu\nu}$ and $2n$
linear differential equations in $\delta\phi_{1},\ldots,\delta\phi_{n}$ and
$\delta\lambda_{1},\ldots,\delta\lambda_{n}$. While, finally, the linearized
equation of trace becomes

$\Box\bar{h}=-\frac{2\kappa^{2}}{B}T^{(0)}\ .$ (45)

###  3.1 The propagation of gravitational waves

The study of only propagation of waves starts from the linearized field
equations in vacuum, that is, in absence of massive objects, because we are
not interested in their generation. Hence, Eqs. (43) and (44) take the form

$\Box\bar{h}_{\mu\nu}=0\ ,$ (46)
$\displaystyle\left(3a_{1}-B\right)\Box\delta\phi_{1}+3\sum_{j=2}^{n}a_{j}\Box\delta\phi_{j}+3\Box\delta\lambda_{1}=0\
,$ $\displaystyle\vdots$ $\displaystyle\Box\delta\phi_{l}-\delta\phi_{l-1}=0\
,$ $\displaystyle\vdots$ (47) $\displaystyle
a_{h}\Box\delta\phi_{1}-\Box\delta\lambda_{h}+\delta\lambda_{h+1}=0\ ,$
$\displaystyle\vdots$ $\displaystyle
a_{n}\Box\delta\phi_{1}-\Box\delta\lambda_{n}=0\ ,$

where $l$ range over $2,\dots,n$ and $h$ over $1,\dots,n-1$. Adopting the
following Fourier decomposition

$h_{\mu\nu}\left(x\right)=\frac{1}{\left(2\pi\right)^{3/2}}\int
d^{3}\mathbf{k}\,\tilde{h}_{\mu\nu}\left(\mathbf{k}\right)e^{ik\cdot x}\ ,$
(48)

and the same decomposition for the other $2n$ scalar fields, we get

$k^{2}\tilde{\bar{h}}_{\mu\nu}=0\ ,$ (49) $\displaystyle
k^{2}\left[\left(3a_{1}-B\right)\widetilde{\delta\phi}_{1}+3\sum_{j=2}^{n}a_{j}\widetilde{\delta\phi}_{j}+3\widetilde{\delta\lambda}_{1}\right]=0\
,$ $\displaystyle\vdots$
$\displaystyle\widetilde{\delta\phi}_{l-1}+k^{2}\widetilde{\delta\phi}_{l}=0\
,$ $\displaystyle\vdots$ (50) $\displaystyle
a_{h}k^{2}\widetilde{\delta\phi}_{1}-k^{2}\widetilde{\delta\lambda}_{h}-\widetilde{\delta\lambda}_{h+1}=0\
,$ $\displaystyle\vdots$ $\displaystyle
k^{2}\left(a_{n}\widetilde{\delta\phi}_{1}-\widetilde{\delta\lambda}_{n}\right)=0\
.$

If we consider a null plane wave $k^{2}=0$, we get

$\widetilde{\delta\phi}_{l-1}=\widetilde{\delta\lambda}_{h+1}=0,\quad\text{with}\quad
l=2,\dots,n\quad\text{and}\quad h=1,\dots,n-1\ ,$ (51)

and

$a_{n}\widetilde{\delta\phi}_{n}+\widetilde{\delta\lambda}_{1}=0\ .$ (52)

Therefore, the linear metric perturbation $h_{\mu\nu}$ Eq. (38) is reduced to
the well-know relation

$\bar{h}_{\mu\nu}=h_{\mu\nu}-\frac{1}{2}\eta_{\mu\nu}h\ .$ (53)

Then, in an suitable gauge that leaves invariant
$\partial_{\mu}\bar{h}^{\mu\nu}=0$, that is with $\Box\xi^{\nu}=0$ under
infinitesimal transformation $x^{\prime\mu}=x^{\mu}+\xi^{\mu}$, the traceless
tensor gravitational wave with $k_{1}^{2}=0$ is

$h_{\mu\nu}^{(k_{1})}\left(x\right)=\frac{1}{\left(2\pi\right)^{3/2}}\int
d^{3}\mathbf{k}\,\widetilde{C}_{\mu\nu}\left(\mathbf{k}\right)e^{ik_{1}\cdot
x}+c.c.\ ,$ (54)

where we have used the four-vector $k_{1}^{\mu}=(\omega_{1},\mathbf{k})$ with
dispersion relation

$\omega_{1}=|\mathbf{k}|\ ,$ (55)

and $\widetilde{C}_{\mu\nu}$ is the transverse and tracefree polarization
tensor in the momentum space, i.e.
$\widetilde{C}^{\sigma}_{\phantom{\sigma}\sigma}=0$ and
$k^{\mu}\widetilde{C}_{\mu\nu}=0$. If we explore a non-null plane wave
$k^{2}\neq 0$, in addition to the trivial solution, we have from Eq. (49)

$\tilde{\bar{h}}_{\mu\nu}\left(\mathbf{k}\right)=0\Rightarrow\tilde{\bar{h}}\left(\mathbf{k}\right)=0\
,$ (56)

and coefficients square matrix $A_{2n}$ of the linear homogeneous system in
$\widetilde{\delta\phi}_{1}$, $\widetilde{\delta\phi}_{l}$,
$\widetilde{\delta\lambda}_{h}$ and $\widetilde{\delta\lambda}_{n}$ Eq. (50)
is given by

$A_{2n}=\left(\begin{array}[]{cccccccccccc}3a_{1}-B&3a_{2}&3a_{3}&\cdots&3a_{n-1}&3a_{n}&3&0&0&\cdots&0&0\\\
1&k^{2}&0&\cdots&0&0&0&0&0&\cdots&0&0\\\
0&1&k^{2}&\cdots&0&0&0&0&0&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&\cdots&1&k^{2}&0&0&0&\cdots&0&0\\\
a_{1}k^{2}&0&0&\cdots&0&0&-k^{2}&-1&0&\cdots&0&0\\\
a_{2}k^{2}&0&0&\cdots&0&0&0&-k^{2}&-1&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
a_{n-1}k^{2}&0&0&\cdots&0&0&0&0&0&\cdots&-k^{2}&-1\\\
a_{n}&0&0&\cdots&0&0&0&0&0&\cdots&0&-1\end{array}\right)\ .$ (57)

Its determinant is

$\lvert
A_{2n}\rvert=\left(-1\right)^{n}k^{2\left(n-1\right)}\biggl{[}\left(6a_{1}-B\right)k^{2\left(n-1\right)}\\\
+\theta(n-2)\sum_{j=2}^{n}6\left(-1\right)^{j+1}a_{j}k^{2\left(n-j\right)}\biggr{]}\
,$ (58)

where $\theta(n-2)$ is Heaviside function that vanishes for $n<2$ and equals
to $1$ otherwise. The square bracket in Eq. (58) is a $(n-1)$-th degree
polynomial in $k^{2}\neq 0$ , that vanishes if

$|A_{2n}|=0\quad\Leftrightarrow\quad\left(6a_{1}-B\right)\prod_{j=2}^{n}\left(k^{2}-k_{j}^{2}\right)=0\
,$ (59)

where $a_{j}$ are such that the equation $\lvert A_{2n}\rvert=0$ admits $n-1$
positive real solutions $k^{2}_{j}\in\mathbb{R}_{+}$. If the next constraint
is fulfilled

$B-6a_{1}=0\ ,$ (60)

the determinant Eq. (58) becomes a $(n-2)$-th degree polynomial in $k^{2}\neq
0$

$|A_{2n}|=6\left(-1\right)^{n}k^{2\left(n-1\right)}\left[\theta(n-2)\sum_{j=2}^{n}\left(-1\right)^{j+1}a_{j}k^{2\left(n-j\right)}\right]\
.$ (61)

Non-zero roots of the polynomial (61) in $k^{2}$ under constraint Eq. (60),
varying $n$ and under suitable conditions for coefficients $a_{j}$, are

$\displaystyle n=1$ $\displaystyle|A_{2}|=0$ $\displaystyle
k^{2}\in\mathbb{R}_{+}\;\text{solutions}$ $\displaystyle\text{there are
$\infty^{1}$ solutions}\ ,$ (62) $\displaystyle n=2$
$\displaystyle|A_{4}|=-6a_{2}k^{2}=0$ $\displaystyle\text{there is no
solutions},$ (63) $\displaystyle n=3$
$\displaystyle|A_{6}|=-6k^{4}\left[-a_{2}k^{2}+a_{3}\right]=0$
$\displaystyle\text{one solution $k_{2}^{2}$}\ ,$ (64) $\displaystyle n=4$
$\displaystyle|A_{8}|=6k^{6}\left[-a_{2}k^{4}+a_{3}k^{2}-a_{4}\right]=0$
$\displaystyle\text{two solutions $k_{2}^{2}$ and $k_{3}^{2}$}\ ,$ (65)
$\displaystyle n$ $\displaystyle|A_{2n}|=0$ there are $n-2$ solutions
$\displaystyle k_{2}^{2},\ldots,k_{n-1}^{2}\ .$ (66)

While for the following other constraint

$B-6a_{1}\neq 0\ ,$ (67)

and $k^{2}=k_{j}^{2}=M_{j}^{2}\neq 0$ with $j=2,\dots,n$, we have $n-1$
solutions for $n\geq 2$ and no solution like $k^{2}\neq 0$ for $n=1$. From the
relation between the traces $\bar{h}$ and $h$ Eq. (40) and the first equation
of the linear system Eq. (50) in $\mathbf{k}$-space we obtain in $k^{2}\neq 0$
case for which both $\tilde{\bar{h}}_{\mu\nu}\left(\mathbf{k}\right)$ and
$\tilde{\bar{h}}\left(\mathbf{k}\right)$ vanish

$\tilde{h}\left(\mathbf{k}\right)=-\frac{4}{3}\widetilde{\delta\phi}_{1}\left(\mathbf{k}\right)\
,$ (68)

and consequently

$\tilde{h}_{\mu\nu}\left(\mathbf{k}\right)=-\frac{\eta_{\mu\nu}}{3}\widetilde{\delta\phi}_{1}\left(\mathbf{k}\right)\
.$ (69)

Then, the $n-1$ non-null wavelike solutions, linked to ${j}$-th mode, are of
the type

$h_{\mu\nu}^{(k_{j})}\left(x\right)=\frac{1}{\left(2\pi\right)^{3/2}}\int
d^{3}\mathbf{k}\,\frac{\eta_{\mu\nu}}{4}\tilde{A}_{j}\left(\mathbf{k}\right)e^{ik_{j}\cdot
x}+c.c.\ ,$ (70)

with the four-vector $k_{j}^{\mu}=(\omega_{j},\mathbf{k})$ where the $j$-th
amplitude in $\mathbf{k}$ space is

$\tilde{A}_{j}\left(\mathbf{k}\right)=\tilde{h}^{(k_{j})}\left(\mathbf{k}\right)\
,$ (71)

for the $j$-th dispersion relation

$\omega_{j}\left(\mathbf{k}\right)=\sqrt{|\mathbf{k}|^{2}+M_{j}^{2}}\ .$ (72)

To show that $k_{j}$-waves are massive, namely the linear perturbation of
scalar fields satisfy an equation like Klein-Gordon, according to the
linearized system Eq. (47), we can rewrite it as equivalent system

$\displaystyle\left(3a_{1}-B\right)\Box\delta\phi_{1}+3\sum_{j=2}^{n}a_{j}\Box\delta\phi_{j}+3\Box\delta\lambda_{1}=0\
,$ $\displaystyle\vdots$
$\displaystyle\Box^{l}\delta\phi_{l}-\Box\delta\phi_{1}=0\ ,$
$\displaystyle\vdots$ (73) $\displaystyle
a_{h}\Box\delta\phi_{1}-\Box\delta\lambda_{h}+\delta\lambda_{h+1}=0\ ,$
$\displaystyle\vdots$
$\displaystyle\sum_{j=1}^{n}a_{j}\Box^{n-j+1}\delta\phi_{1}-\Box^{n}\delta\lambda_{1}=0\
,$

with $l=2,\dots,n-1$ and $h=1,\dots,n-1$. Then, by acting of the differential
operator $\Box^{n-1}$ on the first equation of the system Eq. (73) we get

$\left[\left(6a_{1}-B\right)\Box^{n}+3\sum_{j=2}^{n}a_{j}\left(\Box^{j-1}+\Box^{n-j+1}\right)\right]\delta\phi_{1}=0\
,$ (74)

that can be decomposed as

$\left[\left(6a_{1}-B\right)\prod_{j=1}^{n}\left(\Box+k_{j}^{2}\right)\right]\delta\phi_{1}=0\
,$ (75)

under the condition that $k_{j}^{2}\neq 0$ when $j$ run over $2,\dots,n$, are
roots of $(n-1)$-th degree polynomial Eq. (58) in $k^{2}$ and $k_{1}^{2}=0$.
If we define the new scalar fields $\delta\Phi_{h}$, with
$h\in\\{2,\dots,n\\}$, as [23]

$\delta\Phi_{h}=\prod_{j=1,j\neq
h}^{n}\left(\Box+M_{j}^{2}\right)\delta\phi_{1}\ ,$ (76)

we obtain the Klein-Gordon equation for each scalar fields $\delta\Phi_{h}$
with mass $M_{h}\neq 0$ and $M_{1}=0$, that is

$\left(\Box+M_{h}^{2}\right)\delta\Phi_{h}=0\ ,$ (77)

and accordingly

$\delta\Phi_{h}=\frac{\prod_{j=1,j\neq
h}^{n}\left(M_{j}^{2}-M_{h}^{2}\right)}{(2\pi)^{3/2}}\int
d^{3}\mathbf{k}\widetilde{\delta\phi}_{1}(\mathbf{k})e^{ik_{h}\cdot x}\ .$
(78)

In summary, the gravitational wave of higher order non-local gravity, from
superposition principle can be expressed as

$h_{\mu\nu}\left(x\right)=\frac{1}{\left(2\pi\right)^{3/2}}\int
d^{3}\mathbf{k}\,\widetilde{C}_{\mu\nu}\left(\mathbf{k}\right)e^{ik_{1}\cdot
x}\\\ +\frac{1}{\left(2\pi\right)^{3/2}}\sum_{j=2}^{n}\int
d^{3}\mathbf{k}\,\frac{\eta_{\mu\nu}}{4}\tilde{A}_{j}\left(\mathbf{k}\right)e^{ik_{j}\cdot
x}+c.c.\ ,$ (79)

that is, one massless wave $k_{1}^{2}=0$ and $n-1$ massive waves
$k_{j}^{2}=M_{j}^{2}\neq 0$ under the constraint $B\neq 6a_{1}$ occur. It
should be noted that when $B=6a_{1}$, massive waves are $n-2$, because one
non-null wave is inhibit and the $j$ index in Eq. (79) run from $2$ to $n-1$.
For the study of GWs in higher order local gravity governed by Lagrangian
$R+\sum_{h=0}^{n}a_{h}R\Box^{h}R$, see [24], and for a confront between
general relativity and experiments, see [25, 26, 27, 28].

###  3.2 The degenerate case $n$=1

In this particular case, the gravitational part of the action Eq. (7) takes
the following form [29]

$S[g]=\frac{1}{2\kappa^{2}}\int
d^{4}x\sqrt{-g}\left(R+a_{1}R\Box^{-1}R\right)\ ,$ (80)

The gravitational wave can be derived as solutions of the linearized field
equations of the free action Eq. (80). Thus, according to the system Eqs. (46)
and (47), we get

$\begin{gathered}\Box\bar{h}_{\mu\nu}=0\ ,\\\
\left(3a_{1}-B\right)\Box\delta\phi_{1}+3\Box\delta\lambda_{1}=0\ ,\\\
a_{1}\Box\delta\phi_{1}-\Box\delta\lambda_{1}=0\ ,\end{gathered}$ (81)

that in momentum space, according to decomposition Eq.(48), equations (81)
yield

$\begin{gathered}k^{2}\tilde{\bar{h}}_{\mu\nu}=0\ ,\\\
k^{2}\left[\left(3a_{1}-B\right)\widetilde{\delta\phi}_{1}+3\widetilde{\delta\lambda}_{1}\right]=0\
,\\\
k^{2}\Bigl{(}a_{1}\widetilde{\delta\phi}_{1}-\widetilde{\delta\lambda}_{1}\Bigr{)}=0\
.\end{gathered}$ (82)

When $k^{2}=0$ Eqs. (82) admit the standard plus and times solutions of
general relativity. While, when $k^{2}\neq 0$ the following coefficient
matrix, that is

$A_{2}=\begin{pmatrix}3a_{1}-B&3\\\ a_{1}&-1\end{pmatrix}\ ,$ (83)

has the following determinant

$|A_{2}|=(-1)(6a_{1}-B)\ .$ (84)

It vanish only when $B=6a_{1}$. Then, for $B\neq 6a_{1}$ there are no
$k^{2}\neq 0$ solutions, while for $B=6a_{1}$ we get

$\tilde{\bar{h}}_{\mu\nu}=0\quad\text{and}\quad\widetilde{\delta\lambda}_{1}=a_{1}\widetilde{\delta\phi}_{1}\
.$ (85)

From the gauge Eqs. (38),(39) and (40) and from relations Eqs. (85), it gets

$\tilde{h}=-\frac{4}{3}\widetilde{\delta\phi}_{1}=\tilde{A}\ ,$ (86)
$\tilde{h}_{\mu\nu}=-\frac{1}{3}\eta_{\mu\nu}\widetilde{\delta\phi}_{1}\ ,$
(87)

from which the gravitational wave is obtained

$h_{\mu\nu}\left(x\right)=\frac{1}{\left(2\pi\right)^{3/2}}\int
d^{3}\mathbf{k}_{1}\,\widetilde{C}_{\mu\nu}\left(\mathbf{k}_{1}\right)e^{ik_{1}\cdot
x}\\\ +\frac{1}{\left(2\pi\right)^{3/2}}\int
d^{3}\mathbf{k}\,\frac{\eta_{\mu\nu}}{4}\tilde{A}\left(\mathbf{k}\right)e^{ik\cdot
x}+c.c.\ .$ (88)

That is, when $B=6a_{1}$ it can be regarded as a degenerate case, because
every value $k^{2}\in\mathbb{R}_{+}$ is solutions of system Eq. (82), given
that it is not the solution of any algebraic equation. If we can set
$\phi_{1}^{(0)}=\lambda_{1}^{(0)}=0$, the constraint Eq. (60) becomes

$a_{1}=\frac{1}{6}\ ,$ (89)

that is, a dimensionless constant because also $\Box^{-1}R$ is
dimensionless.The gravitational action Eq. (7) takes the simple form

$S_{g}[g]=\frac{1}{2\kappa^{2}}\int
d^{4}x\sqrt{-g}\biggl{(}R+\frac{1}{6}R\,\Box^{-1}R\biggr{)}\ ,$ (90)

and in units $G=c=1$ in two-dimensional curved spacetime, its non-local part
reproduces the Polyakov effective action [30, 31]

$S_{g}[g]=\frac{1}{96}\int d^{2}x\sqrt{-g}R\,\Box^{-1}R\ .$ (91)

###  3.3 The case $n$=2

In this second case, the gravitational action Eq. (7) takes the following form

$S[g]=\frac{1}{2\kappa^{2}}\int
d^{4}x\sqrt{-g}\left(R+a_{1}R\Box^{-1}R+a_{2}R\Box^{-2}R\right)\ .$ (92)

The differential system in linear approximation Eq. (47) becomes

$\begin{gathered}\Box\bar{h}_{\mu\nu}=0\ ,\\\
\left(3a_{1}-B\right)\Box\delta\phi_{1}++3a_{2}\Box\delta\phi_{2}+3\Box\delta\lambda_{1}=0\
,\\\ \Box\delta\phi_{2}-\delta\phi_{1}=0\ ,\\\
a_{1}\Box\delta\phi_{1}-\Box\delta\lambda_{1}+\delta\lambda_{2}=0\ ,\\\
a_{2}\Box\delta\phi_{1}-\Box\delta\lambda_{2}=0\ ,\end{gathered}$ (93)

that in momentum space, according to decomposition Eq. (48), yields

$\begin{gathered}k^{2}\tilde{\bar{h}}_{\mu\nu}=0\ ,\\\
k^{2}\left[\left(3a_{1}-B\right)\widetilde{\delta\phi}_{1}+3a_{2}\widetilde{\delta\phi}_{2}+3\widetilde{\delta\lambda}_{1}\right]=0\
,\\\ \widetilde{\delta\phi}_{1}+k^{2}\widetilde{\delta\phi}_{2}=0\ ,\\\
a_{1}k^{2}\widetilde{\delta\phi}_{1}-k^{2}\widetilde{\delta\lambda}_{1}-\widetilde{\delta\lambda}_{2}=0\
,\\\
k^{2}\Bigl{(}a_{2}\widetilde{\delta\phi}_{1}-\widetilde{\delta\lambda}_{2}\Bigr{)}=0\
.\end{gathered}$ (94)

For $k^{2}=0$ the Eq. (94) gives back the two standard plus and times
gravitational wave, proper to general relativity. Now let’s see under what
conditions non-null waves like $k^{2}\neq 0$ are obtained. Then, the
coefficient matrix of the last four equations of the Eq. (94) yields

$A_{4}=\begin{pmatrix}3a_{1}-B&3a_{2}&3&0\\\ 1&k^{2}&0&0\\\
a_{1}k^{2}&0&-k^{2}&-1\\\ a_{2}&0&0&-1\end{pmatrix}\ ,$ (95)

whose determinant takes a form as

$|A_{4}|=k^{2}\left[(6a_{1}-B)k^{2}-6a_{2}\right]\ .$ (96)

The first equation of Eq. (94) gives $\bar{h}_{\mu\nu}=0$ and its trace
vanishes $\bar{h}=0$. Thus, if $B=6a_{1}$, there is no non-null solutions
because $|A_{4}|\neq 0$, while for $B\neq 6a_{1}$ the determinant vanishes for
one non-null solution

$k_{2}^{2}=\frac{6a_{2}}{6a_{1}-B}\ ,$ (97)

if

$\frac{6a_{2}}{6a_{1}-B}>0\ .$ (98)

Hence from Eq. (69), we have

$\tilde{h}_{\mu\nu}(\mathbf{k}_{2})=-\frac{1}{3}\eta_{\mu\nu}\tilde{\phi}_{1}(\mathbf{k}_{2})=\frac{\eta_{\mu\nu}}{4}\tilde{A}_{2}(\mathbf{k}_{2})\
.$ (99)

Then, the combination of null and non-null waves $k_{1}$ and $k_{2}$ give us

$h_{\mu\nu}\left(x\right)=\frac{1}{\left(2\pi\right)^{3/2}}\int
d^{3}\mathbf{k}\,\widetilde{C}_{\mu\nu}\left(\mathbf{k}\right)e^{ik_{1}\cdot
x}\\\ +\frac{1}{\left(2\pi\right)^{3/2}}\int
d^{3}\mathbf{k}\,\frac{\eta_{\mu\nu}}{4}\tilde{A}_{2}\left(\mathbf{k}\right)e^{ik_{2}\cdot
x}+c.c.\ .$ (100)

##  4 Polarizations via geodesic deviation equation and helicity

The equation of geodesic deviation provides an effective tool to investigate
the polarizations of gravitational radiation when the wave invests a small
region of spacetime, as a measure of relative acceleration between nearby
geodesic, on which two freely test masses move slowly. Hence, we consider a
wave $h_{\mu\nu}\left(t-v_{g}z\right)$ propagating in $+\hat{z}$ direction,
where $v_{g}$ is the group velocity in units where $c=1$ defined as

$v_{g}=\frac{d\omega}{dk_{z}}\ ,$ (101)

that for the dispersion relation of ours $j$-th modes Eq. (72) becomes

$v_{g_{j}}=\frac{c^{2}k_{z}}{\sqrt{k_{z}^{2}+M_{j}^{2}}}\ .$ (102)

Now, let us start from the geodesic deviation equation for only spatial
components in a locally Lorentz normal coordinate system, that reduces to the
form

$\ddot{x}^{i}=-R^{i}_{\phantom{i}0k0}x^{k}\ ,$ (103)

with the Latin index range over the set $\left\\{1,2,3\right\\}$. The
separation vector is $x^{\mu}$ and $R^{i}_{\phantom{i}0k0}$ are the only six
measurable components, so-called electric components [32]. Inserting the
linearized electric components of the Riemann tensor
$R^{\left(1\right)}_{\phantom{1}i0j0}$ expressed in terms of the metric
perturbation $h_{\mu\nu}$

$R^{\left(1\right)}_{\phantom{1}i0j0}=\frac{1}{2}\left(h_{i0,j0}+h_{j0,i0}-h_{ij,00}-h_{00,ij}\right)\
,$ (104)

into Eq. (103), we find a linear non-homogeneous system of differential
equations

$\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\\\
\qquad\qquad\qquad\quad+\frac{1}{2}\left(h_{02,03}-h_{23,00}\right)y+\frac{1}{2}\left(2h_{03,03}-h_{33,00}-h_{00,33}\right)z\end{cases}\
.$ (105)

Thus, we keep $k_{1}^{\mu}=\left(\omega_{1},0,0,k_{z}\right)$ four-vector
fixed, which from dispersion relation Eq. (55), means to keep $k_{z}$. For a
massless plane wave traveling in $+\hat{z}$ direction, that is, $k_{1}^{2}=0$,
which propagates at speed $c$, Eq. (79) yields

$h^{(k_{1})}_{\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^{(\times)}_{\mu\nu}\right]e^{i\omega_{1}\left(t-z\right)}+c.c.\
,$ (106)

with two polarization tensors

$\epsilon^{(+)}_{\mu\nu}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&0&0&0\\\
0&1&0&0\\\ 0&0&-1&0\\\ 0&0&0&0\end{pmatrix}\ ,$ (107)
$\epsilon^{(\times)}_{\mu\nu}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&0&0&0\\\
0&0&1&0\\\ 0&1&0&0\\\ 0&0&0&0\end{pmatrix}\ ,$ (108)

and angular frequency $\omega_{1}=k_{z}$. Moreover, when $k_{j}^{2}\neq
0$,i.e., for a massive plane wave, propagating in $+\hat{z}$ direction keeping
$k_{j}^{\mu}=\left(\omega_{j},0,0,k_{z}\right)$ fixed, that according to
dispersion relation Eq. (72) means $k_{z}$ fixed, Eq. (79) becomes

$h^{(k_{j})}_{\mu\nu}\left(t,z\right)=\frac{\tilde{A}_{j}\left(k_{z}\right)}{4}\eta_{\mu\nu}e^{i\left(\omega_{j}t-k_{z}z\right)}+c.c.\
,$ (109)

where here the propagation speed of $j$-th mode is less than $c$ and $j$ run
over $\\{2,\dots,n\\}$. In a more compact form, the metric linear perturbation
$h_{\mu\nu}$, traveling in the $+\hat{z}$ direction with $k_{z}$ fixed, is
given by

$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^{(\times)}_{\mu\nu}\right]e^{i\omega_{1}\left(t-z\right)}\\\
+\sum_{j=2}^{n}\tilde{\epsilon}^{\left(s_{j}\right)}_{\mu\nu}\left(k_{z}\right)e^{i\left(\omega_{j}t-k_{z}z\right)}+c.c.\
,$ (110)

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

$\tilde{\epsilon}^{\left(s_{j}\right)}_{\mu\nu}\left(k_{z}\right)=\frac{\eta_{\mu\nu}}{4}\tilde{A}_{j}\left(k_{z}\right)=\left(\epsilon^{(TT)}_{\mu\nu}-\sqrt{2}\epsilon_{\mu\nu}^{(b)}-\epsilon_{\mu\nu}^{(l)}\right)\frac{\tilde{A}_{j}\left(k_{z}\right)}{4}\
,$ (111)

and the three polarization tensors are explicitly the following

$\displaystyle\epsilon^{(TT)}_{\mu\nu}$
$\displaystyle=\begin{pmatrix}1&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\end{pmatrix}\ ,$ $\displaystyle\epsilon^{(b)}_{\mu\nu}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{pmatrix}0&0&0&0\\\ 0&1&0&0\\\
0&0&1&0\\\ 0&0&0&0\end{pmatrix}\ ,$ $\displaystyle\epsilon^{(l)}_{\mu\nu}$
$\displaystyle=\begin{pmatrix}0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
0&0&0&1\end{pmatrix}\ .$ (112)

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\\}$
fulfill the orthonormality relations

$\text{Tr}\left\\{\epsilon^{(a)}\epsilon^{(b)}\right\\}=\epsilon_{\mu\nu}^{(a)}\epsilon^{(b)\mu\nu}=\delta^{ab}\quad\text{with}\quad
a,b\in\left\\{+,\times,TT,b,l\right\\}\ .$ (113)

The helicity of wave Eq. (110) can be derived by performing a rotation of
$\theta$ angle along the $z$ propagation direction, on the polarization
tensors as follows

$\epsilon_{\mu\nu}^{\prime}=R_{\mu}^{\alpha}\epsilon_{\alpha\beta}R_{\beta}^{\nu}\
,$ (114)

where the matrix of rotation is

$R=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&\cos\theta&\sin\theta\\\
0&0&-\sin\theta&\cos\theta\end{pmatrix}\ .$ (115)

Hence the wave has helicity $s$ is if polarization tensor changes as

$\epsilon_{\mu\nu}^{\prime}=e^{\pm is\theta}\epsilon_{\mu\nu}\ .$ (116)

Thus, for the $\omega_{1}$ mode with $k_{1}^{2}=0$, the polarization tensors
$\epsilon_{\mu\nu}^{(+)}$ and $\epsilon_{\mu\nu}^{(\times)}$ transform as

$\epsilon^{\prime}=e^{\pm 2i\theta}\epsilon\ ,$ (117)

that is, we obtain the tensor modes with helicity equal to two. While for
$\omega_{j}$ mode with $k_{j}^{2}\neq 0$, the polarization tensors
$\epsilon_{\mu\nu}^{(s_{j})}$ transform as

$\Bigl{(}\epsilon^{(s_{j})}\Bigr{)}^{\prime}=\epsilon^{(s_{j})}\ ,$ (118)

that is, we obtain $j$ scalar modes with helicity equal to zero. The
polarization tensor of the single scalar $j$-th mode
$\tilde{\epsilon}^{\left(s_{j}\right)}_{\mu\nu}$, is a mixed state made up of
a combination of longitudinal and transverse scalar modes, produced by the
single degree of freedom $\tilde{A}_{j}$, similarly to $f(R)$ gravity which
has three d.o.f.: $\tilde{\epsilon}^{(+)}$, $\tilde{\epsilon}^{(\times)}$ and
$\tilde{A}$ [33]. However, as we will see, the transverse component weighs
more than the longitudinal one if we consider slightly massive waves. Then,
considering that the polarization of the wave is determined only by the
spatial components of the polarization tensor
$\tilde{\epsilon}^{\left(s_{j}\right)}_{\mu\nu}$, we use its restriction to
only spatial components, i.e., $\tilde{\epsilon}^{\left(s_{j}\right)}_{m,n}$,
that give

$\tilde{\epsilon}^{\left(s_{j}\right)}_{m,n}=\Bigl{(}\sqrt{2}\epsilon_{m,n}^{(b)}+\epsilon_{m,n}^{(l)}\Bigl{)}\frac{\tilde{A}_{j}\left(k_{z}\right)}{4}\
,$ (119)

where $(m,n)$ range over $(1,2,3)$. Thus, inserting the generic gravitational
wave $z$-axes propagating at $k_{z}$ fixed into the system Eq. (105), we get
in geometrized units where $c=G=1$

$\left\\{\begin{array}[]{l}\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]e^{i\omega_{1}\left(t-z\right)}\\\
\qquad\qquad\qquad\qquad\qquad\qquad-\frac{1}{8}\sum_{j=2}^{n}\omega^{2}_{j}\tilde{A}_{j}\left(k_{z}\right)xe^{i\left(\omega_{j}t-k_{z}z\right)}+c.c.\\\
\\\
\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]e^{i\omega_{1}\left(t-z\right)}\\\
\qquad\qquad\qquad\qquad\qquad\qquad-\frac{1}{8}\sum_{j=2}^{n}\omega^{2}_{j}\tilde{A}_{j}\left(k_{z}\right)ye^{i\left(\omega_{j}t-k_{z}z\right)}+c.c.\\\
\\\
\ddot{z}(t)=-\frac{1}{8}\sum_{j=2}^{n}M^{2}_{j}\tilde{A}_{j}\left(k_{z}\right)ze^{i\left(\omega_{j}t-k_{z}z\right)}+c.c.\end{array}\right.\
,$ (120)

that can be integrated, assuming that displacements are small, as

$\begin{cases}x(t)=x_{0}+\frac{1}{2}\left[\tilde{\epsilon}^{\left(+\right)}x_{0}+\tilde{\epsilon}^{\left(\times\right)}y_{0}\right]e^{i\omega_{1}\left(t-z\right)}\\\
\qquad\qquad\qquad\qquad+\frac{1}{8}\sum_{j=2}^{n}\tilde{A}_{j}\left(k_{z}\right)x_{0}e^{i\left(\omega_{j}t-k_{z}z\right)}+c.c.\\\
\\\
y(t)=y_{0}+\frac{1}{2}\left[\tilde{\epsilon}^{\left(\times\right)}x_{0}-\tilde{\epsilon}^{\left(+\right)}y_{0}\right]e^{i\omega_{1}\left(t-z\right)}\\\
\qquad\qquad\qquad\qquad+\frac{1}{8}\sum_{j=2}^{n}\tilde{A}_{j}\left(k_{z}\right)y_{0}e^{i\left(\omega_{j}t-k_{z}z\right)}+c.c.\\\
\\\
z(t)=z_{0}+\frac{1}{8\omega^{2}_{j}}M^{2}_{j}\tilde{A}_{j}\left(k_{z}\right)z_{0}e^{i\left(\omega_{j}t-k_{z}z\right)}+c.c.\end{cases}.$
(121)

Hence, the solution Eqs. (121) in the case $k_{1}^{2}=0$ for massless waves
$h_{\mu\nu}^{(k_{1})}$, that is, exactly null plane, associated to angular
frequency $\omega_{1}$, gives

$\begin{cases}x(t)=x_{0}+\frac{1}{2}\left[\tilde{\epsilon}^{\left(+\right)}x_{0}+\tilde{\epsilon}^{\left(\times\right)}y_{0}\right]e^{i\omega_{1}\left(t-z\right)}+c.c.\\\
y(t)=y_{0}+\frac{1}{2}\left[\tilde{\epsilon}^{\left(\times\right)}x_{0}-\tilde{\epsilon}^{\left(+\right)}y_{0}\right]e^{i\omega_{1}\left(t-z\right)}+c.c.\\\
z(t)=z_{0}+c.c.\end{cases}\ ,$ (122)

namely, we recover the two pure, plus and cross, standard transverse tensor
modes predicted by general relativity,.

Otherwise, in the pure case of a massive wave $h^{(k_{j})}_{\mu\nu}$, i.e,
non-null plane waves, with $0\neq
k_{j}^{2}=M^{2}_{j}=\omega^{2}_{j}-k^{2}_{z}$ for the $j$-th mode of angular
frequency $\omega_{j}$, the solution Eq. (121) becomes

$\left\\{\begin{array}[]{l}x(t)=x(0)+\frac{1}{8}\tilde{A}_{j}\left(k_{z}\right)x(0)e^{i\left(\omega_{j}t-k_{z}z\right)}+c.c.\\\
\\\
y(t)=y(0)+\frac{1}{8}\tilde{A}_{j}\left(k_{z}\right)y(0)e^{i\left(\omega_{j}t-k_{z}z\right)}+c.c.\\\
\\\
z(t)=z(0)+\frac{1}{8\omega^{2}_{j}}M^{2}_{j}\tilde{A}_{j}\left(k_{z}\right)z(0)e^{i\left(\omega_{j}t-k_{z}z\right)}+c.c.\end{array}.\right.$
(123)

When we are sufficiently far from the radiation source, we can suppose that
$M^{2}_{j}$ is very small, i.e., that the plane waves are nearly null. Then,
we can expand our quantities with respect to a parameter $\gamma$, which takes
into account the difference in speed between nearly null plane waves with
speed $v_{g}$ and exactly null plane waves with speed $c$ where the parameter
$\gamma$ vanishes. Also, the expansion coefficient $\gamma$ reads [34]

$\gamma=\left(\frac{c}{v_{g}}\right)^{2}-1\ .$ (124)

By using Landau symbols, namely little-$o$ and big-$\mathcal{O}$ notation, we
can keep $k_{z}$ fixed and expand in terms of our parameter $\gamma$ the
following quantities

$\frac{M^{2}_{j}}{k_{z}^{2}}=\frac{1}{k_{z}^{2}}\left[\frac{\omega_{j}^{2}}{c^{2}}-k_{z}^{2}\right]=\left(\frac{\omega_{j}}{ck_{z}}\right)^{2}-1=\gamma\
,$ (125)

and in $c=1$ units we have

$\frac{\omega_{j}}{k_{z}}=\sqrt{1+\gamma}=\left(1+\frac{1}{2}\gamma\right)+o\left(\gamma\right)\
,$ (126)
$e^{i\left(\omega_{j}t-k_{z}z\right)}=e^{ik_{z}\left(t-z\right)}+\mathcal{O}\left(\gamma\right)\
,$ (127)

from which it gets to first order in $\gamma$

$\frac{M^{2}_{j}}{\omega_{j}^{2}}=\gamma+o\left(\gamma\right)\ .$ (128)

Hence, taking into account the previous expansions, the solution Eq. (123) for
the mode $\omega_{j}$ give us

$\left\\{\begin{array}[]{l}x(t)=x(0)+\frac{1}{8}\tilde{A}_{j}\left(k_{z}\right)x(0)e^{ik_{z}\left(t-z\right)}+\mathcal{O}\left(\gamma\right)+c.c.\\\
\\\
y(t)=y(0)+\frac{1}{8}\tilde{A}_{j}\left(k_{z}\right)y(0)e^{ik_{z}\left(t-z\right)}+\mathcal{O}\left(\gamma\right)+c.c.\\\
\\\
z(t)=z(0)+\frac{1}{8}\gamma\tilde{A}_{j}\left(k_{z}\right)z(0)e^{ik_{z}\left(t-z\right)}+\mathcal{O}\left(\gamma^{2}\right)+c.c.\end{array}\
,\right.$ (129)

that is

$\Delta z(t)=o\left(\Delta x(t)\right)\quad\text{and}\quad\Delta
z(t)=o\left(\Delta y(t)\right)\quad\text{for}\quad\gamma\to 0\ .$ (130)

This suggests that longitudinal modes are infinitesimal in higher order than
transverse modes when $\gamma$ tends to zero and so only the breathing tensor
polarization $\epsilon^{(b)}_{m,n}$ survives. Thus, when a gravitational wave
strikes a sphere of particles of radius $r=\sqrt{x^{2}(0)+y^{2}(0)+z^{2}(0)}$,
from the Eq. (129) the radiation will distort it into an ellipsoid described
by

$\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}\
,$ (131)

where the principal axes to zero-th order in $\gamma$ are given by

$\displaystyle\rho_{1}(t)=1+\frac{1}{4}\tilde{A}_{j}\left(k_{z}\right)\cos\left[k_{z}\left(t-z\right)+\phi\right]+\mathcal{O}\left(\gamma\right)\
,$ (132) $\displaystyle\rho_{2}(t)=1+\mathcal{O}\left(\gamma\right)\ .$ (133)

So, to lowest order only $\rho_{1}$ varies between its maximum and minimum
value. This means that the ellipsoid swings only on $xy$-plane between two
circumferences of minimum and maximum radius, reproducing the additional
transverse scalar breathing polarization which has zero helicity to lowest
order in $\gamma$ [35].

According to these considerations, under constraint $B\neq 6a_{1}$, the
gravity $R+\sum_{h=1}^{n}a_{h}R\Box^{-h}R$, linear in $R$ and $\Box^{-k}$
[36], has $n+1$ degrees of freedom: two of these,
$\tilde{\epsilon}^{\left(+\right)}$ and
$\tilde{\epsilon}^{\left(\times\right)}$, generate the standard tensor modes
of general relativity, while each $n-1$ degrees of freedom $\tilde{A}_{j}$
gives rise to a further breathing scalar mode. Instead, for $B=6a_{1}$ one
scalar mode is inhibited for all $n$ greater than two, while in $n=1$ a
degenerate state with continuous infinity of solutions $k^{2}\neq 0$ occurs.
In summary, generally $R+\sum_{h=1}^{n}a_{h}R\Box^{-h}R$ gravity has three
polarizations, $(+)$, $(\times)$ and $(b)$, but $n+1$ modes, namely two
massless 2-helicity transverse tensor modes and $n-1$ massive 0-helicity
scalar mode, all purely transverse to lowest order in $\gamma$, exactly like
$f(R)$ gravity (see for a discussion [33, 37, 38, 39, 40]).

##  5 Polarizations via NP formalism and helicity

Newman-Penrose (NP) formalism is a additional method to analyze polarizations
that works out for massless waves. Even if it is not directly applicable to
massive waves, NP formalism can be applied to slightly massive gravitational
waves generalizing it to waves propagating along nearly null geodesics [34,
41]. Furthermore, the little group $E\left(2\right)$ classification fails for
massive waves but can be recovered to first order in the small expansion
parameter $\gamma$, that is for slightly massive waves. Let us define a local
quasi-orthonormal null tetrad basis $\left(k,l,m,\bar{m}\right)$, as a new
basis [35, 42, 43]

$\displaystyle k$ $\displaystyle=\partial_{t}+\partial_{z}\ ,$ $\displaystyle
l$ $\displaystyle=\frac{1}{2}\left(\partial_{t}-\partial_{z}\right)\ ,$ (134)
$\displaystyle m$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\partial_{x}+i\partial_{y}\right)\ ,$
$\displaystyle\bar{m}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\partial_{x}-i\partial_{y}\right)\ ,$
(135)

which has to fulfill the relations

$\displaystyle k\cdot l=-m\cdot\bar{m}=1\ ,$ $\displaystyle k\cdot k=l\cdot
l=m\cdot m=\bar{m}\cdot\bar{m}=0\ ,$ (136) $\displaystyle k\cdot
m=k\cdot\bar{m}=l\cdot m=l\cdot\bar{m}=0\ ,$

that is, the Minkowski metric tensor $\eta_{\mu\nu}$ of signature $-2$, can be
expressed as

$\eta^{\mu\nu}=2k^{(\mu}l^{\nu)}-2m^{(\mu}\bar{m}^{\nu)}\ .$ (137)

Therefore we can raise and lower the tetrad indices by the metric of the
tetrad $\eta_{ab}$

$\eta_{ab}=\eta^{ab}=\begin{pmatrix}0&1&0&0\\\ 1&0&0&0\\\ 0&0&0&-1\\\
0&0&-1&0\end{pmatrix}\ ,$ (138)

where $(a,b)$ run over $(k,l,m,\bar{m})$. We now can decompose the Riemann
tensor into three irreducible parts, namely Ricci decomposition, related to
Weyl tensor, Ricci tensor and Ricci scalar. The four-dimensional Weyl tensor
$C_{\mu\nu\rho\sigma}$ in coordinate basis, is defined as the trace free part
of the Riemann tensor

$C_{\mu\nu\rho\sigma}=R_{\mu\nu\rho\sigma}-2g_{[\mu|[\rho}R_{\sigma]|\nu]}+\frac{1}{3}g_{\mu[\rho}g_{\sigma]\nu}R\
,$ (139)

that in tetrad form becomes

$C_{abcd}=R_{abcd}-2\eta_{[a|[c}R_{d]|b]}+\frac{1}{3}\eta_{a[c}\eta_{d]b}R\ ,$
(140)

according to tetrad components of the generic tensor $P_{abcd\dots}$ expressed
in terms of the local coordinate basis as

$P_{abcd\cdots}=P_{\alpha,\beta\gamma\delta\cdots}a^{\alpha}b^{\beta}c^{\gamma}d^{\delta}\cdots\
,$ (141)

where $(a,b,c,d,\dots)$ run over $(k,l,m,\bar{m})$. So, can be defined fifteen
Newman-Penrose quantities, expressing the Weyl tensor, the Ricci tensor and
the Ricci scalar in the null tetrad basis Eqs. (134) and (135). The NP
amplitudes are, specifically, the five Weyl-NP $\Psi$ complex scalars
$\\{\Psi_{0},\Psi_{1},\Psi_{2},\Psi_{2},\Psi_{4}\\}$, expressed in tetrad
components of the Weyl tensor as

$\displaystyle\Psi_{0}$ $\displaystyle\equiv-C_{kmkm}\ ,$ (142)
$\displaystyle\Psi_{1}$ $\displaystyle\equiv-C_{klkm}\ ,$
$\displaystyle\Psi_{2}$ $\displaystyle\equiv-C_{km\bar{m}l}\ ,$
$\displaystyle\Psi_{3}$ $\displaystyle\equiv-C_{kl\bar{m}l}\ ,$
$\displaystyle\Psi_{4}$ $\displaystyle\equiv-C_{\bar{m}l\bar{m}l}\ ,$

and the seven Ricci-NP scalars $\Phi$ and $\Lambda$, of which four real
$\\{\Phi_{00},\Phi_{11},\Phi_{22},\Lambda\\}$ and three complex
$\\{\Phi_{10},\Phi_{20},\Phi_{21}\\}$, expressed in tetrad components of Ricci
tensor and Ricci scalar as

$\displaystyle\quad\Phi_{02}\equiv-\frac{1}{2}R_{mm\ ,}$ (143)
$\displaystyle\left\\{\begin{array}[]{l}\Phi_{01}\equiv-\frac{1}{2}R_{km}\\\
\Phi_{12}\equiv-\frac{1}{2}R_{lm}\\\ \end{array}\ ,\right.$
$\displaystyle\left\\{\begin{array}[]{l}\Phi_{00}\equiv-\frac{1}{2}R_{kk}\\\
\Phi_{11}\equiv-\frac{1}{4}\left(R_{kl}+R_{m\bar{m}}\right)\\\
\Phi_{22}\equiv-\frac{1}{2}R_{ll}\end{array}\ ,\right.$
$\displaystyle\left\\{\begin{array}[]{l}\Phi_{10}\equiv-\frac{1}{2}R_{k\bar{m}}=\Phi_{01}^{*}\\\
\Phi_{21}\equiv-\frac{1}{2}R_{l\bar{m}}=\Phi_{21}^{*}\\\ \end{array}\
,\right.$
$\displaystyle\quad\Phi_{20}\equiv-\frac{1}{2}R_{\bar{m}\bar{m}}=\Phi_{02}^{*}\
,$ $\displaystyle\qquad\Lambda\equiv\frac{R}{24}\ .$

For the purpose of expand the slightly massive gravitational radiation in
terms of null plane waves, we define the four-wavevector
$\left(k_{2}^{\prime}\right)^{\mu}$ in units $c=1$ associated to the nearly
null plane wave traveling in $+z$ direction as

$k_{j}^{\mu}=\left(\omega_{j},0,0,k_{z}\right)=\omega_{j}k_{j}^{\prime\mu}\ ,$
(144)

with $j\in\\{2,\ldots,n\\}$ where

$k^{\prime\mu}_{j}=\left(1,0,0,v_{g_{j}}\right)\ .$ (145)

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

$\tilde{u}=t-v_{g}z\ ,$ (146)

that gives us

$\nabla_{\mu}\tilde{u}=\left(k_{j}^{\prime}\right)_{\mu}\ .$ (147)

So, we expand $k_{j}^{\prime}$ with respect to NP tetrads basis as

$k_{j}^{\prime\mu}=\left(1+\gamma_{k}\right)k^{\mu}+\gamma_{l}l^{\mu}+\gamma_{m}m^{\mu}+\bar{\gamma}_{m}\bar{m}^{\mu}\
,$ (148)

where the expansion coefficients $\gamma_{k},\gamma_{l},\gamma_{m}$ are of
same order of the $\gamma$, already previous defined in Eq. (124). Given the
freedom of the observer to orient his reference system, it is possible to
choose orientation in such a way that $k_{j}^{0}=k^{0}$ and $k_{j}^{3}\propto
k^{3}$, where $k^{0}$ is the angular frequency of null wave and $k^{3}$ the
third component of its vector wave. Hence, Eq. (148) gives
$\gamma_{l}=-2\gamma_{k}$ and $\gamma_{m}=0$, and this implies

$k_{j}^{\prime\mu}=k^{\mu}+\gamma_{k}\left(k^{\mu}-2l^{\mu}\right)\ ,$ (149)

or

$k_{j}^{\prime\mu}=\left(1,0,0,1+2\gamma_{k}\right)\ .$ (150)

Then, we note that the parameters $\gamma$ and $\gamma_{k}$ are of the same
order, according to

$\gamma_{k}=-\frac{1}{4}\gamma+\small{o}\left(\gamma\right).$ (151)

The derivatives of Riemann tensor can be expressed as

$R_{\alpha\beta\gamma\delta,\mu}=\frac{\partial
R_{\alpha\beta\gamma\delta}}{\partial\tilde{u}}\nabla_{\mu}\tilde{u}=\left(k_{j}^{\prime}\right)_{\mu}\dot{R}_{\alpha\beta\gamma\delta}\
,$ (152)

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

$\displaystyle R_{\alpha\beta\gamma\delta,k}$
$\displaystyle=-2\gamma_{k}\dot{R}_{\alpha\beta\gamma\delta}\ ,$ (153)
$\displaystyle R_{\alpha\beta\gamma\delta,l}$
$\displaystyle=\left(1+\gamma_{k}\right)\dot{R}_{\alpha\beta\gamma\delta}\ ,$
(154) $\displaystyle R_{\alpha\beta\gamma\delta,m}$ $\displaystyle=0\ ,$ (155)

joined with differential Bianchi identity

$R_{ab[cd,e]}=0\ ,$ (156)

to zeroth order in $\gamma_{k}$, the only non-zero tetrad components of
Riemann tensor $R_{abcd}$ that remain, are terms of form $R_{lplq}$ with
$(p,q)$ range over $(k,m,\bar{m)}$. In the nearly null plane waves
approximation, only four complex NP tetrad components are independent and non-
vanishing to first order in $\gamma_{k}$, that is, from Eq. (140), they are

$\displaystyle\Psi_{2}\left(\tilde{u}\right)$
$\displaystyle=-\frac{1}{6}R_{lklk}+\mathcal{O}\left(\gamma_{k}\right)\ ,$
(157) $\displaystyle\Psi_{3}\left(\tilde{u}\right)$
$\displaystyle=-\frac{1}{2}R_{lkl\bar{m}}+\mathcal{O}\left(\gamma_{k}\right)\
,$ (158) $\displaystyle\Psi_{4}\left(\tilde{u}\right)$
$\displaystyle=-R_{l\bar{m}l\bar{m}}\ ,$ (159)
$\displaystyle\Phi_{22}\left(\tilde{u}\right)$ $\displaystyle=-R_{lml\bar{m}}\
.$ (160)

The four amplitudes $\\{\Psi_{2},\Psi_{3},\Psi_{4},\Phi_{22}\\}$ expressed in
terms of tetrad components of metric perturbation $h_{ab}$, are given by

$\displaystyle\Psi_{2}\left(\tilde{u}\right)$
$\displaystyle=\frac{1}{12}\ddot{h}_{kk}+\mathcal{O}\left(\gamma_{k}\right)\
,$ (161) $\displaystyle\Psi_{3}\left(\tilde{u}\right)$
$\displaystyle=\frac{1}{4}\ddot{h}_{k\bar{m}}+\mathcal{O}\left(\gamma_{k}\right)\
,$ (162) $\displaystyle\Psi_{4}\left(\tilde{u}\right)$
$\displaystyle=\frac{1}{2}\ddot{h}_{\bar{m}\bar{m}}+\mathcal{O}\left(\gamma_{k}\right)\
,$ (163) $\displaystyle\Phi_{22}\left(\tilde{u}\right)$
$\displaystyle=\frac{1}{2}\ddot{h}_{m\bar{m}}+\mathcal{O}\left(\gamma_{k}\right)\
.$ (164)

These NP scalars, according to their behavior under the subgroup of Lorentz
transformations which leaves $\mathbf{k}_{j}$ unchanged, namely the little
group $E(2)$, show the following four helicity values $s$, for each of them

$\displaystyle\Psi_{2}\quad s=0\ ,$ $\displaystyle\Psi_{4\quad}s=2\ ,$ (165)
$\displaystyle\Psi_{3}\quad s=1\ ,$ $\displaystyle\Phi_{22}\quad s=0\ .$ (166)

These NP quantities allow the classification of quasi-Lorentz invariant GWs,
know as Petrov classification. For a gravitational radiation propagating along
$+\hat{z}$ axis, we can express the four NP amplitudes both in terms of the
electric components of the Riemann tensor $R_{i0j0}$ and its linearized
components from Eq. (104). Indeed, considering the identities

$\displaystyle R_{kml\bar{m}}$ $\displaystyle=0$ $\displaystyle\rightarrow$
$\displaystyle\,\,\,R_{1313}+R_{2323}=R_{0101}+R_{0202}\ ,$ (167)
$\displaystyle R_{kmk\bar{m}}$ $\displaystyle=0$ $\displaystyle\rightarrow$
$\displaystyle\,\,\,-2R_{0131}-2R_{0232}=2\left(R_{0101}+R_{0202}\right)\ ,$
(168) $\displaystyle R_{klkm}$ $\displaystyle=0$ $\displaystyle\rightarrow$
$\displaystyle\,\,\,R_{0301}=-R_{0331}\quad\text{and}\quad R_{0302}=-R_{0332}\
,$ (169) $\displaystyle R_{kmlm}$ $\displaystyle=0$ $\displaystyle\rightarrow$
$\displaystyle\,\,\,R_{3132}=R_{0102}\quad\text{and}\quad
R_{3131}-R_{3232}=R_{0101}-R_{0202}\ ,$ (170) $\displaystyle R_{kmkm}$
$\displaystyle=0$ $\displaystyle\rightarrow$
$\displaystyle\,\,\,\begin{cases}2R_{3101}-2R_{3202}=-R_{3131}+R_{3232}-R_{0101}+R_{0202}\\\
R_{0102}=-R_{0132}-R_{0231}-R_{3132}\end{cases},$ (171)

occurs in terms of $R_{i0j0}$ that

$\displaystyle\Psi_{2}\left(\tilde{u}\right)=$
$\displaystyle-\frac{1}{6}R_{0303}+\mathcal{O}\left(\gamma_{k}\right),$ (172)
$\displaystyle\Psi_{3}\left(\tilde{u}\right)=$
$\displaystyle-\frac{1}{2\sqrt{2}}R_{0301}+\frac{1}{2\sqrt{2}}iR_{0302}+\mathcal{O}\left(\gamma_{k}\right),$
(173) $\displaystyle\Psi_{4}\left(\tilde{u}\right)=$
$\displaystyle-\frac{1}{2}R_{0101}+\frac{1}{2}R_{0202}+iR_{0102}+\mathcal{O}\left(\gamma_{k}\right),$
(174) $\displaystyle\Phi_{22}\left(\tilde{u}\right)=$
$\displaystyle-\frac{1}{2}R_{0101}-\frac{1}{2}R_{0202}+\mathcal{O}\left(\gamma_{k}\right),$
(175)

while, in term of metric perturbation, the NP scalars reads

$\displaystyle\Psi_{2}\left(\tilde{u}\right)=$
$\displaystyle-\frac{1}{12}\left(2h_{03,03}-h_{00,33}-h_{33,00}\right)+\mathcal{O}\left(\gamma_{k}\right),$
(176) $\displaystyle\Psi_{3}\left(\tilde{u}\right)=$
$\displaystyle-\frac{1}{4\sqrt{2}}\left(h_{01,03}-h_{13,00}\right)+\frac{1}{4\sqrt{2}}i\left(h_{02,03}-h_{23,00}\right)+\mathcal{O}\left(\gamma_{k}\right),$
(177) $\displaystyle\Psi_{4}\left(\tilde{u}\right)=$
$\displaystyle\frac{1}{4}\left(h_{11,00}-h_{22,00}\right)-2ih_{12,00}+\mathcal{O}\left(\gamma_{k}\right),$
(178) $\displaystyle\Phi_{22}\left(\tilde{u}\right)=$
$\displaystyle\frac{1}{4}\left(h_{11,00}+h_{22,00}\right)+\mathcal{O}\left(\gamma_{k}\right)\
.$ (179)

From Eqs. (110), for a general congruence of non-null geodesics of GWs, we
obtain, for a massless mode $\omega_{1}$ and $n-1$ massive modes $\omega_{j}$,
at $\mathbf{k}$ fixed to lowest order in $\gamma_{k}$, the following
expressions

$\displaystyle\Psi_{2}\left(\tilde{u}\right)$
$\displaystyle=\mathcal{O}\left(\gamma_{k}\right),$ (180)
$\displaystyle\Psi_{3}\left(\tilde{u}\right)$
$\displaystyle=\mathcal{O}\left(\gamma_{k}\right)\ ,$
$\displaystyle\Psi_{4}\left(\tilde{u}\right)$
$\displaystyle=-\omega_{1}^{2}\biggl{[}\Bigl{(}\epsilon^{(+)}\left(\omega_{1}\right)e^{i\omega_{1}\left(t-z\right)}+c.c.\Bigr{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-i\left(\epsilon^{(\times)}\left(\omega_{1}\right)e^{i\omega_{1}\left(t-z\right)}+c.c.\right)\biggr{]}+\mathcal{O}\left(\gamma_{k}\right)\
,$ $\displaystyle\Phi_{22}\left(\tilde{u}\right)$
$\displaystyle=\frac{k_{z}^{2}}{8}\sum_{j=2}^{n}\tilde{A}_{j}\left(k_{z}\right)e^{ik_{z}\left(t-z\right)}+c.c.+\mathcal{O}\left(\gamma_{k}\right)\
,$

according to

$\omega_{j}=k_{z}+\mathcal{O}\left(\gamma_{k}\right)\ ,$ (181)

in units $c=1$. Thereby, from Eqs. (180) follows that, to zeroth order in
$\gamma_{k}$, the $\omega_{1}$ mode has helicity 2 because only $\Psi_{4}\neq
0$ while $\omega_{j}$ modes have helicity 0 because $\Phi_{22}\neq 0$. However
the quasi-Lorentz invariant $E(2)$ class of GWS generated by non-local
gravitational theory $R+\sum_{h=1}^{n}a_{h}R\Box^{-h}R$ is $N_{3}$, because
for a generic non-null wave the NP quantities are $\Psi_{2}=\Psi_{3}=0$ and
$\Psi_{4}\neq 0\neq\Phi_{22}$ to zeroth order in $\gamma_{k}$. This means that
the presence or absence of all modes is observer-independent. The components
of driving-force matrix $S(t)$ can be expressed in terms of the six electrics
components $R_{i0j0}$, as

$S_{ij}(t)\equiv R_{i0j0}\ ,$ (182)

and the six amplitudes with direction $\hat{\mathbf{k}}=\mathbf{e}_{z}$ as

$\displaystyle p_{1}\left(\mathbf{e}_{z},t\right)$
$\displaystyle\equiv\Psi_{2}\ ,$ (183) $\displaystyle
p_{2}\left(\mathbf{e}_{z},t\right)$ $\displaystyle\equiv\text{Re}\;\Psi_{3}\
,$ $\displaystyle p_{3}\left(\mathbf{e}_{z},t\right)$
$\displaystyle\equiv\text{Im}\;\Psi_{3}\ ,$ $\displaystyle
p_{4}\left(\mathbf{e}_{z},t\right)$ $\displaystyle\equiv\text{Re}\;\Psi_{4}\
,$ $\displaystyle p_{5}\left(\mathbf{e}_{z},t\right)$
$\displaystyle\equiv\text{Im}\;\Psi_{4}\ ,$ $\displaystyle
p_{6}\left(\mathbf{e}_{z},t\right)$ $\displaystyle\equiv\Phi_{22}\ .$

In matrix form, $S(t)$ i terms of six polarization matrices
$W_{A}(\mathbf{e}_{z})$ defined by

$.\begin{aligned} W_{1}\left(\mathbf{e}_{z}\right)=&-6\begin{pmatrix}0&0&0\\\
0&0&0\\\ 0&0&1\end{pmatrix}\
,&W_{2}\left(\mathbf{e}_{z}\right)=&-2\sqrt{2}\begin{pmatrix}0&0&1\\\ 0&0&0\\\
1&0&0\end{pmatrix}\ ,\\\
W_{3}\left(\mathbf{e}_{z}\right)=&2\sqrt{2}\begin{pmatrix}0&0&0\\\ 0&0&1\\\
0&1&0\end{pmatrix}\
,&W_{4}\left(\mathbf{e}_{z}\right)=&-\begin{pmatrix}1&0&0\\\ 0&-1&0\\\
0&0&0\end{pmatrix}\ ,\\\
W_{5}\left(\mathbf{e}_{z}\right)=&\begin{pmatrix}0&1&0\\\ 1&0&0\\\
0&0&0\end{pmatrix}\
,&W_{6}\left(\mathbf{e}_{z}\right)=&-\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&0\end{pmatrix}\ \end{aligned}$ (184)

yields

$S\left(t\right)=\sum_{A}p_{A}\left(\mathbf{e}_{z},t\right)W_{A}\left(\mathbf{e}_{z}\right)\
,$ (185)

where the index $A$ ranges over $\\{1,2,3,4,5,6\\}$ [42, 43].

Accordingly, 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 [37, 38, 25, 44]. Then, the six polarization
amplitudes $p_{A}\left(\mathbf{e}_{z},t\right)$ in terms of the NP scalars for
our non-local waves read

$\displaystyle p_{1}^{\left(l\right)}\left(\mathbf{e}_{z},t\right)$
$\displaystyle=\mathcal{O}\left(\gamma_{k}\right)\ ,$ (186) $\displaystyle
p_{2}^{\left(x\right)}\left(\mathbf{e}_{z},t\right)$
$\displaystyle=p_{3}^{\left(y\right)}\left(\mathbf{e}_{z},t\right)=\mathcal{O}\left(\gamma_{k}\right)\
,$ $\displaystyle p_{4}^{\left(+\right)}\left(\mathbf{e}_{z},t\right)$
$\displaystyle=-\omega_{1}^{2}\tilde{\epsilon}^{\left(+\right)}\left(\omega_{1}\right)e^{i\omega_{1}t}+c.c.+\mathcal{O}\left(\gamma_{k}\right)\
,$ $\displaystyle p_{5}^{\left(\times\right)}\left(\mathbf{e}_{z},t\right)$
$\displaystyle=\omega_{1}^{2}\tilde{\epsilon}^{\left(\times\right)}\left(\omega_{1}\right)e^{i\omega_{1}t}+c.c.+\mathcal{O}\left(\gamma_{k}\right)\
,$ $\displaystyle p_{6}^{\left(b\right)}\left(\mathbf{e}_{z},t\right)$
$\displaystyle=\frac{k_{z}^{2}}{8}\sum_{j=2}^{n}\tilde{A}_{j}\left(k_{z}\right)e^{ik_{z}t}+c.c.+\mathcal{O}\left(\gamma_{k}\right)\
.$

To first order in $\gamma_{k}$, from Eqs. (186) only the amplitudes
$p_{4}^{\left(+\right)}$, $p_{5}^{\left(\times\right)}$ and
$p_{6}^{\left(b\right)}$ survive, while others are suppressed. Precisely
$p_{4}^{\left(+\right)}$ and $p_{5}^{\left(\times\right)}$ correspond to
amplitudes of plus and cross modes of frequency $\omega_{1}$ predicted by
general relativity while the amplitude $p_{6}^{\left(b\right)}$ is associated
to the $j$ transverse breathing scalar modes of frequency $\omega_{j}$, always
to lowest order. In a nutshell, the gravitational radiation for non-local
gravity $R+\sum_{h=1}^{n}a_{h}R\Box^{-h}R$ exhibits three polarizations, i.e.
two tensor and, generally, one scalar, all transverse and $n+1$ modes.
Specifically, is always present the standard radiation that is the associated
with two modes, $(+)$ and $(\times)$, massless 2-helicity transverse and
tensor, of frequency $\omega_{1}$, governed by two degrees of freedom (d.o.f.)
$\tilde{\epsilon}^{(+)}(\omega_{1})$ and
$\tilde{\epsilon}^{(\times)}(\omega_{1})$. In addition under certain
conditions are present $n-1$ massive 0-helicity transverse scalar modes of
frequency $\omega_{j}$, namely the breathing mode, governed by $n-1$ d.o.f.
$\tilde{A}_{j}(k_{z})$. It is worth observing that the same approach can be
used also for higher order theories of gravity both in Riemannian,
teleparallel and non-metric framework. For details see [45, 46].

##  6 Tables

A classification of solutions corresponding to null and non-null gravitational
waves in higher order non-local gravity under some constraints, are summarized
in following Tables 1 and 2. In the Table 3 instead, we point out Lagrangians
in which the scalar radiation is prohibited while in Table 4 are summarized
the count of massive scalar waves. The polarizations and helicities of GWs
with their related modes, taking into account the constraints, are showed in
the following Tables 5 and 6. Finally, the main results of this paper on
gravitational waves in higher order non-local gravity are summarized in the
Table 7, where we assumed $\phi_{1}^{(0)}=\lambda_{1}^{(0)}=0$.

Table 1: Non-null wavelike solutions in $\Box^{-n}$-gravity. Gravitational | $B=6a_{1}$ | Solutions | $B\neq 6a_{1}$ | Solutions
---|---|---|---|---
Lagrangian | $k^{2}\neq 0$ |  | $k^{2}\neq 0$ |
$R+a_{1}R\Box^{-1}R$ |  $\infty^{1}$ solutions | $k^{2}\in\mathbb{R}_{+}$ | no solution |
$R+\sum_{h=1}^{2}a_{h}R\Box^{-h}R$ | no solution |  | one solution | $k_{2}^{2}=\frac{6a_{2}}{6a_{1}-B}>0$
$R+\sum_{h=1}^{3}a_{h}R\Box^{-h}R$ | one solution | $k_{2}^{2}=\frac{a_{3}}{a_{2}}>0$ | two solutions | $k_{2,3}^{2}=\frac{\pm\sqrt{3a_{2}^{2}-12a_{1}a_{3}+2a_{3}B}-3a_{2}}{B-6a_{1}}>0$
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$
$R+\sum_{h=1}^{n}a_{h}R\Box^{-h}R$ |  $n-2$ solutions | $k_{2}^{2},\cdots,k_{n-1}^{2}>0$ |  $n-1$ solutions | $k_{2}^{2},\cdots,k_{n}^{2}>0$
Table 2: Null wavelike solutions in $\Box^{-n}$ -gravity.  Gravitational Lagrangian | Solutions |
---|---|---
$R+\sum_{h=1}^{n}a_{h}R\Box^{-h}R$ |  $k_{1}^{2}=0$ for all $n$ |
Table 3: Lagrangians where the scalar radiation is forbidden. Gravitational | Constraint | GWs
---|---|---
Lagrangian |  |
$R+a_{1}R\Box^{-1}R$ | $B\neq 6a_{1}$ |  $(+)$ and $(\times)$ tensor modes
$R+a_{1}R\Box^{-1}R+a_{2}R\Box^{-2}R$ | $B=6a_{1}$ |  $(+)$ and $(\times)$ tensor modes
Table 4: Number of massive transverse scalar modes with helicity $0$ in higher order non-local gravity. Constraint | Number of modes or d.o.f. | Order of $\Box^{-n}$
---|---|---
$B\neq 6a_{1}$ | $n-1$ | $n\geq 1$
$B=6a_{1}$ | $n-2$ | $n\geq 2$
$B=6a_{1}$ | $\infty$ | $n=1$
Table 5: Polarizations and modes of gravitational waves in higher order non-local gravity for $B\neq 6a_{1}$. Gravitational | Conditions | Frequency | Polarization | Type | d.o.f. | Helicity | Mass
---|---|---|---|---|---|---|---
Lagrangian |  |  | modes |  |  |  |
$R+\sum_{h=1}^{n}a_{h}R\Box^{-h}R$0 |  0$n\geq 1$ |  0$\omega_{1}$ | transverse | tensor | 2 | 2 | 0
| $k_{1}^{2}=0$ |  | $(+),(\times)$ |  |  |  |
|  0$n\geq 2$ |  0$\omega_{2},\ldots,\omega_{n}$ | transverse | scalar | $n-1$ | 0 | $M_{2},\ldots,M_{n}$
| $k_{j}^{2}\neq 0$ |  | $(b)$ |  |  |  |
Table 6: Polarizations and modes of gravitational waves in higher order non-local gravity for $B=6a_{1}$. Gravitational | Conditions | Frequency | Polarization | Type | d.o.f. | Helicity | Mass
---|---|---|---|---|---|---|---
Lagrangian |  |  | modes |  |  |  |
$R+\sum_{h=1}^{n}a_{h}R\Box^{-h}R$0 |  0$n\geq 1$ |  0$\omega_{1}$ | transverse | tensor | 2 | 2 | 0
| $k_{1}^{2}=0$ |  | $(+),(\times)$ |  |  |  |
degenerate state 00 |  0$n=1$ |  0$\omega\in\mathbb{R}_{+}$ | transverse | scalar | $\infty$ | 0 | $M\in\mathbb{R}_{+}$
| $k^{2}\neq 0$ |  | $(b)$ |  |  |  |
| $n=2$ | no solutions |  |  |  |  |
| $k^{2}\neq 0$ |  |  |  |  |  |
|  0$n\geq 3$ |  0$\omega_{2},\ldots,\omega_{n-1}$ | transverse | scalar | $n-2$ | 0 | $M_{2},\ldots,M_{n-1}$
| $k_{j}^{2}\neq 0$ |  | $(b)$ |  |  |  |
Table 7: Polarizations and modes of gravitational waves in in higher order non-local gravity.  Gravitational | Constraint | Gravitational wave | Properties
---|---|---|---
Lagrangian |  | solutions |
$R+a_{1}R\Box^{-1}R$ | $a_{1}\neq\frac{1}{6}$ | two $(+)$ and $(\times)$ modes $k_{1}^{2}=0$ | tensor, transverse,massless, $2$-helicity
$n=1$ |  |  |
| $a_{1}=\frac{1}{6}$ | two $(+)$ and $(\times)$ modes $k_{1}^{2}=0$ | tensor, transverse, massless, $2$-helicity
|  |  $\infty$ breathing modes $k^{2}\neq 0$ | scalar, transverse, massive, $0$-helicity
$R+\sum_{h=1}^{2}a_{h}R\Box^{-h}R$ | $a_{1}\neq\frac{1}{6}$ | two $(+)$ and $(\times)$ modes $k_{1}^{2}=0$ | tensor, transverse, massless, $2$-helicity
|  | one breathing mode $k_{2}^{2}\neq 0$ | scalar, transverse, massive, $0$-helicity
$n=2$ |  |  |
| $a_{1}=\frac{1}{6}$ | two $(+)$ and $(\times)$ modes $k_{1}^{2}=0$ | tensor, transverse, massless, $2$-helicity
$R+\sum_{h=1}^{n}a_{h}R\Box^{-h}R$ | $a_{1}\neq\frac{1}{6}$ | two $(+)$ and $(\times)$ modes $k_{1}^{2}=0$ | tensor, transverse, massless, $2$-helicity
|  |  $(n-1)$ breathing modes $k_{j}^{2}\neq 0$ | scalar, transverse, massive, $0$-helicity
$n\geq 3$ |  |  |
| $a_{1}=\frac{1}{6}$ | two $(+)$ and $(\times)$ modes $k_{1}^{2}=0$ | tensor, transverse, massless, $2$-helicity
|  |  $(n-2)$ breathing modes $k_{j}^{2}\neq 0$ | scalar, transverse, massive, $0$-helicity

##  7 Conclusions

The most interesting result of this paper is the presence of transverse scalar
gravitational waves under certain conditions in some orders of gravity, that
being massive propagate at a speed $v_{g}$ below $c$. More specifically, after
localizing the non-local action and deriving the related field equations by a
variational principle, we linearized them weakly perturbing the tensor and
scalar fields, because suppose to study the radiation at very high distances
from the sources that had emitted it. Therefore, the linearized equations in
vacuum in a suitable gauge were resolved and we got non-null gravitational
waves $k^{2}\neq 0$. Indeed, while the tensor radiation predicted by Einstein
is always present in gravity of any order under any constraint, further scalar
radiation is prohibited in $\Box^{-1}$ gravity for $B\neq 6a_{1}$ and in
$\Box^{-2}$ gravity for $B=6a_{1}$, see Table 3. If $n\geq 2$ and $B\neq
6a_{1}$ instead, the non-local theory exhibits $n-1$ wavelike scalar solutions
$k^{2}\neq 0$, where each $\Box^{-1}$ introduces a single massive scalar
solution. While if $n\geq 2$ and $B=6a_{1}$, a single $k^{2}\neq 0$ scalar
mode is suppressed and the non-local model shows $n-2$ wavelike scalar
solutions. It is crucial to point out that under constraint $B=6a_{1}$, the
$\Box^{-1}$ gravity shows a degenerate behavior in the meaning that any
massive scalar waves $k^{2}\neq 0$ is solution of linearized filed equations
in vacuum and under further assumption $\phi_{1}^{(0)}=\Lambda_{1}^{(0)}=0$
which implies $a_{1}=1/6$, the non-local $\Box^{-1}$ action restricted to two-
dimensional spacetime is reduced to the Polyakov effective action. The
polarizations of this radiation were analyzed both by means of geodesic
deviation equation and Newmann-Penrose formalism. Thus, a transverse breathing
polarization was found for each $j$-th modes corresponding to each angular
frequency $\omega_{j}$ to lowest order in $\gamma$, a parameter that takes
into account the difference in speed between the slightly massive wave and the
massless one. The massive waves have helicity equal to zero and belong to
$N_{3}$ class, according to $E(2)$ invariant classification of Petrov.

It would be very interesting in the future to study also the production of
gravitational waves and not only their propagation in higher order non-local
gravity, even in theories non-curvature based such as those torsion or non-
metricity based. Being able to measure the further massive scalar breathing
polarization could be a signature of the viable theory of gravitation.

## Acknowledgments

MC acknowledges the Istituto Nazionale di Fisica Nucleare (INFN) Sez. di
Napoli, Iniziative Specifiche QGSKY, and the Istituto Nazionale di Alta
Matematica (INdAM), gruppo GNFM, for the support.

## References

  * [1] L. Buoninfante, G. Lambiase, L. Petruzziello, _Eur.Phys.J.C_ 81 (2021) 10, 928.
  * [2] L. Buoninfante, G. Lambiase, Y. Miyashita, W. Takebe, M. Yamaguchi, _Phys.Rev.D_ 101 (2020) 8, 084019.
  * [3] L. Buoninfante, G. Lambiase, M. Yamaguchi, _Phys.Rev.D_ 100 (2019) 2, 026019.
  * [4] L. Buoninfante, A. Ghoshal, G. Lambiase, A. Mazumdar, _Phys.Rev.D_ 99 (2019) 4, 044032.
  * [5] L. Buoninfante, A.S. Cornell, G. Harmsen, A.S. Koshelev, G. Lambiase, A. Mazumadra, _Phys.Rev.D_ 98 (2018) 8, 084041 .
  * [6] L Buoninfante, G. Lambiase, A. Mazumdar, _Nucl.Phys.B_ 944 (2019), 114646.
  * [7] L. Buoninfante, A.S. Koshelev, G. Lambiase, A. Mazumdar, _JCAP_ 09 (2018), 034.
  * [8] G. V. Efimov, Non-Local Quantum Theory of the Scalar Field, _Commun. math. Phys._ 5 (1967), 42.
  * [9] B. Mashhoon, _Nonlocal Gravity_ , International Series of Monographs on Physics 167 (Oxford University Press, 2017).
  * [10] B. Mashhoon, Toward a non-local theory of gravitation, _Annalen Phys._ 16 (2007), 57.
  * [11] L.D. Landau and E.M. Lifshitz, _Electrodynamics of Continuous Media_ , Pergamon, Oxford (1960).
  * [12] J.D. Jackson, _Classical Electrodynamics_ 3rd edn, Wiley, Hoboken, NJ. (1999).
  * [13] S. Capozziello and F. Bajardi, Nonlocal gravity cosmology: An overview, _Int. J. Mod. Phys. D_(2021), 2230009. https://doi.org/10.1142/S0218271822300099
  * [14] L. Modesto and L. Rachwał, Nonlocal quantum gravity: A review, _Int. J. Mod. Phys. D_ 26 (2017), 1730020.
  * [15] F. Bouchè, S. Capozziello, V. Salzano, K. Umetsu, Testing non-local gravity by clusters of galaxies, https://doi.org/10.48550/arXiv.2205.03216.
  * [16] G. Calcagni, _Classical and Quantum Cosmology_ , Graduate Texts in Physics Springer, Switzerland (2017).
  * [17] L. Modesto, I.L. Shapiro, Super-renormalizable quantum gravity with complex ghosts, _Phys.Lett. B_ 755 (2016), 279.
  * [18] S. Nojiri and S. D. Odintsov, Modified non-local-F(R) gravity as the key for the inflation and dark energy, _Physics Letters B_ 659 (2008), 821.
  * [19] S. Deser and R. P. Woodard, Nonlocal Cosmology, _Phys. Rev. Lett._ 99 (2007), 111301.
  * [20] S. Deser and R.P. Woodard, Nonlocal cosmology II. Cosmic acceleration without fine tuning or dark energy, _JCAP_ 06 (2019), 034.
  * [21] S. Deser and R.P. Woodard, Observational Viability and Stability of Nonlocal Cosmology, _J. Cosmol. Astropart. Phys._ 1311 (2013), 036.
  * [22] R. V. Wagoner, Scalar-Tensor Theory and Gravitational Waves, _Phys. Rev. D_ 1 (1970), 3209.
  * [23] H. J. Schmidt, Variational derivatives of arbitrarily high order and multi-inflation cosmological models, Class. Quantum Grav. 7 (1990), 1023.
  * [24] S. Capozziello, M. Capriolo and L. Caso, Weak field limit and gravitational waves in higher order gravity, Int. J. Geom. Methods Mod. Phys. 16 No.03 (2019), 1950047.
  * [25] C. M. Will, The Confrontation between General Relativity and Experiment, _Living Rev. Relativity_ 17 (2014), 4.
  * [26] B. P. Abbott et al., LIGO Scientific and Virgo Collaborations, GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, _Phys.Rev.Lett._ 119 (2017), 161101.
  * [27] B. P. Abbott et al., LIGO Scientific and Virgo Collaborations, Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A, _ApJ Lett._ 848 (2017), L13.
  * [28] B. P. Abbott et al., LIGO Scientific and Virgo Collaborations, Search for Tensor, Vector, and Scalar Polarizations in the Stochastic Gravitational-Wave Background, _Phys. Rev. Lett._ 120 (2018), 201102.
  * [29] S. Capozziello and M. Capriolo , Gravitational waves in non-local gravity, Class. Quantum Grav. 38 (2021), 175008.
  * [30] J. Boos, Effects of Non-locality in Gravity and Quantum Theory, Ph.D. Thesis, University of Alberta Edmonton Canada, 2020.
https://doi.org/10.48550/arXiv.2009.10856

  * [31] A. Codello, Polyakov effective action from functional renormalization group equation, _Annals of Physics_ 325 (2010), 1727.
  * [32] N. Straumann, _General Relativity_ , (Springer-Verlag, Dordrecht, 2013).
  * [33] T. Katsuragawa, T. Nakamura, T. Ikeda and S. Capozziello, Gravitational Waves in $F(R)$ Gravity: Scalar Waves and the Chameleon Mechanism, _Phys. Rev. D_ 99 (2019), 124050.
  * [34] C. M. Will, Theory and Experiment in Gravitational Physics (Cambridge University Press, Cambridge, 2018). doi:10.1017/9781316338612
  * [35] E. Poisson and C. M. Will, Gravity Newtonian, Post-Newtonian, Relativistic (Cambridge University Press, Cambridge, 2014). doi:10.1017/CBO9781139507486.
  * [36] S. Capozziello, M. Capriolo, S. Nojiri, Considerations on gravitational waves in higher-order local and non-local gravity, Physics Letters B 810 (2020), 135821.
  * [37] M. E. S. Alves, O. D. Miranda, and J. C. N. de Araujo, Probing the f(R) formalism through gravitational wave polarizations, _Phys. Lett. B_ 679 (2009), 401.
  * [38] C. Bogdanos, S. Capozziello, M. De Laurentis and S. Nesseris, Massive, massless and ghost modes of gravitational waves from higher-order gravity, _Astropart. Phys._ 34 (2010), 236 .
  * [39] S. Capozziello, C. Corda, and M. De Laurentis, Massive gravitational waves from f(R) theories of gravity: Potential detection with LISA, _Phys. Lett. B_ 669 (2008), 255.
  * [40] S. Capozziello, R. Cianci, M. De Laurentis and S. Vignolo, Testing metric-affine $f(R)$-gravity by relic scalar gravitational waves, _Eur. Phys. J. C_ 70 (2010), 341.
  * [41] E. Newman and R. Penrose, An approach to gravitational radiation by a method of spin coefficients, _J. Math. Phys._ 3 (1962), 566.
  * [42] D. M. Eardley, D. L. Lee, and A. P. Lightman, Gravitational-Wave Observations as a Tool for Testing Relativistic Gravity, _Phys. Rev. D_ 8 (1973), 3308.
  * [43] D. M. Eardley, D. L. Lee, A. P. Lightman, R. V. Wagoner and C. M. Will, Gravitational-Wave Observations as a Tool for Testing Relativistic Gravity, _Phys. Rev. Lett._ 30 (1973), 884.
  * [44] M. Maggiore and A. Nicolis, Detection strategies for scalar gravitational waves with interferometers and resonant spheres, _Phys. Rev. D_ 62 (2000), 024004.
  * [45] S. Capozziello, M. Capriolo and L. Caso, Weak field limit and gravitational waves in f(T, B) teleparallel gravity, _Eur. Phys. J. C_ 80 (2020), 156.
  * [46] S. Capozziello, M. Capriolo and L. Caso, Gravitational waves in higher order teleparallel gravity, Class. Quantum Grav. 37 (2020), 235013.


Intrinsic charm in the nucleon and forward production of charm:
a new constrain from IceCube Neutrino Observatory

Rafał Maciuła1$\star$, Victor P. Goncalves2 and Antoni Szczurek1,3

1 Institute of Nuclear Physics, Polish Academy of Sciences, ul. Radzikowskiego
152, PL-31-342 Kraków, Poland

2 Instituto de Física e Matemática, Universidade Federal de Pelotas (UFPel),

Caixa Postal 354, CEP 96010-900, Pelotas, RS, Brazil

3 University of Rzeszów, PL-35-959 Rzeszów, Poland

*<EMAIL_ADDRESS>
Proceedings for the XXVIII International Workshop
on Deep-Inelastic Scattering and Related Subjects,
Stony Brook University, New York, USA, 12-16 April 2021
10.21468/SciPostPhysProc.?

## Abstract

The predictions for the atmospheric neutrino flux at high energies strongly
depend on the contribution of prompt neutrinos, which are determined by the
production of charmed meson in the atmosphere at very forward rapidities. Here
we estimate the related cross sections taking into account the presence of an
intrinsic charm (IC) component in the proton wave function. The impact on the
predictions for the prompt neutrino flux is investigated assuming different
values for the probability to find the IC in the nucleon.

##  1 Introduction

Recent experimental results obtained by the LHC, the Pierre Auger and IceCube
Neutrino Observatories have challenged understanding of many interesting
aspects of Quantum Chromodynamics in the high energy limit. In particular, in
recent years, IceCube measured the astrophysical and atmospheric neutrinos
fluxes at high energies (see e.g. [1]) while different collaborations from the
LHC performed several analyses of the heavy meson production at high energies
and forward rapidities (see e.g. [2]). Those different origin data sets are
strictly interrelated, since the description of the heavy meson production at
the LHC and higher center of mass energies is fundamental to make precise
predictions of the prompt neutrino flux [3], which is expected to dominate the
atmospheric $\nu$ flux for large neutrino energies [4].

This aspect motivates the development of new and/or more precise approaches to
describe the perturbative and nonperturbative regimes of the Quantum
Chromodynamics (QCD) needed to describe the charmed meson production in a
kinematical range beyond that reached in hadronic collisions at the LHC. For
this new kinematical range, some topics are theme of intense debate. An
important question, which motivates the present study, is whether the current
and future IceCube data can put some constraints on the intrinsic charm
concept in the nucleon.

##  2 Formalism

The atmospheric neutrinos are produced in cosmic-ray interactions with nuclei
in Earth’s atmosphere [4]. While at low neutrino energies ($E_{\nu}\lesssim
10^{5}$ GeV), these neutrinos arise from the decay of light mesons (pions and
kaons), and the associated flux is denoted as the conventional atmospheric
neutrino flux [5], for larger energies it is expected that the prompt
atmospheric neutrino flux associated with the decay of hadrons containing
heavy flavours become important [6].

Calculations of the prompt atmospheric neutrino flux at the detector level
depend on the description of the production and decay of the heavy hadrons as
well as the propagation of the associated particles through the atmosphere.
Following our previous studies [3, 7], we estimated the expected prompt
neutrino flux in the detector $\phi_{\nu}$ using the $Z$-moment method [6],
which implies that $\phi_{\nu}$ can be estimated using the geometric
interpolation formula

$\displaystyle\phi_{\nu}=\sum_{H}\frac{\phi_{\nu}^{H,low}\cdot\phi_{\nu}^{H,high}}{\phi_{\nu}^{H,low}+\phi_{\nu}^{H,high}}\,\,.$
(1)

where $H=D^{0},D^{+},D_{s}^{+}$, $\Lambda_{c}$ for charmed hadrons and
${\phi_{\nu}^{H,low}}$ and ${\phi_{\nu}^{H,high}}$ are solutions of a set of
coupled cascade equations for the nucleons, heavy meson and lepton (and their
antiparticles) fluxes in the low- and high-energy ranges, respectively. They
can be expressed in terms of the nucleon-to-hadron ($Z_{NH}$), nucleon-to-
nucleon ($Z_{NN}$), hadron-to-hadron ($Z_{HH}$) and hadron-to-neutrino
($Z_{H\nu}$) $Z$-moments, as follows [6]:

$\displaystyle\phi_{\nu}^{H,low}$ $\displaystyle=$
$\displaystyle\frac{Z_{NH}(E)\,Z_{H\nu}(E)}{1-Z_{NN}(E)}\phi_{N}(E,0)\,,$ (2)
$\displaystyle\phi_{\nu}^{H,high}$ $\displaystyle=$
$\displaystyle\frac{Z_{NH}(E)\,Z_{H\nu}(E)}{1-Z_{NN}(E)}\frac{\ln(\Lambda_{H}/\Lambda_{N})}{1-\Lambda_{N}/\Lambda_{H}}\frac{m_{H}ch_{0}}{E\tau_{H}}f(\theta)\,\phi_{N}(E,0)\,,$
(3)

where $\phi_{N}(E,0)$ is the primary flux of nucleons in the atmosphere,
$m_{H}$ is the decaying particle’s mass, $\tau_{H}$ is the proper lifetime of
the hadron, $h_{0}=6.4$ km, $f(\theta)\approx 1/\cos\theta$ for
$\theta<60^{o}$, and the effective interaction lengths $\Lambda_{i}$ are given
by $\Lambda_{i}=\lambda_{i}/(1-Z_{ii})$, with $\lambda_{i}$ being the
associated interaction length ($i=N,H$). For $Z_{H\nu}$, our treatment of the
semileptonic decay of $D$-hadrons follows closely Ref. [8]. For a detailed
discussion of the cascade equations, see e.g. Ref. [6]. Assuming that the
incident flux can be represented by protons ($N=p$), the charmed hadron
$Z$-moments are given by

$\displaystyle
Z_{pH}(E)=\int_{0}^{1}\frac{dx_{F}}{x_{F}}\frac{\phi_{p}(E/x_{F})}{\phi_{p}(E)}\frac{1}{\sigma_{pA}(E)}\frac{d\sigma_{pA\rightarrow
H}(E/x_{F})}{dx_{F}}\,\,,$ (4)

where $E$ is the energy of the produced particle (charmed meson), $x_{F}$ is
the Feynman variable, $\sigma_{pA}$ is the inelastic proton-Air cross section
and $d\sigma/dx_{F}$ is the differential cross section for the charmed meson
production.

As discussed in Ref. [9], the cross section for charm production at large
forward rapidities, which is the region of interest for estimating the prompt
$\nu_{\mu}$ flux [3], can be expressed as

$\displaystyle d\sigma_{pp\rightarrow charm}=d\sigma_{pp\rightarrow
charm}(gg\rightarrow c\bar{c})+d\sigma_{pp\rightarrow charm}(cg\rightarrow
cg)\,\,,$ (5)

where the first and second terms represent the contributions associated with
the $gg\rightarrow c\bar{c}$ and $cg\rightarrow cg$ mechanisms, with the
corresponding expressions depending on the factorization scheme assumed in the
calculations. In Ref. [9], a detailed comparison between the collinear, hybrid
and $k_{T}$-factorization approaches was performed. In what follows, we will
focus on the hybrid factorization model. In this approach, the differential
cross sections for $gg^{*}\rightarrow c\bar{c}$ and $cg^{*}\rightarrow cg$
mechanisms are given by

$\displaystyle d\sigma_{pp\rightarrow charm}(gg\rightarrow c\bar{c})=\int
dx_{1}\int\frac{dx_{2}}{x_{2}}\int
d^{2}k_{t}\,g(x_{1},\mu^{2})\,{\cal{F}}_{g^{*}}(x_{2},k_{t}^{2},\mu^{2})\,d\hat{\sigma}_{gg^{*}\rightarrow
c\bar{c}}\;,$ (6) $\displaystyle d\sigma_{pp\rightarrow charm}(cg\rightarrow
cg)=\int dx_{1}\int\frac{dx_{2}}{x_{2}}\int
d^{2}k_{t}\,c(x_{1},\mu^{2})\,{\cal{F}}_{g^{*}}(x_{2},k_{t}^{2},\mu^{2})\,d\hat{\sigma}_{cg^{*}\rightarrow
cg}\,\,,$ (7)

where $g(x_{1},\mu^{2})$ and $c(x_{1},\mu^{2})$ are the collinear PDFs in the
projectile, ${\cal{F}}_{g^{*}}(x_{2},k_{t}^{2},\mu^{2})$ is the unintegrated
gluon distribution (gluon uPDF) of the proton target, $\mu^{2}$ is the
factorization scale of the hard process and the subprocesses cross sections
are calculated assuming that the small-$x$ gluon is off mass shell and are
obtained from a gauge invariant tree-level off-shell amplitude. In our
calculations $c(x_{1},\mu^{2})$, similarly $\bar{c}(x_{1},\mu^{2})$, contain
the intrinsic charm component.

As emphasized in Ref. [9], the hybrid model, already at leading-order, takes
into account radiative higher-order corrections associated with extra hard
emissions that are resummed by the gluon uPDF. In the numerical calculations
below the intrinsic charm PDFs are taken at the initial scale $m_{c}=1.3$ GeV,
so the perturbative charm contribution is intentionally not taken into account
when discussing IC contributions.

Considering the $cg^{*}\rightarrow cg$ mechanism one has to deal with the
massless partons (minijets) in the final state. The relevant formalism with
massive partons is not yet available. Therefore it is necessary to regularize
the cross section that has a singularity in the $p_{t}\rightarrow 0$ limit. We
follow here the known prescription adopted in Pythia, where a special
suppression factor is introduced at the cross section level. The form factor
depends on a free parameter $p_{t0}$, which can be fixed using experimental
data for the $D$ meson production in $p+p$ and $p+^{4}He$ collisions at
$\sqrt{s}=38.7$ GeV and 86 GeV, respectively (see e.g. Ref. [10]). In
numerical calculations below we use $p_{t0}=2$ GeV.

The predictions for the charm production strongly depend on modelling of the
partonic content of the proton [9]. In particular, the contribution of the
charm \- initiated process is directly associated with the description of the
extrinsic and intrinsic components (for a recent review see, e.g. Ref. [11]).
Differently from the extrinsic charm quarks/antiquarks that are generated
perturbatively by gluon splitting, the intrinsic one have multiple connections
to the valence quarks of the proton and thus is sensitive to its
nonperturbative structure. The presence of an intrinsic component implies a
large enhancement of the charm distribution at large $x$ ($>0.1$) in
comparison to the extrinsic charm prediction. In recent years, the presence of
an intrinsic charm component have been included in the initial conditions of
the global parton analysis [12], resulting in IC distributions that are
compatible with the world experimental data. However, its existence is still a
subject of intense debate, mainly associated with the amount of intrinsic
charm in the proton wave function, which is directly related to the magnitude
of the probability to find an intrinsic charm or anticharm ($P_{ic}$) in the
nucleon.

In our analysis we will consider the collinear PDFs given by the CT14nnloIC
parametrization [12] from a global analysis assuming that the $x$-dependence
of the intrinsic charm component is described by the BHPS model [13]. Another
important ingredient is the modelling of
${\cal{F}}_{g^{*}}(x_{2},k_{t}^{2},\mu^{2})$. In our analysis here we will use
the uPDF derived using the Kimber-Martin-Ryskin (KMR) prescription [14].

##  3 Numerical results

|
---|---
(a) | (b)
Figure 1: Predictions of the hybrid model for (a) the Feynman $x_{F}$ \-
distributions for charm particles and (b) the prompt neutrino flux (rescaled
by $E_{\nu}^{3}$).

In Fig. 1 (a), we present our predictions for the Feynman $x_{F}$ distribution
of charm particles produced in $pp$ collisions at the atmosphere, considering
an incident proton with an energy of $E_{p}=10^{8}$ GeV and the KMR model for
the uPDF. We present separately the contribution associated with the
$cg\rightarrow cg$ mechanism and the sum of the two mechanisms, denoted by
“cg” and “gg + cg”, respectively. Moreover, we compare the IC predictions,
obtained using the CT14nnloIC parametrization for $P_{ic}=1\%$, with those
obtained disregarding the presence of the intrinsic component (denoted No IC
hereafter). One has that for small $x_{F}(\equiv x_{1}-x_{2})$, the charm
production is dominated by the $gg\rightarrow c\bar{c}$ mechanism, which is
expected since for $x_{F}\approx 0$ and high energies both longitudinal
momentum fractions $x_{i}$ are very small and the proton structure is
dominated by gluons. For the No IC case, the contribution of the
$cg\rightarrow cg$ mechanism is smaller than the gluon fusion one for all
values of $x_{F}$. In contrast, when intrinsic charm is included, the behavior
of the distribution in the intermediate $x_{F}$ range ($0.06\leq x_{F}\leq
0.6$) is strongly modified. Such a behaviour is expected, since for this
kinematical range, the charm production depends on the description of the
partonic content of the incident proton at large values of the Bjorken $x$
variable. The impact on the predictions for the prompt neutrino flux is
presented in Fig. 1 (b). As expected from the analysis performed in Ref. [3],
where we find that the dominant contribution to the neutrino flux comes
typically from $x_{F}$ in the region $0.2<x_{F}<0.5$, one has that the flux is
enhanced by one order of magnitude when intrinsic charm is included. In
agreement with the results presented in Fig. 1 (a), the contribution of the
$cg\rightarrow cg$ mechanism is negligible for the No IC case. However, it
becomes dominant in the IC case, with the normalization of the prompt flux
dependent on the amount of IC.

Figure 2: Comparison between our predictions and the experimental IceCube data
[1] for the atmospheric $\nu_{\mu}$ flux for the KMR uPDFs.

In Fig. 2 we present our results for the atmospheric $\nu_{\mu}$ flux, scaled
by a factor $E_{\nu}^{2}$, which is the sum of the conventional and prompt
contributions. The predictions were obtained considering different values for
$P_{ic}$ in the calculation of the prompt contribution. Moreover, for the
conventional atmospheric neutrino flux we assume the result derived in Ref.
[5]. The resulting predictions are compared with the IceCube data obtained in
Ref. [1] for the zenith-averaged flux of atmospheric neutrinos. One has that
the prompt contribution enhances the flux at large neutrino energies, with the
enhancement being strongly dependent on the magnitude of the $cg\rightarrow
cg$ mechanism. If this mechanism is disregarded, the results represented by
“Conv. + gg” in the figures indicate that the impact of the prompt flux is
small in the current kinematical range probed by IceCube. On the other hand,
the inclusion of the $cg\rightarrow cg$ mechanism implies a large enhancement
of the prompt flux at large $E_{\nu}$, with the associated magnitude being
strongly dependent on the value of $P_{ic}$. Our results for the KMR uPDF,
presented in Fig. 2, indicate that a value of $P_{ic}$ larger than $1.5\%$
implies a prediction for neutrino flux that overestimate the IceCube data at
high energies. This result sets the upper limit for the intrinsic charm
ammount in the nucleon to $P_{ic}\leq 1.5\%$. Surely, future data can be more
restrictive in the acceptable range of values for $P_{ic}$. It was also shown
in our original paper [15] that presence of saturation effects leads to a
slightly larger values for $P_{ic}$ that are not discarded by the current
IceCube data.

##  4 Conclusion

We have investigated the impact of the intrinsic charm component in the hadron
wave function, which carries a large fraction of the hadron momentum, on the
prompt neutrino flux. Our results has indicated that the inclusion of the
$cg\rightarrow cg$ mechanism has a strong effect on the prompt neutrino flux.
In particular, when the IC component is present, such a mechanism determines
the energy dependence of the flux at high energies, with the normalization
dependent on the value assumed for the probability to find the IC in the
proton wave function.

## Acknowledgements

This study was supported by the Polish National Science Center grant
UMO-2018/31/B/ST2/03537 and by the Center for Innovation and Transfer of
Natural Sciences and Engineering Knowledge in Rzeszów and by the Brazilian
funding agencies CNPq, FAPERGS and INCT-FNA (process number 464898/2014-5).

## References

  * [1] M. G. Aartsen et al. [IceCube], Development of a General Analysis and Unfolding Scheme and its Application to Measure the Energy Spectrum of Atmospheric Neutrinos with IceCube, Eur. Phys. J. C 75, no.3, 116 (2015), 10.1140/epjc/s10052-015-3330-z.
  * [2] R. Aaij et al. [LHCb], Measurements of prompt charm production cross-sections in $pp$ collisions at $\sqrt{s}=13$ TeV, JHEP 03, 159 (2016) [erratum: JHEP 09, 013 (2016); erratum: JHEP 05, 074 (2017)], 10.1007/JHEP03(2016)159.
  * [3] V. P. Goncalves, R. Maciuła, R. Pasechnik and A. Szczurek, Mapping the dominant regions of the phase space associated with $c\bar{c}$ production relevant for the prompt atmospheric neutrino flux, Phys. Rev. D 96, no.9, 094026 (2017), 10.1103/PhysRevD.96.094026.
  * [4] M. Ahlers, K. Helbing and C. Pérez de los Heros, Probing Particle Physics with IceCube, Eur. Phys. J. C 78, no.11, 924 (2018), 10.1140/epjc/s10052-018-6369-9.
  * [5] M. Honda, T. Kajita, K. Kasahara, S. Midorikawa and T. Sanuki, Calculation of atmospheric neutrino flux using the interaction model calibrated with atmospheric muon data Phys. Rev. D 75, 043006 (2007), 10.1103/PhysRevD.75.043006.
  * [6] P. Gondolo, G. Ingelman and M. Thunman, Charm production and high-energy atmospheric muon and neutrino fluxes Astropart. Phys. 5, 309-332 (1996), 10.1016/0927-6505(96)00033-3.
  * [7] V. P. Goncalves, R. Maciuła and A. Szczurek, From $D_{s}^{\pm}$ production asymmetry at the LHC to prompt $\nu_{\tau}$ at IceCube, Phys. Lett. B 794, 29-35 (2019), 10.1016/j.physletb.2019.05.026.
  * [8] A. Bhattacharya, et al., Prompt atmospheric neutrino fluxes: perturbative QCD models and nuclear effects, JHEP 11, 167 (2016), 10.1007/JHEP11(2016)167.
  * [9] R. Maciuła and A. Szczurek, Intrinsic charm in the nucleon and charm production at large rapidities in collinear, hybrid and kT-factorization approaches JHEP 10, 135 (2020), 10.1007/JHEP10(2020)135.
  * [10] R. Maciuła and A. Szczurek, Impact of the LHCb $p\\!+\\!^{4}\\!H\\!e$ fixed-target $D^{0}/\bar{D^{0}}$ data on the intrinisic $c\bar{c}$ component in the nucleon, https://arxiv.org/abs/2105.09370.
  * [11] S. J. Brodsky, G. I. Lykasov, A. V. Lipatov and J. Smiesko, Novel Heavy-Quark Physics Phenomena, Prog. Part. Nucl. Phys. 114, 103802 (2020), 10.1016/j.ppnp.2020.103802.
  * [12] T. J. Hou, et al., CT14 Intrinsic Charm Parton Distribution Functions from CTEQ-TEA Global Analysis JHEP 02, 059 (2018), 10.1007/JHEP02(2018)059.
  * [13] S. J. Brodsky, P. Hoyer, C. Peterson and N. Sakai, The Intrinsic Charm of the Proton, Phys. Lett. B 93, 451-455 (1980), 10.1016/0370-2693(80)90364-0.
  * [14] G. Watt, A. D. Martin and M. G. Ryskin, Unintegrated parton distributions and inclusive jet production at HERA, Eur. Phys. J. C 31, 73-89 (2003), 10.1140/epjc/s2003-01320-4.
  * [15] V. P. Goncalves, R. Maciula and A. Szczurek, Impact of intrinsic charm amount in the nucleon and saturation effects on the prompt atmospheric $\nu_{\mu}$ flux for IceCube, https://arxiv.org/abs/2103.05503.