Also at ]Universidade Estadual Paulista (UNESP), Faculdade de Engenharia e
Ciências, Departamento de Física e Química, Guaratinguetá, SP, 12516-410,
Brazil.

# Testing gravity with gauge-invariant polarization states of gravitational
waves:
Theory and pulsar timing sensitivity

Márcio E.S. Alves  [<EMAIL_ADDRESS>Universidade Estadual Paulista
(UNESP), Instituto de Ciência e Tecnologia,
São José dos Campos, SP, 12247-004, Brazil

###### Abstract

The determination of the polarization modes of gravitational waves (GWs) and
their dispersion relations is a crucial task for scrutinizing the viability of
extended theories of gravity. A tool to investigate the polarization states of
GWs is the well-known formalism developed by Eardley, Lee, and Lightman (ELL)
[Phys. Rev. D 8, 3308 (1973)] which uses the Newman-Penrose (NP) coefficients
to determine the polarization content of GWs in metric theories of gravity.
However, if the speed of GWs is smaller than the speed of light, the number of
NP coefficients is greater than the number of polarizations. To overcome this
inconvenience we use the Bardeen formalism to describe the six possible
polarization modes of GWs considering general dispersion relations for the
modes. The definition of a new gauge-invariant quantity enables an unambiguous
description of the scalar longitudinal polarization mode. We apply the
formalism to General Relativity, scalar-tensor theories, $f(R)$-gravity, and a
wide class of quadratic gravity. To obtain a bridge between theory and
experiment, we derive an explicit relation between a physical observable (the
derivative of the frequency shift of an electromagnetic signal), and the
gauge-invariant variables. From this relation, we find an analytical formula
for the pulsar timing rms response to each polarization mode. To estimate the
sensitivity of a single pulsar timing we focus on the case of a dispersion
relation of a massive particle. The sensitivity curves of the scalar
longitudinal and vector polarization modes change significantly depending on
the value of the effective mass. The detection (or absence of detection) of
the polarization modes using the pulsar timing technique has decisive
implications for alternative theories of gravity. Finally, investigating a
cutoff frequency in the pulsar timing band can lead to a more stringent bound
on the graviton mass than that presented by ground-based interferometers.

††preprint: APS/123-QED

##  I Introduction

The gravitational wave (GW) events detected so far by the Advanced LIGO and
Advanced Virgo interferometers have shown their ability to impact our
knowledge of physics and astrophysics. These observations offer a unique
opportunity to test General Relativity (GR) in the dynamical regime.

All the extensions to Einstein’s theory predict modifications to the
conventional GW signal due to one or more of three aspects, namely, changes in
the waveform due to particularities in the generation mechanism, changes in
the propagation due to new dispersion relations or differences in the
interaction of the wave with the background geometry, and the number of
independent polarization states of GWs.

Considering these effects, tests of gravity performed with the data of the
three observing runs of Advanced LIGO and Advanced Virgo have shown that all
the observed events are consistent with GR [1, 2, 3]. However, the planned
increase in the sensitivity of the detectors, the new generations of
interferometers, the pulsar timing technique and the future space-based GW
detectors as LISA will be able to produce stringent tests to GR.

In the case of the polarization states of GWs, a detection indicating the
presence of a polarization mode beyond the usual plus and cross polarizations
would imply a violation of Einstein’s theory. In general, an alternative
theory of gravity in four dimensions can predict up to six polarization modes
of GWs, namely, two tensor, two vector, one scalar transversal, and one scalar
longitudinal [4, 5].

To check the presence or absence of such modes in a specific theory is
appropriate to consider the evaluation of gauge-invariant quantities to
warrant that they are related to truly physical observables. The most common
strategy is to use the formalism developed by Eardley, Lee, and Lightman (ELL)
fifty years ago [4, 5]. Within the ELL framework, the irreducible parts of the
linearized Riemann tensor are related to the polarization states of GWs. These
components can be written on a tetrad basis to obtain the so-called Newman-
Penrose (NP) coefficients.

The NP formalism [6, 7] is a powerful framework to describe the spacetime
properties of a Lorentzian geometry in a coordinate-independent manner. It has
many applications not only for GR but also in the scope of extended theories
of gravity including the study of perturbed spacetimes in such theories (see,
e.g., [8, 9, 10]).

In the ELL framework, the specific case of perturbations around a flat
background is considered, resulting in two real and two complex NP
coefficients describing the polarization states of GWs in any four-dimensional
metric theory of gravity. The ELL formalism has been applied in the scope of
several theories to reveal the polarization properties of GWs (see, e.g., [11,
12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]).

In recent years, however, some criticisms have been raised in the literature
regarding the use of the original ELL formalism in theories that present one
or more massive modes in the linearized regime [27, 28, 29, 30]. This is the
case of a huge class of alternative theories, including the massive version of
the Brans-Dicke theory [13], Horndeski theory [29] and $f(R)$-gravity [12, 17,
27]. The main criticism is related to the fact that if GWs travel at a speed
different from the speed of light, then other NP coefficients, beyond the
original four, would be non-null. Consequently, the original ELL formalism
would be incomplete and could result in misleading conclusions for some
theories. Thus, in this case, to have a complete description of the
polarizations one needs new NP variables. In fact, to describe six
polarizations there are not only four but nine NP coefficients representing
fourteen components of the Riemann tensor [30]. Certainly, some coefficients
could have a greater amplitude when compared to others depending on the
dispersion relation and frequency.

However, there is a gauge-invariant alternate formalism to identify the
polarization modes in a metric theory of gravity. It consists of decomposing
the metric into irreducible components according to their properties under
spatial rotation [(3+1) decomposition] and then constructing gauge-invariant
combinations of the metric perturbations. In cosmological perturbation theory,
such quantities are known as Bardeen variables [31, 32]. Flanagan and Hughes
[33] have used these variables in GR to describe perturbations over a
Minkowski background. Recently, Wagle et al. [22] used these gauge-invariant
variables to study the polarizations of GWs in the context of two theories,
namely, the dynamical Chern-Simons gravity and the Einstein-dilaton-Gauss-
Bonnet gravity. They have found that the ELL formalism and the (3+1)
decomposition lead to the same result in both cases. In the case of $f(R)$
gravity, Moretti et al. [34] have applied the gauge-invariant formalism to
investigate the propagating degrees of freedom within this class of theories.
Using the same formalism, Dong et al. [35] have shown that the number of
polarization modes of GWs depends on the choice of parameters in the
generalized Proca theory. Other works in the literature have used this
framework to describe GWs in alternative theories of gravity (see, e.g., [36,
37, 38]). The formalism of Bardeen variables has the advantage that the same
number of variables are applicable to describe the polarization modes of GWs
to any frequency.

The aim of the present article is twofold. First, we review the Bardeen
formalism and show how it applies in identifying the polarization modes of GWs
for any metric theory of gravity. We discuss the advantages of this formalism
when compared to the EEL formalism. The second aim is to estimate the pulsar
timing sensitivities considering the gauge-invariant variables to evaluate the
response to the polarization modes.

Analyzing the pulsar timing sensitivity, Alves and Tinto have shown that this
technique is significantly more sensitive to non-transverse polarizations than
the usual plus and cross polarizations [39]. To a lesser extent, the same
effect appears in the LISA frequency band [40]. Other authors have noticed the
enhancement in the sensitivity to non-transverse polarization modes [41, 42].
The effect appears if the wavelength of GWs is of the same order or smaller
than the scale size of the detector [39, 40]. Therefore, one does not expect
this effect to appear in detectors operating within the long-wavelength limit
as in ground-based interferometers.

In deriving the pulsar timing sensitivity, previous works consider a specific
gauge choice (one of the most common is the synchronous gauge). Although this
is a usual procedure, the resulting sensitivity curve may not be appropriate
to evaluate the detectability of the polarization modes in all alternative
theories of gravity. This is because a given gauge is not necessarily
applicable in all theories and residual gauge freedom may be present,
resulting in misleading interpretations. Furthermore, the most important
aspect is that since gauge-invariant quantities can be related to truly
physical observables, detector sensitivity is expected to be derived in terms
of them.

In this sense, here we derive in a complete gauge-invariant fashion the pulsar
timing sensitivities to each polarization mode considering the dispersion
relation of a massive particle. We show how the sensitivity can be written in
terms of the Bardeen variables used to describe the polarization modes of GWs
in a wide range of metric theories of gravity.

The article is organized as follows. In Section (II) we give a short overview
of the ELL formalism. In Section (III) we describe the Bardeen formalism and
show how it can be applied in describing the polarization modes of GWs for any
theory of gravity. We apply the formalism to the case of GR, scalar-tensor
theories of gravity, $f(R)$-gravity, and a wide class of quadratic gravity.
The relation between the one-way response to GWs and the gauge-invariant
variables, as well as the pulsar timing sensitivity, are obtained in Section
(IV). Finally, we conclude the article with Section (V). Throughout the
article, we use the metric signature $(-,+,+,+)$ and units such that
$c=\hbar=1$ unless otherwise mentioned.

##  II An overview of the ELL formalism

In the original ELL formalism for the determination of polarization modes of
GWs, Eardley et al. [4, 5] considered GWs propagating in the $+z$ direction at
the speed of light and defined a null complex tetrad
$(\mathbf{k,l,m,\bar{m}})$. This tetrad is related to the Cartesian tetrad
$(\mathbf{e}_{t},\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z})$ by

$\mathbf{k}=\frac{1}{\sqrt{2}}(\mathbf{e}_{t}+\mathbf{e}_{z}),$ (1)
$\mathbf{l}=\frac{1}{\sqrt{2}}(\mathbf{e}_{t}-\mathbf{e}_{z}),$ (2)
$\mathbf{m}=\frac{1}{\sqrt{2}}(\mathbf{e}_{x}+i\mathbf{e}_{y}),$ (3)
$\bar{\mathbf{m}}=\frac{1}{\sqrt{2}}(\mathbf{e}_{x}-i\mathbf{e}_{y}).$ (4)

It is easy to verify that the tetrad vectors obey the relations:

$-\mathbf{k}\cdot\mathbf{l}=\mathbf{m}\cdot\bar{\mathbf{m}}=1,$ (5)
$\mathbf{k}\cdot\mathbf{m}=\mathbf{k}\cdot\bar{\mathbf{m}}=\mathbf{l}\cdot\mathbf{m}=\mathbf{l}\cdot\bar{\mathbf{m}}=0.$
(6)

To denote components of tensors with respect to the null tetrad basis we use
Roman subscripts, that is, $P_{abc\dots}\equiv
P_{\alpha\beta\gamma\dots}a^{\alpha}b^{\beta}c^{\gamma}\dots$, where
$(a,b,c,\dots)$ run over $(\mathbf{k,l,m,\bar{m}})$ and
$(\alpha,\beta,\gamma,\dots)$ run over $(t,x,y,z)$.

The Riemann curvature tensor $R_{\alpha\beta\gamma\delta}$ can be split into
the irreducible parts: the Weyl tensor, the traceless Ricci tensor and the
curvature scalar, whose tetrad components can be named, respectively as
$\Psi_{A}$, $\Phi_{AB}$ and $\Lambda$ following the notation of Newman and
Penrose [6, 7]. In general, in a four-dimensional space, we have ten $\Psi$’s,
nine $\Phi$’s, and one $\Lambda$ which are all algebraically independent.
However, when we restrict ourselves to null plane waves, we find that the
differential and algebraic properties of $R_{\alpha\beta\gamma\delta}$ reduce
the number of independent components to six. In the above tetrad, we can
choose the following quantities to represent these components [4, 5]

$\Psi_{2}=-\frac{1}{6}R_{lklk},$ (7) $\Psi_{3}=-\frac{1}{2}R_{lkl\bar{k}},$
(8) $\Psi_{4}=-R_{l\bar{m}l\bar{m}},$ (9) $\Phi_{22}=-R_{lml\bar{m}}.$ (10)

Notice that since $\Psi_{3}$ and $\Psi_{4}$ are complex, each represents two
independent polarizations. For these six components, three are transverse to
the direction of propagation, with two representing quadrupolar deformations
[${\rm Re}(\Psi_{4})$ and ${\rm Im}(\Psi_{4})$] and one monopolar deformation
($\Phi_{22}$). Three modes are longitudinal, with one an axially symmetric
stretching mode in the propagation direction ($\Psi_{2}$), and one quadrupolar
mode in each one of the two orthogonal planes containing the propagation
direction [${\rm Re}(\Psi_{3})$ and ${\rm Im}(\Psi_{3})$].

The above formalism is still accurate if the speed of GWs is close to the
speed of light. In fact, corrections to the ELL formalism are of the order
$O(\epsilon\delta\mathcal{R})$ to the case of a nearly null GW [43], where
$\epsilon=(c/v_{\rm gw})^{2}-1$, $v_{\rm gw}$ is the speed of GWs and
$\delta\mathcal{R}$ is some component of the perturbed Riemann tensor.
Therefore, considering the current upper bound for the graviton mass from the
observations of binary black hole mergers ($m_{g}\leq 1.27\times
10^{-23}~{}{\rm eV}/c^{2}$) [3], we find $\epsilon\sim 10^{-21}$ for the
frequency 0.1 kHz (considering the dispersion relation of a massive graviton).
For the frequency of 1 mHz, we obtain $\epsilon\sim 10^{-11}$. Thus, in these
cases $O(\epsilon\mathcal{R})$ is several orders of magnitude smaller than the
NP quantities for null waves [which are of the order $O(\mathcal{R})$ ]. We
conclude that such corrections are undetectable in the frequency band of
ground-based and space-based interferometers.

