Gauge-Invariant Gravitational Wave Polarization in Metric f(R) Gravity with Cosmological Implications

1. Introduction

The study of gravitational-wave (GW) polarizations provides a powerful way to distinguish General Relativity (GR) from alternative theories of gravity. In the classic classification by Eardley, Lee, and Lightman using the Newman–Penrose (NP) formalism [1], a general metric theory of gravity may contain up to six possible GW polarization states. In GR only two of these, the familiar tensor “plus” (⊕) and “cross” (⊗) modes, are present, corresponding to the two radiative degrees of freedom (DoF) of the metric. Extensions of GR often introduce additional scalar and/or vectorial modes whose presence modifies the relative displacement of freely-falling test particles.

A particularly well-known example is Brans–Dicke theory [2], in which the scalar field gives rise to an additional transverse breathing mode. More generally, recent analyses using both the NP formalism and the irreducible (3+1) decomposition [3,4] have confirmed that the number of NP polarization states does not necessarily coincide with the number of radiative DoF in a theory. This mismatch appears naturally in scalar–tensor theories and in metric $f\left(R\right)$ gravity, where the Ricci scalar perturbation introduces a new massive scalar propagating mode that obeys a Klein–Gordon equation with an effective mass $m_{\psi }$ [5,6,7]. For massless propagation, this scalar mode produces a single independent polarization: a transverse breathing distortion of a ring of test particles. When the scalar mode is massive ($m_{\psi }\ne 0$), however, the scalar sector induces a mixture of transverse breathing and longitudinal motion along the propagation direction; these two NP amplitudes are not independent but jointly encode a single scalar radiative DoF. Thus, metric $f\left(R\right)$ gravity contains three radiative DoF but may exhibit up to four NP polarization amplitudes in the massive case.

It is important to note recent discussions regarding the interpretation of NP quantities in theories containing massive modes. For gravitational waves whose group velocity differs from the speed of light, some NP components that vanish in GR no longer vanish identically, raising subtleties concerning the mapping between NP scalars and physical polarizations [1]. This occurs because nonluminal propagation renders the wave vector non-null, so the standard NP null tetrad cannot be aligned with the direction of propagation, and the usual decoupling between radiative and non-radiative components breaks down [3,8]. These issues do not invalidate the NP approach but motivate complementary gauge-invariant formalisms.

Several modified gravity models exhibit similar features. In modified Gauss–Bonnet gravity $f\left(G\right)$, for example, tensor waves propagate as in GR while an additional massive scalar mode appears [9]. In massive gravity theories studied via Bardeen variables [10], a normally non-radiative scalar mode becomes dynamical, constituting the helicity-0 component of a massive graviton. By contrast, quadratic theories such as Einstein–Dilaton–Gauss–Bonnet and dynamical Chern–Simons gravity can display the same polarization content as GR in their linearized limit [11]. Extensions involving explicit matter–geometry couplings, including $f(R,T)$ and $f(R,T^{\varphi })$ gravity [12], further illustrate the subtlety of polarization counting. Although $f(R,T)$ and $f\left(R\right)$ gravity share the same far-field polarization structure in vacuum ($T=0$), their source-dependent dynamics differ. The $f(R,T^{\varphi })$ theory, involving an additional scalar $\varphi$, leads to distinguishable polarization patterns even in vacuum. These analyses also emphasize that the scalar mode $\psi$ belonging to $f\left(R\right)$ gravity must be clearly distinguished from any additional matter scalar $\varphi$.

Gravitational radiation in linearized metric $f\left(R\right)$ gravity has been studied in both power-series models

$$f\left(R\right)=\sum _{n=0}^N{a}_n{R}^n,$$

(1)

and in specific subclasses such as the Starobinsky model. Solar-system tests show that metric $f\left(R\right)$ gravity reproduces light deflection with the same post-Newtonian parameter $\Gamma =1$ as GR, although perihelion precession can differ [13]. Studies of waveform phases in extreme-mass-ratio inspirals suggest that deviations from GR may be detectable in some regimes. Polarization analyses performed in both $f\left(R\right)$ and Horndeski theories [8] highlight the challenge of detecting longitudinal scalar modes with laser interferometers, whereas Pulsar Timing Arrays may offer greater sensitivity. The mixed longitudinal–breathing nature of the massive scalar propagating mode has been explicitly confirmed using the geodesic deviation equation [14]. Additional applications of $f\left(R\right)$ gravity include the study of gravitational radiation from white dwarfs with sub- and super-Chandrasekhar masses, where all relevant polarization amplitudes were estimated using Green-function methods [15].

A de Sitter background is particularly well motivated for analyzing the propagation of tensor and scalar modes in metric $f\left(R\right)$ gravity. It provides an excellent approximation to late-time cosmic acceleration driven by dark energy and also captures the quasi-exponential “slow-roll” inflationary phase in the early Universe. (Here “slow-roll” refers to the regime in which the inflaton’s kinetic energy remains small compared to its potential energy, yielding an almost constant Hubble parameter.) Background curvature affects dispersion relations, mode mixing, and asymptotic behavior of gravitational waves, which motivates studying the massive scalar mode and tensor modes directly on de Sitter space [16,17]. Because GR with a cosmological constant supports only two tensor polarizations, de Sitter space provides a clean setting for isolating any additional modes arising from $f\left(R\right)$ gravity and for making cosmological links between inflationary physics and late-time acceleration.

The structure of this paper is as follows. In Section 2 we derive the field equations of metric $f\left(R\right)$ gravity on a de Sitter background and obtain the Klein–Gordon equation for the extra scalar mode with mass $m_{\psi }$. Section 3 develops the perturbation of the Ricci tensor, and the resulting linearized field equations. Section 4 performs the (3+1) irreducible decomposition into scalar, vector, and tensor sectors. In Section 5 we specialize to the model $f\left(R\right)=R+\alpha R^2-2\Lambda$ and present its explicit linearized dynamics. Section 6 identifies the gauge-invariant Bardeen variables and derives the physical polarization content. Finally, in Section 7 we verify these results using the geodesic deviation equation, demonstrating that the obtained polarization modes are physically realized in the relative acceleration of freely falling test particles.

2. Field Equations of Metric $\mathit{f}\left(\mathit{R}\right)$ Gravity on a de Sitter Background

In metric $f\left(R\right)$ gravity, the Einstein–Hilbert Lagrangian density R is replaced by a general function $f\left(R\right)$,

$$S={\displaystyle \frac{1}{2\kappa }}\int d^{4}x\sqrt{-g}f\left(R\right)+S_m[g_{\mu \nu },{\Psi }_m],$$

(2)

where $S_m$ is the matter action, ${\Psi }_m$ collectively denotes the matter fields, and

$$\kappa \equiv 8\pi G=2M_{\mathrm{Pl}}^{-2}$$

(3)

in terms of the reduced Planck mass $M_{\mathrm{Pl}}$ (we use units $c=1$). Varying (2) with respect to the metric and following, e.g., [18], one obtains the metric $f\left(R\right)$ field equations

$$f^{\prime }\left(R\right)R_{\mu \nu }-{\displaystyle \frac{1}{2}}{g}_{\mu \nu }f\left(R\right)+\left(g_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu }\right)f^{\prime }\left(R\right)=\kappa T_{\mu \nu },$$

(4)

where $f^{\prime }\left(R\right)\equiv df/dR$, $\square \equiv g^{\rho \sigma }{\nabla }_{\rho }{\nabla }_{\sigma }$, and

$$T_{\mu \nu }=-{\displaystyle \frac{2}{\sqrt{-g}}}{\displaystyle \frac{\delta \left(\sqrt{-g}{S}_m\right)}{\delta g^{\mu \nu }}}$$

(5)

is the matter stress–energy tensor. In vacuum we set $T_{\mu \nu }=0$ and (4) reduces to

$$f^{\prime }\left(R\right)R_{\mu \nu }-{\displaystyle \frac{1}{2}}{g}_{\mu \nu }f\left(R\right)+\left(g_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu }\right)f^{\prime }\left(R\right)=0.$$

(6)

Taking the trace of (4) yields the scalar (trace) equation

$$3\square f^{\prime }\left(R\right)+R{f}^{\prime }\left(R\right)-2f\left(R\right)=\kappa T,$$

(7)

where $T\equiv g^{\mu \nu }{T}_{\mu \nu }$. In vacuum,

$$3\square f^{\prime }\left(R\right)+R{f}^{\prime }\left(R\right)-2f\left(R\right)=0.$$

(8)

The trace equation will be the starting point for identifying the massive scalar propagating mode (the “scalaron”) and its effective mass.

2.1. $f\left(R\right)$ Gravity, Scalar–Tensor form, and the Chameleon Mechanism

The field equations (4) contain higher derivatives of the metric through $\square f^{\prime }\left(R\right)$ and ${\nabla }_{\mu }{\nabla }_{\nu }{f}^{\prime }\left(R\right)$. A useful way to expose the extra scalar degree of freedom and to analyze screening—namely, the suppression of scalar-mediated fifth forces in high-density environments via an environment-dependent effective mass—is to recast metric $f\left(R\right)$ gravity as a scalar–tensor theory via a conformal transformation (see, e.g., [6,7]).

Define

$$F\left(R\right)\equiv f^{\prime }\left(R\right)>0,F\left(R\right)=exp\left(-{\displaystyle \frac{2\beta \varphi }{M_{\mathrm{Pl}}}}\right),$$

(9)

where $\varphi$ is a scalar field. In metric $f\left(R\right)$ gravity the coupling parameter is fixed to $\beta =1/\sqrt{6}$, reflecting the universal strength of the scalar coupling to matter. We then introduce the conformal transformation

$${\overline{g}}_{\mu \nu }=F\left(R\right)g_{\mu \nu }=exp\left(-{\displaystyle \frac{2\beta \varphi }{M_{\mathrm{Pl}}}}\right)g_{\mu \nu },$$

(10)

which maps the Jordan frame metric $g_{\mu \nu }$ to the Einstein frame metric ${\overline{g}}_{\mu \nu }$. In the Jordan frame, matter is minimally coupled to $g_{\mu \nu }$ and freely falling test particles follow geodesics of $g_{\mu \nu }$. In the Einstein frame, the gravitational sector takes the Einstein–Hilbert form plus a canonical scalar field, while matter acquires a $\varphi$-dependent non-minimal coupling with respect to the Einstein-frame metric ${\overline{g}}_{\mu \nu }$.

In terms of ${\overline{g}}_{\mu \nu }$ and $\varphi$, the action (2) becomes

$$S=\int d^{4}x\sqrt{-\overline{g}}\left[{\displaystyle \frac{M_{\mathrm{Pl}}^2}{2}}\overline{R}-{\displaystyle \frac{1}{2}}{\overline{g}}^{\mu \nu }{\overline{\nabla }}_{\mu }\varphi {\overline{\nabla }}_{n}u\varphi -V\left(\varphi \right)\right]+S_m\left(A^2\left(\varphi \right){\overline{g}}_{\mu \nu },{\Psi }_m\right),$$

(11)

where ${\overline{\nabla }}_{\mu }$ and $\overline{R}$ are the covariant derivative and Ricci scalar of ${\overline{g}}_{\mu \nu }$, and

$$A^2\left(\varphi \right)=F^{-1}\left(R\right)=exp\left(+{\displaystyle \frac{2\beta \varphi }{M_{\mathrm{Pl}}}}\right)$$

(12)

is the conformal factor relating the Einstein and Jordan metrics in the matter sector. The scalar potential is

$$V\left(\varphi \right)={\displaystyle \frac{M_{\mathrm{Pl}}^2}{2}}{\displaystyle \frac{RF\left(R\right)-f\left(R\right)}{F^2\left(R\right)}}={\displaystyle \frac{1}{\kappa }}{\displaystyle \frac{R{f}^{\prime }\left(R\right)-f\left(R\right)}{f^{\prime 2}\left(R\right)}}.$$

(13)

The corresponding Einstein-frame field equations can be written as

$$\begin{array}{cc}\hfill {\overline{G}}_{\mu \nu }& ={\displaystyle \frac{1}{M_{\mathrm{Pl}}^2}}\left(T_{\mu \nu }^{\left(\varphi \right)}+T_{\mu \nu }^{\left(m\right)}\right),\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \overline{\square }\varphi & =V^{\prime }\left(\varphi \right)-{\displaystyle \frac{\beta }{M_{\mathrm{Pl}}}}{T}^{\left(m\right)},\hfill \end{array}$$

(15)

where $T_{\mu \nu }^{\left(\varphi \right)}$ is the scalar-field stress tensor and $T_{\mu \nu }^{\left(m\right)}$ is the Einstein-frame matter tensor. For a spatially homogeneous scalar field $\varphi =\varphi \left(t\right)$ in a spatially flat FRW background, this equation reduces to the standard cosmological Klein–Gordon equation

