The Role of Longitudinal Polarizations in Horndeski and Macroscopic Gravity: Introducing Gravitational Plasmas

Abstract

We discuss some general and relevant features of longitudinal gravitational modes in Horndeski gravity and their interaction with matter media. Adopting a gauge-invariant formulation, we clarify how massive scalar and vector fields can induce additional transverse and longitudinal excitations, resulting in breathing, vector, and longitudinal polarizations. We review, then, the interaction of standard gravitational waves with a molecular medium, outlining the emergence of effective massive gravitons, induced by the net quadrupole moment due to molecule deformation. Finally, we investigate the interaction of the massive mode in Horndeski gravity with a noncollisional medium, showing that Landau damping phenomenon can occur in the gravitational sector as well. That allows us to introduce the concept of “gravitational plasma”, where inertial forces associated with the background field play the role of cold ions in electromagnetic plasma.

1. Introduction

The recent observation of gravitational waves emitted from black hole and neutron star coalescence [1,2,3,4,5,6,7,8,9], besides offering an important confirmation to weak field predictions of general relativity [10,11], opens new perspectives in the investigation of astrophysical and cosmological phenomena [12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27]. It is clear that the choice of the Einstein–Hilbert action for the gravitational dynamics corresponds just to a criterion of simplicity, ensuring the second-order differential character of the field equations (see Lovelock theorem [28,29]).

The study of more general formulations was originally motivated by the expectation that the singularities emerging in general relativity, such as the initial big-bang singularity or the diverging curvature characterizing the center of Schwarzschild black holes, could be solved in the presence of different geometrical terms in the Lagrangian; for instance, with higher-order contributions [30,31,32] or further contractions [33,34,35,36,37,38,39]. Higher-order terms were also prompted by the request of recovering renormalizability at the quantum level, and quadratic curvature invariants were shown to properly deal with the problem [40,41,42,43,44]. On the other hand, since their inclusion would evidently have implied field equations of a higher order than two, the issue of possible Ostrogradsky instabilities led Horndeski to define the most general scalar tensor theory of gravity, endowed with an additional scalar field beside the metric, devoid of pathologies and preserving the second-order nature of the equations of motion [45]. Recently, the motivation for extended theories of gravity has been also driven by the request to provide a satisfactory explanation for the Universe’s dark components, e.g., dark matter and dark energy [46,47,48,49,50,51,52,53,54,55].

It is important to stress that reformulations of the Einstein–Hilbert action can be carried on not only in the direction of a revised Lagrangian, but also according to a metric-affine perspective, where the affine connection is assumed to be an independent geometrical object with respect to the metric, and whose form has to be dynamically derived by evaluating its very equation of motion [56,57,58,59]. These alternative approaches are usually enriched with a nontrivial geometric structure, where curvature is accompanied by torsion and nonmetricity [56,60,61,62,63,64]—they turned out to be capable of tackling singularity problems [65,66,67,68,69,70,71,72,73] or cosmological observations [39,74,75,76,77] from an original standpoint.

Clearly, such new formulations must be reconciled with observational and experimental constraints, and deviations from general relativity in the Solar System are very strict [78,79]. Nonetheless, increasing interest emerged for the possibility to observe modifications in the detected gravitational waves, such as additional polarization modes [80,81,82,83,84,85,86] or anomalous dispersive and dissipative behaviors [87], hopefully originating from extended dynamical paradigms.

Here, we will concentrate our attention on the morphology and propagation of gravitational waves in the context of Horndeski gravity, in the broader sense of scalar and vector couplings [88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103]. In particular, in this framework, we will analyze the polarizations induced in the geodesic deviation equation by the massive scalar or vector modes, discussing to some extent the limit of vanishing mass. Furthermore, we will study the interaction of gravitational waves with a material medium, having the morphology of a molecular array or plasmalike features. In this respect, we will introduce some new subtle concepts such as “Macroscopic Gravity Theory” [104] or “Gravitational Plasma” [105].

In more detail, regarding the investigation of the gravitational modes in Horndeski gravity, we review and partially extend the gauge-invariant formulation of the linear theory [106] in order to unequivocally identify the physical degrees of freedom. Specifically, we arrive to a clear demonstration that scalar and vector fields are responsible for a superposition of transverse and longitudinal excitations, whose properties depend crucially on the massive nature of the additional polarizations. This overview of the gauge-invariant formulation fixes a specific and quite general phenomenological signature for Horndeski theories, and offers remarkable suggestions about the device that could optimize the detection of such anomalous modes.

In investigating the interaction of gravitational waves with matter media, we first review some well-known results about the attenuation of the signal due, for instance, to the dissipation properties of the traveled medium [107,108], to the redshift in an expanding Universe [109,110,111], or to the interaction with a cosmological neutrino background [112,113,114]. Then, we face the question concerning the propagation of waves in a molecular or a plasmalike medium. In the case of matter medium characterized by molecular substructures [104], as it occurs in astrophysical systems such as galaxies, we show that standard gravity acquires an effective massive behavior, resulting in five degrees of freedom. In this scenario, anomalous, massivelike polarizations appear and their features are discussed in some detail. This phenomenon is very similar to what happens to electromagnetic waves interacting with a plasma: the photon acquires an effective mass and a longitudinal polarization is present inside the plasma [115,116]. Analogously, the interaction of the graviton in a molecular medium can be restated in terms of effective massive states, which propagate with subluminal velocity.

Then, we consider a different medium, made up of free-falling particles in the gravitational field, whose perturbations are now described by Horndeski gravity. We speak of “Gravitational Plasma” when we can choose a local inertial frame, somehow analogous to a “neutralizing field”, which in real plasma is typically associated to cold ions. The idea is that massive modes of the Horndeski formulation can interact with such a system analogously to longitudinal electromagnetic waves in real plasmas of ions and electrons. In other words, we attribute the local gravitational field with the neutralizing role of inertial forces and we analyze the propagation of the massive mode in a self-consistent way, as for Langmuir modes in the electromagnetic case [117]. The settling of the Landau damping for the massive mode is then confirmed via a kinetic description of the medium-wave interaction. This result establishes an important parallelism between electromagnetism in plasma and Horndeski gravity for free-falling particles, i.e., noncollisional “astrophysical” or “cosmological” plasma.

The present review offers an interesting tool to deepen the comprehension of the role of longitudinal gravity waves, as predicted by the extended formulation of gravity, can play when interacting with physical media of different typologies. The review is structured as follows: In Section 2, we discuss the propagation of massive scalar and vector fields in Horndeski theories, outlining the different kinds of polarizations carried out by these additional degrees of freedom and the stresses induced on test particles via the geodesic deviation equation. In Section 3, we review some known results about the interaction of gravitational waves with matter and we eventually establish formal and phenomenological analogies between electromagnetic plasmas and linearized gravitational waves in matter. We firstly introduce the concept of the effective massive graviton in the presence of a molecular stress energy tensor; then, we discuss the idea of the gravitational Landau damping for longitudinal scalar modes. Finally, in Section 4, conclusions are drawn.

2. Longitudinal Degrees of Freedom and Gauge-Invariant Technique

Let us focus our analysis on metric perturbations over Minkowski spacetime, i.e.,

$$g_{\mu \nu }={\eta }_{\mu \nu }+h_{\mu \nu },$$

(1)

with ${\eta }_{\mu \nu }=diag(-1,1,1,1)$ and where we are assuming there exists some reference frame in which $|h_{\mu \nu }|\ll 1$ holds. It follows, then, that in order for $g^{\mu \rho }{g}_{\rho \nu }={\delta }_{\nu }^{\mu }+\mathcal{O}\left(h^2\right)$ to remain consistent, the inverse metric must retain the form

$$g^{\mu \nu }={\eta }^{\mu \nu }-h^{\mu \nu }.$$

(2)

It is worth noting that in the context of modified theories of gravity, the set of ground states actually available can be significantly larger than the “trivial” vacuum Minkowski solution. The enhanced dynamics featuring these theoretical settings, indeed, can in principle sustain additional configurations with respect to the globally flat state described by ${\eta }_{\mu \nu }$. However, since the current detection of gravitational waves by means of ground-based interferometers is ultimately not affected in an appreciable way by background curvature (for truly local effects, such as Newtonian noise, see [118]), we can safely restrict our treatment to gravitational field excitations traveling on Minkowksi background. As we are primarily interested in determining the number and type of additional polarizations, we can as a matter of fact ignore any background configuration featured by a nonvanishing curvature, which we expect could instead affect evolutionary properties of the wave perturbation, such as amplitude and frequency. We refer, for instance, to [119,120,121], where the effect of a cosmological constant is discussed in general relativity and Brans–Dicke theories, outlining effects eventually measurable with the pulsar timing array technique [122,123,124], or [11,111,125,126,127,128,129] where propagation on cosmological background are also considered. Given such premises, we decompose the metric perturbation $h_{\mu \nu }$, which we assume to be vanishing for $r\to \infty$ (asymptotically flatness is recovered), in its irreducible components (see [130,131]):

$$\begin{array}{cc}\hfill h_{00}& =2\varphi ,\hfill \\ \hfill h_{0i}& ={\beta }_i+{\partial }_i\gamma ,\hfill \\ \hfill h_{ij}& =h_{ij}^{TT}+{\displaystyle \frac{1}{3}}H{\delta }_{ij}+{\partial }_{(i}{ϵ}_{j)}+\left({\partial }_i{\partial }_j-{\displaystyle \frac{1}{3}}{\delta }_{ij}\mathsf{\Delta }\right)\lambda ,\hfill \end{array}$$

(3)

where ${\delta }_{ij}$ is the Kronecker delta, $\mathsf{\Delta }\equiv {\partial }_i{\partial }^i$ is the flat Laplacian operator, and symmetrization is defined as $A_{\left(ij\right)}\equiv \frac{1}{2}(A_{ij}+A_{ji})$. It is then evident that the splitting (3) is somehow redundant, given that we defined 16 new functions for describing a symmetric tensor of rank two 1. In order to restore the proper number of components, it is thus necessary to introduce the set of constraints

$$\begin{array}{cc}\hfill {\partial }^i{\beta }_i& =0\hfill \\ \hfill {\partial }^i{h}_{ij}^{TT}& =0\hfill \\ \hfill {\eta }^{ij}{h}_{ij}^{TT}& =0\hfill \\ \hfill {\partial }^i{ϵ}_i& =0,\hfill \end{array}$$

(4)