On the other hand, they can become important for lower frequencies. For
instance, in the band of pulsar timing arrays, we find $\epsilon\sim 1$ for
the above mass of the graviton, and by considering frequencies of the order of
nanohertz. In this case, to have a complete description of the polarizations
the ELL formalism needs to be amended by including other NP coefficients.
Considering plane waves propagating at a speed $v_{\rm gw}$, Hyun et al. [30]
have expressed the polarizations of GWs in terms of nine NP scalars, namely,
$\Psi_{0}$, $\Psi_{1}$, $\Psi_{2}$, $\Psi_{3}$, $\Psi_{4}$, $\Phi_{00}$,
$\Phi_{02}$, $\Phi_{22}$, $\Lambda$. Since $\Psi_{0}$, $\Psi_{1}$, $\Psi_{3}$,
$\Psi_{4}$, $\Phi_{02}$ are complex, the nine scalars represent fourteen
components of the Riemann curvature tensor needed to describe six polarization
modes of GWs. One can use the formalism of Bardeen variables described in the
next section to overcome this inconvenience.

##  III Describing the polarization states of GWs within the Bardeen framework

###  III.1 Helicity decomposition and gauge invariant perturbations

In this section, we introduce the helicity decomposition of the metric
perturbation and define the gauge invariant variables that can be computed
from them.

Let us start by expanding the metric around flat space
$g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$ with $|h_{\mu\nu}|\ll 1$. The
perturbation $h_{\mu\nu}$ can be decomposed considering the behavior of its
components under spatial rotations as follows

$h_{00}=2\psi,$ (11) $h_{0i}=\beta_{i}+\partial_{i}\gamma,$ (12)
$h_{ij}=-2\phi\delta_{ij}+\left(\partial_{i}\partial_{j}-\frac{1}{3}\delta_{ij}\nabla^{2}\right)\lambda+\frac{1}{2}(\partial_{i}\epsilon_{j}+\partial_{j}\epsilon_{i})+h_{ij}^{\rm
TT},$ (13)

where the vector and tensor quantities are subject to the following
constraints

$\partial_{i}\beta^{i}=0,~{}~{}\partial_{i}\epsilon^{i}=0,$ (14)
$\partial^{j}h_{ij}^{\rm TT}=0,~{}~{}\delta^{ij}h_{ij}^{\rm TT}=0.$ (15)

Therefore, from the 10 degrees of freedom of $h_{\mu\nu}$ we have four scalar
degrees of freedom $\\{\psi,\phi,\gamma,\lambda\\}$, four degrees of freedom
in the two transverse vectors $\\{\beta_{i},\epsilon_{i}\\}$ and two degrees
of freedom in the transverse-traceless (TT) spatial tensor $h^{\rm TT}_{ij}$.

The gauge transformation of the metric perturbation

$h_{\mu\nu}\rightarrow
h_{\mu\nu}-(\partial_{\mu}\xi_{\nu}+\partial_{\nu}\xi_{\mu}),$ (16)

with $|\partial_{\mu}\xi_{\nu}|$ small preserves $|h_{\mu\nu}|\ll 1$, thus it
is a symmetry of the linearized theory in general. To understand how the
harmonic variables behave under this gauge transformation, notice that the
4-vector $\xi_{\mu}$ can also be decomposed as

$\xi_{0}=A,$ (17) $\xi_{i}=B_{i}+\partial_{i}C,$ (18)

with

$\partial_{i}B^{i}=0.$ (19)

Therefore, considering Eq. (16) and following the symmetry of the
transformations under spatial rotations, we find that the gauge
transformations of the scalar harmonic variables read

$\displaystyle\psi$ $\displaystyle\rightarrow\psi-\dot{A},$
$\displaystyle\phi$ $\displaystyle\rightarrow\phi+\frac{1}{3}\nabla^{2}C,$
$\displaystyle\lambda$ $\displaystyle\rightarrow\lambda-2C,$
$\displaystyle\gamma$ $\displaystyle\rightarrow\gamma-\dot{C}-A.$ (20)

The transformations of the vectors are

$\displaystyle\beta_{i}$ $\displaystyle\rightarrow\beta_{i}-\dot{B_{i}},$
$\displaystyle\epsilon_{i}$ $\displaystyle\rightarrow\epsilon_{i}-2B_{i},$
(21)

while $h_{ij}^{TT}$ is gauge-invariant

$\displaystyle h_{ij}^{TT}\rightarrow h_{ij}^{TT}.$ (22)

The Riemann tensor and, correspondingly, the Einstein tensor are gauge
invariant quantities. Therefore, one possible way of dealing with this gauge
freedom is to impose gauge conditions on the metric perturbations. This is the
usual way in GW physics. Several gauges are possible, some of the most common
gauge choices are the synchronous gauge and the Newtonian gauge. The latter
fixes the gauge uniquely, however, the conditions imposed by the former leave
a residual gauge freedom. This ambiguity implies the existence of unphysical
modes when the gravitational equations are solved. Particularly, this can lead
to an ambiguity in the determination of truly propagating GW modes in
alternative theories of gravity. In the Bardeen words “only gauge-invariant
quantities have any inherent physical meaning” [31].

In this sense, Bardeen [31] has constructed gauge-invariant quantities from
combinations of the above scalar and vector variables to deal with
perturbations in a FLRW background spacetime. In this article, we consider
solely the Minkowski metric for the background. From the transformations
(III.1) we see that we can obtain the following two gauge-invariant scalar
combinations

$\displaystyle\Phi$ $\displaystyle=-\phi-\frac{1}{6}\nabla^{2}\lambda,$ (23)
$\displaystyle\Psi$
$\displaystyle=-\psi+\dot{\gamma}-\frac{1}{2}\ddot{\lambda}.$ (24)

In the same way, from the transformations (III.1) we obtain one gauge-
invariant transverse spatial vector

$\Xi_{i}=\beta_{i}-\frac{1}{2}\dot{\epsilon}_{i},~{}~{}\partial_{i}\Xi^{i}=0.$
(25)

Thus, we have six gauge-invariant degrees of freedom: two scalars, $\Psi$ and
$\Phi$, two degrees of freedom in the spatial vector $\Xi_{i}$ and two degrees
of freedom in the transverse-traceless spatial tensor $h_{ij}^{TT}$. These
gauge-invariant variables are the flat background version of the well-known
Bardeen variables.

We can now write the electric components of the first order perturbed Riemann
tensor using the (3+1) decomposition of the metric perturbations

$\delta
R_{i0j0}=\partial_{i}\partial_{j}\Psi-\ddot{\Phi}\delta_{ij}+\frac{1}{2}(\partial_{i}\dot{\Xi}_{j}+\partial_{j}\dot{\Xi}_{i})-\frac{1}{2}\ddot{h}_{ij}^{\rm
TT}.$ (26)

As expected, the perturbed Riemann tensor depends only on the gauge-invariant
variables $\Phi$, $\Psi$, $\Xi_{i}$, and $h_{ij}^{\rm TT}$. It will also be
useful to know the components of the perturbed Ricci tensor

$\displaystyle\delta R_{00}$ $\displaystyle=\nabla^{2}\Psi-3\ddot{\Phi},$ (27)
$\displaystyle\delta R_{0i}$
$\displaystyle=-2\partial_{i}\dot{\Phi}-\frac{1}{2}\nabla^{2}\Xi_{i},$ (28)
$\displaystyle\delta R_{ij}$
$\displaystyle=-\partial_{i}\partial_{j}(\Phi+\Psi)-\delta_{ij}(-\ddot{\Phi}+\nabla^{2}\Phi)$
$\displaystyle-\frac{1}{2}(\partial_{i}\dot{\Xi}_{j}+\partial_{j}\dot{\Xi}_{i})-\frac{1}{2}\Box
h_{ij}^{\rm TT},$ (29)

and the perturbed curvature scalar

$\delta R=-2\nabla^{2}\Psi+6\ddot{\Phi}-4\nabla^{2}\Phi.$ (30)

From the above expressions we find the components of the perturbed Einstein
tensor

$\delta G_{00}=-2\nabla^{2}\Phi,~{}~{}~{}\delta G_{0i}=\delta R_{0i},$ (31)
$\delta G_{ij}=\delta R_{ij}-\frac{1}{2}\delta_{ij}\delta R.$ (32)

If $\dot{\Psi}\neq 0$ it will also prove useful to define the following gauge-
invariant variable

$\Theta\equiv\eta_{\Psi}^{2}\Psi-\Phi,$ (33)

where

$\eta_{\Psi}\equiv\left|\frac{\nabla\Psi}{\dot{\Psi}}\right|.$ (34)

The physical meaning of $\Theta$ will be clarified in the next section.

###  III.2 Description of the polarization states of GWs with Bardeen
variables

The six degrees of freedom encoded in the four gauge-invariant variables
defined above can be radiative or non-radiative depending on the underlying
theory of gravity. It is well known that the transverse-traceless tensor
$h_{ij}^{\rm TT}$ is the only radiative quantity in the GR theory. Moreover,
for those theories of gravity predicting spin-0 polarization modes, it is
expected that the scalars $\Phi$ and $\Psi$ are related through the field
equations. This issue will be analyzed in the next section through examples in
the context of some theories of gravity. In the present section, we describe
in a general way the six polarization states of GWs employing the Bardeen
variables. To this end, we suppose that all gauge-invariant quantities are
radiative and independent.

Since the variables are radiative, they are functions of the retarded time

$u\equiv t-\frac{\vec{k}\cdot\vec{r}}{\omega},$ (35)

where $\vec{k}$ is the GW wave vector and $\omega$ is the angular frequency.
Consider, for instance, the scalar $\Psi=\Psi(u)$. In terms of the coordinates
$t$ and $\vec{r}$ it obeys

$\nabla\Psi=-\frac{\vec{k}_{\Psi}}{\omega}\dot{\Psi},$ (36)

where $\vec{k}_{\Psi}$ is the wave vector of the $\Psi$ variable.

Using this equation in the definition (34) we identify $\eta_{\Psi}$ as the
dispersion relation

$\eta_{\Psi}(\omega)=\frac{k_{\Psi}(\omega)}{\omega},$ (37)

with $k_{\Psi}=|\vec{k}_{\Psi}|$.

Since each variable can have a different dispersion relation, we express them
as functions of four different retarded times
$u_{A}=t-\eta_{A}\hat{k}\cdot\vec{r}$, and the four dispersion relations are
expressed by the functions $\eta_{A}(\omega)=k_{A}(\omega)/\omega$ with
$A=\Phi$, $\Psi$, $V$ and $T$.

Hence, using the definition of the variable $\Theta$ (33) the electric
components of the perturbed Riemann tensor (26) now read

$\displaystyle\delta R_{0i0j}$
$\displaystyle=\hat{k}_{i}\hat{k}_{j}\Theta^{\prime\prime}-(\delta_{ij}-\hat{k}_{i}\hat{k}_{j})\Phi^{\prime\prime}$
$\displaystyle-\frac{1}{2}\eta_{V}(\hat{k}_{i}\Xi_{j}^{\prime\prime}+\hat{k}_{j}\Xi_{i}^{\prime\prime})-\frac{1}{2}h_{ij}^{{\rm
TT}{\prime\prime}},$ (38)

where primes denote derivatives with respect to the retarded times and
$\hat{k}_{i}$ are components of the unit wave vector. The transverse
conditions can be written as $\hat{k}^{i}\Xi_{i}=0$ and $\hat{k}^{i}h^{\rm
TT}_{ij}=0$. The two tensor polarization states are described by $h_{ij}^{\rm
TT}$ as usual. The two spin-1 polarization states are described by the vector
$\Xi_{i}$. Although $\Xi_{i}$ is transverse to the direction of propagation of
the GW, notice that it enters $\delta R_{0i0j}$ with a term that is
proportional to $\hat{k}_{i}$, which comes from the spatial derivative of
$\Xi_{i}$. Therefore, we arrive at the known result that vector polarization
affects the curvature in the transverse and longitudinal directions.

In the second term in the right-hand-side of Eq. (III.2), we recognize the
quantity multiplying the scalar Bardeen variable $\Phi$ as the projection
operator $P_{ij}=(\delta_{ij}-\hat{k}_{i}\hat{k}_{j})$ which has the property
$\hat{k}^{i}P_{ij}=0$. This operator projects any spatial vector on the
subspace orthogonal to the direction of propagation of the GW. Thus, this term
represents the scalar-transverse polarization mode. Finally, the term
$\hat{k}_{i}\hat{k}_{j}\Theta^{\prime\prime}$ expresses the overall
longitudinal effect from the scalar variables and then $\Theta$ is the degree
of freedom responsible for describing the scalar-longitudinal polarization
mode.

To summarize, the six possible polarization modes of GWs can be described by
the six degrees of freedom present in the gauge-invariant variables
$h_{ij}^{\rm TT}$, $\Xi_{i}$, $\Phi$, and $\Theta$. Since this result is valid
for any dispersion relation, it turns out that the Bardeen formalism is much
simpler than the ELL formalism in the determination of the polarization modes.

###  III.3 Lorentz transformations of the gauge-invariant variables

Although the Bardeen variables are gauge-invariant quantities, generally they
are not Lorentz-invariant. This means that distinct observers connected by
boosts can measure different polarization contents of a given theory of
gravity.