$$\ddot{\varphi }+3H\dot{\varphi }+V^{\prime }\left(\varphi \right)={\displaystyle \frac{\beta }{M_{\mathrm{Pl}}}}{T}^{\left(m\right)},$$

(16)

which in the vacuum limit ($T^{\left(m\right)}=0$) becomes

$$\ddot{\varphi }+3H\dot{\varphi }+V^{\prime }\left(\varphi \right)=0,$$

(17)

the equation solved in inflationary and background cosmological applications [6,19,20]. This explicitly links the early-time inflationary dynamics of the model to its late-time cosmological implications.

• Jordan vs Einstein frame in practice.

The Jordan frame is the one in which matter is minimally coupled and experimental observables (such as test-particle trajectories and detector responses) are most directly interpreted. The Einstein frame is mathematically convenient for analyzing the dynamics of the extra scalar degree of freedom and for discussing stability, screening mechanisms, and cosmological evolution. Physical predictions are frame-independent provided one consistently transforms both the metric and matter variables. In this paper, we perform the gravitational-wave analysis in the Jordan frame (where the Bardeen variables and metric perturbations are defined), while the Einstein-frame description is used only to clarify the scalar–tensor structure and the chameleon mechanism.

• Chameleon mechanism in $f\left(R\right)$ gravity.

The additional scalar degree of freedom in metric $f\left(R\right)$ gravity mediates a universal fifth force [6,7] through its coupling to the trace of the matter stress–energy tensor. The scalar field $\varphi$ in (11) couples universally to matter via $A^2\left(\varphi \right){\overline{g}}_{\mu \nu }$ and can mediate this fifth force unless its effective mass becomes large in high-density environments. The chameleon mechanism exploits the density dependence of the effective potential

$$V_{\mathrm{eff}}\left(\varphi \right)=V\left(\varphi \right)+\rho A\left(\varphi \right),A\left(\varphi \right)=exp\left(+{\displaystyle \frac{\beta \varphi }{M_{\mathrm{Pl}}}}\right),$$

(18)

where $\rho$ is the local matter density. In regions of high density, $V_{\mathrm{eff}}$ develops a minimum at which the effective mass

$$m_{\mathrm{eff}}^2\left(\rho \right)\equiv {\left.{\displaystyle \frac{d^2{V}_{\mathrm{eff}}}{d{\varphi }^2}}\right|}_{{\varphi }_{\min }\left(\rho \right)}$$

(19)

is large, so that the scalar-mediated force is short-ranged and consistent with local tests of gravity. In low-density environments (cosmological scales) the minimum shifts and $m_{\mathrm{eff}}$ can become small enough for the scalar to drive cosmic acceleration or leave imprints on structure formation [21,22,23].

For an $f\left(R\right)$ model to exhibit viable chameleon behavior, the scalar potential $V\left(\varphi \right)$ derived from (13) must satisfy certain conditions in at least part of field space,

$$V^{\prime }\left(\varphi \right)<0,V^{\prime \prime }\left(\varphi \right)>0,V^{\prime \prime \prime }\left(\varphi \right)<0,$$

(20)

which translate into nontrivial constraints on the form of $f\left(R\right)$ and its derivatives [22]. This ensures that the scalar field can be heavy in high-density regions while remaining light enough on cosmological scales to influence late-time acceleration.

2.2. Slow-Roll Inflation and Scalaron Dynamics in the Einstein Frame

The scalar–tensor reformulation of metric $f\left(R\right)$ gravity introduced above provides a natural framework for discussing early-Universe inflation driven by the scalaron. When written in the Einstein frame, the scalar field $\varphi$ obeys the field equation obtained by varying the Einstein-frame action (11),

$$\overline{\square }\varphi =V^{\prime }\left(\varphi \right),$$

(21)

where $V\left(\varphi \right)$ is given in Equation (13). For a spatially homogeneous scalar field $\varphi =\varphi \left(t\right)$ evolving in a spatially flat FLRW background, $d{\overline{s}}^2=-d{t}^2+a^2\left(t\right)d{\overrightarrow{x}}^2$, this equation reduces to the standard cosmological Klein–Gordon equation [6,19,20]

$$\ddot{\varphi }+3H\dot{\varphi }+V^{\prime }\left(\varphi \right)=0,$$

(22)

where $H=\dot{a}/a$. The background expansion is governed by the Friedmann equation

$$3M_{\mathrm{Pl}}^2{H}^2={\displaystyle \frac{1}{2}}{\dot{\varphi }}^2+V\left(\varphi \right).$$

(23)

Inflation occurs when the potential dominates the kinetic energy and the field slowly rolls along $V\left(\varphi \right)$. This is quantified by the slow-roll parameters

$$ϵ\left(\varphi \right)\equiv {\displaystyle \frac{M_{\mathrm{Pl}}^2}{2}}{\left({\displaystyle \frac{V^{\prime }\left(\varphi \right)}{V\left(\varphi \right)}}\right)}^2,\eta \left(\varphi \right)\equiv M_{\mathrm{Pl}}^2{\displaystyle \frac{V^{\prime \prime }\left(\varphi \right)}{V\left(\varphi \right)}},$$

(24)

where inflation requires $ϵ\ll 1$ and $\left|\eta \right|\ll 1$. Under these conditions, Equation (22) reduces to the familiar slow-roll equation

$$3H\dot{\varphi }\simeq -V^{\prime }\left(\varphi \right),$$

(25)

and the Hubble parameter satisfies $H^2\simeq V\left(\varphi \right)/\left(3M_{\mathrm{Pl}}^2\right)$.

For the Starobinsky-type models considered in this work, the potential $V\left(\varphi \right)$ possesses a nearly flat region at large curvature ($R\gg \Lambda$), ensuring that the slow-roll conditions (24) are naturally satisfied. In this regime the scalaron behaves as an inflaton with an effective mass $m_{\varphi }^2\simeq V^{\prime \prime }\left(\varphi \right)$, and the inflationary predictions coincide with those of the well-known $R+\alpha R^2$ model [5,24]. Observationally, the slow-roll phase gives rise to a nearly scale-invariant spectrum of primordial curvature perturbations and a suppressed tensor-to-scalar ratio, in excellent agreement with current CMB constraints.

Although slow-roll inflation operates at curvature scales far above those relevant for present-day gravitational-wave detectors, the same underlying extra scalar degree of freedom of $f\left(R\right)$ gravity governs both regimes. In the inflationary context this degree of freedom is commonly referred to as the scalaron, while at the level of linear perturbations it appears as a propagating massive scalar mode. In particular, the mass of the scalar perturbation at the de Sitter solution,

$$m_{\psi }^2={\displaystyle \frac{1}{3}}\left[{\displaystyle \frac{f^{\prime }\left(R_d\right)}{f^{\prime \prime }\left(R_d\right)}}-R_d\right],$$

(26)

controls the propagation of the scalar polarization of gravitational waves and provides a link between the early-time inflationary dynamics of the model and its late-time cosmological implications.

• de Sitter solutions and cosmological implications.

We are particularly interested in constant-curvature de Sitter solutions and small perturbations around them. For a vacuum constant-curvature background with $R=R_d=\mathrm{const}$ and $T=0$, the trace equation (8) reduces to the algebraic condition

$$R_d{f}^{\prime }\left(R_d\right)-2f\left(R_d\right)=0.$$

(27)

Any function $f\left(R\right)$ that admits a solution of (27) possesses a de Sitter solution with curvature $R_d$. (Anti–de Sitter solutions correspond to constant-curvature solutions with $R_d<0$ and must be analyzed separately.) For the modified Starobinsky model

$$f\left(R\right)=R+\alpha R^2-2\Lambda ,$$

(28)

the de Sitter curvature $R_d$ is determined by

$$R_d(1+2\alpha R_d)-2(R_d+\alpha R_d^2-2\Lambda )=0,$$

(29)

which reduces to $R_d\simeq 4\Lambda$ in the late-time, low-curvature regime $\alpha R_d\ll 1$.1 Thus, constant-curvature solutions in $f\left(R\right)$ gravity provide a unified framework for modeling both early-time inflation and late-time dark-energy–dominated epochs, and they form the natural background for our gravitational-wave polarization analysis.

In the remainder of this paper, we will work in the Jordan frame and treat the extra propagating scalar mode directly in terms of the curvature perturbation $\delta R$. To avoid confusion of notation, we will use:

• $\varphi$ for the Einstein-frame scalar field entering the scalar–tensor and inflationary description in this subsection; and
• $\delta R$ (or, equivalently, a canonically normalized scalar perturbation $\psi$ with mass $m_{\psi }$) for the extra propagating scalar mode that appears in the linearized Jordan-frame field equations and in the Bardeen-variable analysis.

2.3. Metric Perturbations Around a de Sitter Background

We now consider small perturbations around a background solution $g_{\mu \nu }$ which solves the vacuum field equations (6). In particular, we will later specialize to a de Sitter background satisfying (27). The perturbed metric is written as

$$g_{\mu \nu }⟶{\tilde{g}}_{\mu \nu }=g_{\mu \nu }+h_{\mu \nu },{\tilde{g}}^{\mu \nu }=g^{\mu \nu }-h^{\mu \nu }+\mathcal{O}\left(h^2\right),$$

(30)

where indices on $h_{\mu \nu }$ are raised and lowered with the background metric $g_{\mu \nu }$.

To linear order, the curvature quantities and the function $f\left(R\right)$ expand as

$$\begin{array}{cc}\hfill \tilde{R}& =R+\delta R+\mathcal{O}\left(h^2\right),\hfill \\ \hfill \tilde{f}\left(R\right)& =f\left(R\right)+f^{\prime }\left(R\right)\delta R+\mathcal{O}\left(h^2\right),\hfill \\ \hfill {\tilde{f}}^{\prime }\left(R\right)& =f^{\prime }\left(R\right)+f^{\prime \prime }\left(R\right)\delta R+\mathcal{O}\left(h^2\right),\hfill \end{array}$$

(31)

where R is the background Ricci scalar and $\delta R$ is its perturbation.

It is convenient to separate background and perturbed covariant derivatives. Denoting by ${\tilde{\nabla }}_{\mu }$ the covariant derivative associated with ${\tilde{g}}_{\mu \nu }$ and by ${\nabla }_{\mu }$ the one associated with $g_{\mu \nu }$, the difference between them acting on a generic tensor ${T^{{\nu }_1\dots {\nu }_k}}_{{\rho }_1\dots {\rho }_l}$ is (see, e.g., [25])

$${\tilde{\nabla }}_{\mu }{T^{{\nu }_1\dots {\nu }_k}}_{{\rho }_1\dots {\rho }_l}={\nabla }_{\mu }{T^{{\nu }_1\dots {\nu }_k}}_{{\rho }_1\dots {\rho }_l}+\sum _i{C^{{\nu }_i}}_{\mu \sigma }{T^{{\nu }_1\dots \sigma \dots {\nu }_k}}_{{\rho }_1\dots {\rho }_l}-\sum _j{C^{\sigma }}_{\mu {\rho }_j}{T^{{\nu }_1\dots {\nu }_k}}_{{\rho }_1\dots \sigma \dots {\rho }_l},$$

(32)

where the connection difference ${C^{\rho }}_{\mu \nu }$ is

$${C^{\rho }}_{\mu \nu }={\displaystyle \frac{1}{2}}{\tilde{g}}^{\rho \sigma }\left({\nabla }_{\mu }{\tilde{g}}_{\nu \sigma }+{\nabla }_{\nu }{\tilde{g}}_{\mu \sigma }-{\nabla }_{\sigma }{\tilde{g}}_{\mu \nu }\right).$$

(33)

All quantities without tildes refer to the background metric. These relations allow one to express perturbed curvature tensors and the perturbed trace equation in terms of $h_{\mu \nu }$ and $\delta R$.

Varying the vacuum trace equation (8) and keeping terms linear in the perturbations yields

$$\begin{array}{cc}\hfill 0& =\delta \left[\square f^{\prime }\left(R\right)+{\displaystyle \frac{1}{3}}R{f}^{\prime }\left(R\right)-{\displaystyle \frac{2}{3}}f\left(R\right)\right]\hfill \\ \hfill & =\left[\square f^{\prime \prime }\left(R\right)+f^{\prime \prime }\left(R\right)\square -{\displaystyle \frac{1}{3}}\left(f^{\prime }\left(R\right)-R{f}^{\prime \prime }\left(R\right)\right)\right]\delta R\hfill \\ \hfill & -\left[h^{\mu \rho }{\nabla }_{\mu }+{\displaystyle \frac{1}{2}}{g}^{\mu \nu }{g}^{\rho \sigma }\left({\nabla }_{\mu }{h}_{\nu \sigma }+{\nabla }_{\nu }{h}_{\mu \sigma }-{\nabla }_{\sigma }{h}_{\mu \nu }\right)\right]\left(f^{\prime \prime }\left(R\right){\partial }_{\rho }R\right),\hfill \end{array}$$

(34)

where we have used (31) and (32). Equation (34) is valid for a general background.