which, by fixing six additional conditions, allows us to recover the ten independent elements carried by $h_{\mu \nu }$. In general, when matter is also included in the analysis, the metric perturbation contains a combination of gauge and physical degrees of freedom, where the latter encompass, in turn, radiative and nonradiative fields. The first ones represent self-sustained perturbations, which can propagate freely in vacuum, while the second ones are only excited in the presence of matter sources. In general relativity, however, this crucial distinction is partially obscured when dealing with specific gauge choice, for instance, in the Lorentz gauge ${\partial }_{\mu }{h}^{\mu \nu }=0$, since in this case, the equation of motion can be simply rearranged in the form

$$\square h_{\mu \nu }=-2\kappa T_{\mu \nu },$$

(5)

which seems to suggest that all the remaining six components of the metric perturbation could be radiative. Only in vacuum can the residual gauge freedom be further exploited to extract from $h_{\mu \nu }$ the two physical degrees of freedom truly propagating, i.e., the transverse and traceless states encoded in $h_{ij}^{TT}$, corresponding to the well-known cross “×” and plus “+” polarizations. Therefore, given that it seems sensible to reformulate the problem without fixing a priori the system of coordinates—that is to say, in a gauge-invariant way, in order to obtain directly from the equations of the motion the dynamic properties of the theory. We stress, however, that such a statement has to be considered only at the linearized level, since in performing (1), we have already selected a special gauge. Let us consider the infinitesimal coordinate transformation

$$x^{\prime \mu }=x^{\mu }+{\xi }^{\mu }\left(x\right),$$

(6)

which, once we take into account the general transformation rule for the metric

$$g_{\mu \nu }^{\prime }\left(x^{\prime }\right)=\frac{\partial x^{\rho }}{\partial x^{\prime \mu }}\frac{\partial x^{\sigma }}{\partial x^{\prime \nu }}{g}_{\rho \sigma }\left(x\right),$$

(7)

leads to the gauge transformation for the perturbation $h_{\mu \nu }$:

$$h_{\mu \nu }^{\prime }\left(x^{\prime }\right)=h_{\mu \nu }\left(x\right)-2{\partial }_{(\mu }{\xi }_{\nu )}.$$

(8)

It is therefore clear that for (1) to remain consistent, i.e., $|h_{\mu \nu }|\ll 1$, we require $\partial \xi \sim \mathcal{O}\left(h\right)$. In this respect, it is interesting to note that we are not compelled to demand for ${\xi }^{\mu }$ itself to be of the same order of $h_{\mu \nu }$ but only its derivative, which also guarantees that we could consistently solve (6) for $x^{\mu }$. We introduce, then, the combinations 2

$$\begin{array}{cc}\hfill \mathsf{\Phi }& =-\varphi +\dot{\gamma }-{\displaystyle \frac{1}{2}}\ddot{\lambda }\hfill \\ \hfill \mathsf{\Theta }& ={\displaystyle \frac{1}{3}}\left(H-\mathsf{\Delta }\lambda \right)\hfill \\ \hfill {\mathsf{\Xi }}_i& ={\beta }_i-{\displaystyle \frac{1}{2}}{\dot{ϵ}}_i,\hfill \end{array}$$

(9)

where the dot denotes time derivative, which is easy to see to constitute, together with $h_{ij}^{TT}$, the correct number of gauge-invariant quantities. By other words, instead of dealing directly with the entire tensor $h_{\mu \nu }$, we can simply consider the six independent components carried by the two scalars $\mathsf{\Phi }$ and $\mathsf{\Theta }$, the divergenceless vector ${\mathsf{\Xi }}_i$, and the transverse and traceless matrix $h_{ij}^{TT}$. Hence, we proceed by performing an equivalent decomposition on the stress energy tensor, i.e.,

$$\begin{array}{cc}\hfill T_{00}& =\rho ,\hfill \\ \hfill T_{0i}& =S_i+{\partial }_{i}S,\hfill \\ \hfill T_{ij}& ={\sigma }_{ij}+P{\delta }_{ij}+{\partial }_{(i}{\sigma }_{j)}+\left({\partial }_i{\partial }_j-{\displaystyle \frac{1}{3}}{\delta }_{ij}\mathsf{\Delta }\right)\sigma ,\hfill \end{array}$$

(10)

which are accompanied by an analogous set of constraints, namely,

$$\begin{array}{cc}\hfill {\partial }_i{S}^i& =0\hfill \\ \hfill {\partial }^i{\sigma }_{ij}& =0\hfill \\ \hfill {\eta }^{ij}{\sigma }_{ij}& =0\hfill \\ \hfill {\partial }_i{\sigma }^i& =0.\hfill \end{array}$$

(11)

Moreover, by taking into account the linearized conservation law ${\partial }_{\mu }{T}^{\mu \nu }=0$, we can derive an additional set of relations, i.e.,

$$\begin{array}{cc}\hfill \mathsf{\Delta }S& =\dot{\rho },\hfill \\ \hfill \mathsf{\Delta }\sigma & =-{\displaystyle \frac{3}{2}}P+{\displaystyle \frac{3}{2}}\dot{S},\hfill \\ \hfill \mathsf{\Delta }{\sigma }_i& =2{\dot{S}}_i,\hfill \end{array}$$

(12)

which further reduces the number of truly independent components of the stress energy tensor to the set $\{\rho ,P,S_i,{\sigma }_{ij}\}$. Before proceeding to apply the formalism to a concrete scenario of modified gravity, it can be instructive to briefly discuss how it works in standard GR. To this aim, let us expand the to the first order in $h_{\mu \nu }$ the Riemann and Ricci curvature tensor, i.e.,

$$\begin{array}{cc}\hfill R_{\rho \sigma \mu \nu }^{\left(1\right)}& =\frac{1}{2}\left({\partial }_{\sigma }{\partial }_{\mu }{h}_{\rho \nu }+{\partial }_{\rho }{\partial }_{\nu }{h}_{\sigma \mu }-{\partial }_{\sigma }{\partial }_{\nu }{h}_{\rho \mu }-{\partial }_{\rho }{\partial }_{\mu }{h}_{\sigma \nu }\right)\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill R_{\mu \nu }^{\left(1\right)}& =\frac{1}{2}\left({\partial }_{\mu }{\partial }_{\rho }h{}_{\nu }^{\rho }+{\partial }_{\nu }{\partial }_{\rho }h{}_{\mu }^{\rho }-{\partial }_{\mu }{\partial }_{\nu }h-\square h_{\mu \nu }\right),\hfill \end{array}$$

(14)

where the trace h is defined from the Minkowski metric, i.e., $h\equiv {\eta }^{\mu \nu }{h}_{\mu \nu }$, and $\square \equiv {\partial }^{\mu }{\partial }_{\mu }$ denotes the flat d’Alambert operator. It follows from (14) that the Ricci scalar reads as

$$R^{\left(1\right)}={\eta }^{\mu \nu }{R}_{\mu \nu }^{\left(1\right)}={\partial }_{\mu }{\partial }_{\nu }{h}^{\mu \nu }-\square h.$$

(15)

Then, by taking advantage of the decomposition just discussed, we rewrite the Ricci tensor and the Ricci scalar in terms of gauge-invariant quantities, namely,

$$\begin{array}{cc}\hfill R_{00}^{\left(1\right)}& =\mathsf{\Delta }\mathsf{\Phi }-{\displaystyle \frac{3}{2}}\ddot{\mathsf{\Theta }}\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill R_{0i}^{\left(1\right)}& =-{\displaystyle \frac{1}{2}}\mathsf{\Delta }{\mathsf{\Xi }}_i-{\partial }_i\dot{\mathsf{\Theta }}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill R_{ij}^{\left(1\right)}& =-{\partial }_{(i}{\dot{\mathsf{\Xi }}}_{j)}-{\partial }_i{\partial }_j\left(\mathsf{\Phi }+{\displaystyle \frac{1}{2}}\mathsf{\Theta }\right)-{\displaystyle \frac{1}{2}}\square \left({\delta }_{ij}\mathsf{\Theta }+h_{ij}^{TT}\right)\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill R^{\left(1\right)}& =3\ddot{\mathsf{\Theta }}-2\mathsf{\Delta }(\mathsf{\Theta }+\mathsf{\Phi }),\hfill \end{array}$$

(19)

which when plugged into the well-known Einstein equations $G_{\mu \nu }=\kappa T_{\mu \nu }$ and accounting for the conservation law for the stress energy tensor leads us to the system

$$\begin{array}{cc}\hfill \mathsf{\Delta }\mathsf{\Phi }& ={\displaystyle \frac{\kappa }{2}}\left(3P+\rho -3\dot{S}\right)\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \mathsf{\Delta }\mathsf{\Theta }& =-\kappa \rho \hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill \mathsf{\Delta }{\mathsf{\Xi }}_i& =-2\kappa S_i\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill \square h_{ij}^{TT}& =-2\kappa {\sigma }_{ij}.\hfill \end{array}$$

(23)

We see that only the transverse and traceless part $h_{ij}^{TT}$ satisfies a wave equation, while the other components are just gauge degrees of freedom with no radiative behavior as they obey Poisson (Laplace in vacuum) equations. Now, when we deal with alternative theories of gravity, we expect in the linearized Einstein equations an additional contribution depending on combinations of the curvature invariants and possibly some new fields we introduced in our theoretical setting, such as scalar or vector degrees of freedom coupled in different ways to the metric. The idea behind the gauge-invariant formulation is then to seek for a suitable set of variables whose equations could be still rearranged in the form (20)–(23), provided some redefinition of the gauge-invariant components be performed. To be more specific, let us assume our theory be only endowed with an extra scalar degree along with a novel vector field, whose equations of motion can be generically described at the linearized level as

$${\mathcal{E}}^{\phi }=0,{\mathcal{E}}_{\mu }^A=0,$$

(24)

where ${\mathcal{E}}^{\phi },{\mathcal{E}}_{\mu }^A$ enclose differential operators at most order two (we are considering in the broadest sense, the scenario of Horndeski and generalized Proca theories [88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103]), depending on the d’Alembert operator and containing eventually terms stemming from a potential in the Lagrangian. In this case, it seems reasonable to define a new set of gauge-invariant variables given by