A finite Lorentz transformation can be built up from a sequence of
infinitesimal Lorentz transformations. Therefore, if a quantity is invariant
under infinitesimal boosts, it is also invariant under finite boosts. Thus, to
evaluate the behavior of the gauge-invariant variables under boosts let us
consider an infinitesimal Lorentz transformation $x^{\mu}\rightarrow
x^{\prime\mu}={\Lambda^{\mu}}_{\nu}x^{\nu}$ with
${\Lambda^{\mu}}_{\nu}={\delta^{\mu}}_{\nu}+{\omega^{\mu}}_{\nu}$. The metric
perturbation changes according to

$h_{\mu\nu}^{\prime}=h_{\mu\nu}+\delta h_{\mu\nu},$ (39)

where

$\delta
h_{\mu\nu}={\omega_{\mu}}^{\rho}h_{\rho\nu}+{\omega_{\nu}}^{\rho}h_{\mu\rho}.$
(40)

The non-null components of the quantity $\delta h_{\mu\nu}$ are small compared
with $h_{\mu\nu}$, however when considering a sequence of infinitesimal
Lorentz transformations the transformed $h^{\prime}_{\mu\nu}$ can be
significantly different from $h_{\mu\nu}$.

Using this transformation and remembering the decomposition of the metric
perturbation in terms of the harmonic variables given by Eqs. (11), (12) and
(13) one can find the change of each variable due to infinitesimal Lorentz
transformations. If one restricts to boosts $\omega_{ij}=0$ and we have the
following expressions for the change in the gauge-invariant variables
[44]111Notice that our definition of the variable $\Psi$ differs from the
definition given by the Ref. [44] by a minus sign, i.e., one should change
$\Psi\rightarrow-\Psi$ to compare the equations presented in both articles.

$\delta\Phi=\frac{1}{2}{\omega_{0}}^{i}\Xi_{i},$ (41)
$\delta\Psi=-2{\omega_{0}}^{i}\nabla^{-2}\partial_{i}(\dot{\Phi}-\dot{\Psi})+{\omega_{0}}^{i}\Xi_{i}-\frac{3}{2}{\omega_{0}}^{i}\nabla^{-2}\ddot{\Xi_{i}},$
(42) $\displaystyle\delta\Xi_{i}$
$\displaystyle={\omega_{0}}^{j}\nabla^{-2}{\Big{[}}\Box
h_{ij}^{TT}-(\partial_{i}\dot{\Xi}_{j}+\partial_{j}\dot{\Xi}_{i})$
$\displaystyle-2(\partial_{i}\partial_{j}-\delta_{ij}\nabla^{2})(\Phi-\Psi){\Big{]}},$
(43) $\displaystyle\delta h_{ij}^{TT}$
$\displaystyle=\omega_{0i}\Xi_{j}+\omega_{0j}\Xi_{i}-\delta_{ij}{\omega_{0}}^{k}\Xi_{k}$
$\displaystyle+{\omega_{0}}^{k}\nabla^{-2}{\Big{[}}\partial_{i}\partial_{j}\Xi_{k}-\partial_{i}\partial_{k}\Xi_{j}-\partial_{j}\partial_{k}\Xi_{i}$
$\displaystyle-(\partial_{i}\dot{h}_{jk}^{TT}+\partial_{j}\dot{h}_{ik}^{TT}){\Big{]}}.$
(44)

Therefore, a consequence of the decomposition scheme is that the Bardeen
variables transform among themselves under boosts.

As we will show in the next section, analyzing the linearized vacuum field
equations of a theory of gravity, one can find the governing equations for the
Bardeen variables and possible relations between the scalars $\Phi$ and
$\Psi$. Although the variables are not Lorentz invariant, in general, it is
interesting to notice some particular cases:

  1. 1.

$h_{ij}^{TT}\neq 0,~{}\Box h_{ij}^{TT}=0,~{}\Xi_{i}=0,~{}\Phi=\Psi$. The
variables $\Phi$, $\Psi$ and $\Xi_{i}$ are Lorentz invariant.

  2. 2.

$h_{ij}^{TT}\neq 0,~{}\Box h_{ij}^{TT}\neq 0,~{}\Xi_{i}=0,~{}\Phi=\Psi$. The
variables $\Phi$ and $\Psi$ are Lorentz invariant.

  3. 3.

$h_{ij}^{TT}\neq 0,~{}\Xi_{i}=0,~{}\Phi\neq\Psi$. The variable $\Phi$ is the
only Lorentz invariant.

As a consequence of cases 1 and 2, the variable $\Theta$ is also Lorentz
invariant, and thus all the Lorentz observers measure the same scalar-
longitudinal mode. The scalar-transversal polarization mode is Lorentz
invariant in the three cases. On the other hand, the vector polarization mode
is invariant only in the special case 1 for which it is null for all Lorentz
observers. If $\Xi_{i}\neq 0$ and/or $h_{ij}^{TT}\neq 0$, the vector and
tensor modes are not Lorentz invariant in general.

###  III.4 Polarization states of GWs in some theories of gravity

The procedure of determination of the polarization modes in an alternative
theory of gravity starts, as usual, with the linearization of the vacuum field
equations of the theory. The equations for metric perturbations should be
written in terms of the Bardeen variables. Finally, from these equations, it
will be possible to determine which variable represents a truly radiative
mode, the number of independent degrees of freedom, and each dispersion
relation related to the polarization modes. In this section, we illustrate
this procedure by evaluating the polarization modes of GWs for General
Relativity, scalar-tensor theories of gravity, and $f(R)$-gravity.

####  III.4.1 General Relativity

Let us consider the Einstein-Hilbert action in the absence of matter sources
with a vanishing cosmological constant

$I=\frac{1}{16\pi G}\int d^{4}x\sqrt{-g}R.$ (45)

The Einstein equations for vacuum is obtained if we vary this action with
respect to the metric $g_{\mu\nu}$. They are

$G_{\mu\nu}=0.$ (46)

Let us expand the metric about the Minkowski spacetime
$g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$ and write the Einstein equations to
first order in $h_{\mu\nu}$. If we use the gauge-invariant quantities
described previously we obtain the components of the Einstein tensor as given
by the Eqs. (31) and (32). From the 00 component, we obtain

$\nabla^{2}\Phi=0.$ (47)

With this result and using Eq. (30) in the trace equation
$g^{\mu\nu}G_{\mu\nu}=-R=0$ we have

$\ddot{\Phi}-\frac{1}{3}\nabla^{2}\Psi=0.$ (48)

From Eq. (47) we conclude that in the absence of matter

$\Phi=0,$ (49)

and with this result in the Eq. (48) we have $\nabla^{2}\Psi=0$ and then

$\Psi=0.$ (50)

Therefore, the two gauge-invariant scalars are non-propagating degrees of
freedom in Einstein theory and vanish in the absence of matter fields. If we
use this result in the $0i$ components of the Eq. (46) along with Eq. (31) we
obtain a similar result for the vector modes

$\nabla^{2}\Xi_{i}=0~{}\Rightarrow~{}\Xi_{i}=0.$ (51)

Applying the above results and Eq. (32) in the spatial components of the field
equations $G_{ij}=0$ we find the following equation for gauge-invariant tensor
perturbation

$\square h_{ij}^{TT}=0.$ (52)

Therefore, we see that only the tensor degrees of freedom are radiative since
they obey a wave equation. In the absence of matter $\Phi=\Psi=0$ and
$\Xi_{i}=0$, we have only two polarization states represented by
$h_{ij}^{TT}$, which are the usual $+$ and $\times$ polarizations. Since the
tensor modes propagate at the speed of light we have $\eta_{T}=1$ for GR. From
the Lorentz transformations of the Bardeen variables given by Eqs. (41), (42),
(III.3) and (III.3), we see that the absence of scalar and vector polarization
modes is a Lorentz invariant statement in the GR case.

####  III.4.2 Scalar-tensor theories of gravity

For simplicity, we restrict to scalar-tensor theories of gravity whose action
can be written in the following form in the absence of matter [46, 47]

$I=\frac{1}{16\pi}\int d^{4}x\sqrt{-g}\left[\varphi
R-\frac{\varpi(\varphi)}{\varphi}\nabla^{\mu}\varphi\nabla_{\mu}\varphi+V(\varphi)\right].$
(53)

However, as will be clear at the end of this subsection, the results we will
find are valid for more general scenarios encompassed by the Horndeski theory.

In the theory described by (53), gravity is mediated not only by the metric
but also by a scalar field $\varphi$, $\varpi(\varphi)$ is a coupling function
and $V(\varphi)$ is a generic scalar field potential. Varying this action with
respect to the metric and scalar field we obtain

$\displaystyle G_{\mu\nu}$
$\displaystyle=\frac{1}{2}\varphi^{-1}V(\varphi)g_{\mu\nu}+\varpi(\varphi)\varphi^{-2}{\Big{(}}\nabla_{\mu}\varphi\nabla_{\nu}\varphi$
$\displaystyle-\frac{1}{2}g_{\mu\nu}\nabla_{\alpha}\varphi\nabla^{\alpha}\varphi{\Big{)}}+\varphi^{-1}(\nabla_{\mu}\nabla_{\nu}\varphi-
g_{\mu\nu}\square\varphi),$ (54)

and

$\square\varphi+\frac{\varphi
V^{\prime}(\varphi)-2V(\varphi)}{3+2\varpi(\varphi)}=-\frac{\varpi^{\prime}(\varphi)\nabla_{\alpha}\varphi\nabla^{\alpha}\varphi}{3+2\varpi(\varphi)},$
(55)

where a prime denotes derivative with respect to $\varphi$.

In the weak-field limit, we can expand the metric about the Minkowski
background as in the GR case, while the scalar field is expanded as

$\varphi=\varphi_{0}+\delta\varphi,~{}~{}\delta\varphi\ll\varphi_{0},$ (56)

where $\varphi_{0}$ is the asymptotic value of $\varphi$ far from the system
that is generating GWs. Expanding the potential and the coupling function
about $\varphi_{0}$ up to the second order we have

$\displaystyle V(\varphi)$
$\displaystyle=V(\varphi_{0})+V^{\prime}(\varphi_{0})\delta\varphi+\frac{1}{2}V^{\prime\prime}(\varphi_{0})\delta\varphi^{2}+{\mathcal{O}}(\delta\varphi^{3}),$
(57) $\displaystyle\varpi(\varphi)$
$\displaystyle=\varpi(\varphi_{0})+\varpi^{\prime}(\varphi_{0})\delta\varphi+\frac{1}{2}\varpi^{\prime\prime}(\varphi_{0})\delta\varphi^{2}+{\mathcal{O}}(\delta\varphi^{3}).$
(58)

Far from the system we assume that the spacetime is asymptotically Minkowski.
Thus, imposing that the background Einstein tensor vanishes in the Eq.
(III.4.2), we are lead to $V(\varphi_{0})=V^{\prime}(\varphi_{0})=0$. In this
limit, from (55) we obtain the equation for the first-order perturbation of
the scalar field

$\left(\square-m^{2}\right)\delta\varphi=0,$ (59)

where now $\square$ is the D’Alembertian operator for the Minkowski background
metric and the mass of the scalar field is defined by

$m^{2}\equiv-\frac{\varphi_{0}V^{\prime\prime}(\varphi_{0})}{3+2\varpi_{0}},$
(60)

where $\varpi_{0}=\varpi(\varphi_{0})$. The first-order perturbation of Eqs.
(III.4.2) become

$\delta
G_{\mu\nu}+\left(\eta_{\mu\nu}\square-\partial_{\mu}\partial_{\nu}\right)\frac{\delta\varphi}{\varphi_{0}}=0,$
(61)

where $\delta G_{\mu\nu}$ is the linearized Einstein tensor. Replacing the Eq.
(31) in the 00 component of this equation we find

$\nabla^{2}\left(\Phi+\frac{1}{2}\frac{\delta\varphi}{\varphi_{0}}\right)=0.$
(62)

Hence, in the absence of matter, we conclude that the gauge-invariant scalar
$\Phi$ is proportional to the scalar field perturbation

$\Phi=-\frac{1}{2}\frac{\delta\varphi}{\varphi_{0}}.$ (63)

Moreover, using the trace of Eqs. (61), the Eq. (59) and the perturbed
expression of the Ricci scalar (30), it is easy to show that

$\Psi=\Phi=-\frac{1}{2}\frac{\delta\varphi}{\varphi_{0}}.$ (64)

The $0i$ components of Eqs. (61) together with Eq. (31) lead to the following
equation for the gauge-invariant vector perturbation

$\nabla^{2}\Xi_{i}=0,$ (65)

which, in the absence of matter, gives

$\Xi_{i}=0.$ (66)

Finally, using Eqs. (32), (64) and (66), the spatial components of Eqs. (61)
result

$\square h_{ij}^{TT}=0.$ (67)

Therefore, we conclude that in the scalar-tensor theory of gravity, there are
three radiative degrees of freedom. The two tensor degrees of freedom obey the
wave equation (67) in the same manner as in the Einstein theory. They
propagate at the speed of light and, therefore, $\eta_{T}=1$. On the other
hand, there is a scalar degree of freedom $\Phi=\Psi$ which obeys a Klein-
Gordon type equation

$\left(\square-m^{2}\right)\Phi=0.$ (68)

This equation has a solution $\Phi\propto{e}^{ik_{\alpha}x^{\alpha}}$ with the
wave 4-vector $k^{\alpha}\equiv(\omega,\vec{k})$ respecting the dispersion
relation

