Gravitational Wave Propagation in Starobinsky Inflationary Model

In recent decades, modified gravity theories have emerged as a significant alternative for explaining cosmological phenomena tied to fundamental physics [1, 2]. These theories are motivated by the need to expand our understanding of gravity to address problems that General Relativity (GR) cannot solve without invoking dark matter and dark energy, such as the accelerated expansion of the universe [3, 4] and the dynamics of large-scale structures [5, 6]. Among these, $f(R)$ theories generalize the Einstein-Hilbert Lagrangian by introducing nonlinear terms in the Ricci scalar $R$, offering new perspectives on the curvature effects of spacetime across various gravitational regimes [7, 8, 9]. Viable models within $f(R)$ gravity theories [10, 11, 12], show promise in addressing both cosmological and astrophysical phenomena beyond standard GR [13, 14, 15, 16]. These models introduce a broader range of gravitational behaviors that remain consistent with observations while avoiding many of the issues found in alternative theories [5, 17, 18, 19]. Consequently, $f(R)$ gravity provides a compelling framework for exploring modifications to GR, with applications ranging from large-scale cosmic expansion to localized strong gravitational fields, such as black holes, where they predict modified horizon structures and unique stability properties [20, 21, 22].
In this sense, Starobinsky’s model [23], $f(R)=R+\alpha R^{2}$, with $\alpha$ being a parameter; initially proposed within the framework of cosmic inflation, stands as one of the earliest alternatives to scalar-field inflation models [24, 25]. This model is characterized by the inclusion of an $R^{2}$ term, introducing new degrees of freedom in the field equations, resulting in additional modes that can impact the early Universe’s behavior [26, 24], particularly during inflation, where quantum effects may play a significant role [27, 28], as well as in producing perturbations that affect the evolution of large-scale structures such as galaxies and galaxy clusters.
With the successful detection of gravitational waves [29], the framework of $f(R)$ gravity could gain significant attention [30, 31, 32, 33], since important works have been developed in order to explore how these modified theories could influence wave propagation, polarization modes, and potential deviations from predictions in GR [34, 35, 36, 37, 38, 39]. Hopefully the next-generation detectors will allow to observe a wider range of astrophysical events, providing further evidence of gravitational waves and their properties [40, 41]. These advancements may reveal distinctive features in gravitational wave signals that could indicate the influence of $f(R)$ modifications, potentially setting them apart from standard relativistic predictions.
Considering this context, we focus this work on linearizing the field equations of $f(R)$ gravity in the Starobinsky model within the weak-field approximation, analyzing the propagation of perturbations under this modified gravity theory to identify distinctive differences between this model and GR.

This paper is organized as follows: in Section II, we linearize the field equations of the $f(R)$ theory for the model $R+R^{2}/(6m^{2})$ in the weak-field regime, expressing the field equations in terms of the perturbation $\bar{h}_{\mu\nu}$ relative to the background Minkowski metric. In Section III, we solve the trace equation using Green’s functions and auxiliary fields. In Section IV, we calculate the components of $\bar{h}_{\mu\nu}$, considering both the massive contribution induced by the quadratic term and the propagation effects of the field. Then, in Section V, we calculate the perturbation found for a binary star system, and its differences compared to GR. Finally, in Section VI, we discuss and analyze the results.

The $f(R)$ theory is an extension of GR, reformulating the Einstein-Hilbert action in terms of a nonlinear function of the Ricci curvature scalar $R$

where $S_{M}$ is the contribution from matter and energy, and $\kappa=8\pi G$. The field equations governing the dynamics of the metric tensor $g_{\mu\nu}$, relate the spacetime geometry to the distribution of matter and energy, and are derived by evaluating the critical points of the action $S$. In the metric formalism are written as

where $f=f(R)$, $F=F(R)=f^{\prime}(R)$, $F_{;\mu\nu}=\nabla_{\mu}\nabla_{\nu}F$ is the covariant derivative, and with the definition of the D’Alembertian operator $\Box=\nabla_{\sigma}\nabla^{\sigma}=g^{\sigma\rho}\nabla_{\sigma}\nabla_{\rho}$, $F_{;\sigma}^{;\sigma}=\Box F$. Moreover, the energy momentum tensor is defined as

Equations (2) incorporate the additional degrees of freedom introduced by the function $f(R)$, resulting in a framework where the curvature itself acts as a dynamical field, allowing the theory to be more flexible in addressing a broader range of gravitational phenomena. The trace of field equations (2) is obtained by multiplying the equations by the metric tensor

where $T=T_{\mu\nu}g^{\mu\nu}$. Now, we will take into account the Starobinsky quadratic model [23]

where the parameter $m$ is identified with the inflaton mass [24, 25]. With this function of $f(R)$, field equations take the form

where we have used the Einstein tensor