$$\begin{array}{cc}\hfill & {\mathsf{\Phi }}^{\left(\phi \right)}\equiv \mathsf{\Phi }+a\phi \hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill & {\mathsf{\Theta }}^{\left(\phi \right)}\equiv \mathsf{\Theta }+b\phi \hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & {\mathsf{\Xi }}_i^{\left(A\right)}\equiv {\mathsf{\Xi }}_i+c{A}_i,\hfill \end{array}$$

(27)

where $a,b,c$ are coefficients related to the structure of the model, such that the following hold:

$$\begin{array}{cc}\hfill \mathsf{\Delta }{\mathsf{\Phi }}^{\left(\phi \right)}& ={\displaystyle \frac{\kappa }{2}}\left(3P+\rho -3\dot{S}\right)\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill \mathsf{\Delta }{\mathsf{\Theta }}^{\left(\phi \right)}& =-\kappa \rho \hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill \mathsf{\Delta }{\mathsf{\Xi }}_i^{\left(A\right)}& =-2\kappa S_i\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill \square h_{ij}^{TT}& =-2\kappa {\sigma }_{ij}.\hfill \end{array}$$

(31)

Thus, even if ${\mathsf{\Phi }}^{\left(\phi \right)},{\mathsf{\Theta }}^{\left(\phi \right)},{\mathsf{\Xi }}_i^{\left(A\right)}$ are still nonradiative, they now actually contain dynamical degrees of freedom according to (24) and this can be demonstrated to be capable of exciting additional polarizations of the gravitational wave spectrum. Let us consider, indeed, a sphere of test masses and let us analyze the effect due to the crossing of a gravitational wave, i.e., let us study the linearized geodesic deviation equation

$$\frac{d^2\delta X^i}{d{t}^2}\simeq -R{}_{0j0}^i{X}_0^j,$$

(32)

where $\delta X^i$ is the displacement induced by the gravitational wave with respect to the initial position $X_0^i$ and we select the comoving frame where the four-velocity is normalized as $U^{\mu }=(-1,0,0,0)$. Straight computation shows that the Riemann components can be locally expressed in terms of gauge-invariant quantities as

$$R_{i0j0}=-{\displaystyle \frac{1}{2}}{\ddot{h}}_{ij}^{TT}+{\partial }_i{\partial }_j\mathsf{\Phi }-{\displaystyle \frac{1}{2}}{\delta }_{ij}\ddot{\mathsf{\Theta }}+{\partial }_{(i}{\dot{\mathsf{\Xi }}}_{j)}.$$

(33)

It is clear that in GR, where the only dynamical part is provided by the purely tensor modes $h_{ij}^{TT}$, we just recover the standard plus and cross polarizations. In this case, however, $\mathsf{\Theta }$, $\mathsf{\Phi }$, and ${\mathsf{\Xi }}_i$ are not static degrees and the Poisson Equations (28)–(30) actually hold for the modified set of gauge-invariant variables ${\mathsf{\Phi }}^{\left(\phi \right)},{\mathsf{\Theta }}^{\left(\phi \right)}$, and ${\mathsf{\Xi }}_i^{\left(A\right)}$. This implies that we can rearrange (33) in the more convenient form

$$R_{i0j0}=\underset{\mathrm{Static}\mathrm{part}}{\underbrace{{\partial }_i{\partial }_j{\mathsf{\Phi }}^{\left(\phi \right)}-{\displaystyle \frac{1}{2}}{\delta }_{ij}{\ddot{\mathsf{\Theta }}}^{\left(\phi \right)}+{\partial }_{(i}{\dot{\mathsf{\Xi }}}_{j)}^{\left(A\right)}}}+\underset{\mathrm{Radiative}\mathrm{part}}{\underbrace{\left(\frac{b}{2}{\delta }_{ij}{\partial }_t^2-a{\partial }_i{\partial }_j\right)\phi -c{\partial }_{(i}{\dot{A}}_{j)}-{\displaystyle \frac{1}{2}}{\ddot{h}}_{ij}^{TT}}},$$

(34)

so that we expect, in general, that the evolution of $\phi$ and A could give rise to novel phenomenological effects besides the well-established tensorial strains. To be more specific, let us focus on the radiative part of the gravitational wave, which we assume to propagate in all its components along the z axes and let us switch off for the sake of clarity the standard tensor modes. The geodesic deviation equation then reads as

$$\begin{array}{cc}\hfill & \frac{d^2\delta X}{d{t}^2}\simeq -\frac{b}{2}{\partial }_t^2\phi X_0+\frac{c}{2}{\partial }_z{\dot{A}}_x{Z}_0\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill & \frac{d^2\delta Y}{d{t}^2}\simeq -\frac{b}{2}{\partial }_t^2\phi Y_0+\frac{c}{2}{\partial }_z{\dot{A}}_y{Z}_0\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill & \frac{d^2\delta Z}{d{t}^2}\simeq \left[\left(-\frac{b}{2}{\partial }_t^2+a{\partial }_z^2\right)\phi +c{\partial }_z{\dot{A}}_z\right]Z_0+\frac{c}{2}\left({\partial }_z{\dot{A}}_x{X}_0+{\partial }_z{\dot{A}}_y{Y}_0\right).\hfill \end{array}$$

(37)

We see that the scalar field $\phi$ is now responsible for a breathing mode in the transverse plane $xy$, i.e., a conformal polarization with no preferred directions in the plane, together with a longitudinal polarization along the direction of propagation of the gravitational wave. The vector field $A_{\mu }$, instead, induces cross polarizations in the planes $xz$ and $yz$, commonly denoted as vector-x and vector-y polarizations, and an additional contribution to the longitudinal mode. We stress, however, that this very type of excitation depends crucially on the mass terms eventually present in (24) (for details concerning the way of restoring the $U\left(1\right)$ gauge-invariance of $A_{\mu }\left(x\right)$ in the presence of a mass term via the Stueckelberg trick, see [100,132,133,134,135]). If we assume, for instance, that (24) have the simple form of Klein–Gordon, respectively, Proca equations, i.e.,

$$(\square -M_{\phi }^2)\phi =0,(\square -M_A^2)A_{\mu }=0,$$

(38)

we see that in the limit $M_A\to 0$, the transversality condition on the vector $A_{\mu }$ compels us to set $A_z=0$, so that its contribution simply reduces to the vector modes. When we perform the same limit for the scalar field, instead, the perturbation along the z axes cancels out only if the identity $2a=b$ holds, so that the massless condition cannot guarantee for itself the absence of a longitudinal mode. We emphasize, however, that strictly speaking, this holds even for the vector polarizations, since both the vector modes are always endowed with strains along the direction of propagation of the gravitational perturbation. We conclude, therefore, that longitudinality, in the broader sense just discussed, is quite a common feature whenever additional polarizations are taken into account. Now, for future convenience (see Section 3), let us consider by way of example a generic Horndeski model with an additional scalar degree of freedom [45,86], described by the action

$$S=\frac{1}{2\kappa }\int d^{4}x\sqrt{-g}\sum _{i=2}^5{L}_i,$$

(39)

where the terms $L_i$ are given by

$$\begin{array}{cc}\hfill & L_2=K(\phi ,X)\hfill \\ \hfill & L_3=-G_3(\phi ,X)\square \phi \hfill \\ \hfill & L_4=G_4(\phi ,X)R+G_{4,X}\left({(\square \phi )}^2-{\phi }_{\mu \nu }{\phi }^{\mu \nu }\right)\hfill \\ \hfill & L_5=G_5(\phi ,X)G_{\mu \nu }{\phi }^{\mu \nu }+\frac{1}{6}{G}_{5,X}\left({(\square \phi )}^3-3\square \phi {\phi }_{\mu \nu }{\phi }^{\mu \nu }+2\phi {}_{\nu }^{\mu }\phi {}_{\rho }^{\nu }\phi {}_{\mu }^{\rho }\right),\hfill \end{array}$$

(40)

and we introduce the notation

$$X\equiv -\frac{1}{2}{\nabla }_{\mu }\phi {\nabla }^{\mu }\phi ,{\phi }_{\mu \nu }\equiv {\nabla }_{\mu }{\nabla }_{\nu }\phi .$$

(41)

Then, by expanding all the fields and taking into account perturbations at the first order, we can obtain the system of linearized equations (in vacuum):

$$\begin{array}{cc}\hfill & G_{\mu \nu }^{\left(1\right)}-\frac{G_{4,\phi }\left(0\right)}{G_4\left(0\right)}\left({\partial }_{\mu }{\partial }_{\nu }-{\eta }_{\mu \nu }\square \right)\varphi =0\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill & \left(\square -M^2\right)\varphi =0,\hfill \end{array}$$

(43)

with the effective mass of the scalar mode given by

$$M^2=-\frac{K_{,\phi \phi }\left(0\right)}{K_{,X}\left(0\right)-2G_{3,\phi }\left(0\right)+\frac{3G_{4,\phi }^2\left(0\right)}{G_4\left(0\right)}},$$

(44)

where the functions are all evaluated in background values, $\phi ={\phi }_0,X=0$. We stress that the observation of the multimessenger signals associated with GW170817-GRB170817A [4,136] put severe constraints on gravitational waves speed propagation in vacuum, restricting the form of the terms actually feasible of (29) (see [22,137,138,139] for more details and implications in cosmological settings). With respect to [105], we are not considering the coupling with the matter (see the discussion in Section 3 or [86]), since we are only interested in the phenomenology of gravitational wave polarizations on the sphere of test masses, i.e., we are neglecting at all the backreaction due to the displacement of the medium particles. It is important to note that the form of (42) and (43) encompass a large variety of modified theories of gravity, such as metric $f\left(R\right)$ models [85,106] or generalized hybrid metric-Palatini gravity [140,141], even if in this very last case, we deal with two additional scalar fields whose interaction can lead to interesting phenomena such as beatings (see [140] for more details). Now, a bit of manipulation reveals that the coefficients $a,b$ take the simple form

$$2a=b=\frac{G_{4,\phi }\left(0\right)}{G_4\left(0\right)},$$

(45)

so that (35)–(37) simplifies to

$$\begin{array}{cc}\hfill & \frac{d^2\delta X}{d{t}^2}\simeq -\frac{b}{2}{\partial }_t^2\phi X_0\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill & \frac{d^2\delta Y}{d{t}^2}\simeq -\frac{b}{2}{\partial }_t^2\phi Y_0\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill & \frac{d^2\delta Z}{d{t}^2}\simeq -\frac{b}{2}\left({\partial }_t^2-{\partial }_z^2\right)\phi Z_0,\hfill \end{array}$$