$\omega^{2}=k^{2}+m^{2},$ (69)

and, therefore, the function $\eta_{A}(\omega)=k_{A}(\omega)/\omega$ is given
by

$\eta_{\Phi}(\omega)=\eta_{\Psi}(\omega)=\sqrt{1-\left(\frac{m}{\omega}\right)^{2}}.$
(70)

Notice that this is a propagating mode provided that $\omega>m$. Hence, $m$ is
a cutoff frequency for the massive scalar degree of freedom.

Evaluating the scalar longitudinal gauge-invariant variable defined by Eq.
(33) we find

$\Theta=-\left(\frac{m}{\omega}\right)^{2}\Phi.$ (71)

Thus, there is a non-zero contribution to the Riemann tensor in the direction
of propagation of the scalar GW. The meaning of this result is that although
one has only one scalar degree of freedom, it generates the effects of the two
scalar polarization states, the scalar transversal and the scalar longitudinal
modes. If $m=0\Rightarrow\Theta=0$, the longitudinal effect vanishes and one
restores the result of the original massless Brans-Dicke theory for which
there is only the scalar transversal polarization mode [4, 5].

Notice that GWs in the scalar-tensor theories of gravity enter in case 1
discussed in Section III.3. Therefore the variables $\Phi$, $\Theta$ and
$\Xi_{i}$ are Lorentz invariant quantities.

Although in the present derivation we have considered the action (53), the
results are valid for the Horndeski theory [48], which is the most general
scalar-tensor theory of gravity with second-order equations of motion. This is
because our results depend essentially on the weak field equations (59) and
(61). The linearized field equations of Horndeski theory acquire exactly these
forms, with a redefinition of the mass $m$ (see Eqs. (17) and (18) of the Ref.
[29]). In Ref. [29], there is a similar discussion about the polarization
states in scalar-tensor theories, though a different approach has been used.
However, if the Palatini formalism is used to study the Horndeski theory, the
number of scalar polarization modes can vary depending on the chosen
parameters [25].

For purposes that will be clear in the following sections, notice that from
the trace of Eq (61) and Eq. (59) we obtain $\delta
R=3m^{2}\delta\varphi/\varphi_{0}$. Thus, provided $m\neq 0$ the linearized
field equations in scalar-tensor theory can be regarded as

$\delta
G_{\mu\nu}+\frac{1}{3m^{2}}\left(\eta_{\mu\nu}\Box-\partial_{\mu}\partial_{\nu}\right)\delta
R=0,$ (72)

with trace given by

$\left(\Box-m^{2}\right)\delta R=0.$ (73)

In this way, the scalar field does not appear explicitly in the equations
except in the definition of the mass.

####  III.4.3 $f(R)$-gravity

The action for the $f(R)$-gravity is defined as an extension of the Einstein-
Hilbert action which, in the absence of matter, has the following form

$I=\frac{1}{16\pi G}\int d^{4}x\sqrt{-g}f(R),$ (74)

where $f(R)$ is an arbitrary function of the Ricci scalar. If we vary this
action with respect to the metric $g_{\mu\nu}$ we obtain the vacuum field
equations

$f^{\prime}R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}f-\nabla_{\mu}\nabla_{\nu}f^{\prime}+g_{\mu\nu}\square
f^{\prime}=0,$ (75)

where in this section we use a prime to denote the derivative with respect to
$R$. Additionally, the trace of Eq. (75) gives

$\square f^{\prime}+\frac{Rf^{\prime}-2f}{3}=0.$ (76)

It is well known that the $f(R)$ gravity is equivalent to a scalar-tensor
theory of gravity [49]. Therefore, we expect the same results for the
polarization modes as obtained in the previous section. We show this
equivalence by directly solving the equations in the weak-field approximation.

First, notice that Minkowski is not a vacuum solution of the theory.
Therefore, different from Einstein’s gravity, to study vacuum GWs in
$f(R)$-gravity we should expand the metric around a non-flat background metric
[50]

$g_{\mu\nu}=g_{\mu\nu}^{\rm(b)}+h_{\mu\nu},$ (77)

where $g_{\mu\nu}^{\rm(b)}$ is a background metric with constant curvature (de
Sitter or anti-de Sitter). In this sense, the perturbed Ricci scalar and the
perturbed function $f(R)$ become

$\displaystyle R$ $\displaystyle=R_{\rm b}+\delta R+{\mathcal{O}}(h^{2}),$
(78) $\displaystyle f(R)$ $\displaystyle=f(R_{\rm b})+f^{\prime}(R_{\rm
b})\delta R+\frac{1}{2}f^{\prime\prime}(R_{b})\delta
R^{2}+{\mathcal{O}}(h^{3}),$ (79)

where $R_{\rm b}$ is the constant background curvature scalar.

With this expansion in the Eq. (76) we obtain

$\left(\square-m^{2}\right)\delta R=0,$ (80)

where

$m^{2}\equiv\frac{1}{3}\left(\frac{f^{\prime}_{\rm b}}{f^{\prime\prime}_{\rm
b}}-R_{\rm b}\right),$ (81)

and $f_{\rm b}=f(R_{\rm b})$. Notice that now all the covariant derivatives
are evaluated using the background metric.

Moreover, from the Eq. (75) we obtain the following equation for the
perturbation of the Ricci tensor [50]

$\displaystyle\delta R_{\mu\nu}$
$\displaystyle+\left(\frac{f^{\prime\prime}_{\rm b}}{f^{\prime}_{\rm
b}}R_{\mu\nu}^{\rm b}-\frac{1}{2}g_{\mu\nu}^{\rm b}\right)\delta
R-\frac{1}{2}\frac{f_{\rm b}}{f^{\prime}_{\rm b}}\delta g_{\mu\nu}$
$\displaystyle+\frac{f^{\prime\prime}_{\rm b}}{f^{\prime}_{\rm
b}}\left(g_{\mu\nu}^{\rm(b)}\square-\nabla_{\mu}\nabla_{\nu}\right)\delta
R=0.$ (82)

At the scale size of the GW detectors, one can assume a nearly Minkowski
background metric $g_{\mu\nu}^{\rm(b)}\approx\eta_{\mu\nu}$ and
$R_{\mu\nu}^{\rm b}\approx 0$. Let us assume $f(R)$ models for which $f(R_{\rm
b})\approx 0$ at this limit, but in general $f_{\rm b}^{\prime}\neq 0$ and
$f_{\rm b}^{\prime\prime}\neq 0$. This is the case of some $f(R)$ models which
are viable alternatives to explain the accelerated expansion of the Universe
[50]. In this limit the d’Alembertian operator in Eq. (80) is
$\square=\eta^{\mu\nu}\partial_{\mu}\partial_{\nu}$, and Eq. (III.4.3)
simplifies to

$\delta
G_{\mu\nu}+\frac{1}{3m^{2}}\left(\eta_{\mu\nu}\square-\partial_{\mu}\partial_{\nu}\right)\delta
R=0,$ (83)

which is identical to Eq. 72. Therefore, we can follow the same procedure as
in the case of scalar-tensor theories to find the equations for the gauge-
invariant variables in the $f(R)$-gravity

$\left(\square-m^{2}\right)\Phi=0,$ (84) $\nabla^{2}\Xi_{i}=0,$ (85) $\square
h_{ij}^{TT}=0,$ (86)

where

$\Psi=\Phi=-\frac{R}{6m^{2}},$ (87)

and again we conclude that

$\Xi_{i}=0.$ (88)

Thus, as in scalar-tensor theories of gravity, we conclude that $f(R)$-gravity
presents three propagating degrees of freedom. The usual two tensor modes
propagating at the speed of light and one scalar degree of freedom with
transversal and longitudinal behavior. Such a conclusion was found previously
by Moretti et al. [34], where the authors have also used gauge-invariant
variables to describe the polarization modes, though with some slight
differences when compared with our derivation. An analog discussion about the
degrees of freedom and the polarization states of GWs in $f(R)$-gravity can be
found, e.g., in Ref. [27] where the Lorentz gauge has been used. Regarding
Lorentz transformations, GWs in the scope of $f(R)$-gravity are also within
case 1 of Section III.3.

####  III.4.4 A class of quadratic theories of gravity

Now, let us apply the gauge-invariant formalism to a wide class of alternative
theories of gravity considered by Yunes and Stein [51] in the scope of
nonspinning black holes and the parametrized post-Einsteinian framework for
GWs. The modified Einstein-Hilbert action is given by

$\displaystyle I=$ $\displaystyle\int d^{4}x\sqrt{-g}\Big{\\{}\frac{R}{16\pi
G}+\tilde{\alpha}_{1}f_{1}(\varphi)R^{2}+\tilde{\alpha}_{2}f_{2}(\varphi)R_{\mu\nu}R^{\mu\nu}$
$\displaystyle+\tilde{\alpha}_{3}f_{3}(\varphi)R_{\lambda\mu\nu\kappa}R^{\lambda\mu\nu\kappa}+\tilde{\alpha}_{4}f_{4}(\varphi){R_{\lambda\mu\nu\kappa}}^{\ast}R^{\lambda\mu\nu\kappa}$
$\displaystyle-\frac{\beta}{2}[\nabla_{\mu}\varphi\nabla^{\mu}\varphi+V(\varphi)]\Big{\\}},$
(89)

where ${R_{\lambda\mu\nu\kappa}}^{\ast}$ is the dual of the Riemann tensor,
$\varphi$ is a scalar-field, and $(\tilde{\alpha}_{i},\beta)$ are coupling
constants. Theories described by actions of this type can be motivated by the
low-energy expansions of string theory. Notice that the above action differs
from that considered in [38].

Since we have already discussed the role of the scalar fields in the
polarization modes of GWs, let us consider only the case of constant
couplings, and the scalar field is absent ($\beta=0$). Therefore we can set
$f_{i}(\varphi)=1$, and from the variation of the action (III.4.4) we obtain
the field equations

$G_{\mu\nu}+\alpha_{1}\mathcal{H}_{\mu\nu}+\alpha_{2}\mathcal{I}_{\mu\nu}+\alpha_{3}\mathcal{J}_{\mu\nu}=0,$
(90)

where $\alpha_{i}=16\pi G\tilde{\alpha}_{i}$, and

$\displaystyle\mathcal{H}_{\mu\nu}\equiv~{}$ $\displaystyle
2R_{\mu\nu}R-\frac{1}{2}g_{\mu\nu}R^{2}-2\nabla_{\mu}\nabla_{\nu}R+2g_{\mu\nu}\Box
R,$ (91) $\displaystyle\mathcal{I}_{\mu\nu}\equiv~{}$ $\displaystyle\Box
R_{\mu\nu}+2R_{\mu\lambda\nu\kappa}R^{\lambda\kappa}-\frac{1}{2}g_{\mu\nu}R_{\lambda\kappa}R^{\lambda\kappa}$
$\displaystyle+\frac{1}{2}g_{\mu\nu}\Box R-\nabla_{\mu}\nabla_{\nu}R,$ (92)
$\displaystyle\mathcal{J}_{\mu\nu}\equiv~{}$ $\displaystyle
8R^{\lambda\kappa}R_{\mu\lambda\nu\kappa}-2g_{\mu\nu}R^{\lambda\kappa}R_{\lambda\kappa}+4\Box
R_{\mu\nu}$
$\displaystyle-2RR_{\mu\nu}+\frac{1}{2}g_{\mu\nu}R^{2}-2\nabla_{\mu}\nabla_{\nu}R.$
(93)

The trace of Eq. (90) is

$2(3\alpha_{1}+\alpha_{2}+\alpha_{3})\Box R-R=0.$ (94)

Notice that in the case of $f(R)$-gravity, a Klein-Gordon type equation was
found for the first-order perturbation of the Ricci scalar [see Eq. (80)]. On
the other hand, the above equation for quadratic gravity is valid for the
Ricci scalar in general, and not only for its perturbation. Nonetheless, in
the present article, we are concerned only with the first-order perturbation
of the Eq. 94 considering some background metric

$2(3\alpha_{1}+\alpha_{2}+\alpha_{3})\Box\delta R-\delta R=0.$ (95)

To obtain the perturbation equations we can start considering a non-flat
background metric, and then make the flat approximation in the region of a
detector as we did for $f(R)$-gravity. However, in the present case, we
consider a Minkowski background metric from the beginning for simplicity. As
we saw in the previous subsection, the difference between the two approaches
is a redefinition of the mass parameters of the theory (see Eq. 81).

Thus, with a Minkowski background metric, the perturbation of the Eq. (90)
reads

$\displaystyle(\alpha_{2}$ $\displaystyle+4\alpha_{3})\Box\delta
G_{\mu\nu}+\delta G_{\mu\nu}$
$\displaystyle+(2\alpha_{1}+\alpha_{2}+2\alpha_{3})(\eta_{\mu\nu}\Box-\partial_{\mu}\partial_{\nu})\delta
R=0.$ (96)

It is worth noticing that the above equation is a generalization of the Eq.
(83) including a term with a D’Alembertian operator applied to $\delta
G_{\mu\nu}$. Notice that within the present framework, it is possible to
identify four subclasses of theories from quadratic gravity, not only two as
shown by Tachinami et al. [38]. In what follows we show that these subclasses
can differ in the number of independent radiative degrees of freedom,
polarizations, and the expressions for the scalar polarizations.

  * •