Considering minimal deviations, $h_{\mu\nu}$, from flat spacetime, that is, small fluctuations around the background Minkowski metric, $\eta_{\mu\nu}=\text{diag}(-1,1,1,1)$, so that the metric tensor can be expressed as

such that $h_{\mu\nu}\ll 1$. In this weak-field approximation, the quadratic terms are not taking into account in equations (6), therefore they are simplified to

and using the definition of the trace reverse tensor

through which it is possible to simplify the Einstein tensor, that is

and by choosing the Lorentz gauge condition

we arrive to the linearized field equations

that is, in the formulation of $f(R)$ gravity using the Starobinsky model, the additional term $1/(3m^{2})$ acts as a correction that introduces mass effects associated with the extra degrees of freedom, particularly with the scalar field incorporated through the $R^{2}$ term. This expression is more complex than in GR, and its deviation reflects the additional modes due to the effective mass $m$. However, in the limit $m\to\infty$ (so that $f(R)\to R$ in Eq. (5)), the additional term vanishes, and the standard wave equation of GR is recovered. Taking the trace of Eq. (13) we obtain

or using the D’Alembertian operator in flat space $\Box=-\partial_{t}^{2}+\nabla^{2}$,

thus, field equations (13) are rewritten as

and returning to the initial field $h$

Equation (15) describes the evolution of a field $\bar{h}(x^{\sigma})$ under the influence of a source $T(x^{\sigma})$, with a mass term $m^{2}$. Due to its structure, it can be divided into two stages, each with a distinct physical meaning and mediated by a fictitious field. To see this, suppose $H(x^{\sigma})$ is a fictitious field generated by the real source $T(x^{\sigma})$; it can be regarded as a perturbation or an “auxiliary wave” that propagates as a massless field:

thus, by Eq. (15), we can think of $H(x^{\sigma})$ as a fictional source that produces $\bar{h}(x^{\sigma})$ through the inhomogeneous Klein-Gordon equation

where we see that the propagation of the perturbation is modified by the presence of the mass term. Interpreting Eq (15) in two stages, reflects the fact that modified gravity can introduce interactions or auxiliary fields that mediate the effective gravitational response.
To solve these coupled differential equations, Green’s functions can be used, as they hold significant physical meaning regarding the propagation of perturbations in spacetime. Suppose $G_{1}$ is the solution of

where $\Box_{x}$ denotes derivative with respect coordinates $x^{\sigma}$ and $\delta^{(4)}(x^{\sigma}-z^{\sigma})$ is the four dimensional Dirac delta function. Therefore, we can express $\bar{H}(x^{\sigma})$ as the sum of all signals to the past generated by the source $T$ in the point $z^{\sigma}$, i.e.

and equivalently, if $G_{2}(x^{\sigma}-y^{\sigma})$ satisfies

we can express the field $\bar{h}$ at the point $x^{\sigma}$, as

or in terms of the source $T$,

in other words, the function $G_{1}$ represents the spacetime response to the source $T$ located at the point $z^{\sigma}$. It describes the propagation of the perturbation from $z^{\sigma}$ to the point $y^{\sigma}$, generating the field $H$ there. In turn, $G_{2}$ propagates the influence of $H$ both on and inside the light cone (at subluminal speeds) to the point $x^{\sigma}$, where the field $\bar{h}$ is measured.
Note that at the limit $m\to\infty$ in Eq. (19), $\bar{h}(x^{\sigma})=H(x^{\sigma})$, and Eq. (18) takes the same form as in the GR case. Equivalently,

by virtue of the finiteness of $T$. From a physical perspective, in this limit, the propagation of the field becomes constrained, and the Green’s function is concentrated in an infinitesimally small region around $y^{\sigma}$. Mathematically, this is expressed as $G_{2}(x^{\sigma}-y^{\sigma})\to-\delta^{(4)}(x^{\sigma}-y^{\sigma})$, so that the influence of the massive field is localized at the point $y^{\sigma}$ and rapidly diminishes outside of it.
The solution of Eq. (20) is the well known retarded Green function

where the theta function, $\theta(x)=1$ only for $x>0$, $s_{xz}^{2}=t_{xz}^{2}-r_{xz}^{2}$, $t_{xz}=x^{0}-z^{0}$ and $r_{xz}^{2}=\sum_{i=1}^{3}(x^{i}-z^{i})^{2}$. To ensure the propagator is non-zero in Minkowski spacetime, the signal generated at $z^{\sigma}$ must lie in the past of the event $x^{\sigma}$, with a positive proper time interval between them, thereby upholding causality. Likewise the solution For Eq. (22) is

where

$J_{1}(x)$ is the Bessel function of the first kind, and $x^{0}>y^{0}>z^{0}$, from which