(48)

which, under the hypothesis of a monochromatic wave perturbation with planar wave front, described by the dispersion relation ${\omega }^2=k^2+M^2$, can be further rearranged as

$$\begin{array}{cc}\hfill & \frac{d^2\delta X}{d{t}^2}\simeq \frac{b{\omega }^2}{2}\phi X_0\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill & \frac{d^2\delta Y}{d{t}^2}\simeq \frac{b{\omega }^2}{2}\phi Y_0\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill & \frac{d^2\delta Z}{d{t}^2}\simeq \frac{b{M}^2}{2}\phi Z_0.\hfill \end{array}$$

(51)

In this form, some peculiar characteristic are easily displayed. Firstly, we note as in the limit $M\to 0$, the transversality is correctly recovered, while the requirement of having no additional perturbation at all implies $b=0$, which when reformulated in terms of the function $G_4(\phi ,X)$ leads us to consider only minimal coupling of the scalar field with the Ricci scalar. Secondly, we observe for $b\ne 0$ that the two stresses are always in phase, with the amplitude of the breathing mode greater than the amplitude of the longitudinal one.

3. Gravitational Waves in Matter

In this section, we address the broad and debated issue of the propagation of gravitational waves within matter. The problem has a long history in the literature and the great number of works on the theme can be divided into two main categories, depending on the approach to the characterization of the material medium. In fact, the latter can be described either by continuous fields or by a large number of particles, whose statistics are depicted by a distribution function. In the first case, we deal with a hydrodynamic approach, and the dynamics can be traced, in principle, with fully analytic methods, whilst in the second, the problem is treated by following kinetic theory prescriptions, resulting in an intrinsically statistic picture. Let us enumerate the most important results achieved by following the hydrodynamic approximation. In [107,108,142,143], the amplitude damping of gravitational perturbations either on an expanding or generic background was analyzed. The main result is that damping occurs only when the propagation in a dissipative fluid, characterized by a definite viscosity, is considered. The same issue is addressed also in [144], where the thorny issue of the gauge freedom pertaining to the wave throughout the propagation in the medium, strictly connected to the count of the physical and fictitious degrees of freedom, is carefully examined. The authors demonstrate that there exists two physical polarization states, as in vacuum, and confirm again the previous results regarding the characteristic time of absorption. An interesting phenomenological application following from these findings is given in [145], where the authors illustrate that observations of gravitational waves signals through LIGO–Virgo interferometers can constrain the viscosity of the cosmological fluid. In [146], a modified dispersion relation and extra modes of polarizations are shown to appear when the medium considered is a spherical cloud (perfect fluid in Schwarzschild metric); the case of pressureless matter (dust) is studied in [147] through the application of a WKB approximation. The authors derive two different dispersion relations: the first is shown to be simply degenerated and corresponding to tensorial perturbations (gravitational waves), whilst the second is characterized by a double degeneracy, describing the dispersion of both density and vorticity perturbations, i.e., novel degrees of freedom arising when the propagation of perturbations in matter is considered.

Further, in the case of kinetic approach to the description of the material medium, a brief review of the main results achieved can be provided. For instance, in [148,149] is considered the interaction of transverse traceless gravitational waves with noncollisional particles on a Minkowski background. In both works, it is recognized that Landau damping [150] arises only if the wave phase velocity is subluminal throughout the propagation within the medium. Moreover, it is shown that such a condition is fulfilled only by including anisotropies in the distribution function describing the equilibrium configuration. The case of collisionless particles on a globally flat FLRW background is considered in [109,151,152]. In all three works, it is assumed that the equilibrium configuration of the medium is described by a Jüttner distribution [153] and the appearance of extra modes of oscillations, corresponding to scalar and vector perturbations, is highlighted. The results in [151,152], derived under the assumption of frequency of the wave much larger than the Hubble parameter, show that no damping is expected. On the contrary, the analysis in [109], pursued without the aforementioned assumption, proves that vector perturbations were strongly suppressed in the early Universe (the damping rate is found to be proportional to ${\eta }^{-4}$, being $\eta$ the comoving time). The role of collisions is examined in [154,155]. Particularly, in [154] is quantified the damping of gravitational waves coming from scattering processes between elementary particles. In [155] is described the arising of a peculiar mechanism: on one hand, collisions macroscopically act as viscosity, causing the appearance of damping, but on the other hand, they tend to suppress anisotropies, which—as already mentioned—are a source of damping. Therefore, the late-time evolution is substantially undamped and only primordial perturbations are significantly suppressed. The case of the interaction of cosmological gravitational waves with neutrinos is considered in [110,112,113,114]. Specifically, in [110,114], neutrinos are assumed to be already decoupled from the cosmological bath; therefore, collisions are ignored. The maximum damping is calculated for perturbations that enter the horizon during the radiation-dominated era and equals roughly $35\%$, independently from cosmological parameters and frequencies of the radiation. Collisions are included in [112] in order to describe the interaction of primordial gravitational waves with neutrinos that are still in equilibrium with the cosmological bath. It is found that the amount of absorbed radiation is proportional to the abundance of neutrinos and is independent from the frequency considered or the details of the collisions. The transition from collisional to collisionless regimes is addressed in [113], where it is shown that, during the transition, the damping is frequency-dependent. In particular, the appearance of a characteristic spectral signature in the region of the nHz is outlined. An analytic treatment of the integro-differential equation obtained in [110] is provided in [156], where asymptotic solutions are evaluated in terms of Bessel functions. In [157], a term that gives account for the spin of the particles is included in Vlasov equation. However, explicit calculations show that the spin contribution does not couple with traceless transverse tensorial metric perturbations. Lastly, in [111] is considered the case of the interaction of gravitational waves with cold dark matter. It is shown that the dominant effect at astrophysical scale is constituted by a frequency-dependent modification in the propagation speed of the wave. Even though it is outlined that, in principle, the power spectrum of primordial gravitational waves carries signatures of the interaction with dark matter, the effects are found to be small, due to the highly nonrelativistic nature of the medium traversed.

A third and far less-pursued approach to the description of the material medium is the so-called molecular picture. In this case, the matter content is imagined as a set of pointlike particles that can be grouped into molecules. The external microscopic field, i.e., the gravitational wave, alters the structure of the molecule, generating an induced quadrupole moment. At a macroscopic level, the average quadrupole moment of the whole medium is then responsible for peculiar modifications in the gravitational wave dynamics, as for instance dispersion or damping. This line of research begins with the work from Szekeres [158], in which is provided the definition of the quadrupole polarization tensor of the molecule in terms of the separation vectors between point-particles and centers of mass of the molecules. In addition to this, it is found that when the Weyl tensor is assumed to be the fundamental dynamical variable, dispersion occurs and no extra polarizations are excited by the propagation within matter. A prosecution of this work is provided in [159], with a particular focus on the static limit, i.e., the modified Newtonian potential affected by the contribution of the average quadrupole moment of the medium. Lastly, in [160] is outlined the possibility of the settling of an amplitude damping, when the molecules are treated on an anti-de-Sitter background.

In the next sections, we will illustrate two works regarding the propagation of gravitational waves in matter. In the first [104], it is shown that gravitational waves from general relativity are endowed with five physical degrees of freedom when propagating in a medium of massive spherical molecules. The analysis is pursued in the theoretical framework of macroscopic gravity, and the constitutive relation between the molecule quadrupole tensor and the external gravitational wave is derived from a simple phenomenological model. The number of radiating degrees of freedom and the dispersion suffered by the macroscopic radiation allow us to establish an effective analogy with the propagation of massive gravitational waves, as described by the Fierz–Pauli theory [161]—the linear approximation of the full theory of massive gravity, which has been introduced only in recent years [162]. In the second work [105], the authors analyzed the interaction between gravitational waves from Horndeski theories and a noncollisional medium of massive particles. The kinetic treatment enforced outlines the presence of Landau damping for the scalar massive mode, when an inequality between theory parameters and physical quantities describing the medium is satisfied.

3.1. Macroscopic Gravity

Longitudinal stresses induced by gravitational waves are not a specific trait of alternative theories of gravity, but can arise also when the propagation of gravitational waves from general relativity within matter is addressed. Indeed, it is a well-known fact that the TT-gauge, which manifests the transversality of the two physical degrees of freedom $h_{+}$ and $h_{\times }$, is achievable only in vacuum. When the general case is considered, instead, the transverse nature of gravitational waves is not ensured and the number of physical degrees of freedom must be deduced from the equation of motion. With respect to the discussion in Section 2, where it was shown that in the matter, we have to deal with spurious gauge degrees of freedom as well; here, we refer to the description of the interaction with the traversed medium in a self-consistent way. In other words, we are interested in expressing the perturbation in the stress-energy tensor due to the presence of the wave through a constitutive relation, i.e., in terms of the components of the wave itself. The resulting equation of motion is closed in the metric perturbation and, in general, exhibits a reduced gauge freedom with respect to the vacuum case, which can lead to fundamental differences in the number of physical degrees of freedom. In this case, the typical Poisson equations ruling the scalar and vector components of the metric perturbation are drastically affected, resulting in novel radiative modes besides the tensor polarizations. We shall consider, in particular, the case of gravitational waves from general relativity traveling in a material medium of pointlike particles grouped into molecules. In [158], it is demonstrated that, under covariant averaging, it is possible to extract the molecular moments from the general matter distribution. The first-order expression of the macroscopic field equation is then

$$〈G_{\mu \nu }〉=\kappa \left(T_{\mu \nu }^{\left(f\right)}+\frac{1}{2}{Q}_{\mu \rho \nu \sigma ,}^{\rho \sigma }\right),$$

(52)

where $〈·〉$ is the covariant averaging applied to the microscopic field equation, $T_{\mu \nu }^{\left(f\right)}$ is the free stress-energy tensor of the pressureless dust composed by the molecules’ centers of mass and $Q_{\mu \rho \nu \sigma }$ is the quadrupole tensor, describing the molecule structure. An expression of the latter is also provided in [158], given in terms of the separation vectors between each pointlike particle and the center of mass of the relative molecule. The external (microscopic) gravitational wave is responsible for a perturbation in the trajectories of the pointlike particles, which can be calculated in terms of the amplitude of the wave itself from the geodesic deviation equation. When the relation is inserted into the expression of the quadrupole tensor, the set of constitutive relations