Subclass 1: $(3\alpha_{1}+\alpha_{2}+\alpha_{3})>0$,
$(\alpha_{2}+4\alpha_{3})<0$

In this subclass, we have the following equations from Eq. (95) and the 00
component of Eq. (III.4.4)

$\left(\Box-m^{2}\right)\delta R=0,$ (97)

and

$\left(\Box-M^{2}\right)\Phi=-\gamma\delta R,$ (98)

where

$\displaystyle m^{2}\equiv$
$\displaystyle\frac{1}{2(3\alpha_{1}+\alpha_{2}+\alpha_{3})},~{}~{}M^{2}\equiv-\frac{1}{(\alpha_{2}+4\alpha_{3})},$
$\displaystyle\gamma\equiv$
$\displaystyle\frac{2\alpha_{1}+\alpha_{2}+2\alpha_{3}}{2(\alpha_{2}+4\alpha_{3})}=\frac{1}{6}\left[1-\left(\frac{M}{m}\right)^{2}\right].$
(99)

From the $0i$ components of Eq. (III.4.4) and the equations for the scalars we
obtain the equation for the gauge-invariant vector

$\left(\Box-M^{2}\right)\Xi_{i}=0.$ (100)

Finally, from the $ij$ components of Eq. (III.4.4) we obtain

$\partial_{i}\partial_{j}\left(\Box\Psi-M^{2}\Psi+\gamma
R\right)+\frac{1}{2}\Box\left(\Box h_{ij}^{TT}-M^{2}h_{ij}^{TT}\right)=0.$
(101)

Taking the divergence of the above equation and remembering that
$\partial^{i}h_{ij}^{TT}=0$, we obtain

$\left(\Box-M^{2}\right)\Psi=-\gamma\delta R,$ (102)

and comparing with the Eq. (98) we conclude that

$\Psi=\Phi.$ (103)

Using this identity in the definition of the perturbed curvature scalar (30),
the equations (97) and (98) result

$\Box\left(\Box-m^{2}\right)\Phi=0.$ (104)

Finally, we obtain the equations for the tensor gauge-invariant quantities by
using the equation (102) back in Eq. (101)

$\Box\left(\Box h_{ij}^{TT}-M^{2}h_{ij}^{TT}\right)=0.$ (105)

Therefore, we conclude that all the gauge-invariant variables describe the
propagating modes of GWs. From Eq. (100) we have two propagating vector
degrees of freedom with mass $M$. For the scalar sector, we can obtain two
independent propagating solutions from Eq. (104), one describes a massless
scalar and the other a massive scalar with mass $m$. A similar conclusion can
be drawn for the tensor sector from Eq. (105). Now, we have two massless
degrees of freedom and two degrees of freedom of a massive tensor field with
mass $M$. Thus, the quadratic gravity in this general case presents 8
propagating degrees of freedom.

All six polarization modes of GWs are present in this case. Only the massive
solution of the scalar $\Phi$ contributes to the scalar longitudinal
polarization described by the variable $\Theta$. On the other hand, both
solutions contribute to the scalar transversal mode. Similarly, the massless
and the massive solutions of the tensor variable describe the $+$ and $\times$
polarizations simultaneously. Finally, since $\Xi_{i}$ are two non-null
propagating modes, none of the variables are Lorentz invariant, as seen from
the Eqs. (41), (42),(III.3), and (III.3).

As a remark notice that, for the scalar sector the solution of a system of
equations similar to Eqs. (97) and (98) in the presence of matter was found by
Vilhena et al. [52].

  * •

Subclass 2: $(2\alpha_{1}+\alpha_{2}+2\alpha_{3})=0$

This is a subclass of theories whose solution can be obtained from the
previous general solution. From Eq. (III.4.4) notice that $\gamma=0$ and so
there is only one mass scale

$m^{2}=M^{2}=\frac{1}{2(\alpha_{1}-\alpha_{3})}.$ (106)

Using $\gamma=0$ in Eq. (98) and (102) we find that the scalars respect the
same Klein-Gordon type equation as well as $\delta R$ [see Eq. (97)]. On the
other hand, we can not conclude that $\Psi=\Phi$ as in the previous case, but
they are related through Eq. (30). If we further assume that the scalars,
$\Psi$, $\Phi$, and $\delta R$ are functions of the retarded time, the Eq.
(30) leads to

$\delta R=2\omega^{2}\left[\eta^{2}\Psi+(2\eta^{2}-3)\Phi\right],$ (107)

where $\eta=\sqrt{1-(m/\omega)^{2}}$.

Thus, this subclass presents only one independent scalar degree of freedom. It
is a propagating mode provided $\alpha_{1}>\alpha_{3}$. The number of degrees
of freedom in the vector and tensor sectors is the same as before since they
respect the same equations. We conclude that seven degrees of freedom describe
the six polarization modes of GWs.

  * •

Subclass 3: $(3\alpha_{1}+\alpha_{2}+\alpha_{3})=0$

Within this condition, the Eq. (95) reduces to

$\delta R=0,$ (108)

and the Eq. (III.4.4) simplifies to

$(\alpha_{2}+4\alpha_{3})\Box\delta G_{\mu\nu}+\delta G_{\mu\nu}=0.$ (109)

From these equations, we obtain $\Phi$ and $\Psi$ satisfying a Klein-Gordon
equation, each with the same mass

$M^{2}=\frac{1}{3(\alpha_{1}-\alpha_{3})}.$ (110)

The relation between the two scalars can be found after combining Eqs. (30)
and (108). Assuming they are oscillatory functions of the retarded time, we
find

$\Psi=\left(\frac{3}{\eta^{2}}-2\right)\Phi,$ (111)

where $\eta=\sqrt{1-(M/\omega)^{2}}$, and the scalar longitudinal polarization
is described by

$\Theta=2(1-\eta^{2})\Phi=2\left(\frac{M}{\omega}\right)^{2}\Phi.$ (112)

The two vector modes also have mass $M$ and the tensor modes are described by
two massless and two massive degrees of freedom with the mass $M$. Thus, again
we have seven independent degrees of freedom describing the six polarization
states of GWs. As before, we must have $\alpha_{1}>\alpha_{3}$ to the massive
modes represent propagating modes.

  * •

Subclass 4: $\alpha_{2}=-4\alpha_{3}$

In this subclass, the term with a D’Alembertian operator applied to the
Einstein tensor vanishes in Eqs. (III.4.4) and we obtain

$\delta
G_{\mu\nu}+(2\alpha_{1}+\alpha_{2}+2\alpha_{3})(\eta_{\mu\nu}\Box-\partial_{\mu}\partial_{\nu})\delta
R=0,$ (113)

which is identical in form to the Eq. (61) of the scalar-tensor theory and to
the Eq. (83) of $f(R)$-gravity. Therefore, we obtain the same results as in
those cases. We have only one independent scalar degree of freedom, and the
scalars respect the identity

$\Psi=\Phi=-\frac{\delta R}{6m^{2}},$ (114)

with the mass $m$ given by the definition (III.4.4). In the present case, it
simplifies to

$m^{2}=\frac{1}{6(\alpha_{1}-\alpha_{3})}.$ (115)

Again, the scalar degree of freedom is a propagating mode provided
$\alpha_{1}>\alpha_{3}$.

The vector components vanish identically and we have only the two massless
propagating degrees of freedom in the tensor sector. Theories within this
subclass present four polarizations, the scalar transversal described by
$\Phi$, the scalar longitudinal $\Theta=-(m/\omega)^{2}\Phi$, and the two
tensor polarizations propagating at the speed of light. Furthermore, the
variables $\Phi$, $\Theta$, and $\Xi_{i}$ are Lorentz invariant quantities
since this subclass of theories entered in case 1 discussed in Section III.3.

####  III.4.5 Towards a general parametrization

The linearized field equations of the theories studied in the preceding
subsections are subclasses of the following equations

$\displaystyle A_{1}$ $\displaystyle\Box\delta G_{\mu\nu}+\delta G_{\mu\nu}$
$\displaystyle+\frac{1}{3}\left(\eta_{\mu\nu}\Box-\partial_{\mu}\partial_{\nu}\right)\left[(A_{1}+A_{2})\delta
R+3B\frac{\delta\varphi}{\varphi_{0}}\right]=0,$ (116)

where $A_{1}$ and $A_{2}$ are constants with dimension of square length and
$B$ is a dimensionless constant. The above equation must be supplemented by

$\Box\delta\varphi=0,$ (117)

if there is a coupling with a massless scalar field ($B\neq 0$). On the other
hand, the Eq. (III.4.5) with $B=0$ and $A_{1}=0$ is enough to describe the
case of non-minimal coupling with a massive scalar field. This was shown at
the end of subsection III.4.2.

Thus, comparing Eq. (III.4.5) with those found for scalar-tensor theory,
$f(R)$-gravity and quadratic gravity [see Eqs. (61), (72), (83), and
(III.4.4)], we see that it is possible to reduce the problem of determining
the number of radiative degrees of freedom and the polarization content of a
given theory to the problem of finding three parameters, namely, $A_{1}$,
$A_{2}$, and $B$. This is true if a theory of gravity has linearized field
equations given by (III.4.5), which encompasses a wide variety of theories as
summarized in Table 1. This table shows that these three parameters are
closely related to the gauge-invariant variables. Determining their values or
relations between them implies a theory’s polarization content and the number
of independent radiative degrees of freedom.

####  III.4.6 Number of radiative degrees of freedom versus number of
polarization modes

In the previous subsections, we noticed a distinction between the number of
radiative degrees of freedom and the number of polarization modes of GWs in
gravity theories. Discussions and some criticisms appear in several recent
works, e.g., [27] and [29]. In the language of [4, 5], the number of
polarization modes corresponds to how GWs interact with a sphere of test
particles. It has nothing to do with the number of independent dynamical
degrees of freedom of the linearized theory, which can be smaller or bigger
than the number of polarization modes. In the case of massless scalar GWs, for
instance, relative acceleration between the particles in the sphere is
observed in the direction orthogonal to the wave vector $\vec{k}$. On the
other hand, for a massive scalar mode, relative accelerations are also
generated for particles located in the wave’s propagation direction. In the
latter case, we say the theory presents two scalar polarization modes, though
these modes are not independent. Therefore, GWs in the $f(R)$-gravity and the
scalar-tensor theories, for instance, present four polarization modes if
$m\neq 0$, with three independent radiative degrees of freedom. For these
theories, the number of polarization modes and independent degrees of freedom
agrees only if $m=0$. In the present work, we have defined a gauge-invariant
variable $\Theta$ in Eq. (33) which enables an unambiguous determination of
the existence of the scalar longitudinal GW mode. For other theories studied
in this article, the number of radiative degrees of freedom and the number of
polarization modes are depicted in Table 1. We see that the number of
independent radiative degrees of freedom and the number of polarization modes
of a theory depends essentially on the parameters $A_{1},~{}A_{2}$, and $B$.

Table 1: Here we show how the number of independent radiative degrees of freedom (d.o.f.) and non-null polarization modes of GWs depend on the parameters of Eq. (III.4.5). | # Independent | $h_{ij}^{TT}$ | $\Xi_{i}$ | $\Phi$ | $\Theta$ | Examples of theories
---|---|---|---|---|---|---
| radiative d.o.f. |  |  |  |  |
$A_{1}=A_{2}=B=0$ | 2 | ✓ | 0 | 0 | 0 | GR
$A_{1}=A_{2}=0,~{}B=1$ | 3 | ✓ | 0 | ✓ | 0 | Brans-Dicke
$A_{1}=0,~{}A_{2}>0,~{}B=0$ | 3 | ✓ | 0 | ✓ | ✓ | Horndeski, $f(R)$-gravity,
|  |  |  |  |  | subclass 4 of quadratic gravity
$A_{1}<0,~{}A_{2}>0,~{}B=0$ | 8 | ✓ | ✓ | ✓ | ✓ | subclass 1 of quadratic gravity
$A_{1}=-A_{2},~{}A_{2}>0,~{}B=0$ | 7 | ✓ | ✓ | ✓ | ✓ | subclass 2 of quadratic gravity
$A_{1}<0,~{}A_{2}=0,~{}B=0$ | 7 | ✓ | ✓ | ✓ | ✓ | subclass 3 of quadratic gravity

##  IV Pulsar timing sensitivity

###  IV.1 Gauge-invariance and physical observables

Once we have defined the polarization modes of GWs in terms of gauge-invariant
variables, the sensitivities of the GW detectors are necessarily connected to
such variables. This is because gauge-invariant quantities express truly
physical observables.

In the present work, we focus on the sensitivity of the pulsar timing
technique to each polarization mode of GWs. Previous works consider the
synchronous gauge in evaluating the pulsar timing sensitivity. However, using
this gauge leaves residual gauge freedom and, apart from this, could not be
appropriate to use such a gauge in all metric theories of gravity. Moreover,
any GW waveform originating from a compact binary system, for instance, should
be expressed using gauge-invariant variables to ensure it is a physical
quantity. To evaluate the detectability of such a signal, it should be
compared against a sensitivity curve which was also evaluated using the same
gauge-invariant quantities.