For the de Sitter backgrounds of interest in this work, the Ricci scalar is constant,

$$R=R_d=\mathrm{const},$$

(35)

so that ${\nabla }_{\mu }R=0$ and $f^{\prime }\left(R_d\right)$, $f^{\prime \prime }\left(R_d\right)$ are constants. In this case the second line of (34) vanishes, $\square f^{\prime \prime }\left(R_d\right)=0$, and we obtain the simplified scalar perturbation equation

$$\left[\square -{\displaystyle \frac{1}{3}}\left({\displaystyle \frac{f^{\prime }\left(R_d\right)}{f^{\prime \prime }\left(R_d\right)}}-R_d\right)\right]\delta R=0,$$

(36)

or, equivalently,

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

(37)

with

$$m_{\psi }^2={\displaystyle \frac{1}{3}}\left[{\displaystyle \frac{f^{\prime }\left(R_d\right)}{f^{\prime \prime }\left(R_d\right)}}-R_d\right].$$

(38)

Here $m_{\psi }$ is the effective mass of the scalar propagating mode associated with the curvature perturbation $\delta R$ in the de Sitter background. Equation (37) is a Klein–Gordon equation for $\delta R$ and describes the propagation of a massive scalar mode in addition to the usual tensor modes of general relativity. In later sections we will relate $\delta R$ to a gauge-invariant Bardeen combination and denote the corresponding massive scalar propagating mode by $\psi$.

• Cosmological interpretation of $m_{\psi }$.

The mass scale $m_{\psi }^{-1}$ determines the range of the scalar-mediated interaction and the characteristic dispersion of the scalar polarization of gravitational waves in a de Sitter background. On sub-horizon scales with $k\gg a{m}_{\psi }$, the scalar mode behaves effectively massless and can, in principle, contribute to additional polarization signatures. On scales $k\ll a{m}_{\psi }$ the mode is strongly suppressed, consistent with local gravity constraints. For the modified Starobinsky model (28), one finds $m_{\psi }^2\approx 1/\left(6\alpha \right)$ in the high-curvature regime, linking the mass scale to both early-time inflationary dynamics and late-time modifications of gravitational-wave propagation in cosmology.

The detailed decomposition of the metric perturbations into scalar, vector, and tensor Bardeen variables, and the identification of the corresponding polarization modes, will be carried out in the following sections.

3. Perturbations of the Ricci Tensor $\mathit{\delta }{\mathit{R}}_{\mathit{\mu }\mathit{\nu }}$ and Scalar Dynamics in $\mathit{f}\left(\mathit{R}\right)$ Gravity

The evolution of cosmological perturbations in $f\left(R\right)$ gravity influences both the expansion history of the Universe and the propagation of gravitational waves across different cosmological epochs [26]. To understand how the gravitational field responds to small deviations from a background metric—whether a cosmological FRW background, a black hole spacetime, or, as in this work, a de Sitter background—it is necessary to compute the perturbation of curvature quantities. Since the Ricci tensor enters directly in the field equations, its perturbation represents the leading-order correction to the spacetime curvature and is essential for identifying the massive scalar propagating mode present in $f\left(R\right)$ theories.

Furthermore, gauge transformations in $f\left(R\right)$ gravity are complicated by the presence of higher derivatives of R. A fully gauge-invariant description of perturbations therefore requires determining how the scalar curvature perturbation $\delta R$ interacts with metric perturbations $h_{\mu \nu }$ through $\delta R_{\mu \nu }$. This provides a crucial intermediate step on the way to constructing the gauge-invariant Bardeen potentials in later sections.

To obtain the perturbed field equation, we expand the metric as $g_{\mu \nu }=g_{\mu \nu }^d+h_{\mu \nu }$ and linearize each term of the $f\left(R\right)$ field equation (6). Using $\delta f=f^{\prime }\left(R_d\right)\delta R$ and $\delta f^{\prime }=f^{\prime \prime }\left(R_d\right)\delta R$, and recalling that $R_d$ is constant, the variation of covariant derivative terms such as ${\nabla }_{\mu }{\nabla }_{\nu }{f}^{\prime }\left(R\right)$ must be treated carefully. In general,

$$\delta \left({\nabla }_{\mu }{\nabla }_{\nu }{f}^{\prime }\right)={\nabla }_{\mu }{\nabla }_{\nu }\delta f^{\prime }-\left(\delta {\Gamma }_{\mu \nu }^{\lambda }\right){\nabla }_{\lambda }{f}^{\prime }.$$

(39)

However, on a constant-curvature de Sitter background one has

$${\nabla }_{\lambda }{f}^{\prime }\left(R_d\right)=0,$$

(40)

so the connection-variation term vanishes identically. As a result, at linear order

$$\delta \left({\nabla }_{\mu }{\nabla }_{\nu }{f}^{\prime }\right)={\nabla }_{\mu }{\nabla }_{\nu }\delta f^{\prime }.$$

(41)

After accounting for this simplification, the variation of $(g_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu })f^{\prime }$ yields the operator $(g_{\mu \nu }^d\square -{\nabla }_{\mu }{\nabla }_{\nu })\delta f^{\prime }$.

Perturbing the vacuum field equation (6) using a de Sitter background ($R=R_d=\mathrm{const}$) yields the linearized equation

$$f^{\prime }\left(R_d\right)\delta R_{\mu \nu }+R_{\mu \nu }^d\delta f^{\prime }\left(R_d\right)-{\displaystyle \frac{1}{2}}\left[g_{\mu \nu }^d\delta f\left(R_d\right)+h_{\mu \nu }f\left(R_d\right)\right]+\left(g_{\mu \nu }^d\square -{\nabla }_{\mu }{\nabla }_{\nu }\right)\delta f^{\prime }\left(R_d\right)=0,$$

(42)

where $g_{\mu \nu }^d$ is the background de Sitter metric.

To obtain (42), we decompose the metric as

$$g_{\mu \nu }=g_{\mu \nu }^d+h_{\mu \nu },$$

(43)

with indices on $h_{\mu \nu }$ raised and lowered using $g_{\mu \nu }^d$. Since de Sitter space is a constant-curvature solution of metric $f\left(R\right)$ gravity, the background satisfies

$${\nabla }_{\mu }{f}^{\prime }\left(R_d\right)=0,$$

(44)

which follows directly from the constancy of $R_d$. This condition reflects the fact that $f^{\prime }\left(R_d\right)$ is a nonzero constant fixed by the chosen $f\left(R\right)$ model and represents the effective gravitational coupling on the background.

The nonperturbed trace equation (8) evaluated on a constant-curvature vacuum background yields the algebraic de Sitter condition

$$R_d{f}^{\prime }\left(R_d\right)=2f\left(R_d\right),$$

(45)

which determines the allowed background curvature $R_d$ of the $f\left(R\right)$ theory. Substituting (45) into the background field equation (6) and using $R=R_d$ and $T_{\mu \nu }=0$ gives

$$R_{\mu \nu }^d{f}^{\prime }\left(R_d\right)-{\displaystyle \frac{1}{2}}{g}_{\mu \nu }^{d}f\left(R_d\right)=0,$$

(46)

which implies

$$R_{\mu \nu }^d={\displaystyle \frac{1}{2}}{\displaystyle \frac{f\left(R_d\right)}{f^{\prime }\left(R_d\right)}}{g}_{\mu \nu }^d.$$

(47)

Using (45) once more yields

$$R_{\mu \nu }^d={\displaystyle \frac{1}{4}}{R}_d{g}_{\mu \nu }^d,$$

(48)

showing that the background is an Einstein space and, in fact, maximally symmetric.

For a spatially flat FLRW spacetime, the Ricci scalar is $R=6(2H^2+\dot{H})$. In de Sitter space $\dot{H}=0$, so $R_d=12H_d^2$ and $H_d=\sqrt{R_d/12}$.

A de Sitter Universe undergoes exponential expansion,

$$a\left(t\right)\propto e^{H_{d}t},H_d=\sqrt{R_d/12},$$

(49)

and this constant-curvature solution will serve as the background for our perturbative analysis.

Dividing (42) by $f^{\prime }\left(R_d\right)$ and substituting $\delta f=f_d^{\prime }\delta R$ and $\delta f^{\prime }=f_d^{\prime \prime }\delta R$ gives

$$\begin{array}{c}\hfill \delta R_{\mu \nu }+\left({\displaystyle \frac{f_d^{\prime \prime }}{f_d^{\prime }}}{R}_{\mu \nu }^d-{\displaystyle \frac{1}{2}}{g}_{\mu \nu }^d\right)\delta R-{\displaystyle \frac{1}{2}}{\displaystyle \frac{f_d}{f_d^{\prime }}}{h}_{\mu \nu }+{\displaystyle \frac{f_d^{\prime \prime }}{f_d^{\prime }}}\left(g_{\mu \nu }^d\square -{\nabla }_{\mu }{\nabla }_{\nu }\right)\delta R=0.\end{array}$$

(50)

At this stage, it is important to clarify the fate of the algebraic term $-{\displaystyle \frac{1}{2}}{\displaystyle \frac{f_d}{f_d^{\prime }}}{h}_{\mu \nu }$ appearing in Equation (50). On a constant–curvature de Sitter background, the trace condition (45) implies

$${\displaystyle \frac{f_d}{f_d^{\prime }}}={\displaystyle \frac{R_d}{2}},$$

(51)

so that this contribution may be written as $-(R_d/4)h_{\mu \nu }$. This term is proportional to the background curvature scale and contains no derivatives. For gravitational waves of wavelength $\lambda \ll H_d^{-1}$, corresponding to the local inertial (short–wavelength) limit relevant for detector-scale propagation, such curvature-suppressed algebraic terms do not contribute to the dynamical wave equation. They may therefore be consistently neglected, or equivalently absorbed into the background de Sitter curvature. With this understanding, the linearized field equation reduces to Equation (54).

To evaluate (50), we recall the general linearized Ricci tensor [25,27]

$$\delta R_{\mu \nu }={\displaystyle \frac{1}{2}}\left(-\square h_{\mu \nu }-{\nabla }_{\mu }{\nabla }_{\nu }h+{\nabla }^{\rho }{\nabla }_{\mu }{h}_{\rho \nu }+{\nabla }^{\rho }{\nabla }_{\nu }{h}_{\rho \mu }\right)+\mathcal{O}\left(h^2\right),$$

(52)

where $h\equiv {h^{\rho }}_{\rho }$.

For gravitational-wave propagation at detector scales, the wavelength of the perturbation is much shorter than the de Sitter curvature radius $H_d^{-1}$. One may therefore work in the local inertial frame of the background, in which

$$g_{\mu \nu }\approx {\eta }_{\mu \nu },R_{\mu \nu }\approx 0,$$

(53)

while still retaining the nonzero constants $f_d^{\prime }$ and $f_d^{\prime \prime }$. In this limit, curvature-suppressed algebraic terms proportional to $h_{\mu \nu }$ do not contribute to the dynamical propagation of high-frequency gravitational waves.

Under this approximation, (50) reduces to

$$\delta R_{\mu \nu }-{\displaystyle \frac{1}{2}}{\eta }_{\mu \nu }\delta R+{\displaystyle \frac{f_d^{\prime \prime }}{f_d^{\prime }}}({\eta }_{\mu \nu }\square -{\partial }_{\mu }{\partial }_{\nu })\delta R=0.$$

(54)

The linearized scalar curvature is

$$\delta R\equiv R\left[h\right],$$

(55)

and the linearized Einstein tensor reduces to

$$\delta G_{\mu \nu }=\delta R_{\mu \nu }-{\displaystyle \frac{1}{2}}{\eta }_{\mu \nu }\delta R\left[h\right].$$

(56)

Substituting (55) and (56) into (54) yields the perturbed field equation

$$\delta G_{\mu \nu }+{\displaystyle \frac{1}{3m_{\psi }^2}}({\eta }_{\mu \nu }\square -{\partial }_{\mu }{\partial }_{\nu })\delta R=0,$$

(57)

where the effective scalar mass $m_{\psi }$ arises from the trace of the linearized field equations and is given by

$$m_{\psi }^2={\displaystyle \frac{1}{3}}\left({\displaystyle \frac{f_d^{\prime }}{f_d^{\prime \prime }}}-R_d\right).$$

(58)

In the Minkowski limit $(R_d\to 0)$, this reduces to

$$m_{\psi }^2={\displaystyle \frac{1}{3}}{\displaystyle \frac{f_d^{\prime }}{f_d^{\prime \prime }}}.$$

(59)

Once the effective scalar mass $m_{\psi }^2$ is identified, its physical meaning becomes transparent by considering a Fourier (plane-wave) decomposition of the gauge-invariant variables. In a constant-curvature background, linear perturbations admit the ansatz $X(t,\mathbf{x})=X_k{e}^{-i(\omega t-\mathbf{k}·\mathbf{x})}$, which diagonalizes the spatial Laplacian. Substituting this into the Klein–Gordon–type equation $(\square -m_{\psi }^2)\delta R=0$ yields the dispersion relation ${\omega }^2=k^2+m_{\psi }^2$. Thus the extra scalar mode propagates as a massive mode, in contrast to the transverse tensor polarizations, which remain massless. This plane-wave form clarifies how the additional scalar polarization arises in metric $f\left(R\right)$ gravity.