$$\begin{array}{cc}\hfill Q_{i0j0}& =\frac{1}{2}{ϵ}_g{h}_{ij,00}\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill Q_{0ijk}& =0\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill Q_{ijkl}& =0\hfill \end{array}$$

(55)

is obtained. The parameter ${ϵ}_g$ is the gravitational dielectric constant, quantifying the amount of induced quadrupole moment with respect to the strength of the external wave. Its value is determined by physical parameters characterizing the molecule

$${ϵ}_g={\displaystyle \frac{MN{L}^2}{4{\omega }_0^2}},$$

(56)

with M and L are the molecule mass and radius, respectively; N is the number of molecules per unit volume; and ${\omega }_0$ is the proper frequency of the molecule, defined in terms of the molecule mass density $\rho$ as

$${\omega }_0=\sqrt{{\displaystyle \frac{\kappa \rho }{6}}}.$$

(57)

Let us describe in detail the form of the macroscopic field equation descending from the adoption of the set of constitutive relations (53)–(55). The first aspect that has to be remarked is that the macroscopic field Equation (52) is endowed with a global diffeomorphism invariance, due to the covariant nature of the quadrupole tensor. Therefore, it is possible to require that the averaged metric perturbation $〈{\overline{h}}_{\mu \nu }〉$—which is nothing more than the usual trace-reversed matrix $〈{\overline{h}}_{\mu \nu }〉=〈h_{\mu \nu }〉-\frac{1}{2}{\eta }_{\mu \nu }〈h〉$—must satisfy the Hilbert gauge condition

$${\partial }^{\mu }〈{\overline{h}}_{\mu \nu }〉=0.$$

(58)

As a result, the ten free components pertaining to the general rank two symmetric tensor $〈{\overline{h}}_{\mu \nu }〉$ are reduced to six. In particular, by choosing the reference frame in order to have the direction of propagation of the perturbation along the z axis and writing the components of the perturbation in the Fourier space, the following matrix is obtained

$${\overline{h}}_{\mu \nu }=\left(\begin{array}{cccc}{\left(\frac{k}{\omega }\right)}^2{\overline{h}}_{33}& -\frac{k}{\omega }{\overline{h}}_{13}& -\frac{k}{\omega }{\overline{h}}_{23}& -\frac{k}{\omega }{\overline{h}}_{33}\\ ⋮& {\overline{h}}_{11}& {\overline{h}}_{12}& {\overline{h}}_{13}\\ ⋮& \cdots & {\overline{h}}_{22}& {\overline{h}}_{23}\\ ⋮& \cdots & \cdots & {\overline{h}}_{33}\end{array}\right),$$

(59)

where the symbol $〈·〉$ is dropped, here and in the following, for the sake of convenience. An ulterior gauge transformation that could further reduce the number of free components of the metric perturbation is not, in general, achievable. Indeed, it is a well-known result that, once Hilbert gauge fixing has been enforced, one can cancel out a certain number of components if, and only if, they result to be a solution of d’Alembert equation. Therefore, the analysis of the equations of motion is the way to reveal the number of physical degrees of freedom of the theory. Let us display the $\left(00\right)$, $\left(0i\right)$, and $\left(ij\right)$ components of the field equation, which read, respectively,

$$\begin{array}{cc}\hfill \square {\overline{h}}_{33}& =-4\pi G{ϵ}_g{\partial }_0^4{\overline{h}}_{33}+2\pi G{ϵ}_g{\displaystyle \frac{{\omega }^2}{k^2}}\mathsf{\Delta }{\partial }_0^2\overline{h}\hfill \\ \hfill \square {\overline{h}}_{3i}& =-4\pi G{ϵ}_g{\partial }_0^4{\overline{h}}_{3i}-2\pi G{ϵ}_g{\displaystyle \frac{\omega }{k}}{\partial }_0^3{\partial }_i\overline{h}\hfill \\ \hfill \square {\overline{h}}_{ij}& =-4\pi G{ϵ}_g{\partial }_0^4{\overline{h}}_{ij}+2\pi G{ϵ}_g{\eta }_{ij}{\partial }_0^4\overline{h}.\hfill \end{array}$$

(60)

From these, it is possible to calculate a wave equation for the trace $\overline{h}$, resulting in

$$\square \overline{h}=-2\pi G{ϵ}_g\left({\displaystyle \frac{k^2}{{\omega }^2}}-1\right){\partial }_0^4\overline{h}.$$

(61)

The descending dispersion relations $\omega =\omega \left(k\right)$ are easily calculated

$$\omega \left(k\right)=\pm k\omega \left(k\right)=\pm \sqrt{2}m,$$

(62)

where it has been introduced the constant wavenumber $m={\left(4\pi G{ϵ}_g\right)}^{-\frac{1}{2}}$. This finding shows that the trace can be written as the superposition of a solution of d’Alembert equation, which, as noted before, can be canceled out with a gauge transformation that preserves Hilbert gauge, and a local oscillation, devoid of propagative features (in other words, the group velocity $\frac{d\omega }{dk}$ is strictly null for the second branch). Hence, we are free to set $\overline{h}=0$, bringing the count of degrees of freedom to five. After this fixing, the remaining free components are solutions of the equation

$$\square {\overline{h}}_{\mu \nu }=-{\displaystyle \frac{1}{m^2}}{\partial }_0^4{\overline{h}}_{\mu \nu }.$$

(63)

It is immediate to notice that the fourth-order wave equation obtained for the five macroscopic degrees of freedom is not manifestly covariant: the reason for this fact is that the constitutive relations (53)–(55) have been calculated in a definite reference frame, i.e., the one that is comoving with the center of mass of the molecule. As is well-known, any higher-derivative theory can be affected by the Ostrogradsky instability, resulting in an energy spectrum unbounded from below. A rigorous proof of the global stability of the macroscopic theory should be carried out by showing the singularity of the connected effective Lagrangian; however, in our case, the latter cannot be defined, due to the lack of covariance of the equation of motion. Therefore, we limit our analysis to a phenomenological level and we calculate the dispersion relation coming from (63), which reads

$${\omega }_{\pm }^2\left(k\right)=-{\displaystyle \frac{m^2}{2}}\pm \sqrt{{\displaystyle \frac{m^4}{4}}+m^2{k}^2}.$$

(64)

As it can be easily inferred, the plus-signed branch of the dispersion relation implies pure dispersion of the signal, whereas the minus-signed is related to damping and enhancement phenomena. However, by performing the limit in which the material medium is removed ${ϵ}_g\to 0$ (which is equivalent to $m\to \infty$), one finds that only ${\omega }_{+}\left(k\right)$ goes back to be a vacuum dispersion relation $\omega \left(k\right)=\pm k$. On the contrary, the minus-signed branch has a divergent behavior in the same limit: this suspiciously nonphysical feature allows us to consider the solutions characterized by ${\omega }_{-}\left(k\right)$ as fictitious. Hence, in the following, we will focus only on purely dispersive solutions. We calculate the following group velocity

$$v_g\left(k\right)={\displaystyle \frac{m^{2}k}{2\sqrt{\left(\sqrt{m^2{k}^2+{\displaystyle \frac{m^4}{4}}}-{\displaystyle \frac{m^2}{2}}\right)\left(m^2{k}^2+{\displaystyle \frac{m^4}{4}}\right)}}},$$

(65)

which, in the long-wavelengths limit $k\ll m$, reduces to the approximate form

$$v_g\left(k\right)\simeq 1-{\displaystyle \frac{3k^2}{2m^2}}.$$

(66)

We outline that the group velocity shows the good behavior of being subluminal for all wavenumbers. We are now interested in describing the polarization content of the five macroscopic degrees of freedom: we proceed to exploit the condition $\overline{h}=0$ via gauge fixing

$$\begin{array}{cc}\hfill {\overline{h}}_{11}& ={\overline{h}}_{+}+{\overline{h}}^{*}\hfill \\ \hfill {\overline{h}}_{22}& =-{\overline{h}}_{+}+{\overline{h}}^{*}\hfill \\ \hfill {\overline{h}}_{33}& ={\displaystyle \frac{2{\omega }^2}{k^2-{\omega }^2}}{\overline{h}}^{*}.\hfill \end{array}$$

(67)

As seen in the previous section, we study the linearized geodesic deviation equation in the comoving frame

$$\frac{d^2\delta X^i}{d{t}^2}=-R{}_{0j0}^i{X}_0^j$$

(68)

separately for each free component of ${\overline{h}}_{\mu \nu }$. First, we recover the standard plus and cross tensorial polarizations carried by the components named ${\overline{h}}_{+}$ and ${\overline{h}}_{12}$. Then, we display the sets of equations of motion, which are obtained when the components ${\overline{h}}_{13}$ and ${\overline{h}}_{23}$ are kept non-null:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d^2\delta X}{d{t}^2}}& ={\displaystyle \frac{k^2-{\omega }^2}{2}}{\overline{h}}_{13}{Z}_0\hfill \\ \hfill {\displaystyle \frac{d^2\delta Y}{d{t}^2}}& =0\hfill \\ \hfill {\displaystyle \frac{d^2\delta Z}{d{t}^2}}& ={\displaystyle \frac{k^2-{\omega }^2}{2}}{\overline{h}}_{13}{X}_0\hfill \end{array}$$

(69)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d^2\delta X}{d{t}^2}}& =0\hfill \\ \hfill {\displaystyle \frac{d^2\delta Y}{d{t}^2}}& ={\displaystyle \frac{k^2-{\omega }^2}{2}}{\overline{h}}_{23}{Z}_0\hfill \\ \hfill {\displaystyle \frac{d^2\delta Z}{d{t}^2}}& ={\displaystyle \frac{k^2-{\omega }^2}{2}}{\overline{h}}_{23}{Y}_0.\hfill \end{array}$$

(70)

We recognize that the displacements induced by these components on a sphere of test particles are equivalent to vector-x and vector-y polarizations, respectively (see Section 2). Eventually, when we plug the component ${\overline{h}}^{*}$ into the geodesic deviation equation, we obtain the following

$$\begin{array}{cc}\hfill {\displaystyle \frac{d^2\delta X}{d{t}^2}}& =-{\displaystyle \frac{{\omega }^2}{2}}{\overline{h}}^{*}{X}_0\hfill \\ \hfill {\displaystyle \frac{d^2\delta Y}{d{t}^2}}& =-{\displaystyle \frac{{\omega }^2}{2}}{\overline{h}}^{*}{Y}_0\hfill \\ \hfill {\displaystyle \frac{d^2\delta Z}{d{t}^2}}& =-\left(k^2-{\omega }^2\right){\overline{h}}^{*}{Z}_0.\hfill \end{array}$$