Taking into account these aspects, we estimate the sensitivity considering the
gauge-invariant variables to evaluate the pulsar timing response to the
polarization modes. The elementary observable for interferometric detectors
and the pulsar timing technique is the “one-way” fractional frequency shift
$y(t)=[\nu(t)-\nu_{0}]/\nu_{0}$, where $\nu(t)$ is the frequency of an
electromagnetic signal at the time of reception $t$ and $\nu_{0}$ is the
unperturbed frequency. To achieve our goals we need a relation between $y$ and
the gauge-invariant variables. However, the majority of the derivations
appearing in the literature are gauge-dependent. A full gauge-invariant
derivation of $y$ was obtained by Koop and Finn [53]. Their derivation is
quite general and includes the GW effect and all possible contributions from
the background curvature (e.g., Rømer delay, aberration, Shapiro time delay,
and other effects appear naturally). Furthermore, no assumptions were made
about the size of the detectors compared to the GW wavelength or on the
dispersion relation of GWs.

For the present article, it is enough to consider the special case of a
Minkowski background. If in addition, we consider that the source and the
receiver of the electromagnetic signal are at rest in the same global Lorentz
frame, the equation for $y$ can be written as [53]

$\frac{dy(t)}{dt}=-\int_{0}^{\lambda_{R}}\delta{R}_{0i0j}n^{i}n^{j}d\lambda,$
(118)

where $n^{i}$ is the spatial unit vector in the direction of the link between
the source and the receiver, and $\lambda$ is the photon’s affine parameter
along its unperturbed trajectory. Therefore, the GW contribution to the time
derivative of the frequency shift $y$ is given by the projection of the
Riemann tensor integrated along the unperturbed null geodesic linking the
source and the receiver of the electromagnetic signal. Recently, Błaut [54]
found the same result although using a quite different approach.

Notice that in the Koop and Finn derivation, there is no specification of the
field equations of the underlying theory of gravity. The validity of their
derivation lies in using the Riemann tensor as the fundamental quantity to
describe the spacetime geometry and the geodesic deviation equation has the
form as it appears in GR. Therefore, Eq. (118) is valid for all four-
dimensional metric theories of gravity with these properties as, for instance,
those theories presented in Section III. Thus, this equation is appropriate
for obtaining the sensitivity of interferometers and of the pulsar timing
technique to the polarization modes of GWs in alternative theories of gravity
in a gauge-invariant fashion.

To evaluate the one-way response let us consider the emitter of a light signal
located at point 1 at a distance $L$ from the receiver. The receiver is
located at point 2 at the origin of the coordinate system. The trajectory of
the light signal can be parametrized as

$t=t_{2}-(L-\lambda),~{}~{}\vec{r}=(L-\lambda)\hat{n},$ (119)

with $\lambda\in[0,~{}L]$; $t_{1}$ is the time of emission, $t_{2}=t_{1}+L$ is
the time of reception and $\hat{n}$ is the unit vector pointing from 2 to 1.
Within this parametrization, the retarded times are given by

$u_{A}=t_{2}-(1+\eta_{A}\mu)(L-\lambda),~{}~{}\lambda\in[0,~{}L],$ (120)

where we have defined $\mu\equiv\hat{k}\cdot\hat{n}$.

Now, changing the variable of integration to the retarded time $u_{A}$, the
integration along the unperturbed trajectory of the light signal in Eq. (118)
can be performed. Using the Eq. (III.2) in (118) we find

$\displaystyle\frac{dy}{dt}=$
$\displaystyle-\left(\frac{\mu^{2}}{1+\eta_{S}\mu}\right)\left[\Theta^{\prime}(t)-\Theta^{\prime}\left(t-(1+\eta_{S}\mu)L\right)\right]$
$\displaystyle+\left(\frac{1-\mu^{2}}{1+\eta_{S}\mu}\right)\left[\Phi^{\prime}(t)-\Phi^{\prime}\left(t-(1+\eta_{S}\mu)L\right)\right]$
$\displaystyle+\left(\frac{\eta_{V}\mu~{}n^{i}}{1+\eta_{V}\mu}\right)\left[\Xi_{i}^{\prime}(t)-\Xi_{i}^{\prime}\left(t-(1+\eta_{V}\mu)L\right)\right]$
$\displaystyle+\frac{1}{2}\left(\frac{n^{i}n^{j}}{1+\eta_{T}\mu}\right)\left[h_{ij}^{{\rm
TT}\prime}(t)-h_{ij}^{{\rm TT}\prime}\left(t-(1+\eta_{T}\mu)L\right)\right].$
(121)

In the final expression, we have replaced the time of reception
$t_{2}\rightarrow t$, a prime denotes derivative with respect to the retarded
time and, for simplicity, we have considered
$\eta_{\Psi}=\eta_{\Phi}=\eta_{S}$. In the synchronous gauge and for
$\eta_{A}=1$, the above equation coincides with the frequency shift derived in
Refs. [40, 39].

Therefore, we have obtained an explicit relation between a physical observable
(the time derivative of the frequency shift) with gauge-invariant quantities
that describe the six possible polarization states of GWs. Furthermore, we
have not made any hypotheses regarding the four dispersion relations except
the equality of the dispersion relations of the scalar modes. The result above
applies to pulsar timing, spacecraft Doppler tracking, and ground-based and
space-based interferometric GW detectors. Here we specialize in the case of
the pulsar timing technique.

###  IV.2 Pulsar timing responses and sensitivities

To derive the pulsar timing response, let us consider the Earth located at the
origin of a Cartesian system of coordinates with unit vectors
$(\hat{i},\hat{j},\hat{k})$ oriented in the $x$, $y$ and $z$ directions
respectively. The GW wave vector is in the direction of $\hat{k}$ and the
vector $L\hat{n}$ locates the Pulsar, where $\hat{n}$ is the unit vector
pointing from Earth to the Pulsar. The Pulsar emits electromagnetic signals
continuously which are detected on Earth. With the help of Eq. (IV.1), we can
obtain the Fourier transform of the frequency shift induced by GWs on the
signal emitted by the Pulsar. Using the property of the time derivative of the
Fourier transform and the relation between the time $t$ and the retarded time
we find

$\tilde{y}(f)=\tilde{y}_{SL}(f)+\tilde{y}_{ST}(f)+\tilde{y}_{V}(f)+\tilde{y}_{T}(f),$
(122)

where the Fourier transforms of the induced frequency shifts due to each
gauge-invariant variable are given by

$\displaystyle\tilde{y}_{SL}(f)$
$\displaystyle\equiv-\left(\frac{\mu^{2}}{1+\eta_{S}\mu}\right)H_{SL}(f)\left[1-e^{i2\pi
fL(1+\eta_{S}\mu)}\right],$ (123) $\displaystyle\tilde{y}_{ST}(f)$
$\displaystyle\equiv\left(\frac{1-\mu^{2}}{1+\eta_{S}\mu}\right)H_{ST}(f)\left[1-e^{i2\pi
fL(1+\eta_{S}\mu)}\right],$ (124) $\displaystyle\tilde{y}_{V}(f)$
$\displaystyle\equiv\left(\frac{\eta_{V}\mu~{}n^{i}}{1+\eta_{V}\mu}\right)\tilde{\Xi}_{i}(f)\left[1-e^{i2\pi
fL(1+\eta_{V}\mu)}\right],$ (125) $\displaystyle\tilde{y}_{T}(f)$
$\displaystyle\equiv\frac{1}{2}\left(\frac{n^{i}n^{j}}{1+\eta_{T}\mu}\right)\tilde{h}_{ij}^{TT}(f)\left[1-e^{i2\pi
fL(1+\eta_{T}\mu)}\right],$ (126)

where $H_{SL}(f)\equiv\tilde{\Theta}(f)$ and $H_{ST}(f)\equiv\tilde{\Phi}(f)$
are the frequency-dependent wave amplitude for the scalar longitudinal and
scalar transversal polarizations respectively.

For the polarization mode ‘$A$’ we define the angular response $R_{A}$ of a
single pulsar timing as $R^{2}_{A}\equiv|\tilde{y}_{A}(f)|^{2}/H_{A}^{2}(f)$.
It follows

$R_{SL}^{2}=2\left(\frac{\mu^{2}}{1+\eta_{S}\mu}\right)^{2}\Big{[}1-\cos\big{(}2\pi
fL(1+\eta_{S}\mu)\big{)}\Big{]},$ (127)

and

$R_{ST}^{2}=2\left(\frac{1-\mu^{2}}{1+\eta_{S}\mu}\right)^{2}\Big{[}1-\cos\big{(}2\pi
fL(1+\eta_{S}\mu)\big{)}\Big{]}.$ (128)

In the case of vector and tensor polarization modes, we assume an elliptically
polarized wave and then average over the polarizations to find the response.
For an elliptically polarized vector GW, we have

$\tilde{\Xi}_{i}(f)=H_{V}(f)\left(e^{i\varphi}\sin\Gamma{\epsilon}^{(1)}_{i}+\cos\Gamma{\epsilon}^{(2)}_{i}\right),$
(129)

where $H_{V}(f)$ is the vector wave amplitude, and ${\epsilon}^{(1)}_{i}$ and
${\epsilon}^{(2)}_{i}$ are two orthogonal unit polarization vectors and both
are orthogonal to $\hat{k}$. We use the polarization angles $(\varphi,\Gamma)$
to characterize elliptically polarized waves. Two particular cases are
linearly and circularly polarized waves. The former can be obtained by
choosing $\varphi=0$ representing a vector wave $\mathbf{\Xi}$ linearly
polarized making an angle $\Gamma$ with $\hat{\bm{\epsilon}}^{(2)}$.
Circularly polarized waves are obtained by using $\Gamma=\frac{\pi}{4}$ (the
vector wave has the same amplitude in both directions) and
$\varphi=\pm\frac{\pi}{2}$ (right circularly polarized wave for the plus sign
and left circularly polarized wave for the minus sign). The polarization of a
GW depends essentially on the generating mechanism. In what follows, we
evaluate the response as an average over the polarization angles in the
intervals $\varphi\in[0,~{}2\pi]$ and $\Gamma\in[0,~{}\pi]$.

An analogous expression can be written for the usual elliptically polarized
tensor GWs

$\tilde{h}_{ij}^{TT}(f)=H_{T}(f)\left(e^{i\varphi}\sin\Gamma{\varepsilon}^{+}_{ij}+\cos\Gamma{\varepsilon}^{\times}_{ij}\right),$
(130)

where we can use the pair of orthogonal vectors
$(\hat{\bm{\epsilon}}^{(1)},\hat{\bm{\epsilon}}^{(2)})$ to define the two
polarization tensors

$\displaystyle{\varepsilon}^{+}_{ij}$
$\displaystyle={\epsilon}^{(1)}_{i}{\epsilon}^{(1)}_{j}-{\epsilon}^{(2)}_{i}{\epsilon}^{(2)}_{j},$
(131) $\displaystyle{\varepsilon}^{\times}_{ij}$
$\displaystyle={\epsilon}^{(1)}_{i}{\epsilon}^{(2)}_{j}+{\epsilon}^{(2)}_{i}{\epsilon}^{(1)}_{j}.$
(132)

In our Cartesian coordinate system we chose $\hat{\bm{\epsilon}}^{(1)}$ and
$\hat{\bm{\epsilon}}^{(2)}$ to coincide with $\hat{i}$ and $\hat{j}$
respectively. Moreover, let us consider the usual spherical coordinates
$(\theta,\phi)$ associated with the vector $L\hat{n}$ that locates the Pulsar.
Notice that in this coordinate system $\mu=\hat{n}\cdot\hat{k}=\cos\theta$.
Then, in the case of vector and tensor waves, we can perform an average over
the polarization angles $(\varphi,\Gamma)$ to find the angular pulsar timing
response

$\displaystyle R_{V}^{2}$
$\displaystyle=\frac{1}{3}\left(\frac{\eta_{V}\mu\sqrt{1-\mu^{2}}}{1+\eta_{V}\mu}\right)^{2}(1+\cos^{2}\phi)$
$\displaystyle\times\Big{[}1-\cos\big{(}2\pi fL(1+\eta_{V}\mu)\big{)}\Big{]},$
(133)

and

$\displaystyle R_{T}^{2}$
$\displaystyle=\frac{1}{12}\left(\frac{1-\mu^{2}}{1+\eta_{T}\mu}\right)^{2}(1+\cos^{2}2\phi)$
$\displaystyle\times\Big{[}1-\cos\big{(}2\pi fL(1+\eta_{T}\mu)\big{)}\Big{]},$
(134)

where we have included a factor $1/2$ in both expressions since there are two
vector polarizations in the first case and two tensor polarizations in the
second case.

In the present article we are interested in the sensitivity to single-source
GW signals averaged over the sky and polarization states. The most promising
GW sources in the pulsar timing band are supermassive binary black holes (with
masses in the range $10^{7}-10^{10}~{}{\rm M}_{\odot}$) hosted in the center
of galaxies. The sensitivity is defined by $\sqrt{S_{y}(f)B}/R^{\rm rms}_{A}$,
where $S_{y}(f)$ is the one-sided power spectral density of the noise
affecting the relative frequency shift of pulsar timing, and $B$ is the
bandwidth. Here, we assume $B$ corresponding to an integration time of 10
years ($B=1{\rm cycle}/10~{}{\rm years}$), and a signal-to-noise ratio SNR =
1. The quantity $R^{\rm rms}_{A}$ is the rms of the response evaluated by
performing an average over sources uniformly distributed over the celestial
sphere. Thus, we find the following formula for the pulsar timing rms response

$\displaystyle R^{\rm rms}_{A}(f)=\frac{2\pi fL}{\sqrt{\eta_{A}}}$
$\displaystyle\left[\sum_{n=1}^{5}\frac{a_{n}^{A}(f)}{(2\pi
fL)^{n}}I^{A}_{n}(f)\right]^{1/2},$ (135)