It is worth noting here that Eq. (15) can be decomposed from a different perspective: suppose the field $\bar{h}(x^{\sigma})$ is generated by a fictitious source $K(x^{\sigma})$ through the wave equation

where the source $K$ is produced by the real source $T(x^{\sigma})$ according to

which, again, is the inhomogeneous Klein-Gordon equation. However, when expressing $T$ in integral form using the delta function

and for Eq. (21), we found the relation between the fields $H$ and $K$,

this equation indicates that the two fictitious fields are not independent; rather, the physical field $\bar{h}$ produced by the real source $T$ can be interpreted as being mediated by the coupled fields $K$ and $H$. In this scheme, $T$ generates $H$, $H$ generates $K$, and $K$ in turn produces $\bar{h}$.
Eq. (31) can be written as

and from Eq. (30),

since the physics described by the physical field $\bar{h}$ must be invariant in both approaches, the last integral must be equivalent to (24). Therefore, the integrands can be equated

which can be showed from the definitions of the Green functions. This equality indicates a commutation symmetry between the Green’s functions $G_{1}$ and $G_{2}$, suggesting that when applied sequentially for field propagation, the order in which they are used does not affect the final result. With this, we arrive to the important relation

this expression describes how the propagation of signals from point $z^{\sigma}$ to $x^{\sigma}$ depends on a linear combination of the functions $G_{1}$ and $G_{2}$, where the mass $m$ controls the types propagation. It shows that the propagation of the perturbation includes the influence of both a massless and a massive field. Note that we could also have arrived at this integral from Eq. (22),

With Eq. (37) we can write the trace Eq. (35) as

where it can be seen the first term that corresponds to GR, and a second term that disappears when $m\to\infty$. However, the trace can be written in a simpler way by using the bulk term of $G_{2}$, since

therefore

note that the integrand has an oscillatory behavior due to the Bessel function, whose frequency increases with $m$, as does its amplitude, however as the interval $s_{xz}$ increases, the intensity of the perturbation decreases.

Now that we have established the dependence of $\bar{h}(x^{\sigma})$ on $T(x^{\sigma})$ and the function $g(x^{\sigma}-z^{\sigma})$, we can solve Eq. (13) for the field $\bar{h}_{\mu\nu}$, as a combination of derivatives and integrals of the stress-energy tensor, along with the Green’s functions. In this sense, note that

so that Eq. (41) takes the form

and

by substituting these two results into Eq. (13), we obtain a wave equation for the tensor $\bar{h}_{\mu\nu}$

where the derivatives are respect to $x^{\sigma}$. Therefore we can interpret that the perturbation is produced by a modified, or effective, stress-energy tensor

Using the Green’s function $G_{1}$ again, we can write the solution to the differential equation (45), that is

where it should be noted that the derivatives are now with respect to $y^{\sigma}$. Likewise, it is observed that at $m\to\infty$, the terms associated with the trace $T$ disappear, and the case of GR is obtained.
Since $\mathcal{T}_{00}$ incorporates the effects of the modified theory, including the mass terms, we can express the quadrupole moment tensor as

This approach is advantageous because it allows us to integrate directly over the effective energy density, revealing how the source responds under the influence of the additional terms. Therefore, from (46)

or equivalently

this expression shows how the effective energy density is the sum of the direct component of the source, $T_{00}$, plus a correction that fades as $m$ increases, but is responsible for the finite-range effects. In the same way, the spacelike components of the perturbation can be written as

where the retarded time is defined as $t_{r}=t-r$, that is, the time in which the wave is observed after being generated by the system.

As an illustrative example, let us consider a binary star system, as described in [42]. This system is an ideal candidate for calculating the components of the quadrupole moment tensor (48), as well as the perturbation (51), allowing us to observe the additional terms that arise in the modified theory and compare them with the predictions of GR. Let the positions of two stars of mass $M$, at a distance $r$ from the point of observation, orbiting each other due to their mutual gravitational attraction, at a distance $R$ from the system’s center of mass, and with orbital angular frequency $\Omega$, be

and

from which

and the trace

where

however the energy density can be expressed as follows

where

As $m$ increases, the Green’s function $G_{2}(x^{\sigma}-z^{\sigma})$ becomes concentrated around $x^{\sigma}=z^{\sigma}$, restricting the propagation of the source’s effects to its origin point. Consequently, $\Box G_{2}(x^{\sigma}-z^{\sigma})$ tends to zero outside of this neighborhood and can be considered negligible. Thus, it is reasonable to approximate

in the large-mass regime. In this way, the components of the quadrupole moment tensor are

where

and the gravitational perturbation is

and

Therefore, it is observed that the factor $\xi$ introduces modifications to the amplitude of the perturbation due to the mass of the field. When $m\to\infty$, $\xi\to 1$, recovering the relativistic result. However, the same scenario arises when $R=\Omega^{-1}$, in which case $\xi=1$.
It is important to note that in this kind of systems the angular frequency and the orbital radius are closely related. For example, from a classical perspective