(71)

As it can be deduced from a direct confrontation with (46)–(48), the component ${\overline{h}}^{*}$ is responsible for the simultaneous excitation of a longitudinal polarization superposed to a breathing in the transverse plane. Again, by taking into account the dispersion relation (64), it follows that the two stresses are always in phase, with amplitude of the breathing greater than the amplitude of the longitudinal.

The results collected in this section suggest that it is then possible to build a bridge between gravitational and electromagnetic radiation propagating in matter. Indeed, the five physical degrees of freedom of the macroscopic wave are related to the viable set of states pertaining to a massive spin 2 boson, so it is legitimate to define the concept of effective mass of the graviton, fully analogous with the effective mass acquired by photons when crossing a plasma [115,116]. In the following, we will further exploit this comparison, by introducing the concepts of Landau damping for scalar gravitational waves and neutralizing background for gravitational plasmas.

3.2. Landau Damping for Gravitational Scalar Waves

As previously stated, there is a great number of works in the literature that have investigated the possibility of Landau damping for tensorial gravitational waves, propagating either on a static flat background or over an expanding one. To summarize, it has been found that Landau damping can occur only by including anisotropies in the equilibrium configuration of the medium or by considering a non-null collision term in the Boltzmann equation governing the distribution function time evolution. A second important result achieved by many works on the theme is that a subluminal phase velocity for the signal within the medium results to be a necessary condition for the emergence of the damping phenomenon. With great generality, it has been demonstrated that tensorial massless perturbations on Minkowski cannot satisfy such a request. Is it now important to clarify the writers’ opinion on a technical and conceptual controversy contained in a number of works on the theme: it has been argued that results obtained under the assumption of a flat Minkowski background are not able to satisfactorily describe realistic settings for the fact that, in principle, one should obtain the background metric as the solution of the unperturbed Vlasov problem, i.e., it should be determined by the equilibrium configuration of the medium. Our point of view regarding this issue is based on an analogy with plasma physics. In this field, it is recognized that Landau damping and other collective phenomena are due to the presence of a long-range interaction affecting the great number of electrons contained within a Debye length. This latter concept, in turn, emerges from the screening of the Coulomb electrons’ potential caused by the ions’ distribution. It is therefore evident that, in order to translate the physics of Landau damping in the gravitational sector, it must be provided a mechanism that could mimic the presence of a particle species with negative mass. We argue that it is precisely the possibility of setting a local inertial frame, covering the proximity of a chosen space-time event, which gives the aforementioned phenomenological resemblance between electromagnetic plasma physics and its gravitational counterpart. Indeed, in gravity, one can always make the Christoffel symbols relative to the background almost vanish in a certain neighborhood of a space-time-point, canceling out—for all the particles contained in that region—the mean field generated by the whole medium. We further point out how it is only in such conditions that the assumption of isotropy of the equilibrium configuration of the particles, which is common to many works on the theme, is sufficiently well-grounded, given that any non-null mean field would introduce a privileged direction and inhomogeneities. Hence, in order to study the possibility of Landau damping for gravitational waves, the assumption of a Minkowski background is in all parts acceptable, as long as the considered wavelength of the radiation is smaller than the size of the space-time region to which the treatment can be enforced. Therefore, a sufficient condition for the applicability of this conceptual structure is to look at gravitational perturbations much smaller than the characteristic background scale. As mentioned above, the other fundamental ingredient for the settling of the damping phenomenon is the subluminal phase velocity for modes within the medium. As it emerges from the literature and the previous sections of this work, a subluminal speed of propagation is often related to some sort of massive behavior of the gravitational ripples—be it effective, i.e., emerging at a macroscopic level from the interaction with the medium, or intrinsic, i.e., caused by the peculiar structure of the theory. Hence, it is intriguing to look at the great number of theories of modified gravity that admit the presence of massive degrees of freedom, in order to find out if Landau damping constitutes a viable channel for energy exchange between gravitational radiation and matter.

In this section, we present the analysis of the propagation, in a medium composed by massive particles, of gravitational waves for Horndeski theories, under the assumption of negligible collisions. The general setting of Horndeski theories together with the form of the linearized field equations have been provided in Section 2. Here, we rewrite them in the presence of a nonvanishing stress energy tensor, namely,

$$\begin{array}{cc}\hfill & G_{\mu \nu }^{\left(1\right)}-\frac{G_{4,\phi }\left(0\right)}{G_4\left(0\right)}\left({\partial }_{\mu }{\partial }_{\nu }-{\eta }_{\mu \nu }\square \right)\varphi ={\kappa }^{″}{T}_{\mu \nu }^{\left(1\right)}\hfill \end{array}$$

(72)

$$\begin{array}{cc}\hfill & \left(\square -M^2\right)\varphi ={\kappa }^{\prime }{T}^{\left(1\right)},\hfill \end{array}$$

(73)

where the source terms $T_{\mu \nu }^{\left(1\right)}$ and $T^{\left(1\right)}$ are meant to be perturbations of order $\mathcal{O}\left(h\right)$ induced by the wave itself, whereas the coupling constants ${\kappa }^{\prime }$ and ${\kappa }^{″}$ are expressed in terms of the theory parameters as

$$\begin{array}{cc}\hfill & {\kappa }^{\prime }=\frac{G_4\left(0\right)\kappa }{G_4\left(0\right)\left(K_{,X}\left(0\right)-2G_{3,\phi }\left(0\right)\right)+3G_{4,\phi }^2\left(0\right)}\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill & {\kappa }^{″}=\frac{\kappa }{G_4\left(0\right)}.\hfill \end{array}$$

(75)

In order to separate the scalar degree of freedom from the others, it is useful to introduced the generalized trace-reversed metric perturbation

$${\overline{h}}_{\mu \nu }\equiv h_{\mu \nu }-\frac{1}{2}{\eta }_{\mu \nu }\left(h+2b\varphi \right),$$

(76)

where b is the coefficient introduced in Section 2.

In this manner, it is possible to recast (72) in the familiar form (5), and we can impose on the novel metric variable ${\overline{h}}_{\mu \nu }$ the Lorentz condition, i.e., ${\partial }_{\mu }{\overline{h}}^{\mu \nu }=0$. In the following, we will focus on the purely tensor part of ${\overline{h}}_{\mu \nu }$, and we will neglect its vector component. The scalar part is instead now encoded in the field $\varphi$, which, as we will show, naturally decouples from the tensor perturbations.

The material medium in which gravitational waves are propagating is described by the distribution function $f(\overrightarrow{x},\overrightarrow{p},t)$, whose time evolution is governed by the Vlasov equation

$${\displaystyle \frac{Df}{dt}}={\displaystyle \frac{\partial f}{\partial t}}+{\displaystyle \frac{d{x}^i}{dt}}{\displaystyle \frac{\partial f}{\partial x^i}}+{\displaystyle \frac{d{p}_i}{dt}}{\displaystyle \frac{\partial f}{\partial p_i}}=0$$

(77)

and in terms of which the stress-energy tensor is defined

$$T_{\mu \nu }^{\left(1\right)}={\displaystyle \frac{1}{\sqrt{-g}}}\int d^{3}p{\displaystyle \frac{p_{\mu }{p}_{\nu }}{p^0}}f(\overrightarrow{x},\overrightarrow{p},t).$$

(78)

It must be stressed that the distribution function is assumed to be dependent from the covariant components of the momentum $p_i$ and that the volume element in the space of momenta is $d^{3}p=d{p}_{1}d{p}_{2}d{p}_3$. The external force $\frac{d{p}_i}{dt}$, driving the particles out of the equilibrium configuration, is calculated from the geodesic equation

$${\displaystyle \frac{d{p}^{\mu }}{dt}}+{\Gamma }_{\alpha \beta }^{\mu }{\displaystyle \frac{p^{\alpha }{p}^{\beta }}{p^0}}=0,$$

(79)

where ${\Gamma }_{\alpha \beta }^{\mu }$ are the Christoffel symbols, expressed up to first order in the metric perturbation; $p^0\equiv \sqrt{m^2+g^{ij}{p}_i{p}_j}$ is the particle energy; and m is the particle mass. It results in

$${\displaystyle \frac{d{p}_i}{dt}}={\displaystyle \frac{1}{2p^0}}\left(p_k{p}_l{\displaystyle \frac{\partial {\overline{h}}_{kl}}{\partial x^i}}+b{m}^2{\displaystyle \frac{\partial \varphi }{\partial x^i}}\right).$$

(80)

We assume that gravitational waves start to interact with the matter medium at the fiducial time $t=0$, so that for any negative time, the distribution function is an isotropic solution of the unperturbed equation $f_0\left(p\right)$, where p is the Euclidean modulus of the particle momentum, i.e., $p\equiv \sqrt{{\delta }^{ij}{p}_i{p}_j}$. Then, in order to assure the continuity of the distribution function, we impose that at $t=0$, it is simply given by $f(\overrightarrow{x},\overrightarrow{p},0)=f_0\left(\sqrt{g^{ij}(\overrightarrow{x},0)p_i{p}_j}\right)$.

At first order in the metric perturbation, instead, we have

$$f(\overrightarrow{x},\overrightarrow{p},0)=f_0\left(p\right)-{\displaystyle \frac{f_0^{\prime }\left(p\right)}{2}}\left({\displaystyle \frac{p_i{p}_j}{p}}{\overline{h}}_{ij}(\overrightarrow{x},0)-bp\varphi (\overrightarrow{x},0)\right),$$

(81)

where $f_0^{\prime }\left(p\right)\equiv \frac{d{f}_0}{dp}$. It is then recognizable at this stage a fundamental discrepancy with respect to the electromagnetic analogue: at the initial time, the distribution function acquires a non-null contribution from the discontinuous wave-front, which perturbs the material medium. Then, at any positive time, the presence of the gravitational wave causes an alteration in the particles distribution, which we represent as $\delta f(\overrightarrow{x},\overrightarrow{p},t)$. The linear treatment of the Landau damping phenomenon requires that the order of magnitude of the perturbation of the distribution function with respect to the background values coincides with that of the gravitational wave, i.e., $\frac{\delta f}{f_0}=\mathcal{O}\left(h\right)$. Hence, neglecting all quadratic terms, the linearized Vlasov equation for $\delta f(\overrightarrow{x},\overrightarrow{p},t)$ is obtained