where, for a given polarization $A$, $I_{n}^{A}(f)$ is a set of five
elementary integrals

$I_{n}^{A}(f)=-\int_{2\pi fL(1-\eta_{A})}^{2\pi fL(1+\eta_{A})}x^{n-3}(\cos
x-1)dx,$ (136)

where $n=1,2,3,4,5$. Therefore, the rms of the responses differ only in the
dispersion relations $\eta_{A}$ and in the coefficients $a_{n}^{A}$ given in
the Appendix. The latter can be functions of the frequency if the speed of
propagation of GWs is different from the speed of light except in the case of
the scalar longitudinal polarization for which the coefficients $a_{n}^{SL}$
are independent of frequency for any speed.

Notice that in evaluating the response we have not considered any specific
form for the dispersion relation. Therefore, the analytical expression (135)
is a general result. To evaluate the effect of the dispersion relation on the
GW response, henceforth we consider that each mode $A$ has an effective mass
$m_{A}$ which results in $\eta_{A}(f)=\sqrt{1-(m_{A}/2\pi f)^{2}}$. This
dispersion relation is valid for a wide range of metric theories of gravity as
we have shown in Section III.4. As we have verified, the masses $m_{A}$ depend
on the specific parameters of each theory. In Fig. 1 we show the pulsar timing
rms response for the scalar longitudinal, scalar transversal, vector, and
tensor polarization modes for a typical pulsar distance $L=1~{}{\rm kpc}$.

Figure 1: We show the rms of the pulsar timing response for all the gauge-
invariant variables of the polarization modes. For all cases, we have used the
massive dispersion relation $\eta_{A}(f)=\sqrt{1-(m_{A}/2\pi f)^{2}}$. The
shaded region indicates the frequency range of the pulsar timing technique.
Notice that if the mass is about that of the upper bound of the LIGO detector
($1.27\times 10^{-23}~{}{\rm eV}$) we have remarkable effects in this range.

Our estimated sensitivity is based on the noise model discussed in [55]. It is
assumed that timing fluctuations due to intergalactic and interplanetary
plasma can be adequately calibrated, the intrinsic pulsar rotational noise is
negligible, and the pulse profile is stable. Under these assumptions, the
spectrum of the noise is given by

$S_{y}(f)=[4.0\times 10^{-31}f^{-1}+3.41\times 10^{-8}f^{2}]{\rm Hz}^{-1},$
(137)

where the lower part of the frequency band ($f<3\times 10^{-8}$ Hz) is limited
by the ground clock noise. For higher frequencies, the dominant noise is a
white timing noise due to an uncertainty of $100$ nsec in the time of arrival
of a pulse.

The resulting pulsar timing sensitivities to the polarization states are shown
in Fig. 2. It shows the strength of a sinusoidal gravitational wave required
to achieve a SNR = 1 over an integration time of 10 years.

###  IV.3 Interpretation of the sensitivities

Notice that the sensitivity to the scalar longitudinal mode is some orders of
magnitude better than the sensitivities of other polarizations, and the
sensitivity to the vector modes can be up to five times better than that of
the tensor mode. The response decreases as the wavelength of the GWs is of the
order or larger than the distance from Earth to the Pulsar (long-wavelength
limit). If $m_{A}=0$ this happens for a tiny frequency, far beyond the pulsar
timing frequency band $(10^{-9}-10^{-6}~{}{\rm Hz})$ (see Fig. 1). On the
other hand, for a non-null mass, a fast decrease in the response can occur in
this band as the technique approaches the long-wavelength limit. The cutoff
frequency for which the response vanishes is related to the mass by

$f_{c}=\left(\frac{m}{m_{\rm up}}\right)3.07\times 10^{-9}~{}{\rm Hz},$ (138)

where we have considered the upper bound on the graviton mass imposed by LIGO,
$m_{\rm up}=1.27\times 10^{-23}~{}{\rm eV}/c^{2}$ [3], as a fiducial mass.
Obviously, the effective mass of the vector and scalar polarizations do not
need to respect this upper bound since it was derived from detections of the
tensor modes.

(a) Tensor

(b) Vector

(c) Scalar transversal

(d) Scalar longitudinal

Figure 2: Sensitivity curves expressed as the strength of a sinusoidal
gravitational wave required to achieve a SNR = 1 over an integration time of
10 years. The sensitivities are expressed in terms of the square root of the
power spectral density (PSD). It is assumed that the pulsar is at a distance
of 1 kpc from Earth. We show the effect of the dispersion relation of massive
GWs. The sensitivity to the massless tensor mode is shown in all figures (gray
curve). Notice that as the mass approaches the upper bound of the LIGO
detector ($1.27\times 10^{-23}~{}{\rm eV}$) we have a remarkable change in the
shape of the sensitivity curves mainly for vector and scalar longitudinal
polarizations. If the mass is of this order, the cutoff frequency $f_{c}$ is
in the pulsar timing band [see Eq. (138)]. For lower frequencies, GWs can not
be detected. Therefore, the evidence of a cutoff frequency or even the
evidence that such a cutoff is not on the pulsar timing frequency band can
lead to a more stringent bound of the effective mass of the graviton than that
presented by ground-based interferometers.

In Fig. 3 we show, as an example, the angular response (i.e., the frequency-
dependent antenna pattern) given by Eqs. (127), (128), (IV.2) and (IV.2) at
the frequency $f=3.1~{}{\rm nHz}$ and $\phi=0$ considering a single Pulsar. We
have chosen this frequency as example for two reasons. First because
supermassive binary black holes have higher GW strain for lower frequencies.
The second reason is that this frequency is higher but close to the cutoff
frequency obtained for $m=m_{\rm up}$. Therefore, if the mass has this value,
GWs with a frequency of $f=3.1~{}{\rm nHz}$, for instance, approach the long-
wavelength regime.

(a) Tensor

(b) Vector

(c) Scalar transversal

(d) Scalar longitudinal

Figure 3: The pulsar timing angular response for the tensor, vector, scalar
transversal and scalar longitudinal polarization modes using a massive
dispersion relation $\eta_{A}(f)=\sqrt{1-(m_{A}/2\pi f)^{2}}$ for each case.
We have considered a typical distance of $L=1$ kpc and an angle $\phi=0$.

When the GW approaches the long-wavelength regime, the symmetry of the
response around $\theta=\pi/2$ is restored for all the polarization modes.
This case is shown in red in Fig. 3 for each polarization. The behavior of the
response, in this case, is similar to that of ground-based interferometers. At
$\theta=\pi/2$ the tensor polarization has the maximum response and the
response vanishes for the scalar longitudinal and vector modes. On the other
hand, the response for the scalar transversal mode is identical in form to the
response for tensor polarization. In the same figure, we notice the
oscillations in the response which comes from the square brackets in the Eqs.
(127), (128), (IV.2) and (IV.2). In the present case, the angles for which
this term vanishes for a given frequency $f$ are given by

$\cos\theta_{n}=\frac{1}{\eta_{A}}\left(\frac{n}{fL}-1\right),~{}~{}n=0,1,2,\dots$
(139)

As we mentioned earlier, for the massless case the pulsar timing is out of the
long-wavelength regime for the entire frequency range. In this case, we can
notice an asymmetry of the response of GWs propagating in the parallel
directions of the electromagnetic signal ($\pi/2<\theta<\pi$) with respect to
GWs propagating in the antiparallel directions ($0<\theta<\pi/2$). For GWs
traveling in parallel directions, the response can be some orders of magnitude
higher than those traveling in antiparallel directions. This effect occurs for
tensor, vector, and scalar modes. However, for the scalar longitudinal and
vector modes, one can notice a remarkable enhancement of the response. This
enhancement effect has been noticed for the first time by the present author
and a collaborator [40, 39]. It is associated with the longitudinal behavior
of the mentioned polarization modes and with the relative direction of the GW
wave vector $\vec{k}$ with respect to the direction of propagation of the
electromagnetic signal emitted by the Pulsar.

In the present scope, the physical origin of the enhancement effect in the
response of longitudinal polarizations can be understood in light of Eq.
(118). First of all, remember that in this case, the Riemann curvature tensor
has components not only transverse to the direction of the propagation of the
GW, but also in the longitudinal direction. Since the Riemann tensor is a
function of the retarded time $u$, and $u$ depends on $\theta$, light rays
coming from different directions ‘see’ the curvature generated by the GW
differently. Consider that we are out of the long-wavelength regime. The light
rays traveling in the opposite directions of the GWs pass through several
maxima and minima of the curvature, which makes their frequency change
continuously. Since the final frequency shift measured at Earth is an
integrated effect of the curvature, the result can be zero for some
directions. On the other hand, those light rays propagating parallel or almost
parallel to the GW experience fewer oscillations of the curvature. In this
case, the final effect can be a higher frequency shift when compared with the
anti-parallel case. This is because the average curvature is higher generating
an increase in the response as $\theta\rightarrow\pi$. When one approaches the
long-wavelength regime, the light signals originating from different
directions experience fewer curvature oscillations and the curvature effect in
the frequency shift becomes symmetric. Finally, in the long-wavelength regime,
the curvature oscillations cannot be noticed at all in a one-way light travel.
In this situation, we have the usual frequency-independent antenna patterns of
ground-based interferometers.

The same argument applies in explaining the asymmetry of the transversal
polarizations (scalar transversal and tensor) out of the long-wavelength
regime. But in this case, the curvature goes to zero as one approaches
$\theta=0$ or $\theta=\pi$ suppressing the enhancement effect for
$\theta\rightarrow\pi$.

In the Fig. 2, we notice that the graviton mass has a remarkable effect on the
sensitivity curves as it approaches the upper bound of the LIGO detector
($1.27\times 10^{-23}~{}{\rm eV}$). For tensor and scalar transversal
polarizations, the predominant effect is a limit in the sensitivity
established by the cutoff frequency $f_{c}$ given by the relation (138). On
the other hand, for vector and scalar longitudinal polarization modes, we have
a significant change in the shape of the sensitivity curve including a change
in the frequency of maximum sensitivity. The sensitivity curves for massive
gravitons are indistinguishable from that of the massless case if the
effective mass is two orders of magnitude smaller than that of the LIGO upper
bound in the case of vector and scalar longitudinal polarizations. Whereas for
the transversal polarizations, it is enough that the graviton mass is one
order of magnitude smaller than $m_{\rm up}$. Remember that $m_{\rm up}$ was
obtained from observations of the tensor mode. This means that, in principle,
the effective mass of the vector and scalar polarizations can be greater than
$m_{\rm up}$. If $m$ is about three orders of magnitude higher than $m_{\rm
up}$ these polarizations would be undetectable in the pulsar timing frequency
band.

##  V Conclusion

We have shown that the Bardeen framework enables a clear description of the
six polarization modes of GWs even if each mode has a general dispersion
relation. The response given by Eq. (IV.1) shows an explicit relation between
a physical observable (the derivative of the frequency shift) and the gauge-
invariant variables. Therefore, this relation means we have a bridge between
theory and experiment, avoiding possible ambiguities of gauge choice. A new
gauge-invariant variable was introduced [see Eq. (33)] aiming for an
unambiguous description of the scalar longitudinal polarization mode.

In the case of a single pulsar timing, we obtained an analytical formula for
the rms response [see Eq. (135)] which is valid for any dispersion relation.
In the case of a dispersion relation of a massive particle, we have seen that
it has a significant impact on the pulsar timing sensitivity to scalar
longitudinal and vector GWs. Remarkably, the effects of the mass on the pulsar
timing sensitivity are particularly noticeable if it is of the order of the
LIGO’s upper bound for the graviton mass ($m_{\rm up}$). If the mass is two
orders of magnitude smaller than $m_{\rm up}$, the sensitivity curves are
indistinguishable from the massless case. On the other hand, in the case of
the scalar transversal and the tensor polarization modes, it is enough that
the mass is one order of magnitude smaller than $m_{\rm up}$ to disregard its
effects on the sensitivity. With a dispersion relation of massive particles
the main physical effect in the case of pure transversal modes is a limitation
in the detectability of these modes established by a cutoff frequency that
depends on the mass. Notice that the effects on the sensitivity appear in the
case of pulsar timing because the cutoff frequency we have considered lies in
the pulsar timing frequency band. But, in principle, the cutoff frequency can
be higher than the pulsar timing band in the case of vector and scalar
polarizations. If this happens, such modes would be undetectable by pulsar
timing experiments. In other words, the absence of detection does not imply
that extra polarization states beyond the tensor polarization do not exist. In
the future, we plan to analyze other dispersion relations of GWs appearing in
the literature to check their implications on the pulsar timing sensitivity.

The detection (or absence of detection) of the polarization modes using the
pulsar timing technique has decisive implications for alternative theories of
gravity. Consider, for instance, the case of the theories studied in Section
III.4 for which the tensor mode is massless and the scalar modes can be
massive. Suppose that the scalar mode has a mass of about that of the LIGO
upper bound, therefore for frequencies approaching the cutoff $f_{c}\sim
3\times 10^{-9}$ Hz the sensitivity of the scalar longitudinal polarization
becomes worse than that of the tensor modes. Below this frequency, the scalar
modes could not be detected (or even be produced!). Thus, suppose we are
looking for GWs only in a frequency band below $f_{c}$, and we detect only
tensor polarizations. We could be led to the wrong conclusion that the scalar
modes do not exist. On the other hand, if there is a cutoff frequency for the
scalar modes, but not for the tensor modes, this could corroborate the scalar-
tensor theories of gravity or $f(R)$-gravity. Moreover, this would lead to a
bound on the mass of the scalar mode.