and

which implies that for $R=GM/4$, $\xi=1$, thus, if $R$ is greater than this value, the amplitude of the wave decreases, whereas if $R$ is smaller, the amplitude increases, which could yield measurable differences from GR. Thus, if we are interested in detecting these corrections, compact and high-frequency binary systems would be the most promising candidates.

This work explores the linearized field equations in the Starobinsky $R+R^{2}/(6m^{2})$ model, deriving an effective matter-energy distribution for perturbations through an auxiliary field, $H(x^{\sigma})$, that links the trace equation (15), the Klein-Gordon equation (19), and wave equation (18). This two-step propagation reflects how $H(x^{\sigma})$ mediate the gravitational effect in an indirect way. Equations (18) and (19) were solved using Green’s functions, (26) and (27). It should be noted that $G_{1}(x^{\sigma}-z^{\sigma})$ models the direct propagation of the perturbation at light speed from the source, while $G_{2}(x^{\sigma}-z^{\sigma})$ allows for the additional propagation of the field $H(x^{\sigma})$ with massive characteristics, creating both light-speed and slower components in the gravitational response measured at $x^{\sigma}$. For this massive field, the Green’s function, $G_{2}$ shows exponential decay: the larger the mass, the more localized the field’s effect near the source, in a similar way to the Yukawa potential: a massive field’s influence decays with distance, reflecting a finite interaction range. In the limit $m\to\infty$, $G_{2}$ approaches a delta-like behavior, confining the field entirely to the source point, in contrast to the infinite range of a massless field.
We found the trace $\bar{h}(x^{\sigma})$, Eq. (35), which could be simplified due to Eq. (36). This symmetry can be understood in terms of the commutativity in combining these solutions for the propagation of perturbations between the points $z^{\sigma}\to y^{\sigma}\to x^{\sigma}$. This property could be useful in coupled field theories, where propagation and interaction between points in spacetime can be decomposed into multiple steps.
The solution found for the trace $\bar{h}(x^{\sigma})$, Eq. (41), is written in terms of the of the Bessel function, suggesting that the propagation of $\bar{h}$ exhibits an oscillatory and finite-range nature. This is because the mass term restricts the propagation to a scale defined by $1/m$. At greater distances, the Bessel function decays, reflecting the effect of the mass in limiting the influence of the source. Additionally, the term $J_{1}(ms_{xz})/s_{xz}$ has a decay controlled by $s_{xz}$, indicating that the intensity of the perturbation decreases with the distance between $x^{\sigma}$ and $z^{\sigma}$. As the mass increases, the term $J_{1}(ms_{xz})/s_{xz}$ narrows and concentrates near $s_{xz}=0$, while its amplitude asymptotically grows toward infinity, exhibiting a delta-like behavior in the limit $m\to\infty$. This transition of the function $g(x^{\sigma})$ toward $G_{1}(x^{\sigma})$ in this limit reflects that, as the massive field becomes extremely heavy, its influence on the effective field becomes confined to the source point.
Once the trace was replaced in the field equations, it could be observed that the field $\bar{h}(x^{\sigma})$ was defined by a wave equation, (45), with an effective stress-energy tensor, Eq. (46). This modified source term incorporates contributions not only from the actual energy-momentum tensor $T_{\mu\nu}$ but also from additional terms mediated by the Green’s functions and differential operators, which capture the effects of the massive field components on the propagation of $\bar{h}_{\mu\nu}$. This effective source modifies the propagation characteristics of the wave equation, reflecting the impact of mass terms and other higher-order effects on the field.
We found the solution for the tensor $\bar{h}_{\mu\nu}$, Eq. (47). The terms within the integral, represent both massless and massive field contributions: $T_{\mu\nu}$ corresponds to the massless component, while the terms associated with the trace $T$ are modulated by the mass $m$. Consequently, the effect of $T$ on $\bar{h}_{\mu\nu}$, can reach over long distances for small masses or becoming locally confined for large masses. In fact, as $m\to\infty$, the trace terms vanish, recovering the case of GR.
To calculate the quadrupole moment tensor and compare it with the GR case, in the large-mass regime, the approximation (59), to energy-density was applied, using the fact that $\Box G\approx 0$. In this scenario, the $\mathcal{T}_{00}$ component is dominated by the standard energy density $T_{00}$ along with a correction term that depends on the time derivative of the trace, which also decreases as $m$ increases.
Subsequently we calculated the quadrupole moment and perturbation tensor for a binary star system. This results suggest that high-frequency binary systems could provide measurable deviations from GR due to the modified gravity effects, which may perhaps be observable with the next generation of gravitational wave observatories, offering new insights into the nature of gravity and its effects on astrophysical phenomena.