$$\frac{\partial \delta f}{\partial t}+\frac{p^m}{p^0}\frac{\partial \delta f}{\partial x^m}-\frac{f_0^{\prime }\left(p\right)}{2p}\left(p_i{p}_j\frac{\partial {\overline{h}}_{ij}}{\partial t}-\alpha p^2\frac{\partial \varphi }{\partial t}-b{p}^0{p}^m\frac{\partial \varphi }{\partial x^m}\right)=0.$$

(82)

The system composed by (5) (restricted to the solely spatial indices), (73) and (82) results are closed once the source terms for the metric perturbations are expressed in terms of $\delta f(\overrightarrow{x},\overrightarrow{p},t)$, i.e.,

$$\begin{array}{cc}\hfill & T^{\left(1\right)}=-m^2\int d^{3}p{\displaystyle \frac{\delta f(\overrightarrow{x},\overrightarrow{p},t)}{p^0}}\hfill \end{array}$$

(83)

$$\begin{array}{cc}\hfill & T_{ij}^{\left(1\right)}=\int d^{3}p{\displaystyle \frac{p_i{p}_j}{p^0}}\delta f(\overrightarrow{x},\overrightarrow{p},t).\hfill \end{array}$$

(84)

Now, we search for plane wave solutions of this differential problem: particularly, we set the z axis of our reference frame to be coincident with the direction of propagation of the perturbations, so that the field’s dependence from spatial coordinates is restricted to the sole z. Then, as it is usually done in electromagnetic plasma theory, we perform a Fourier transform of parameter k on the spatial coordinate, together with a Laplace transform of parameter s on the time coordinate t. In the Fourier–Laplace space, the differential problem is transformed into an algebraic one, so that solutions for the distribution function perturbation

$$\delta f^{(k,s)}\left(\overrightarrow{p}\right)=\frac{\frac{f_0^{\prime }\left(p\right)}{2p}\left(p_i{p}_j\left(s{\overline{h}}_{ij}^{(k,s)}-{\overline{h}}_{ij}^{\left(k\right)}\left(0\right)\right)-b\left(p^{2}s+ik{p}_3{p}^0\right){\varphi }^{(k,s)}+b{p}^2{\varphi }^{\left(k\right)}\left(0\right)\right)}{s+ik\frac{p_3}{p^0}}$$

(85)

and for the metric fields

$$\begin{array}{cc}\hfill & {\varphi }^{(k,s)}=\frac{s{\varphi }^{\left(k\right)}\left(0\right)+m^2{\kappa }^{\prime }\int d^{3}p\frac{\delta f^{(k,s)}\left(\overrightarrow{p}\right)}{p^0}}{s^2+k^2+M^2}\hfill \end{array}$$

(86)

$$\begin{array}{cc}\hfill & {\overline{h}}_{ij}^{(k,s)}=\frac{s{\overline{h}}_{ij}^{\left(k\right)}\left(0\right)+2{\kappa }^{″}\int d^{3}p\frac{p_i{p}_j}{p^0}\delta f^{(k,s)}\left(\overrightarrow{p}\right)}{s^2+k^2}\hfill \end{array}$$

(87)

are readily obtained. The Fourier and Fourier–Laplace projections are denoted by ${\varphi }^{\left(k\right)}\left(t\right)$ and ${\varphi }^{(k,s)}$, respectively, and analogously for the other fields. Moreover, without loss of generality, we have set ${\dot{\varphi }}^{\left(k\right)}\left(0\right)$, ${\dot{\overline{h}}}_{ij}^{\left(k\right)}\left(0\right)$, and $\delta f^{\left(k\right)}\left(0\right)$ vanishing. Once the expression for the distribution function perturbation (85) is inserted into (86) and (87), it can be shown that scalar and tensor perturbations within the medium are naturally decoupled, i.e.,

$$\begin{array}{cc}\hfill & {\varphi }^{(k,\omega )}={\displaystyle \frac{\left(-i\omega +ib\pi m^2{\kappa }^{\prime }\int d\rho d{p}_3\rho {\displaystyle \frac{f_0^{\prime }\left(p\right)p}{p^0\omega -k{p}_3}}\right){\varphi }^{\left(k\right)}\left(0\right)}{(k^2+M^2-{\omega }^2){ϵ}^{\left(\varphi \right)}(k,\omega )}}\hfill \end{array}$$

(88)

$$\begin{array}{cc}\hfill & {\overline{h}}_{ij}^{(k,\omega )}={\displaystyle \frac{\left(-i\omega -\frac{i\pi {\kappa }^{\prime \prime }}{2}\int d\rho d{p}_3{\rho }^5{\displaystyle \frac{f_0^{\prime }\left(p\right)}{p(p^0\omega -k{p}_3)}}\right){\overline{h}}_{ij}^{\left(k\right)}\left(0\right)}{(k^2-{\omega }^2){ϵ}^{\left(h\right)}(k,\omega )}},\hfill \end{array}$$

(89)

where the integrals have been converted into cylindrical coordinates by means of the substitution $\left(p_1,p_2\right)\to \left(\rho ,\theta \right)$, with ${\rho }^2=p_1^2+p_2^2$ and $tan\left(\theta \right)=\frac{p_2}{p_1}$, and we have redefined $s\equiv -i\omega$. We have also introduced the complex dielectric functions

$${ϵ}^{\left(\varphi \right)}(k,\omega )=1+{\displaystyle \frac{b\pi m^2{\kappa }^{\prime }}{k^2+M^2-{\omega }^2}}\int d\rho d{p}_3\rho {\displaystyle \frac{f_0^{\prime }\left(p\right)}{p}}{\displaystyle \frac{p^2\omega -k{p}^0{p}_3}{p^0\omega -k{p}_3}}$$

(90)

$${ϵ}^{\left(h\right)}(k,\omega )=1-{\displaystyle \frac{\pi {\kappa }^{\prime }}{2(k^2-{\omega }^2)}}\int d\rho d{p}_3{\rho }^5{\displaystyle \frac{f_0^{\prime }\left(p\right)}{p}}{\displaystyle \frac{\omega }{p^0\omega -k{p}_3}},$$

(91)

whose roots are connected with the allowed oscillation modes within the medium. Indeed, as is well-known from the Laplace transform theory, the inverse Laplace transform is dominated, for long times, by the contribution stemming from the poles in the complex plane of (86) and (87), i.e., the zeros of the dielectric functions given above. Specifically, in order to obtain the full spectrum of solutions, one should find, for each value of the wavenumber k, the intersection between the curves ${\omega }_r\left(k\right)$ and ${\omega }_i\left(k\right)$ along which the dielectric function vanishes 3. However, when the weak damping scenario $|{\omega }_r|\gg |{\omega }_i|$ is investigated, approximate solutions for the dispersion relation ${\omega }_r\left(k\right)$ and the damping coefficient ${\omega }_i\left(k\right)$ are easily obtained from

$${ϵ}_r^{(\varphi ,h)}(k,{\omega }_r)=0$$

(92)

and

$${\omega }_i=-{\left.\frac{{ϵ}_i^{(\varphi ,h)}(k,\omega )}{\frac{\partial {ϵ}_r^{(\varphi ,h)}(k,\omega )}{\partial \omega }}\right|}_{\omega ={\omega }_r}.$$

(93)

In particular, we observe that a null imaginary part of the dielectric function corresponds to a purely real angular frequency, leading to the absence of any damping effect.

On the other hand, by the inspection of (90) and (91), we realize that the only imaginary contribution comes from the integration around the Landau pole ${\left(p^0\omega -k{p}_3\right)}^{-1}$. Thus, it is clearly demonstrated that a necessary condition for the presence of damping is a subluminal phase velocity $v_p\equiv \frac{{\omega }_r}{k}<1$; so, that Landau pole is included in the path of integration. Now, we specialize our analysis to the case in which the equilibrium configuration of the medium is described by a Jüttner distribution

$$f_0\left(p\right)={\displaystyle \frac{n}{4\pi m^2\mathsf{\Theta }{K}_2\left(\frac{m}{\mathsf{\Theta }}\right)}}{e}^{-\frac{\sqrt{m^2+p^2}}{\mathsf{\Theta }}},$$

(94)

where n indicates the number of particles per unit volume, $\mathsf{\Theta }$ represents the temperature of the medium expressed in units of the Boltzmann constant $k_B$, and $K_2\left(·\right)$ is the modified Bessel function of the second kind. We proceed by inserting the Jüttner distribution into the expression of the dielectric functions (90) and (91): we assume that the phase velocity is much greater than the thermal velocity of the medium; hence, $\frac{\omega }{k}\gg \frac{p_3}{p^0}$ for most of the particles. Under this hypothesis, it is possible to expand the denominators in convergent power series and integrate term by term. Equating to zero, the expressions obtained yield the dispersion relations for the gravitational radiation within the medium. For the tensor modes, it is found that

$${\omega }_r^2\left(k\right)=\frac{k^2+12{\omega }_h^2\frac{x-\gamma \left(x\right)}{x^2}+\sqrt{{\left(k^2+12{\omega }_h^2\frac{x-\gamma \left(x\right)}{x^2}\right)}^2+48k^2{\omega }_h^2\frac{\gamma \left(x\right)}{x^2}}}{2},$$

(95)

with $x\equiv \frac{m}{\mathsf{\Theta }}$ being the ratio between rest and thermal energy, and $\gamma \left(x\right)\equiv \frac{K_1\left(x\right)}{K_2\left(x\right)}$ and ${\omega }_h^2\equiv \frac{{\kappa }^{″}nm}{6}$ being the squared proper frequency of the medium for tensor modes (see (57) for a comparison). A simple calculation shows that the phase velocity for tensor gravitational waves is always greater than unity throughout the propagation in the material medium, hence, resulting in a null imaginary part for the angular frequency and in the absence of damping, in agreement with the literature on the theme. Conversely, repeating the same steps for the dielectric function relative to the scalar mode yields

$${\omega }_r^2\left(k\right)={\displaystyle \frac{k^2+M^2-9\gamma {\omega }_0^2+\sqrt{{\left(k^2+M^2-9\gamma {\omega }_0^2\right)}^2+12\gamma {\omega }_0^2{k}^2}}{2}},$$

(96)

where, in this case, ${\omega }_0$ represents the proper frequency of the medium for scalar waves, i.e., ${\omega }_0^2\equiv \frac{b{\kappa }^{\prime }nm}{6}$. Now, in order to have a subluminal phase velocity, it results that the following inequality must hold