These expressions show that the scalar curvature perturbation $\delta R$ propagates as a massive scalar field on the de Sitter background. The mass $m_{\psi }$ controls the range and dispersion of the scalar gravitational-wave mode and depends explicitly on the background curvature $R_d$. Since $R_d$ changes across cosmological epochs, the behavior of the scalar mode encodes information about the underlying cosmic expansion and offers potential observational signatures beyond the standard tensor modes.

4. 3+1 Decomposition and the Scalar, Vector, Tensor Modes of $\mathit{f}\left(\mathit{R}\right)$ Gravity

In this section we analyze the scalar, vector, and tensor perturbations of the metric using the standard scalar–vector–tensor (SVT) decomposition. As discussed in Section 2, gravitational waves detected at astrophysical scales propagate on a spacetime whose curvature radius is much larger than their wavelength. Therefore, for the purpose of the 3+1 decomposition we work in the local Minkowski limit of the de Sitter background,

$$g_{\mu \nu }\approx {\eta }_{\mu \nu },$$

(60)

while retaining the constant background quantities $f_d^{\prime }$, $f_d^{\prime \prime }$, and the mass $m_{\psi }$ of massive scalar propagating mode. Throughout this section we adopt the metric signature $(+---)$. For the scalar sector we work in the longitudinal (Newtonian) gauge, in which the metric perturbations are encoded in the gauge-invariant Bardeen potentials $\Phi$ and $\Psi$.

4.1. Scalar Mode

The 00 component of the linearized Einstein tensor in the longitudinal gauge is

$$\delta G_{00}=-2{\nabla }^2\Phi .$$

(61)

The 00 component of the perturbed $f\left(R\right)$ field equation (57) then takes the form

$$\delta G_{00}-{\displaystyle \frac{1}{3m_{\psi }^2}}{\nabla }^2\delta R=0.$$

(62)

Substituting (61) into (62) yields

$$\begin{array}{c}\hfill -2{\nabla }^2\Phi -{\displaystyle \frac{1}{3m_{\psi }^2}}{\nabla }^2\delta R=0,\end{array}$$

(63)

so that

$$2\Phi +{\displaystyle \frac{\delta R}{3m_{\psi }^2}}=0,$$

(64)

and therefore

$$\Phi =-{\displaystyle \frac{\delta R}{6m_{\psi }^2}}.$$

(65)

The traceless spatial components ($i\ne j$) yield the anisotropy constraint

$$\Phi -\Psi =\Pi ,$$

(66)

where $\Pi$ denotes the (gauge-invariant) anisotropic stress. In vacuum, $\Pi =0$, and the two Bardeen potentials coincide,

$$\Phi =\Psi =-{\displaystyle \frac{\delta R}{6m_{\psi }^2}}.$$

(67)

Thus the curvature perturbation $\delta R$ directly sources the scalar Bardeen potentials even in the absence of matter anisotropic stress, producing the massive scalar (breathing/longitudinal) gravitational-wave polarization predicted in $f\left(R\right)$ gravity. In a cosmological context, this relation links the extra gravitational-wave polarization to the scalar sector of cosmological perturbations and to the scalar mass $m_{\psi }$ on a de Sitter background.

It is worth noting that the expression $\Phi =\Psi =-\delta R/\left(6m_{\psi }^2\right)$ and the mass parameter $m_{\psi }^2$ used in this local 3+1 analysis correspond to the short-wavelength, locally Minkowskian limit of a de Sitter background. In this regime the curvature radius $H^{-1}$ is much larger than the gravitational-wave wavelength, and terms proportional to the background curvature $R_d$ are negligible. Consequently, the scalar mass reduces to $m_{\psi }^2\simeq \frac{1}{3}{f}_d^{\prime }/f_d^{\prime \prime }=1/\left(6\alpha \right)$. In a fully global de Sitter treatment, however, the effective mass contains an additional curvature contribution and takes the form $m_{\psi }^2=\frac{1}{3}(f^{\prime }\left(R_d\right)/f^{\prime \prime }\left(R_d\right)-R_d)$. Thus the local vacuum 3+1 decomposition used here captures the correct propagation physics for gravitational waves measured in a local inertial frame, while the global de Sitter mass governs the long-wavelength, cosmological evolution of the scalar mode.

4.2. Vector Modes

Vector perturbations appear in the $0i$ components of the metric as divergence-free vectors. The gauge-invariant combination is [28,29]

$$V_i\equiv S_i-F_i^{\prime },{\partial }^i{V}_i=0,$$

(68)

where $S_i$ enters the $g_{0i}$ component and $F_i$ the vector part of the spatial metric.

In the vector sector all scalar perturbations vanish, so in particular

$$\delta R=0.$$

(69)

The linearized Ricci tensor reduces to [3,4]

$$\delta R_{0i}=-{\displaystyle \frac{1}{2}}{\nabla }^2{V}_i.$$

(70)

The perturbed field equations in vacuum imply

$$\delta G_{0i}=\delta R_{0i}=0,$$

(71)

and hence

$${\nabla }^2{V}_i=0.$$

(72)

Under localized boundary conditions this yields $V_i=0$. As in general relativity, no vector modes propagate in vacuum metric $f\left(R\right)$ gravity.

4.3. Tensor Modes

We now examine the transverse–traceless (TT) tensor perturbations $h_{ij}^{TT}$. The spatial components of the perturbed $f\left(R\right)$ field equations take the form

$$\delta G_{ij}+{\displaystyle \frac{1}{3m_{\psi }^2}}({\delta }_{ij}\square -{\partial }_i{\partial }_j)\delta R=0.$$

(73)

The SVT decomposition of the perturbed Ricci tensor is [3,4,28,30]

$$\delta R_{ij}=-{\partial }_i{\partial }_j(\Phi +\Psi )-{\delta }_{ij}(-\ddot{\Phi }+{\nabla }^2\Phi )-{\displaystyle \frac{1}{2}}({\partial }_i{\dot{V}}_j+{\partial }_j{\dot{V}}_i)-{\displaystyle \frac{1}{2}}\square h_{ij}^{TT}.$$

(74)

In the pure tensor sector,

$$\Phi =\Psi =0,V_i=0,\delta R=0.$$

Equation (74) reduces to

$$\delta R_{ij}=-{\displaystyle \frac{1}{2}}\square h_{ij}^{TT}.$$

(75)

The term $({\delta }_{ij}\square -{\partial }_i{\partial }_j)\delta R$ in (73) has no TT projection and therefore drops out. Substituting into (73) yields

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

(76)

Thus, the tensor modes in metric $f\left(R\right)$ gravity propagate exactly as in GR: they satisfy the standard wave equation, travel at the speed of light, and possess only two transverse-traceless polarization states. All deviations from GR in gravitational-wave propagation therefore originate exclusively from the massive scalar propagating mode.

5. Analyzing a Specific $\mathit{f}\left(\mathit{R}\right)$ Model

We now specialize the general discussion of Secs. Section 2 and Section 3 to a concrete and widely studied model. One of the simplest and most successful choices is the Starobinsky model

$$f\left(R\right)=R+\alpha R^2+O\left(R^3\right),$$

(77)

which provides a purely geometric mechanism for early-Universe inflation and introduces an additional scalar degree of freedom through the higher-curvature term [5,24]. In this framework inflation is driven by the $R^2$ correction itself, rather than by an independent inflationary field, with the parameter $\alpha$ setting the characteristic inflationary scale and controlling the amplitude of primordial fluctuations [31,32].

For the truncated model (77), the first and second derivatives of $f\left(R\right)$ are

$$f^{\prime }\left(R\right)\approx 1+2\alpha R,$$

(78)

$$f^{\prime \prime }\left(R\right)\approx 2\alpha .$$

(79)

The vacuum field equation (6) takes the form

$$f^{\prime }\left(R\right)R_{\mu \nu }-{\displaystyle \frac{1}{2}}{g}_{\mu \nu }f\left(R\right)+\left(g_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu }\right)f^{\prime }\left(R\right)=0,$$

(80)

with the corresponding trace equation

$$3\square f^{\prime }\left(R\right)+R{f}^{\prime }\left(R\right)-2f\left(R\right)=0.$$

(81)

5.1. Constant-Curvature Backgrounds and the Need for a Modified Starobinsky Model

We now seek constant-curvature vacuum solutions characterized by

$$g_{\mu \nu }=g_{\mu \nu }^d,R=R_d=\mathrm{const},{\nabla }_{\mu }{R}_d=0.$$

(82)

For such backgrounds one has

$$\square f^{\prime }\left(R_d\right)=0,$$

(83)

and the trace equation (81) reduces to the algebraic condition

$$R_d{f}^{\prime }\left(R_d\right)-2f\left(R_d\right)=0.$$

(84)

Substituting the Starobinsky form (77) into (84) yields the unique solution

$$R_d=0.$$

(85)

Thus, within the pure Starobinsky model (77), Minkowski spacetime is the only constant-curvature vacuum solution. In particular, there is no nontrivial de Sitter background with $R_d>0$ that could describe an exponentially expanding late-time Universe.

From a cosmological standpoint this limitation motivates extending the model to include a vacuum-energy contribution. We therefore adopt the modified Starobinsky model

$$f\left(R\right)=R+\alpha R^2-2\Lambda ,$$

(86)

which supplements the inflationary $R+\alpha R^2$ sector with a cosmological constant term. In this model the quadratic curvature term governs early-time inflation [5], while the constant contribution $-2\Lambda$ drives late-time accelerated expansion [7,33].

For the modified model, the constant-curvature condition (84) yields

$$R_d=4\Lambda .$$

(87)

This solution defines a de Sitter background on which we will linearize the field equations. As shown in Section 3, a constant-$R_d$ background of this type satisfies

$$R_{\mu \nu }{|}_{R_d}={\displaystyle \frac{1}{4}}{R}_d{g}_{\mu \nu },$$

(88)

so the spacetime is an Einstein space, characterized by $R_{\mu \nu }\propto g_{\mu \nu }$. This geometric notion should be distinguished from the Einstein frame discussed in Section 2.1, which is obtained from the Jordan frame by a conformal transformation.

From a cosmological perspective, the modified Starobinsky model thus provides a unified description of the background dynamics: at high curvature the $\alpha R^2$ term drives inflation, while at low curvature the $-2\Lambda$ term yields a late-time de Sitter phase with curvature $R_d=4\Lambda$. This de Sitter solution serves as the background for the perturbative and gravitational-wave analyses developed in the following sections.

5.2. Einstein-Frame Potential and Stability Around the de Sitter Point

In Section 2.1 we reviewed the scalar–tensor (Einstein-frame) representation of metric $f\left(R\right)$ gravity. For the modified Starobinsky model

$$f\left(R\right)=R+\alpha R^2-2\Lambda ,$$

(89)

the Einstein-frame scalar potential is

$$V\left(\varphi \right)={\displaystyle \frac{M_{\mathrm{Pl}}^2}{2}}{\displaystyle \frac{\alpha R^2+2\Lambda }{{(1+2\alpha R)}^2}},$$

(90)

where the scalar field $\varphi$ is related to the curvature through the conformal relation

$$f^{\prime }\left(R\right)=1+2\alpha R=exp\left(-{\displaystyle \frac{2\beta \varphi }{M_{\mathrm{Pl}}}}\right),$$

(91)

with $\beta$ defined as in Section 2.1. Equation (91) implicitly defines $R=R\left(\varphi \right)$.

• Chain rule and curvature derivatives.

Using (91), the derivative of R with respect to $\varphi$ is

$${\displaystyle \frac{dR}{d\varphi }}=-{\displaystyle \frac{2\beta }{M_{\mathrm{Pl}}}}{\displaystyle \frac{f^{\prime }\left(R\right)}{f^{\prime \prime }\left(R\right)}}=-{\displaystyle \frac{\beta }{\alpha M_{\mathrm{Pl}}}}(1+2\alpha R),$$

(92)

where we used $f^{\prime \prime }\left(R\right)=2\alpha$. All derivatives of $V\left(\varphi \right)$ follow from repeated application of the chain rule $d/d\varphi =(dR/d\varphi )d/dR.$

• First derivative.

Differentiating (90) with respect to R gives

$${\displaystyle \frac{dV}{dR}}={\displaystyle \frac{M_{\mathrm{Pl}}^2}{2}}{\displaystyle \frac{(1+2\alpha R)(R-4\Lambda )}{{(1+2\alpha R)}^3}}.$$

(93)

Using (92), the first derivative of the potential becomes

$$V^{\prime }\left(\varphi \right)=-{\displaystyle \frac{\beta M_{\mathrm{Pl}}}{2\alpha }}{\displaystyle \frac{R-4\Lambda }{(1+2\alpha R)}}.$$