We have seen that the pulsar timing sensitivity to the scalar longitudinal
mode is some orders of magnitude better than the sensitivity to tensor modes.
However, depending on the theory of gravity this could not be an advantage for
detecting this mode. In the case of the theories we have analyzed, the
amplitude of the scalar longitudinal mode is related to the amplitude of the
scalar transversal mode through a factor $(m/\omega)^{2}$ [see, for instance,
Eq. (71)]. Therefore, if $m$ is much smaller than the smallest detectable
frequency of pulsar timing, the scalar-longitudinal mode can become
undetectable even if the scalar transversal mode is detected. Obviously, these
results apply to scalar-tensor theories of gravity, to $f(R)$-gravity and also
to some subclasses of the quadratic gravity. Other theories may have a
different relation between $\Theta$, $\Phi$, the mass, and the frequency,
leading to different conclusions.

Our analysis shows that the evidence of a cutoff frequency for any
polarization or that such a cutoff is not in the pulsar timing band can lead
to a more stringent bound on the graviton mass than that presented by ground-
based interferometers.

Pulsar timing detection presents a great opportunity to test gravity by
imposing bounds on the polarization modes of GWs. However, to impose such
bounds, it is necessary to compute the expected GW strain in the pulsar timing
band. The most promising source of GWs in the frequency band
$10^{-9}-10^{-6}~{}{\rm Hz}$ is supermassive binary black holes with masses in
the range $10^{7}-10^{10}~{}{\rm M}_{\odot}$. This could be individual sources
or an incoherent superposition of the cosmic population of such systems
forming a stochastic background. The GW strains corresponding to each
polarization mode generated by supermassive black holes depend on the details
of each specific theory. The evaluation of them is out of the scope of the
present article since we are interested in a general formalism for both
describing the polarization modes and their corresponding pulsar timing
sensitivity.

Finally, to evaluate the sensitivity to a stochastic background of
astrophysical or cosmological origin, it is necessary to consider an array of
pulsars and the correlation functions between them, that are distinct for each
polarization mode. The derivation of such sensitivity with gauge-invariant
variables is the subject of a forthcoming publication.

###### Acknowledgements.

The author thanks Dr. Massimo Tinto for helpful discussions and encouragement,
and Lívia R. Alves for continuous encouragement during the development of this
work.

## Appendix

Here we give the frequency-dependent quantities that appear in Eq. (135).

For the scalar-longitudinal response

$\displaystyle a^{SL}_{1}$ $\displaystyle=1,$ (140) $\displaystyle a^{SL}_{2}$
$\displaystyle=-4,$ (141) $\displaystyle a^{SL}_{3}$ $\displaystyle=6,$ (142)
$\displaystyle a^{SL}_{4}$ $\displaystyle=-4,$ (143) $\displaystyle
a^{SL}_{5}$ $\displaystyle=1.$ (144)

For the scalar-transversal response

$\displaystyle a^{ST}_{1}$
$\displaystyle=\left[1-\frac{1}{\eta_{S}^{2}(f)}\right]^{2},$ (145)
$\displaystyle a^{ST}_{2}$
$\displaystyle=\frac{4}{\eta_{S}^{2}(f)}\left[1-\frac{1}{\eta_{S}^{2}(f)}\right],$
(146) $\displaystyle a^{ST}_{3}$
$\displaystyle=\frac{2}{\eta_{S}^{2}(f)}\left[\frac{3}{\eta_{S}^{2}(f)}-1\right],$
(147) $\displaystyle a^{ST}_{4}$ $\displaystyle=-\frac{4}{\eta_{S}^{4}(f)},$
(148) $\displaystyle a^{ST}_{5}$ $\displaystyle=\frac{1}{\eta_{S}^{4}(f)}.$
(149)

For the vector response

$\displaystyle a^{V}_{1}$
$\displaystyle=\frac{1}{4}\left[1-\frac{1}{\eta_{V}^{2}(f)}\right],$ (150)
$\displaystyle a^{V}_{2}$
$\displaystyle=\frac{1}{2}\left[\frac{2}{\eta_{V}^{2}(f)}-1\right],$ (151)
$\displaystyle a^{V}_{3}$
$\displaystyle=\frac{1}{4}\left[1-\frac{6}{\eta_{V}^{2}(f)}\right],$ (152)
$\displaystyle a^{V}_{4}$ $\displaystyle=\frac{1}{\eta_{V}^{2}(f)},$ (153)
$\displaystyle a^{V}_{5}$ $\displaystyle=-\frac{1}{4\eta_{V}^{2}(f)}.$ (154)

For the tensor response

$\displaystyle a^{T}_{1}$
$\displaystyle=\frac{1}{16}\left[1-\frac{1}{\eta_{T}^{2}(f)}\right]^{2},$
(155) $\displaystyle a^{T}_{2}$
$\displaystyle=\frac{1}{4\eta_{T}^{2}(f)}\left[1-\frac{1}{\eta_{T}^{2}(f)}\right],$
(156) $\displaystyle a^{T}_{3}$
$\displaystyle=\frac{1}{8\eta_{T}^{2}(f)}\left[\frac{3}{\eta_{T}^{2}(f)}-1\right]$
(157) $\displaystyle a^{T}_{4}$ $\displaystyle=-\frac{1}{4\eta_{T}^{4}(f)},$
(158) $\displaystyle a^{T}_{5}$ $\displaystyle=\frac{1}{16\eta_{T}^{4}(f)}.$
(159)

## References

## References

  * Abbott  _et al._ [2019] B. P. Abbott _et al._ (LIGO Scientific Collaboration, Virgo Collaboration), Phys. Rev. D 100, 104036 (2019), arXiv: 1903.04467 .
  * Abbott  _et al._ [2021] R. Abbott _et al._ (LIGO Scientific Collaboration, Virgo Collaboration), Phys. Rev. D 103, 122002 (2021), arXiv: 2010.14529 .
  * Abbott  _et al._ [2022] R. Abbott _et al._ (LIGO Scientific Collaboration, Virgo Colaboration, KAGRA Collaboration), (2022), arXiv: 2112.06861 .
  * Eardley and Lightman [1973] D. M. Eardley and A. P. Lightman, Phys. Rev. D 8, 3308 (1973).
  * Eardley  _et al._ [1973] D. M. Eardley, D. L. Lee, A. P. Lightman, R. V. Wagoner, and C. M. Will, Phys. Rev. Lett. 30, 884 (1973).
  * Newman and Penrose [1962a] E. Newman and R. Penrose, J. Math. Phys. 3, 566 (1962a).
  * Newman and Penrose [1962b] E. Newman and R. Penrose, J. Math. Phys. 4, 998 (1962b).
  * Suvorov [2019] A. G. Suvorov, Phys. Rev D 99, 124026 (2019).
  * Švarc, R. and Pravdová, A. and Miškovský, D. [2023] Švarc, R. and Pravdová, A. and Miškovský, D., Phys. Rev D 107, 024036 (2023).
  * Li  _et al._ [2023] D. Li, P. Wagle, Y. Chen, and N. Yunes, Phys. Rev X 13, 021029 (2023).
  * de Paula  _et al._ [2004] W. L. S. de Paula, O. D. Miranda, and R. M. Marinho, Class. Quantum Gravity 21, 4595 (2004).
  * Alves  _et al._ [2009] M. E. S. Alves, O. D. Miranda, and J. C. N. de Araujo, Phys. Lett. B 679, 401 (2009).
  * Alves  _et al._ [2010] M. E. S. Alves, O. D. Miranda, and J. C. N. de Araujo, Class. Quantum Gravity 27, 145010 (2010).
  * Hohmann [2012] M. Hohmann, Phys. Rev. D 85, 084024 (2012).
  * Myung and Moon [2014] Y. S. Myung and T. Moon, J. Cosmol. Astropart. Phys. 10, 043 (2014).
  * Alves  _et al._ [2016] M. E. S. Alves, P. H. R. S. Moraes, J. C. N. de Araujo, and M. Malheiro, Phys. Rev. D 94, 024032 (2016).
  * Sharif and Siddiqa [2017] M. Sharif and A. Siddiqa, Astroph. and Space Science 362, 226 (2017).
  * Bertolami  _et al._ [2018] O. Bertolami, C. Gomes, and F. S. N. Lobo, The Europ. Physical Journal C 78, 303 (2018).
  * Abedi and Capozziello [2018] H. Abedi and S. Capozziello, The Europ. Physical Journal C 78, 474 (2018).
  * Mebarki [2019] N. Mebarki, in  _Journal of Physics Conference Series_, Journal of Physics Conference Series, Vol. 1269 (2019) p. 012015.
  * Toniato [2019] J. D. Toniato, The Europ. Physical Journal C 79, 680 (2019).
  * Wagle  _et al._ [2019] P. Wagle, A. Saffer, and N. Yunes, Phys. Rev. D 100, 124007 (2019).
  * Haghshenas and Azizi [2020] M. Haghshenas and T. Azizi, International Journal of Modern Physics D 29, 2050004 (2020).
  * Gogoi and Goswami [2020] D. J. Gogoi and U. D. Goswami, The European Physical Journal C 80, 1101 (2020).
  * Dong and Liu [2022] Y. Dong and Y. Liu, Phys. Rev. D 105, 064035 (2022).
  * Liu  _et al._ [2022] Y. Liu, Y. Dong, and Y. Liu, (2022), arXiv: 2206.15333 .
  * Liang  _et al._ [2017] D. Liang, Y. Gong, S. Hou, and Y. Liu, Phys. Rev. D 95, 104034 (2017).
  * Gong and Hou [2018] Y. Gong and S. Hou, Universe 4, 85 (2018).
  * Hou  _et al._ [2018] S. Hou, Y. Gong, , and Y. Liu, The European Physical Journal C 78, 378 (2018).
  * Hyun  _et al._ [2019] Y.-H. Hyun, Y. Kim, and S. Lee, Phys. Rev. D 99, 124002 (2019).
  * Bardeen [1980] J. M. Bardeen, Phys. Rev. D 22, 1882 (1980).
  * Mukhanov  _et al._ [1992] V. F. Mukhanov, H. Feldman, and R. H. Brandenberger, Phys. Rept. 215, 203 (1992).
  * Flanagan and Hughes [2005] É.É. Flanagan and S. A. Hughes, New J. Phys 7, 204 (2005).
  * Moretti  _et al._ [2019] F. Moretti, F. Bombacigno, and G. Montani, Phys. Rev. D 100, 084014 (2019).
  * Dong  _et al._ [2023] Y. Dong, Y. Liu, and Y. Liu, (2023), arXiv: 2305.12516 .
  * Gong  _et al._ [2018a] Y. Gong, S. Hou, D. Liang, and E. Papantonopoulos, Phys. Rev D 97, 084040 (2018a).
  * Gong  _et al._ [2018b] Y. Gong, S. Hou, E. Papantonopoulos, and D. Tzortzis, Phys. Rev D 98, 104017 (2018b).
  * Tachinami  _et al._ [2021] T. Tachinami, S. Tonosaki, and Y. Sendouda, Phys. Rev D 103, 104037 (2021).
  * Alves and Tinto [2011] M. E. S. Alves and M. Tinto, Phys. Rev D 83, 123529 (2011).
  * Tinto and Alves [2010] M. Tinto and M. E. S. Alves, Phys. Rev D 82, 122003 (2010).
  * Chamberlin and Siemens [2012] S. J. Chamberlin and X. Siemens, Phys. Rev D 85, 082001 (2012).
  * Błaut [2012] A. Błaut, Phys. Rev D 85, 043005 (2012).
  * Will [2018] C. M. Will,  _Theory and Experiment in Gravitational Physics_ (Cambridge University Press, 2018).
  * Jaccard  _et al._ [2013] M. Jaccard, M. Maggiore, and E. Mitsou, Phys. Rev D 87, 044017 (2013).
  * Note [1] Notice that our definition of the variable $\Psi$ differs from the definition given by the Ref. [44] by a minus sign, i.e., one should change $\Psi\rightarrow-\Psi$ to compare the equations presented in both articles.
  * Wagoner [1970] R. Wagoner, Phys. Rev. D 1, 3209 (1970).
  * Bergmann [1968] P. Bergmann, Int. J. Theor. Phys. 1, 25 (1968).
  * Horndeski [1974] G. Horndeski, Int. J. Theor. Phys. 10, 363 (1974).
  * De Felice and Tsujikawa [2010] A. De Felice and S. Tsujikawa, Living Rev. Relativity 13, 3 (2010).
  * Yang  _et al._ [2011] L. Yang, C.-C. Li, and C.-Q. Geng, JCAP 08, 029.
  * Yunes and Stein [2011] N. Yunes and L. C. Stein, Phys. Rev D 83, 104002 (2011).
  * Vilhena  _et al._ [2021] S. G. Vilhena, L. G. Medeiros, and R. R. Cuzinatto, Phys. Rev D 104, 084061 (2021).
  * Koop and Finn [2014] M. J. Koop and L. S. Finn, Phys. Rev. D 90, 062002 (2014).
  * Błaut [2019] A. Błaut, Class. Quantum Gravity 36, 055004 (2019).
  * Jenet  _et al._ [2011] F. Jenet, J. Armstrong, and M. Tinto, Phys. Rev D 83, 081301 (2011).