$$M^2<6\gamma {\omega }_0^2.$$

(97)

This condition represents a distinctive trait of the gravitational Landau damping, which is absent in the electromagnetic counterpart. By relating the mass of the scalar wave, which depends ultimately on the parameters of the gravitational theory taken into account, to quantities describing physical features of the ensemble of particles, it selects media that are able to damp massive modes from noninteracting ones. We point out, moreover, that the constraints on the speed of propagation of gravitational waves derived in [4,136] are translated in our case to a maximum allowed value for the mass M of the scalar degree of freedom. However, such a result does not represent a limitation for our findings, considering that for decreasing values of the mass of the scalar mode, the set of material media able to induce damping becomes actually larger, as it can be easily observed from (97). In addition to this, it is worth outlining that the limit of vanishing mass does not correspond to general relativity, since, in this case, we actually deal with a Horndeski theory endowed with a purely transverse scalar wave, which acts as a breathing in the orthogonal plane. In this respect, the pursued analysis remains strictly valid: for such a massless scalar mode, indeed, Landau damping always occurs, as (97) is identically satisfied by any medium. Another interesting aspect of inequality (97) is disclosed if we apply our model to cosmology. Indeed, it can be observed that the quantity $\gamma {\omega }_0^2$ is always decreasing for decreasing redshift. Hence, if we assume that the mass of the scalar gravitational wave M does not change with time, then for each specific model of material medium there exists a single redshift (at least negative) below which inequality (97) is no longer satisfied, so that the medium considered becomes unable to damp the scalar radiation. We evaluate the imaginary part of the dielectric function (90) by exploiting the residue theorem; then, we obtain, from (93), the following form of the damping rate

$${\omega }_i\left(k\right)=-\frac{\pi x}{4k{K}_1\left(x\right)}\frac{{\omega }_r^4(k^2+M^2-{\omega }_r^2)e^{-\frac{x}{\sqrt{1-\frac{{\omega }_r^2}{k^2}}}}}{3{\omega }_r^4-2{\omega }_r^2+(k^2+M^2)}.$$

(98)

When the dispersion relation (96) is inserted into this expression it can be easily shown that ${\omega }_i$ results negative for any value of the wavenumber. This is certainly not a surprising fact, as we know from the electromagnetic Landau damping theory that the sign of the imaginary part of the angular frequency coincides with the sign of the derivative of the background distribution function considered at a velocity equal to the phase velocity of the signal. Therefore, given that we have pursued our analysis by assuming a phase velocity of the wave in the fast-tail region of the distribution, an always-negative damping rate is naturally obtained. However, the aforementioned ansatz is not feasible when highly relativistic media are considered: in fact, it can be shown that the Jüttner distribution is characterized by a thermal velocity very close to the speed of light when the value of the parameter x drops below unity. In such cases, it is required a numerical treatment of the integral contained in the expression of the dielectric function. A simple integration with rectangles and the expansion of the expression obtained in terms of the parameter $\delta \equiv 1-\frac{{\omega }_r}{k}$, assumed to be small and positive up to third order, yields a family of curves 4 ${\overline{\omega }}_i={\overline{\omega }}_i\left(\overline{k};x,\overline{m}\right)$ shows relevant similarities with (98), such as, for instance, being always negative with a minimum located roughly around $\overline{k}\simeq \overline{M}$. We report the graphics of the curves obtained by varying the parameter x while keeping the normalized mass fixed $\overline{M}=1$. In particular, in Figure 1, it can be noticed that the value of the maximum damping grows with x until the value $x\simeq 5$ is reached. Then, as is clearly shown in Figure 2, the effect rapidly vanishes as soon as x becomes larger than a few tens. However, in both figures, we observe that damping occurs only in a narrow region of $\overline{k}$: in fact, the damping coefficient results are negligible when wavenumbers slightly larger than $\overline{M}$ are considered. This finding suggests some sort of resemblance with spectral lines, giving rise to the introduction of the concept of gravitational-wave spectroscopy. Indeed, by estimating the amount of damping suffered by perturbations of a certain size and by assuming a definite model of modified gravity, it is possible, in principle, to measure the physical properties of the medium traversed by the gravitational radiation. Conversely, by taking into account a precise model of material medium—i.e., by setting the mass of the particles, their volume density, and the temperature—it is possible to give a quantitative estimate of the damping suffered by scalar waves of a certain spatial size. For instance, considering the interaction between the cosmological dark matter medium and scalar radiation at a redshift $z\approx 2000$, we calculate from our formulae an estimated maximum damping of roughly ${\omega }_i\lesssim {10}^{-16}\mathrm{Hz}$ (see [105] for details). It is interesting to notice that in [163], a similar value for the damping rate was directly inferred from observational data, i.e., ${\omega }_i\lesssim {10}^{-17}\mathrm{Hz}$. The interpretation of such a measurement requires, in the context of general relativity, the introduction of a nonzero viscosity for the medium traversed (see Section 3), so that the measurement of the damping allows the indirect estimate of the viscosity of the Universe. On the other hand, given the value of damping rate obtained in our model, the same data analysis can be reinterpreted within the broader theoretical scenario of modified theories of gravity as the result of kinetic damping suffered by scalar radiation in the interaction with viable models of dark matter. Another important occurrence caused by the possibility of Landau damping for scalar gravitational waves is that the interaction between radiation and material medium tends, on average, to heat up the latter. Then, if we consider scalar modes propagating into the cosmological dark matter medium in the early Universe, we expect that the energy exchange allowed by the phenomenon we have described should result in a global heating of dark matter. This latter fact could have significant observational implications, such as a modification in the time for a structure formation through Jeans instability. However, the linear treatment enforced in our analysis is not sufficient to rigorously quantify the energy exchange rate between radiation and material medium because, in this case, the role of the dominant nonlinear terms in Vlasov equation becomes relevant and the theoretical analysis must be deeply revised.

4. Concluding Remarks

In this review, we performed an interesting analogy between properties of photons in electromagnetic plasmas and longitudinal modes of gravitational waves in modified theories of gravity. In particular, we showed how longitudinal excitations are, in general, related to massive states of the graviton, both when they are carried by additional degrees of freedom explicitly introduced in the Lagrangian and when they emerge at the effective level as the result of some averaging procedure of the matter source. This motivated us to trace a comparison with the behavior of photons in ordinary plasmas, where it is a well-established fact that the interaction of electromagnetic waves with the ions can lead to the appearance of longitudinal modes, attributable to the presence of effective massive photons in the medium.

Firstly, we offered a rather exhaustive discussion about the nature and the propagation of massive scalar and vector modes in Horndeski gravity by presenting a gauge-invariant formulation; this permitted us to unequivocally identify the radiating degrees of freedom, thus, allowing us to study the phenomenological signature of the polarizations via the geodesic deviation equation. In this respect, besides the ordinary tensor modes, we clearly identified scalar and the vector components, responsible for breathing, longitudinal, and vector polarizations. This very general conclusion provided a precise hint towards the direction that the new generation of interferometer devices could be pushed in order to have a chance to see modified gravity physics.

We focused, then, on the interaction between the gravitational radiation and the traveled medium, and determined three different scenarios characterized by peculiar phenomenology. We first discussed the case of a standard gravitational wave propagating in the presence of matter, with the aim of elucidating how the damping of the amplitude can settle down due to the dissipation properties of the medium, or to the redshift associated to the Universe expansion (particular attention was dedicated to the gravitational wave–neutrino thermal background interaction). In this respect, we outlined that an ordinary gravitational wave, propagating inside an homogeneous and isotropic medium, is not affected by damping mechanism in a tangible manner. We revised the subtle question concerning the propagation of ordinary gravitational waves within molecular matter, where, due to the presence of bounded substructures, we are forced to abandon a continuum representation for the medium. It was immediate to recognize that such a situation is often well-implemented in astrophysical and cosmological systems, in which different spatial scales of gravitationally bounded systems are naturally present. This approach was summarized by the construction of a Macroscopic Theory of Gravity, according to which the deformation of the molecular medium can affect the gravitational radiation and give rise at the effective level to a massive graviton. Indeed, in close analogy with the propagation of electromagnetic waves within a real plasma, the molecule backreaction can be described by the emergence of additional modes able to induce longitudinal stresses on test particles. However, it should be remarked that these extra polarizations are detectable only at the macroscopic spatial scale characterizing the averaging process that defines the macroscopic theory.

Finally, we discussed the most recent and maybe innovative phenomenon, concerning the Landau damping suffered by the massive scalar mode in generic Horndeski theories. We revised in detail how, even when the matter medium in which the wave propagates is homogeneous and ideal, the scalar mode is suppressed, outlining that the settling of this very specific mechanism depends crucially on a phenomenological inequality relating to the thermodynamic properties of the medium and the mass of the scalar mode. The analysis followed a kinetic approach in which a noncollisional medium, described by the Vlasov equation, is dynamically coupled to the scalar degree of freedom of Horndeski gravity, and the resulting interaction is expressed in terms of the (gravitational) dielectric properties of the medium (we just considered a linear backreaction of the wave on the particles).

This result opens a new interesting perspective on the generality of such a Landau damping within the context of modified gravity approaches, where further scalar and vector degrees of freedom are taken into account. Indeed, as the analysis carried out in Section 3.2 clearly shows, the longitudinal character of a specific mode is not a necessary condition for the settling of the damping, which rather turns out to depend on the peculiar dispersion relation of the mode in the medium traversed. Landau damping occurs only if the latter allows for a subluminal phase velocity of the considered radiating degree of freedom, at least in some regions of the wavelengths spectrum. In this respect, an intriguing possibility is offered by further extensions of Horndeski theories, such as the teleparallel formulation [164,165] or the Palatini generalization [166], exhibiting a large variety of additional polarizations, both massive and massless, which could be, in principle, subject to kinetic damping phenomena.

We conclude by stressing that the achievement of this very result, i.e., the gravitational Landau damping, was based on a very intriguing point of view, according to which the concept of gravitational plasma can be introduced. A crucial step in this direction must be considered the identification of inertial forces as the neutralizing background, analogous to the role played by cold ions in a electromagnetic plasma. This idea, implicitly formulated in [167], permitted [105] to study the self-consistent oscillations inside a gravitational medium—that is, the gravitational equivalent to Langmuir waves in standard plasmas.