(94)

Thus $V^{\prime }\left(\varphi \right)=0$ precisely at

$$R=R_d=4\Lambda ,$$

(95)

corresponding to the de Sitter background.

• Second derivative.

Applying the chain rule once more yields

$$V^{\prime \prime }\left(\varphi \right)={\displaystyle \frac{{\beta }^2}{\alpha }}{\displaystyle \frac{1-2\alpha R+16\alpha \Lambda }{{(1+2\alpha R)}^2}}.$$

(96)

Evaluating this at the de Sitter point gives

$$V^{\prime \prime }\left({\varphi }_d\right)={\displaystyle \frac{{\beta }^2}{\alpha (1+8\alpha \Lambda )}}>0(\alpha >0),$$

(97)

showing that the de Sitter configuration corresponds to a local minimum of the Einstein-frame potential and is therefore linearly stable.

• Third derivative.

For completeness, the third derivative of the potential is

$$V^{\prime \prime \prime }\left(\varphi \right)={\displaystyle \frac{2{\beta }^3}{\alpha M_{\mathrm{Pl}}}}{\displaystyle \frac{32\alpha \Lambda -2\alpha R+3}{{(1+2\alpha R)}^2}},$$

(98)

which is nonvanishing and controls the leading self-interactions of the scalar mode around the de Sitter minimum. Equation (98) is fully consistent with the explicit expression obtained by direct differentiation in the Einstein frame.

• Stability interpretation.

The conditions

$$V^{\prime }\left({\varphi }_d\right)=0,V^{\prime \prime }\left({\varphi }_d\right)>0$$

(99)

establish that the modified Starobinsky model admits a stable de Sitter vacuum solution in metric $f\left(R\right)$ gravity. The corresponding scalar degree of freedom has positive mass squared, in agreement with the perturbative analysis of Section 3.

The full chameleon mechanism discussed in Section 2.1 requires including the matter coupling through the conformal factor $A\left(\varphi \right)$ and analyzing the density dependence of the effective potential $V_{\mathrm{eff}}\left(\varphi \right)=V\left(\varphi \right)+\rho A\left(\varphi \right)$. For the present discussion, it is sufficient to note that the vacuum Einstein-frame potential $V\left(\varphi \right)$ derived from the modified Starobinsky model admits a stable de Sitter minimum within metric $f\left(R\right)$ gravity. The resulting density-dependent scalar mass and screening behavior are encoded in the same Einstein-frame structure already introduced in Section 2.1.

5.3. Trace Perturbations and the Scalar Mass

We now revisit the trace equation in the Jordan frame in order to extract the explicit mass of the massive scalar propagating mode for our specific $f\left(R\right)$ model and to confirm consistency with the general result obtained in Section 3. We decompose the curvature scalar as

$$R=R_d+\delta R,$$

(100)

where $R_d$ is the constant-curvature de Sitter background.

The trace of the vacuum field equations,

$$3\square f^{\prime }\left(R\right)+R{f}^{\prime }\left(R\right)-2f\left(R\right)=0,$$

(101)

may be linearized about the de Sitter background. Using $\delta f^{\prime }=f^{\prime \prime }\left(R_d\right)\delta R$ and $\delta f=f^{\prime }\left(R_d\right)\delta R$, and retaining terms to first order in $\delta R$, the trace equation reduces to a Klein–Gordon equation for the scalar curvature perturbation,

$$\square \delta R-{\displaystyle \frac{1}{3}}\left({\displaystyle \frac{f^{\prime }\left(R_d\right)}{f^{\prime \prime }\left(R_d\right)}}-R_d\right)\delta R=0.$$

(102)

This form makes explicit that $\delta R$ propagates as a massive scalar field, with effective mass

$$m_{\psi }^2={\displaystyle \frac{1}{3}}\left({\displaystyle \frac{f^{\prime }\left(R_d\right)}{f^{\prime \prime }\left(R_d\right)}}-R_d\right),$$

(103)

in agreement with the general expression derived earlier from the linearized field equations.

We now specialize to the modified Starobinsky model,

$$f\left(R\right)=R+\alpha R^2-2\Lambda .$$

(104)

For this choice one finds

$$f^{\prime }\left(R_d\right)=1+2\alpha R_d,f^{\prime \prime }\left(R_d\right)=2\alpha .$$

(105)

Substituting these expressions into Equation (102) yields

$$\left(\square -{\displaystyle \frac{1}{6\alpha }}\right)\delta R=0,$$

(106)

so that the scalar curvature perturbation satisfies a Klein–Gordon equation with mass

$$m_{\psi }^2={\displaystyle \frac{1}{6\alpha }},$$

(107)

independent of the background curvature $R_d$.

For the modified Starobinsky model, the de Sitter background curvature is fixed by the cosmological constant through

$$R_d=4\Lambda ,$$

(108)

while the mass of the additional scalar degree of freedom is controlled entirely by the quadratic coupling $\alpha$,

$$m_{\psi }^{-1}=\sqrt{6\alpha }.$$

(109)

From a cosmological perspective, $\alpha$ determines the range and dispersion scale of the scalar polarization of gravitational waves, whereas $\Lambda$ fixes the asymptotic de Sitter curvature. This clean separation of roles will be important when we discuss the propagation of the scalar mode and its potential observational signatures in de Sitter cosmology.

6. SVT Decomposition of the Perturbed Ricci Tensor in Metric $\mathit{f}\left(\mathit{R}\right)$ Gravity

In this section we revisit the $(3+1)$ decomposition in the presence of matter sources. Instead of starting from the vacuum perturbed field equation (57), we now consider the linearized field equations of the modified Starobinsky model in a nearly Minkowski background, including the stress–energy tensor $T_{\mu \nu }$:

$$\delta G_{\mu \nu }+2\alpha \left({\eta }_{\mu \nu }\square -{\partial }_{\mu }{\partial }_{\nu }\right)\delta R=\kappa T_{\mu \nu },$$

(110)

where $\alpha >0$ is the $R^2$ coupling, $\delta R$ is the scalar curvature perturbation, and ${\eta }_{\mu \nu }$ is the background Minkowski metric. Using the relation

$$m_{\psi }^2={\displaystyle \frac{1}{6\alpha }},$$

(111)

the term proportional to $\alpha$ can also be written as $(1/3m_{\psi }^2)({\eta }_{\mu \nu }\square -{\partial }_{\mu }{\partial }_{\nu })\delta R$, in agreement with the vacuum analysis.

The Klein–Gordon equation for the massive extra scalar mode (106) in vacuum generalizes in the presence of matter to

$$(\square -m_{\psi }^2)\delta R=m_{\psi }^2\kappa T,$$

(112)

or equivalently

$$\square \delta R=m_{\psi }^2(\delta R+\kappa T),$$

(113)

where $T\equiv {\eta }^{\mu \nu }{T}_{\mu \nu }$ is the trace of the stress–energy tensor. Equation (112) shows that $\delta R$ behaves as a massive scalar field (the massive scalar propagating mode) sourced by the trace T; in the limit $T\to 0$ we recover the vacuum equation.

Using the definition of the Einstein tensor in the flat background,

$$\delta G_{\mu \nu }=\delta R_{\mu \nu }-{\displaystyle \frac{1}{2}}{\eta }_{\mu \nu }\delta R,$$

(114)

and eliminating $\square \delta R$ via (113), the linearized field equation (110) may be written as

$$\delta R_{\mu \nu }-{\displaystyle \frac{1}{2}}{\eta }_{\mu \nu }\delta R+{\displaystyle \frac{1}{3}}{\eta }_{\mu \nu }(\delta R+\kappa T)-2\alpha {\partial }_{\mu }{\partial }_{\nu }\delta R=\kappa T_{\mu \nu },$$

(115)

or, equivalently,

$$\delta R_{\mu \nu }-\left({\displaystyle \frac{1}{6}}{\eta }_{\mu \nu }+2\alpha {\partial }_{\mu }{\partial }_{\nu }\right)\delta R=\kappa \left[T_{\mu \nu }-{\displaystyle \frac{1}{3}}{\eta }_{\mu \nu }T\right].$$

(116)

Equations (115) and (116) are the starting point for the $(3+1)$ decomposition with matter: the left-hand side contains the usual Ricci-tensor perturbation corrected by the massive scalar propagating mode $\delta R$, while the right-hand side involves the traceless combination $T_{\mu \nu }-{\eta }_{\mu \nu }T/3$.

6.1. Irreducible SVT Decomposition of the Metric and Matter

Following [4,10], the SVT decomposition of the metric perturbation $h_{\mu \nu }$ in a nearly Minkowski background reads

$$\begin{array}{cc}\hfill h_{00}& =2\psi ,\hfill \end{array}$$

(117)

$$\begin{array}{cc}\hfill h_{0i}& ={\beta }_i+{\partial }_i\gamma ,\hfill \end{array}$$

(118)

$$\begin{array}{cc}\hfill h_{ij}& =-2\varphi {\delta }_{ij}+\left({\partial }_i{\partial }_j-{\displaystyle \frac{1}{3}}{\delta }_{ij}{\nabla }^2\right)\lambda +{\displaystyle \frac{1}{2}}\left({\partial }_i{ϵ}_j+{\partial }_j{ϵ}_i\right)+h_{ij}^{TT},\hfill \end{array}$$

(119)

where we have defined the new quantities $\psi ,{\beta }_i,\gamma ,\varphi ,{ϵ}_i,\lambda ,h_{ij}^{TT}$ with the assumption that $h_{\mu \nu }\to 0$ as $r\to \infty$. The transverse and traceless conditions are

$$\begin{array}{cc}\hfill {\partial }^i{\beta }_i& =0,\hfill \end{array}$$

(120)

$$\begin{array}{cc}\hfill {\partial }^i{ϵ}_i& =0,\hfill \end{array}$$

(121)

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

(122)

$$\begin{array}{cc}\hfill {\delta }^{ij}{h}_{ij}^{TT}& =0.\hfill \end{array}$$

(123)

Both in [4] and in our earlier Bardeen-variable work, it has been shown how the variables $\psi ,{\beta }_i,\gamma ,\varphi ,{ϵ}_i,\lambda ,h_{ij}^{TT}$ transform under a gauge transformation generated by ${\xi }^a$ with ${\xi }^a\to 0$ as $r\to \infty$. Such transformations are parametrized as

$${\xi }^a=({\xi }^0,{\xi }^i)=(A,B^i+{\partial }^{i}C),$$

(124)

with ${\partial }_i{B}^i=0$. Following the same procedure as in [4], one obtains the gauge-invariant scalar and vector combinations

$$\begin{array}{cc}\hfill \Phi & \equiv -\varphi -{\displaystyle \frac{1}{6}}{\nabla }^2\lambda ,\hfill \end{array}$$

(125)

$$\begin{array}{cc}\hfill \Psi & \equiv -\psi +\dot{\gamma }-{\displaystyle \frac{1}{2}}\ddot{\lambda },\hfill \end{array}$$

(126)

$$\begin{array}{cc}\hfill V_i& \equiv {\beta }_i-{\displaystyle \frac{1}{2}}{\dot{ϵ}}_i,{\partial }_i{V}^i=0.\hfill \end{array}$$

(127)

The tensor perturbation $h_{ij}^{TT}$ is already gauge invariant.

We can perform a similar SVT decomposition of the matter stress–energy tensor $T_{\mu \nu }$ on the right-hand side of the field equations. We write

$$\begin{array}{cc}\hfill T_{00}& =\rho ,\hfill \end{array}$$

(128)

$$\begin{array}{cc}\hfill T_{0i}& =S_i+{\partial }_{i}S,\hfill \end{array}$$

(129)

$$\begin{array}{cc}\hfill T_{ij}& =-P{\delta }_{ij}+{\sigma }_{ij}+({\partial }_i{\sigma }_j+{\partial }_j{\sigma }_i)+\left({\partial }_i{\partial }_j-{\displaystyle \frac{1}{3}}{\delta }_{ij}{\nabla }^2\right)\sigma ,\hfill \end{array}$$

(130)

where $\rho$, S, $S_i$, P, $\sigma$, ${\sigma }_i$, and ${\sigma }_{ij}$ are new scalar, vector, and tensor quantities with the constraints

$$\begin{array}{cc}\hfill {\partial }_i{S}^i& =0,\hfill \end{array}$$

(131)

$$\begin{array}{cc}\hfill {\partial }^i{\sigma }_i& =0,\hfill \end{array}$$

(132)

$$\begin{array}{cc}\hfill {\partial }^i{\sigma }_{ij}& =0,\hfill \end{array}$$

(133)

$$\begin{array}{cc}\hfill {\delta }^{ij}{\sigma }_{ij}& =0,\hfill \end{array}$$

(134)

along with boundary conditions $S\to 0,{\sigma }_i\to 0,\sigma \to 0,{\nabla }^2\sigma \to 0$ as $r\to \infty$ (spatial infinity). The overall minus sign in the isotropic part $-P{\delta }_{ij}$ in (130) will be tracked explicitly in the relations obtained from stress–energy conservation below.

The conservation law

$${\partial }^{\mu }{T}_{\mu \nu }=0,$$

(135)

determines relations between $\rho$, S, $S_i$, P, $\sigma$, ${\sigma }_i$, and ${\sigma }_{ij}$. In the nearly Minkowski background used throughout this section we have ${\partial }^{\mu }=(-{\partial }_0,{\partial }^i)$, so (135) reads

$$-{\partial }_0{T}_{0\nu }+{\partial }^i{T}_{i\nu }=0.$$

(136)

As a useful special case (and for later physical interpretation), we note that a perfect fluid at rest in Minkowski space has

$$T_{\mu \nu }=\left[\begin{array}{cccc}\rho & 0& 0& 0\\ 0& -P& 0& 0\\ 0& 0& -P& 0\\ 0& 0& 0& -P\end{array}\right],$$

(137)

so that $T_{0i}=0$ and $T_{ij}=-P{\delta }_{ij}$. Comparing with (129)–(130), this corresponds to vanishing momentum and anisotropic-stress components ($S_i=0$, $S=0$, ${\sigma }_i=0$, $\sigma =0$, ${\sigma }_{ij}=0$), while retaining the scalars $\rho$ and P. In particular, the trace is

$$T=\rho -3P.$$

(138)

The $\nu =0$ component of (135) gives

$$\dot{\rho }={\partial }^i{T}_{i0}={\partial }^i{T}_{0i}={\nabla }^{2}S,$$

(139)

where we have used (129) and the constraint (131). For $\nu =i$, it is convenient to separate the two pieces entering $-{\partial }_0{T}_{0i}+{\partial }^j{T}_{ji}=0$. First, taking a spatial derivative of (130) yields

$${\partial }^j{T}_{ji}=-{\partial }_{i}P+{\displaystyle \frac{2}{3}}{\nabla }^2{\partial }_i\sigma +{\nabla }^2{\sigma }_i+{\partial }^j{\sigma }_{ji},$$

(140)

and the constraints (132)–(133) imply ${\partial }^j{\sigma }_{ji}=0$. Second, taking a time derivative of (129) gives

$${\partial }_0{T}_{0i}={\dot{S}}_i+{\partial }_i\dot{S}.$$

(141)

Stress–energy conservation for $\nu =i$ then combines (140) and (141) as

$${\partial }_i\left(-P-\dot{S}+{\displaystyle \frac{2}{3}}{\nabla }^2\sigma \right)-{\dot{S}}_i+{\nabla }^2{\sigma }_i=0,$$

(142)

where we used ${\partial }^j{\sigma }_{ji}=0$.

Taking one more spatial derivative of (142) and applying the constraints on ${\sigma }_i$ and $S_i$, we obtain

$${\nabla }^2\left[-P-\dot{S}+{\displaystyle \frac{2}{3}}{\nabla }^2\sigma \right]=0.$$

(143)

Applying the boundary condition at spatial infinity for S, P, and ${\nabla }^2\sigma$ (which also guarantees the uniqueness of the decomposition), we conclude that

$$-P-\dot{S}+{\displaystyle \frac{2}{3}}{\nabla }^2\sigma =0.$$

(144)

Inserting this condition into (142) gives

$${\nabla }^2{\sigma }^i={\dot{S}}^i.$$

(145)

Equation (144) can also be rewritten as

$${\nabla }^2\sigma ={\displaystyle \frac{3}{2}}\dot{S}-{\displaystyle \frac{3}{2}}P.$$

(146)

Equations (139), (145), and (146) are the required set of differential equations that relate the newly defined irreducible matter variables $\rho$, $S_i$, S, P, $\sigma$, ${\sigma }_i$, and ${\sigma }_{ij}$. These results match Eq. (23) of [34], up to differences in notation.

6.2. (3+1) Decomposition of the Perturbed Ricci Tensor

In $f\left(R\right)$ gravity the field equations contain higher-order derivatives of the metric through their dependence on the Ricci scalar R. Unlike in General Relativity, where the Einstein tensor $G_{\mu \nu }$ alone determines the dynamics, $f\left(R\right)$ theories introduce a massive scalar propagating mode, associated with the scalar curvature perturbation $\delta R$. To fully understand how this scalar mode interacts with the usual scalar, vector, and tensor components of the metric perturbation, it is necessary to go beyond the standard metric decomposition and analyze the perturbation of the Ricci tensor $\delta R_{\mu \nu }$ itself.

By expressing $\delta R_{\mu \nu }$ in terms of the gauge-invariant Bardeen variables and the scalar curvature perturbation, we obtain a set of decoupled differential equations that reveal how each mode behaves in the presence of matter. This decomposition provides a more complete and transparent description of the linearized dynamics in $f\left(R\right)$ gravity and is particularly useful for identifying modifications to gravitational-wave propagation and structure formation due to the extra scalar mode.

Following the approach in [4], the components of the Einstein tensor $G_{\mu \nu }$ were decomposed in GR to obtain a set of differential equations for the perturbation of the Einstein tensor in terms of the Bardeen variables. In the case of $f\left(R\right)$ gravity, we instead decompose the perturbed Ricci tensor $\delta R_{\mu \nu }$ in terms of the Bardeen variables, to obtain a new set of differential equations. The components of $\delta R_{\mu \nu }$ are

$$\begin{array}{cc}\hfill \delta R_{00}& ={\nabla }^2\Psi -{\displaystyle \frac{3}{2}}\ddot{\Phi },\hfill \end{array}$$

(147)

$$\begin{array}{cc}\hfill \delta R_{0i}& =-{\displaystyle \frac{1}{2}}{\nabla }^2{V}_i-{\partial }_i\dot{\Phi },\hfill \end{array}$$

(148)

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

(149)

In terms of these Bardeen variables $\Phi ,\Psi ,V_i$, the field equation in the form of Equation (116) can be recast into a set of differential equations, each corresponding to a component of $\delta R_{\mu \nu }$. For example, the 00 component of (116) takes the form

$$\delta R_{00}-{\displaystyle \frac{1}{6m_{\psi }^2}}\left(m_{\psi }^2{\eta }_{00}+2{\partial }_0{\partial }_0\right)\delta R=\kappa \left(T_{00}-{\displaystyle \frac{1}{3}}{\eta }_{00}T\right).$$

(150)

Substituting the expression for $\delta R_{00}$ from Equation (147), using ${\eta }_{00}=-1$ and $T=\rho -3P$, we obtain

$${\nabla }^2\Psi -{\displaystyle \frac{3}{2}}\ddot{\Phi }+{\displaystyle \frac{1}{6}}\delta R-{\displaystyle \frac{1}{3m_{\psi }^2}}\ddot{\delta R}=\kappa \left({\displaystyle \frac{4}{3}}\rho -P\right),$$

(151)

which gives the corresponding differential equation for the 00 component of $\delta R_{\mu \nu }$.

Next we consider the differential equations corresponding to the $0i$ component of the perturbation of the Ricci tensor. Equation (116) gives

$$\delta R_{0i}-{\displaystyle \frac{1}{6m_{\psi }^2}}\left(m_{\psi }^2{\eta }_{0i}+2{\partial }_0{\partial }_i\right)\delta R=\kappa \left(T_{0i}-{\displaystyle \frac{1}{3}}{\eta }_{0i}T\right).$$

(152)

Substituting the expression (148) for $\delta R_{0i}$ in Equation (152), Equation (129) for $T_{0i}$, and using ${\eta }_{0i}=0$, we obtain

$$-{\displaystyle \frac{1}{2}}{\nabla }^2{V}_i-{\partial }_i\dot{\Phi }-{\displaystyle \frac{1}{3m_{\psi }^2}}{\partial }_0{\partial }_i\delta R=\kappa (S_i+{\partial }_{i}S).$$

(153)

At spatial infinity ($r\to \infty$), we impose $S\to 0$, so that ${\partial }_{i}S\to 0$, and similarly ${\partial }_i\dot{\Phi }\to 0$ and ${\partial }_0{\partial }_i\delta R\to 0$. Under these conditions Equation (153) reduces to

$${\nabla }^2{V}_i=-2\kappa S_i.$$

(154)

Separating the longitudinal and transverse parts of (153) and comparing the coefficients of ${\partial }_i$ yields

$$\dot{\Phi }+{\displaystyle \frac{1}{3m_{\psi }^2}}\dot{\delta R}=-\kappa S.$$

(155)

Equations (154) and (155) are the differential equations based on the $0i$ component of $\delta R_{\mu \nu }$ in terms of the Bardeen variables and the massive scalar propagating mode.

Finally, we consider the $ij$ component of the perturbation of the Ricci tensor, which is

$$\delta R_{ij}-{\displaystyle \frac{1}{6m_{\psi }^2}}\left(m_{\psi }^2{\eta }_{ij}+2{\partial }_i{\partial }_j\right)\delta R=\kappa \left(T_{ij}-{\displaystyle \frac{1}{3}}{\eta }_{ij}T\right).$$

(156)

Substituting the expression (149) for $\delta R_{ij}$ and Equation (130) for $T_{ij}$ into Equation (156), and equating coefficients of the independent SVT pieces, we obtain

$$\begin{array}{cc}\hfill -{\partial }_{(i}{\dot{ϵ}}_{j)}& =\kappa {\partial }_{(i}{\sigma }_{j)},\hfill \end{array}$$

(157)

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

(158)

$$\begin{array}{cc}\hfill -{\displaystyle \frac{1}{2}}\square {\delta }_{ij}\Phi -{\displaystyle \frac{1}{3}}{\partial }_i{\partial }_j\delta R& =\kappa \left(-{\displaystyle \frac{1}{3}}{\delta }_{ij}{\nabla }^2\sigma -{\displaystyle \frac{1}{3}}{\eta }_{ij}T\right).\hfill \end{array}$$

(159)

These equations imply

$$\begin{array}{cc}\hfill {\dot{V}}_i& =\kappa {\sigma }_i,\hfill \end{array}$$

(160)

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

(161)

$$\begin{array}{cc}\hfill \Psi +{\displaystyle \frac{1}{2}}\Phi +{\displaystyle \frac{1}{3m_{\psi }^2}}\delta R& =-\kappa \sigma ,\hfill \end{array}$$

(162)

and, using $T=\rho -3P$,

$$\square \Phi +{\displaystyle \frac{1}{3}}\delta R={\displaystyle \frac{2}{3}}\kappa \left({\nabla }^2\sigma -\rho \right).$$

(163)

Equations (160)–(163) are the set of differential equations corresponding to the $ij$ component of $\delta R_{\mu \nu }$ in terms of the Bardeen variables and the matter SVT variables. These results are consistent with those derived in [34], up to differences in notation.

6.3. Cosmological Interpretation of the SVT Equations with the Extra Scalar Degree of Freedom

The system of equations (151), (154), (155), and (160)–(163) allows a direct physical interpretation in cosmology once the background is promoted from Minkowski to a slowly varying FLRW or de Sitter spacetime.

The $\delta R_{00}$ equation (151) is a modified Poisson-type equation: the gravitational potential $\Psi$ is sourced not only by the energy density $\rho$, but also by pressure P, time derivatives of the potential $\Phi$, and the dynamics of the scalar curvature perturbation $\delta R$[35,36]. In GR, the corresponding equation at linear order would involve essentially the Laplacian of $\Psi$ sourced only by $\rho$, with no extra scalar degree of freedom contribution. This modification leads to a scale- and time-dependent effective gravitational coupling, which directly affects the growth of cosmological structure and can be constrained by large-scale structure and weak-lensing surveys.

The $\delta R_{0i}$ sector separates into a transverse (vector) part and a longitudinal (scalar) part. The transverse part, Equation (154) together with Equation (160), shows that vector perturbations $V_i$ are sourced by the transverse momentum density $S_i$ and anisotropic stress ${\sigma }_i$, just as in GR. Thus $f\left(R\right)$ gravity does not introduce new propagating vector modes at linear order. The longitudinal scalar equation (155), however, contains the time derivative of the additional scalar degree of freedom $\dot{\delta R}$, modifying the time evolution of $\Phi$ relative to GR. The time dependence of the gravitational potentials is directly probed by the integrated Sachs–Wolfe (ISW) effect and cross-correlations of CMB maps with large-scale structure.

The $\delta R_{ij}$ equations show that the tensor sector, Equation (161), obeys a wave equation structurally identical to that of GR, but with a source from anisotropic stress. In $f\left(R\right)$ gravity, the background scalar degree of freedom and the modified expansion history can nevertheless change the amplitude damping and effective propagation of gravitational waves over cosmological distances, providing an additional channel to test modifications of gravity with standard sirens.

Finally, the scalar sector of the perturbed field equations provides a direct window into one of the characteristic phenomenological signatures of modified gravity. In linear cosmological perturbation theory, scalar metric perturbations are described by the gauge-invariant Bardeen potentials $\Phi$ and $\Psi$, which coincide in General Relativity in the absence of matter anisotropic stress. Their inequality, $\Phi \ne \Psi$, is commonly referred to as gravitational slip and signals a departure from GR caused either by imperfect fluids or by additional gravitational degrees of freedom [36,37].

In metric $f\left(R\right)$ gravity, the scalar part of the $\delta R_{ij}$ equations, Equations (162) and (163), reveals that gravitational slip arises generically even when the matter anisotropic stress vanishes ($\sigma =0$, equivalently $\Pi =0$). In this case, the difference between the two scalar potentials is instead sourced by the scalar curvature perturbation $\delta R$, reflecting the presence of the propagating scalar mode. This modification of the relation between $\Phi$ and $\Psi$ is a robust signature of $f\left(R\right)$ models and can be observationally constrained through joint analyses of galaxy clustering, redshift-space distortions, and weak gravitational lensing [36,38]. The full SVT decomposition of $\delta R_{\mu \nu }$ thus provides a unified framework for linking the gauge-invariant scalar dynamics of the theory to observable effects in both gravitational-wave physics and cosmology.

7. Geodesic Deviation Method to Find the Polarization Content

The geodesic deviation equation relates the Riemann curvature tensor to the relative acceleration of neighboring geodesics and therefore provides a direct probe of gravitational-wave polarizations in a given theory of gravity [27,39]. In this section we use the geodesic deviation equations to identify the polarization modes of gravitational waves in our specific metric $f\left(R\right)$ model,

$$f\left(R\right)=R+\alpha R^2-2\Lambda ,$$

(164)

for which the scalar curvature perturbation $\delta R$ obeys the massive Klein–Gordon equation

$$\left(\square -m_{\psi }^2\right)\delta R=0,m_{\psi }^2={\displaystyle \frac{1}{6\alpha }},$$

(165)

on a de Sitter background. The scalar perturbation $\delta R$ corresponds to the extra scalar degree of freedom, in addition to the usual tensor modes of GR.

We first work in the local Minkowski patch of the de Sitter background, which is appropriate for interferometric detectors whose size is much smaller than the background curvature radius. We then show how the same polarization structure appears when the calculation is formulated fully on a de Sitter FRW background.

7.1. Local Minkowski Patch of de Sitter

The general geodesic deviation equation is

$${\displaystyle \frac{D^2{\xi }^{\mu }}{d{\tau }^2}}=-{R^{\mu }}_{\alpha \nu \beta }{\xi }^{\nu }{\displaystyle \frac{d{x}^{\alpha }}{d\tau }}{\displaystyle \frac{d{x}^{\beta }}{d\tau }},$$

(166)

where ${\xi }^{\mu }$ is the separation vector between neighboring geodesics and $\tau$ is proper time. For gravitational-wave detectors we work in the weak-field, slow-motion limit: the detector is at rest in the chosen coordinates and far from the source, so

$${\displaystyle \frac{d{x}^i}{d\tau }}\ll {\displaystyle \frac{d{x}^0}{d\tau }}\Rightarrow {\displaystyle \frac{d{x}^{\mu }}{d\tau }}\approx (1,0,0,0),$$

(167)

and we can identify proper time with coordinate time,

$$\tau \simeq t.$$

(168)

In this regime the covariant derivatives in (166) reduce to ordinary time derivatives, and the spatial components of the geodesic deviation equation become

$${\ddot{\xi }}_i=-R_{i0j0}{\xi }^j,$$

(169)

where overdots denote derivatives with respect to t.

In linearized gravity, the Riemann tensor is

$$R_{\mu \nu \rho \sigma }={\displaystyle \frac{1}{2}}\left({\partial }_{\rho }{\partial }_{\nu }{h}_{\mu \sigma }+{\partial }_{\sigma }{\partial }_{\mu }{h}_{\nu \rho }-{\partial }_{\sigma }{\partial }_{\nu }{h}_{\mu \rho }-{\partial }_{\rho }{\partial }_{\mu }{h}_{\nu \sigma }\right),$$

(170)

where $h_{\mu \nu }$ is the metric perturbation on the local Minkowski background ${\eta }_{\mu \nu }$. We decompose $h_{\mu \nu }$ into a transverse-traceless tensor part $h_{\mu \nu }^{TT}$ and a scalar part associated with the scalar curvature perturbation $\delta R$.

For the scalar mode, in a convenient gauge compatible with the Newtonian (longitudinal) gauge used in Section 4, the scalar perturbation can be chosen proportional to the background metric:

$$h_{\mu \nu }^{\left(s\right)}=C\delta R{\eta }_{\mu \nu },$$

(171)

where C is an overall constant that only rescales the amplitude and does not affect the polarization pattern. For simplicity we set $C=1$ below.

Explicitly,

$$\begin{array}{cc}\hfill h_{00}^{\left(s\right)}& =-\delta R,\hfill \end{array}$$

(172)

$$\begin{array}{cc}\hfill h_{0i}^{\left(s\right)}& =0,\hfill \end{array}$$

(173)

$$\begin{array}{cc}\hfill h_{ij}^{\left(s\right)}& =\delta R{\delta }_{ij}.\hfill \end{array}$$

(174)

Substituting into the expression for the Riemann tensor and focusing on $R_{i0j0}$, we find

$$R_{i0j0}=-{\displaystyle \frac{1}{2}}\left({\partial }_0{\partial }_0{h}_{ij}+{\partial }_i{\partial }_j{h}_{00}\right)=-{\displaystyle \frac{1}{2}}\left({\delta }_{ij}\ddot{\delta R}-{\partial }_i{\partial }_j\delta R\right).$$

(175)

Now consider a scalar wave propagating along the $+z$ direction,

$$\delta R=\delta R(t,z).$$

(176)

In the transverse directions x and y,

$${\partial }_x\delta R={\partial }_y\delta R=0,$$

(177)

so that

$$R_{x0x0}=R_{y0y0}=-{\displaystyle \frac{1}{2}}\ddot{\delta R}.$$

(178)

In the longitudinal direction,

$$R_{z0z0}=-{\displaystyle \frac{1}{2}}\left(\ddot{\delta R}-{\partial }_z^2\delta R\right).$$

(179)

Using the massive Klein–Gordon equation for the extra scalar mode in the local Minkowski patch,

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

(180)

which implies

$$\ddot{\delta R}={\partial }_z^2\delta R-m_{\psi }^2\delta R,$$

(181)

we obtain

$$\begin{array}{cc}\hfill R_{x0x0}& =R_{y0y0}=-{\displaystyle \frac{1}{2}}\left({\partial }_z^2\delta R-m_{\psi }^2\delta R\right),\hfill \end{array}$$

(182)

$$\begin{array}{cc}\hfill R_{z0z0}& =+{\displaystyle \frac{1}{2}}{m}_{\psi }^2\delta R.\hfill \end{array}$$

(183)

For a monochromatic plane wave

$$\delta R(t,z)=A{e}^{i(k_{z}z-\omega t)},{\omega }^2=k_z^2+m_{\psi }^2,$$

(184)

the tidal components become

$$\begin{array}{cc}\hfill R_{x0x0}& =R_{y0y0}=-{\displaystyle \frac{1}{2}}{\omega }^2\delta R,\hfill \end{array}$$

(185)

$$\begin{array}{cc}\hfill R_{z0z0}& =+{\displaystyle \frac{1}{2}}{m}_{\psi }^2\delta R.\hfill \end{array}$$

(186)

The geodesic deviation equations

$${\ddot{\xi }}_i=-R_{i0j0}{\xi }^j$$

(187)

then give

$$\begin{array}{cc}\hfill \ddot{X}& =-{\displaystyle \frac{1}{2}}{\omega }^2\delta RX,\hfill \end{array}$$

(188)

$$\begin{array}{cc}\hfill \ddot{Y}& =-{\displaystyle \frac{1}{2}}{\omega }^2\delta RY,\hfill \end{array}$$

(189)

$$\begin{array}{cc}\hfill \ddot{Z}& =-{\displaystyle \frac{1}{2}}{m}_{\psi }^2\delta RZ.\hfill \end{array}$$

(190)

These equations show that the extra scalar degree of freedom induces both a transverse breathing mode (in the X and Y directions) and a longitudinal mode (in the Z direction). This is precisely the expected polarization content for a massive scalar mode.

In pure GR, where only the transverse-traceless tensor $h_{ij}^{TT}$ is present, $\delta R=0$ and the scalar-induced contributions vanish; only the familiar ⊕ and ⊗ tensor modes remain. In metric $f\left(R\right)$ gravity, the nonzero $\delta R$ generates additional breathing and longitudinal polarizations on top of the GR tensor modes.

7.2. Geodesic Deviation on a de Sitter FRW Background

We now sketch how the same polarization structure arises when the calculation is performed directly on the de Sitter background without passing explicitly to a Minkowski patch. In spatially flat FRW coordinates, the de Sitter metric can be written as

$$d{s}^2=-d{t}^2+a^2\left(t\right){\delta }_{ij}d{x}^{i}d{x}^j,a\left(t\right)=e^{H_{d}t},$$

(191)

with constant Hubble parameter $H_d=\sqrt{R_d/12}$ in four dimensions.

We consider small perturbations around this background in Newtonian gauge. Restricting initially to the scalar sector, the perturbed metric takes the form

$$d{s}^2=-(1+2\Phi )d{t}^2+a^2\left(t\right)(1-2\Psi ){\delta }_{ij}d{x}^{i}d{x}^j,$$

(192)

where equality of the Bardeen potentials $\Phi =\Psi$ holds for metric $f\left(R\right)$ gravity on a de Sitter background (Section 4), with

$$\Phi =\Psi =-{\displaystyle \frac{\delta R}{6m_{\psi }^2}}.$$

(193)

Thus the additional scalar degree of freedom is directly encoded in both the temporal and isotropic spatial perturbations of the metric.

To relate the geodesic deviation equation to the detector frame, it is convenient to introduce an orthonormal tetrad adapted to a comoving observer,

$$\begin{array}{cc}\hfill {e_{\hat{0}}}^{\mu }& =(1,0,0,0),\hfill \end{array}$$

(194)

$$\begin{array}{cc}\hfill {e_{\hat{i}}}^{\mu }& ={\displaystyle \frac{1}{a\left(t\right)}}{{\delta }_i}^{\mu },\hfill \end{array}$$

(195)

so that physical (proper) spatial separations are measured with hatted indices. In this orthonormal frame the geodesic deviation equation takes the form

$${\displaystyle \frac{d^2{\xi }^{\hat{i}}}{d{t}^2}}=-{R^{\hat{i}}}_{\hat{0}\hat{j}\hat{0}}{\xi }^{\hat{j}}.$$

(196)

The tidal tensor splits naturally into a background de Sitter contribution and a perturbation induced by scalar and tensor modes,

$${R^{\hat{i}}}_{\hat{0}\hat{j}\hat{0}}={R^{\hat{i}}}_{\hat{0}\hat{j}\hat{0}}{|}_{\mathrm{dS}}+\delta {R^{\hat{i}}}_{\hat{0}\hat{j}\hat{0}}.$$

(197)

For the spatially flat de Sitter background, the nonvanishing Christoffel symbols are ${{\Gamma }^0}_{ij}=a\dot{a}{\delta }_{ij}$ and ${{\Gamma }^i}_{0j}=H_d{{\delta }^i}_j,$ which follow directly from the FRW line element. From these, the coordinate-basis Riemann component relevant for geodesic deviation is

$${R^i}_{0j0}={\displaystyle \frac{\ddot{a}}{a}}{{\delta }^i}_j.$$

(198)

For exact de Sitter expansion $a\left(t\right)=e^{H_{d}t}$, one has $\ddot{a}/a=H_d^2$, yielding

$${R^i}_{0j0}{|}_{\mathrm{dS}}=H_d^2{{\delta }^i}_j.$$

(199)

Projecting onto the orthonormal tetrad, ${R^{\hat{i}}}_{\hat{0}\hat{j}\hat{0}}={e^{\hat{i}}}_k{e_{\hat{j}}}^{\ell }{R^k}_{0\ell 0},$ the factors of $a\left(t\right)$ from the tetrads cancel those implicit in the metric, leaving

$${R^{\hat{i}}}_{\hat{0}\hat{j}\hat{0}}{|}_{\mathrm{dS}}=H_d^2{{\delta }^i}_j.$$

(200)

We now include perturbations. Restoring both scalar and tensor modes, the perturbed FRW metric may be written as

$$d{s}^2=-(1+2\Phi )d{t}^2+a^2\left(t\right)\left[(1-2\Psi ){\delta }_{ij}+h_{ij}^{\mathrm{TT}}(t,\mathbf{x})\right]d{x}^{i}d{x}^j,$$

(201)

where $h_{ij}^{\mathrm{TT}}$ denotes the transverse–traceless tensor perturbation. Introducing the Minkowski metric ${\eta }_{\mu \nu }=\mathrm{diag}(-1,1,1,1)$, all perturbations can be collected into a single tensor

$$h_{00}=2\Phi ,h_{0i}=0,h_{ij}=2\Psi {\delta }_{ij}+h_{ij}^{\mathrm{TT}},$$

(202)

so that the metric assumes the conformal form

$$g_{\mu \nu }(t,\mathbf{x})=a^2\left(t\right)\left[{\eta }_{\mu \nu }+h_{\mu \nu }(t,\mathbf{x})\right].$$

(203)

Expanding the Riemann tensor to first order, $R_{\mu \nu \rho \sigma }=R_{\mu \nu \rho \sigma }^{\left(0\right)}+\delta R_{\mu \nu \rho \sigma }\left[h\right],$ the linearized part depends only on derivatives of $h_{\mu \nu }$. Because the conformal factor multiplies both background and perturbation, one finds

$$\delta R_{\mu \nu \rho \sigma }\left[g\right]=a^2\left(t\right)\delta R_{\mu \nu \rho \sigma }^{\left(\mathrm{Mink}\right)}\left[h\right].$$

(204)

Raising indices and projecting onto the orthonormal tetrad yields

$$\delta {R^{\hat{i}}}_{\hat{0}\hat{j}\hat{0}}={\displaystyle \frac{1}{a^2\left(t\right)}}\delta R_{i0j0}^{\left(\mathrm{Mink}\right)},$$

(205)

where $\delta R_{i0j0}^{\left(\mathrm{Mink}\right)}$ is precisely the tidal matrix obtained in Subsection 7.1.

Therefore, for the massive scalar propagating mode, we can directly carry over the Minkowski result, with careful attention to the overall sign,

$$\begin{array}{cc}\hfill \delta {R^{\hat{x}}}_{\hat{0}\hat{x}\hat{0}}& =\delta {R^{\hat{y}}}_{\hat{0}\hat{y}\hat{0}}=-{\displaystyle \frac{1}{2a^2}}{\omega }^2\delta R,\hfill \end{array}$$

(206)

$$\begin{array}{cc}\hfill \delta {R^{\hat{z}}}_{\hat{0}\hat{z}\hat{0}}& =+{\displaystyle \frac{1}{2a^2}}{m}_{\psi }^2\delta R,\hfill \end{array}$$

(207)

for a monochromatic mode, and similarly for a generic wave packet using the Klein–Gordon equation (180). The overall factor $1/a^2\left(t\right)$ dilutes the tidal amplitude due to cosmic expansion, while leaving the polarization pattern unchanged.

The geodesic deviation equations for physical separations ${\xi }^{\hat{i}}$ are therefore

$$\begin{array}{cc}\hfill {\ddot{\xi }}^{\hat{X}}& =-\left[H_d^2+\delta {R^{\hat{x}}}_{\hat{0}\hat{x}\hat{0}}\right]{\xi }^{\hat{X}},\hfill \end{array}$$

(208)

$$\begin{array}{cc}\hfill {\ddot{\xi }}^{\hat{Y}}& =-\left[H_d^2+\delta {R^{\hat{y}}}_{\hat{0}\hat{y}\hat{0}}\right]{\xi }^{\hat{Y}},\hfill \end{array}$$

(209)

$$\begin{array}{cc}\hfill {\ddot{\xi }}^{\hat{Z}}& =-\left[H_d^2+\delta {R^{\hat{z}}}_{\hat{0}\hat{z}\hat{0}}\right]{\xi }^{\hat{Z}}.\hfill \end{array}$$

(210)

The background term produces the isotropic de Sitter expansion, while the wave-induced part reproduces the same transverse breathing and longitudinal pattern as in the local Minkowski analysis. Thus, cosmological expansion modifies amplitudes but does not change the polarization content.

7.3. Polarization Classification via $R_{i0j0}$

The geodesic deviation equations derived above show explicitly that the additional scalar degree of freedom in our $f\left(R\right)$ model produces both breathing and longitudinal motion of test particles. For completeness, we now review a more formal method to classify the polarization modes using the components of $R_{i0j0}$, following [1,40].

In a local inertial (Minkowski) patch of the spacetime, the perturbed metric may be written in terms of scalar, vector, and tensor perturbations as

$$d{s}^2=-(1+2\Phi )d{t}^2+2E_{i}dtd{x}^i+\left[(1-2\Psi ){\delta }_{ij}+h_{ij}^{TT}\right]d{x}^{i}d{x}^j,$$

(211)

where $\Phi$ and $\Psi$ are the scalar Bardeen potentials (with $\Phi =\Psi$ in the present context), $E_i$ encodes the vector (shear) perturbations, and $h_{ij}^{TT}$ is the transverse–traceless tensor mode. To linear order in the perturbations, the Riemann tensor components entering the geodesic deviation equation are

$$R_{i0j0}={\partial }_i{\partial }_j\Psi -{\displaystyle \frac{1}{2}}{\delta }_{ij}{\partial }_0^2\Phi -{\displaystyle \frac{1}{2}}\left({\partial }_0{\partial }_i{E}_j+{\partial }_0{\partial }_j{E}_i\right)-{\displaystyle \frac{1}{2}}{\partial }_0^2{h}_{ij}^{TT}.$$

(212)

The six possible GW polarization modes can be encoded by writing the tidal tensor $R_{i0j0}$ as a symmetric $3\times 3$ matrix,

$$R_{i0j0}=\left[\begin{array}{ccc}{P}_4+P_6& P_5& P_2\\ P_5& -P_4+P_6& P_3\\ P_2& P_3& P_1\end{array}\right],$$

(213)

where $P_1,\dots ,P_6$ are the six independent polarization amplitudes (scalar longitudinal, two vector modes, two tensor modes, and scalar breathing). They correspond to the six standard polarization patterns shown in Figure 1.2

For a plane wave propagating along the z direction, comparison of (212) with the matrix form (213) yields

$$\begin{array}{cc}\hfill P_1& ={\partial }_3{\partial }_3\Psi -{\displaystyle \frac{1}{2}}{\partial }_0{\partial }_0\Phi ,\hfill \end{array}$$

(214)

$$\begin{array}{cc}\hfill P_2& ={\displaystyle \frac{1}{2}}{\partial }_0{\partial }_3{E}_1,\hfill \end{array}$$

(215)

$$\begin{array}{cc}\hfill P_3& ={\displaystyle \frac{1}{2}}{\partial }_0{\partial }_3{E}_2,\hfill \end{array}$$

(216)

$$\begin{array}{cc}\hfill P_4& =-{\displaystyle \frac{1}{2}}{\partial }_0{\partial }_0{h}_{11}^{TT},\hfill \end{array}$$

(217)

$$\begin{array}{cc}\hfill P_5& =-{\displaystyle \frac{1}{2}}{\partial }_0{\partial }_0{h}_{12}^{TT},\hfill \end{array}$$

(218)

$$\begin{array}{cc}\hfill P_6& =-{\displaystyle \frac{1}{2}}{\partial }_0{\partial }_0\Phi .\hfill \end{array}$$

(219)

Here $E_1$ and $E_2$ encode the vector (shear) polarizations, $h_{ij}^{TT}$ represents the usual ⊕ and ⊗ tensor modes, and $\Phi ,\Psi$ are the scalar Bardeen potentials.

In metric $f\left(R\right)$ gravity we have generic vector perturbations $V_i=0$ in vacuum, so the vector modes are absent and $P_2=P_3=0$. The tensor modes $P_4$ and $P_5$ coincide with those of GR and correspond to the ⊕ and ⊗ polarizations. The remaining scalar modes are encoded in $P_1$ (longitudinal mode, involving both $\Phi$ and $\Psi$) and $P_6$ (breathing mode, involving only $\Phi$). Because $\Phi =\Psi \ne 0$ in our model and are related to the additional scalar degree of freedom via

$$\Phi =\Psi =-{\displaystyle \frac{\delta R}{6m_{\psi }^2}},$$

(220)

both $P_1$ and $P_6$ are nonzero, confirming that the model exhibits a mixed longitudinal and breathing scalar polarization in addition to the two tensor polarizations.

In summary, the geodesic deviation analysis—both in the local Minkowski patch and on the full de Sitter background—shows that the metric $f\left(R\right)$ model $f\left(R\right)=R+\alpha R^2-2\Lambda$ supports:

• two massless tensor modes $(\oplus ,\otimes )$, identical to those of GR;
• one massive scalar mode (the massive scalar propagating mode), which decomposes into a transverse breathing polarization and a longitudinal polarization along the propagation direction.

This pattern agrees with the general expectation for metric $f\left(R\right)$ gravity and provides the polarization content against which current and future GW observations can test this class of models.

8. Conclusion and Future Outlook

In this work we developed a unified and fully gauge-invariant analysis of gravitational-wave polarizations in metric $f\left(R\right)$ gravity, with particular emphasis on the modified Starobinsky model $f\left(R\right)=R+\alpha R^2-2\Lambda .$ Working on a constant-curvature de Sitter background, we reformulated the linearized field equations in terms of Bardeen gauge-invariant variables and the scalar curvature perturbation $\delta R$, thereby making the massive scalar propagating mode manifest. By deriving the Klein–Gordon equation for $\delta R$ directly from the perturbed trace equation, we verified that the scalar mode behaves as a massive propagating field with mass $m_{\psi }^2=1/\left(6\alpha \right)$ on the de Sitter background. This establishes the scalar curvature perturbation $\delta R$ as the source of the additional breathing and longitudinal polarizations absent in General Relativity.

We complemented the Bardeen-variable analysis with a full $(3+1)$ decomposition of the perturbed Ricci tensor, including the presence of matter sources. This approach revealed explicitly how scalar, vector, and tensor perturbations enter the modified field equations and how the scalar sector departs from its GR behavior. In particular, the decomposition demonstrated that (i) the vector sector remains nondynamical and identical to that of GR, (ii) the tensor sector continues to satisfy the standard transverse–traceless wave equation, and (iii) all modifications are encoded in the scalar sector through the dynamical curvature perturbation $\delta R$. The resulting coupled equations for $\Phi$, $\Psi$, and $\delta R$ illustrate the origin of the gravitational slip, modified Poisson equation, and scale-dependent evolution of cosmological perturbations characteristic of $f\left(R\right)$ models.

A complementary geodesic-deviation analysis was carried out in both the local-Minkowski patch of de Sitter spacetime and in the fully covariant de Sitter background. In both cases, the tidal tensor $R_{i0j0}$ depends on the scalar curvature perturbation $\delta R$ and yields the characteristic polarization pattern: two tensor modes (⊕ and ⊗), a breathing mode, and a longitudinal mode. This agrees with the general classification of metric theories admitting up to six polarizations and verifies, by two independent methods, that metric $f\left(R\right)$ gravity predicts exactly three observable polarization sectors: two tensor and one massive scalar.

From a cosmological perspective, the mass $m_{\psi }$ of the massive scalar propagating mode sets the transition scale between GR-like behavior at high wavenumbers and modified gravity effects on large scales. Because the same massive scalar propagating mode controls the late-time background evolution, the growth rate of structure, and the propagation of gravitational waves, future multi-probe observations—combining large-scale structure, weak lensing, CMB anisotropies, pulsar-timing arrays, and gravitational-wave observatories—provide a coherent program for testing the viability of $f\left(R\right)$ gravity on both astrophysical and cosmological scales.

Future Outlook

Several natural extensions follow from the framework developed here:

• Beyond de Sitter backgrounds: The methods employed here can be generalized to slowly evolving FLRW backgrounds, permitting a direct link between gravitational-wave propagation and the time dependence of the mass $m_{\psi }$ of the massive scalar propagating mode in realistic cosmologies.
• Mode mixing and GW propagation: A next step is the study of mode mixing between the tensor and scalar sectors, including amplitude damping and potential dispersion effects in late-time, low-density environments.
• Constraints from forthcoming surveys: Current and future missions (Euclid, LSST, SKA, LISA, pulsar-timing arrays) will significantly improve constraints on gravitational slip, the mass of the scalar mode, and the scale-dependent growth of cosmological perturbations. The gauge-invariant formalism presented here is well suited for connecting theoretical predictions with these upcoming datasets.
• Extension to broader modified-gravity families: The techniques developed in this paper—decomposition of the perturbed Ricci tensor, isolation of the massive scalar propagating mode, fully covariant GW polarization extraction, and geodesic-deviation analysis—can be applied to more general higher-curvature theories such as $f\left(G\right)$ gravity, scalar–tensor Horndeski theories, and Einstein–dilaton–Gauss–Bonnet models.

Overall, the combination of gauge-invariant SVT analysis, Ricci-tensor decomposition, and geodesic deviation provides a robust framework for identifying and interpreting the polarization content of gravitational waves in metric $f\left(R\right)$ gravity. This establishes a consistent pathway for future observational tests capable of distinguishing GR from its simplest and most theoretically motivated extensions.